Medical Imaging Interaction Toolkit  2018.4.99-08619e4f
Medical Imaging Interaction Toolkit
itkHessianMatrixEigenvalueImageFilter.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #ifndef ITKHESSIANMATRIXEIGENVALUEIMAGEFILTER_CPP
14 #define ITKHESSIANMATRIXEIGENVALUEIMAGEFILTER_CPP
15 
17 #include <itkImageRegionIterator.h>
18 #include <vigra/tensorutilities.hxx>
19 #include <vigra/convolution.hxx>
20 #include <mitkCLUtil.h>
21 
22 
23 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
25 {
26  Superclass::GenerateOutputInformation();
27  this->GetOutput(0)->SetDirection(this->GetInput()->GetDirection());
28  this->GetOutput(0)->SetSpacing(this->GetInput()->GetSpacing());
29  this->GetOutput(0)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
30  this->GetOutput(0)->Allocate();
31 
32  this->GetOutput(1)->SetDirection(this->GetInput()->GetDirection());
33  this->GetOutput(1)->SetSpacing(this->GetInput()->GetSpacing());
34  this->GetOutput(1)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
35  this->GetOutput(1)->Allocate();
36 
37  this->GetOutput(2)->SetDirection(this->GetInput()->GetDirection());
38  this->GetOutput(2)->SetSpacing(this->GetInput()->GetSpacing());
39  this->GetOutput(2)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
40  this->GetOutput(2)->Allocate();
41 }
42 
43 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
45 {
46 
47  typedef typename TInputImageType::PixelType InputPixelType;
48 
49  typename TInputImageType::RegionType region = this->GetInput()->GetLargestPossibleRegion();
50  unsigned int xdim = region.GetSize(0);
51  unsigned int ydim = region.GetSize(1);
52  unsigned int zdim = region.GetSize(2);
53 
54  typename TInputImageType::Pointer maske_input = TInputImageType::New();
55  maske_input->SetDirection(this->GetInput()->GetDirection());
56  maske_input->SetSpacing(this->GetInput()->GetSpacing());
57  maske_input->SetRegions(this->GetInput()->GetLargestPossibleRegion());
58  maske_input->Allocate();
59 
60  itk::ImageRegionConstIterator<TInputImageType> iit(this->GetInput(), this->GetInput()->GetLargestPossibleRegion());
61  itk::ImageRegionIterator<TInputImageType> miit(maske_input, maske_input->GetLargestPossibleRegion());
62  itk::ImageRegionConstIterator<TMaskImageType> mit(m_ImageMask, m_ImageMask->GetLargestPossibleRegion());
63 
64  while(!miit.IsAtEnd())
65  {
66  if(mit.Value()!=0)
67  miit.Set(iit.Value());
68  else
69  miit.Set(0);
70 
71  ++miit;
72  ++mit;
73  ++iit;
74  }
75 
76 
77  vigra::Shape3 shape(xdim,ydim,zdim);
78  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > input_image_view(
79  shape ,
80  maske_input->GetBufferPointer());
81 
82 
83  vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > hessian_image(shape);
84  vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image(shape);
85 
86  for(unsigned int i = 0 ; i < zdim; ++i )
87  {
88  vigra::Shape2 slice_shape(xdim,ydim);
89  vigra::MultiArrayView<2, InputPixelType, vigra::StridedArrayTag > image_slice(
90  slice_shape,
91  input_image_view.data()+ (i*xdim*ydim));
92 
93  vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > hessian_slice(
94  slice_shape,
95  hessian_image.data() + (i*xdim*ydim));
96 
97  // vigra::hessianMatrixOfGaussian(input_image_slice_view, hessian_image_slice_view, sigma);
98 
99  // vigra::hessianMatrixOfGaussian(image_slice,hessian_slice,sigma);
100  vigra::hessianMatrixOfGaussian(image_slice,
101  hessian_slice.bindElementChannel(0),
102  hessian_slice.bindElementChannel(1),
103  hessian_slice.bindElementChannel(2),
104  m_Sigma);
105 
106  vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image_slice_view(
107  slice_shape,
108  eigenvalues_image.data() + (i*xdim*ydim));
109 
110  vigra::tensorEigenRepresentation(hessian_slice, eigenvalues_image_slice_view);
111  }
112 
113  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev1_image = eigenvalues_image.bindElementChannel(0);
114  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev2_image = eigenvalues_image.bindElementChannel(1);
115  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev3_image = eigenvalues_image.bindElementChannel(2);
116 
117 
118  for(unsigned int x = 0 ; x < xdim; ++x)
119  for(unsigned int y = 0 ; y < ydim; ++y)
120  for(unsigned int z = 0 ; z < zdim; ++z)
121  {
122  itk::Index<3> indx = {{x,y,z}};
123  this->GetOutput(0)->operator [](indx) = ev1_image(x,y,z);
124  this->GetOutput(1)->operator [](indx) = ev2_image(x,y,z);
125  this->GetOutput(2)->operator [](indx) = ev3_image(x,y,z);
126  }
127 
128 
129 }
130 
131 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
133 {
134  this->m_ImageMask = maskimage;
135 }
136 
137 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
139 {
140  this->SetNumberOfIndexedOutputs(3);
141  this->SetNumberOfIndexedInputs(1);
142 
143  this->SetNthOutput( 0, this->MakeOutput(0) );
144  this->SetNthOutput( 1, this->MakeOutput(1) );
145  this->SetNthOutput( 2, this->MakeOutput(2) );
146 
147 }
148 
149 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
151 {
152 
153 }
154 
155 
156 
157 #endif