Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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