17 #ifndef ITKSTRUCTURETENSOREIGENVALUEIMAGEFILTER_CPP
18 #define ITKSTRUCTURETENSOREIGENVALUEIMAGEFILTER_CPP
21 #include <itkImageRegionIterator.h>
22 #include <vigra/tensorutilities.hxx>
23 #include <vigra/convolution.hxx>
26 template<
class TInputImageType,
class TOutputImageType,
class TMaskImageType>
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();
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();
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();
46 template<
class TInputImageType,
class TOutputImageType,
class TMaskImageType>
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);
61 vigra::Shape3 shape(xdim,ydim,zdim);
62 vigra::MultiArrayView<3, InputPixelType, vigra::StridedArrayTag > input_image_view(
64 this->GetInput()->GetBufferPointer());
66 vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > structuretensor_image(shape);
67 vigra::MultiArray<3, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image(shape);
70 for(
int i = 0 ; i < zdim; ++i )
72 vigra::Shape2 slice_shape(xdim,ydim);
73 vigra::MultiArrayView<2, InputPixelType, vigra::StridedArrayTag > input_image_slice_view(
75 input_image_view.data()+ (i*xdim*ydim));
77 vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > structuretensor_image_slice_view(
79 structuretensor_image.data() + (i*xdim*ydim));
81 vigra::structureTensor(input_image_slice_view, structuretensor_image_slice_view, m_InnerScale, m_OuterScale);
83 vigra::MultiArrayView<2, vigra::TinyVector<InputPixelType, 3> > eigenvalues_image_slice_view(
85 eigenvalues_image.data() + (i*xdim*ydim));
86 vigra::tensorEigenRepresentation(structuretensor_image_slice_view, eigenvalues_image_slice_view);
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);
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)
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);
106 template<
class TInputImageType,
class TOutputImageType,
class TMaskImageType>
109 this->m_ImageMask = maskimage;
112 template<
class TInputImageType,
class TOutputImageType,
class TMaskImageType>
115 this->SetNumberOfIndexedOutputs(3);
116 this->SetNumberOfIndexedInputs(1);
118 this->SetNthOutput( 0, this->MakeOutput(0) );
119 this->SetNthOutput( 1, this->MakeOutput(1) );
120 this->SetNthOutput( 2, this->MakeOutput(2) );
124 template<
class TInputImageType,
class TOutputImageType,
class TMaskImageType>
void SetImageMask(TMaskImageType *maskimage)