Medical Imaging Interaction Toolkit  2016.11.0
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,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 #ifndef ITKHESSIANMATRIXEIGENVALUEIMAGEFILTER_CPP
18 #define ITKHESSIANMATRIXEIGENVALUEIMAGEFILTER_CPP
19 
21 #include <itkImageRegionIterator.h>
22 #include <vigra/tensorutilities.hxx>
23 #include <vigra/convolution.hxx>
24 #include <mitkCLUtil.h>
25 
26 
27 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
29 {
30  Superclass::GenerateOutputInformation();
31  this->GetOutput(0)->SetDirection(this->GetInput()->GetDirection());
32  this->GetOutput(0)->SetSpacing(this->GetInput()->GetSpacing());
33  this->GetOutput(0)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
34  this->GetOutput(0)->Allocate();
35 
36  this->GetOutput(1)->SetDirection(this->GetInput()->GetDirection());
37  this->GetOutput(1)->SetSpacing(this->GetInput()->GetSpacing());
38  this->GetOutput(1)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
39  this->GetOutput(1)->Allocate();
40 
41  this->GetOutput(2)->SetDirection(this->GetInput()->GetDirection());
42  this->GetOutput(2)->SetSpacing(this->GetInput()->GetSpacing());
43  this->GetOutput(2)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
44  this->GetOutput(2)->Allocate();
45 }
46 
47 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
49 {
50 
52 
53  typename TInputImageType::RegionType region = this->GetInput()->GetLargestPossibleRegion();
54  unsigned int xdim = region.GetSize(0);
55  unsigned int ydim = region.GetSize(1);
56  unsigned int zdim = region.GetSize(2);
57 
58  typename TInputImageType::Pointer maske_input = TInputImageType::New();
59  maske_input->SetDirection(this->GetInput()->GetDirection());
60  maske_input->SetSpacing(this->GetInput()->GetSpacing());
61  maske_input->SetRegions(this->GetInput()->GetLargestPossibleRegion());
62  maske_input->Allocate();
63 
64  itk::ImageRegionConstIterator<TInputImageType> iit(this->GetInput(), this->GetInput()->GetLargestPossibleRegion());
65  itk::ImageRegionIterator<TInputImageType> miit(maske_input, maske_input->GetLargestPossibleRegion());
66  itk::ImageRegionConstIterator<TMaskImageType> mit(m_ImageMask, m_ImageMask->GetLargestPossibleRegion());
67 
68  while(!miit.IsAtEnd())
69  {
70  if(mit.Value()!=0)
71  miit.Set(iit.Value());
72  else
73  miit.Set(0);
74 
75  ++miit;
76  ++mit;
77  ++iit;
78  }
79 
80 
81  vigra::Shape3 shape(xdim,ydim,zdim);
82  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > input_image_view(
83  shape ,
84  maske_input->GetBufferPointer());
85 
86 
87  vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > hessian_image(shape);
88  vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image(shape);
89 
90  for(unsigned int i = 0 ; i < zdim; ++i )
91  {
92  vigra::Shape2 slice_shape(xdim,ydim);
93  vigra::MultiArrayView<2, InputPixelType, vigra::StridedArrayTag > image_slice(
94  slice_shape,
95  input_image_view.data()+ (i*xdim*ydim));
96 
97  vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > hessian_slice(
98  slice_shape,
99  hessian_image.data() + (i*xdim*ydim));
100 
101  // vigra::hessianMatrixOfGaussian(input_image_slice_view, hessian_image_slice_view, sigma);
102 
103  // vigra::hessianMatrixOfGaussian(image_slice,hessian_slice,sigma);
104  vigra::hessianMatrixOfGaussian(image_slice,
105  hessian_slice.bindElementChannel(0),
106  hessian_slice.bindElementChannel(1),
107  hessian_slice.bindElementChannel(2),
108  m_Sigma);
109 
110  vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image_slice_view(
111  slice_shape,
112  eigenvalues_image.data() + (i*xdim*ydim));
113 
114  vigra::tensorEigenRepresentation(hessian_slice, eigenvalues_image_slice_view);
115  }
116 
117  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev1_image = eigenvalues_image.bindElementChannel(0);
118  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev2_image = eigenvalues_image.bindElementChannel(1);
119  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev3_image = eigenvalues_image.bindElementChannel(2);
120 
121 
122  for(unsigned int x = 0 ; x < xdim; ++x)
123  for(unsigned int y = 0 ; y < ydim; ++y)
124  for(unsigned int z = 0 ; z < zdim; ++z)
125  {
126  itk::Index<3> indx = {{x,y,z}};
127  this->GetOutput(0)->operator [](indx) = ev1_image(x,y,z);
128  this->GetOutput(1)->operator [](indx) = ev2_image(x,y,z);
129  this->GetOutput(2)->operator [](indx) = ev3_image(x,y,z);
130  }
131 
132 
133 }
134 
135 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
137 {
138  this->m_ImageMask = maskimage;
139 }
140 
141 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
143 {
144  this->SetNumberOfIndexedOutputs(3);
145  this->SetNumberOfIndexedInputs(1);
146 
147  this->SetNthOutput( 0, this->MakeOutput(0) );
148  this->SetNthOutput( 1, this->MakeOutput(1) );
149  this->SetNthOutput( 2, this->MakeOutput(2) );
150 
151 }
152 
153 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
155 {
156 
157 }
158 
159 
160 
161 #endif
itk::SmartPointer< Self > Pointer
unsigned short PixelType
double InputPixelType
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.