Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkStructureTensorEigenvalueImageFilter.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 ITKSTRUCTURETENSOREIGENVALUEIMAGEFILTER_CPP
18 #define ITKSTRUCTURETENSOREIGENVALUEIMAGEFILTER_CPP
19 
21 #include <itkImageRegionIterator.h>
22 #include <vigra/tensorutilities.hxx>
23 #include <vigra/convolution.hxx>
24 
25 
26 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
28 {
29  Superclass::GenerateOutputInformation();
30  this->GetOutput(0)->SetDirection(this->GetInput()->GetDirection());
31  this->GetOutput(0)->SetSpacing(this->GetInput()->GetSpacing());
32  this->GetOutput(0)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
33  this->GetOutput(0)->Allocate();
34 
35  this->GetOutput(1)->SetDirection(this->GetInput()->GetDirection());
36  this->GetOutput(1)->SetSpacing(this->GetInput()->GetSpacing());
37  this->GetOutput(1)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
38  this->GetOutput(1)->Allocate();
39 
40  this->GetOutput(2)->SetDirection(this->GetInput()->GetDirection());
41  this->GetOutput(2)->SetSpacing(this->GetInput()->GetSpacing());
42  this->GetOutput(2)->SetRegions(this->GetInput()->GetLargestPossibleRegion());
43  this->GetOutput(2)->Allocate();
44 }
45 
46 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
48 {
49 
50 
51 
53 
54  typename TInputImageType::RegionType region = this->GetInput()->GetLargestPossibleRegion();
55  unsigned int xdim = region.GetSize(0);
56  unsigned int ydim = region.GetSize(1);
57  unsigned int zdim = region.GetSize(2);
58 
59 
60 
61  vigra::Shape3 shape(xdim,ydim,zdim);
62  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > input_image_view(
63  shape ,
64  this->GetInput()->GetBufferPointer());
65 
66  vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > structuretensor_image(shape);
67  vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image(shape);
68 
69 
70  for(int i = 0 ; i < zdim; ++i )
71  {
72  vigra::Shape2 slice_shape(xdim,ydim);
73  vigra::MultiArrayView<2, InputPixelType, vigra::StridedArrayTag > input_image_slice_view(
74  slice_shape,
75  input_image_view.data()+ (i*xdim*ydim));
76 
77  vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > structuretensor_image_slice_view(
78  slice_shape,
79  structuretensor_image.data() + (i*xdim*ydim));
80 
81  vigra::structureTensor(input_image_slice_view, structuretensor_image_slice_view, m_InnerScale, m_OuterScale);
82 
83  vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image_slice_view(
84  slice_shape,
85  eigenvalues_image.data() + (i*xdim*ydim));
86  vigra::tensorEigenRepresentation(structuretensor_image_slice_view, eigenvalues_image_slice_view);
87  }
88 
89  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev1_image = eigenvalues_image.bindElementChannel(0);
90  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev2_image = eigenvalues_image.bindElementChannel(1);
91  vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > ev3_image = eigenvalues_image.bindElementChannel(2);
92 
93  for(unsigned int x = 0 ; x < xdim; ++x)
94  for(unsigned int y = 0 ; y < ydim; ++y)
95  for(unsigned int z = 0 ; z < zdim; ++z)
96  {
97  itk::Index<3> indx = {{x,y,z}};
98  this->GetOutput(0)->operator [](indx) = ev1_image(x,y,z);
99  this->GetOutput(1)->operator [](indx) = ev2_image(x,y,z);
100  this->GetOutput(2)->operator [](indx) = ev3_image(x,y,z);
101  }
102 
103 
104 }
105 
106 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
108 {
109  this->m_ImageMask = maskimage;
110 }
111 
112 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
114 {
115  this->SetNumberOfIndexedOutputs(3);
116  this->SetNumberOfIndexedInputs(1);
117 
118  this->SetNthOutput( 0, this->MakeOutput(0) );
119  this->SetNthOutput( 1, this->MakeOutput(1) );
120  this->SetNthOutput( 2, this->MakeOutput(2) );
121 
122 }
123 
124 template< class TInputImageType, class TOutputImageType, class TMaskImageType>
126 {
127 
128 }
129 
130 
131 
132 #endif
unsigned short PixelType
double InputPixelType