1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 #ifndef __itkDiffusionTensorPrincipalDirectionImageFilter_txx
18 #define __itkDiffusionTensorPrincipalDirectionImageFilter_txx
24 #include "itkDiffusionTensorPrincipalDirectionImageFilter.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionConstIteratorWithIndex.h"
27 #include "itkImageRegionIterator.h"
29 #include "vnl/vnl_vector.h"
30 #include <boost/progress.hpp>
32 #include <vtkSmartPointer.h>
33 #include <vtkPolyData.h>
34 #include <vtkCellArray.h>
35 #include <vtkPoints.h>
36 #include <vtkPolyLine.h>
38 #define _USE_MATH_DEFINES
43 //#define QBALL_RECON_PI M_PI
45 template< class TTensorPixelType, class TPDPixelType>
46 DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
48 ::DiffusionTensorPrincipalDirectionImageFilter()
49 : m_NormalizeVectors(true)
50 , m_MaxEigenvalue(0.0)
52 this->SetNumberOfRequiredInputs( 1 );
55 template< class TTensorPixelType,
57 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
59 ::BeforeThreadedGenerateData()
61 typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
62 Vector<double,3> spacing = inputImagePointer->GetSpacing();
63 mitk::Point3D origin = inputImagePointer->GetOrigin();
64 itk::Matrix<double, 3, 3> direction = inputImagePointer->GetDirection();
65 ImageRegion<3> imageRegion = inputImagePointer->GetLargestPossibleRegion();
67 if (m_MaskImage.IsNull())
69 m_MaskImage = ItkUcharImgType::New();
70 m_MaskImage->SetSpacing( spacing );
71 m_MaskImage->SetOrigin( origin );
72 m_MaskImage->SetDirection( direction );
73 m_MaskImage->SetRegions( imageRegion );
74 m_MaskImage->Allocate();
75 m_MaskImage->FillBuffer(1);
77 m_NumDirectionsImage = ItkUcharImgType::New();
78 m_NumDirectionsImage->SetSpacing( spacing );
79 m_NumDirectionsImage->SetOrigin( origin );
80 m_NumDirectionsImage->SetDirection( direction );
81 m_NumDirectionsImage->SetRegions( imageRegion );
82 m_NumDirectionsImage->Allocate();
83 m_NumDirectionsImage->FillBuffer(0);
85 itk::Vector< TPDPixelType, 3 > nullVec; nullVec.Fill(0.0);
86 typename OutputImageType::Pointer outputImage = OutputImageType::New();
87 outputImage->SetSpacing( spacing );
88 outputImage->SetOrigin( origin );
89 outputImage->SetDirection( direction );
90 outputImage->SetRegions( imageRegion );
91 outputImage->Allocate();
92 outputImage->FillBuffer(nullVec);
93 this->SetNthOutput(0, outputImage);
96 template< class TTensorPixelType,
98 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
100 ::AfterThreadedGenerateData()
102 vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
103 vtkSmartPointer<vtkPoints> m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
105 typename OutputImageType::Pointer directionImage = static_cast< OutputImageType* >( this->ProcessObject::GetPrimaryOutput() );
106 ImageRegionConstIterator< OutputImageType > it(directionImage, directionImage->GetLargestPossibleRegion() );
108 mitk::Vector3D spacing = directionImage->GetSpacing();
109 double minSpacing = spacing[0];
110 if (spacing[1]<minSpacing)
111 minSpacing = spacing[1];
112 if (spacing[2]<minSpacing)
113 minSpacing = spacing[2];
115 while( !it.IsAtEnd() )
117 typename OutputImageType::IndexType index = it.GetIndex();
118 if (m_MaskImage->GetPixel(index)==0)
124 itk::Vector< float, 3 > pixel = directionImage->GetPixel(index);
125 DirectionType dir; dir[0] = pixel[0]; dir[1] = pixel[1]; dir[2] = pixel[2];
127 if (!m_NormalizeVectors && m_MaxEigenvalue>0)
128 dir /= m_MaxEigenvalue;
130 vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
131 itk::ContinuousIndex<double, 3> center;
132 center[0] = index[0];
133 center[1] = index[1];
134 center[2] = index[2];
135 itk::Point<double> worldCenter;
136 directionImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter );
138 itk::Point<double> worldStart;
139 worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing;
140 worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing;
141 worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing;
142 vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer());
143 container->GetPointIds()->InsertNextId(id);
144 itk::Point<double> worldEnd;
145 worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing;
146 worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing;
147 worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing;
148 id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer());
149 container->GetPointIds()->InsertNextId(id);
150 m_VtkCellArray->InsertNextCell(container);
155 vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
156 directionsPolyData->SetPoints(m_VtkPoints);
157 directionsPolyData->SetLines(m_VtkCellArray);
158 m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData);
161 template< class TTensorPixelType,
163 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
165 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType )
168 typedef itk::DiffusionTensor3D<TTensorPixelType> TensorType;
169 typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
170 typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
172 typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
174 ImageRegionIterator< OutputImageType > outIt(outputImage, outputRegionForThread);
175 InputIteratorType inIt(inputImagePointer, outputRegionForThread );
176 ImageRegionIterator< ItkUcharImgType > nIt(m_NumDirectionsImage, outputRegionForThread );
178 while( !inIt.IsAtEnd() )
180 typename InputImageType::IndexType index = inIt.GetIndex();
181 if (m_MaskImage->GetPixel(index)==0)
189 typename InputImageType::PixelType b = inIt.Get();
190 TensorType tensor = b.GetDataPointer();
192 typename OutputImageType::PixelType dir;
193 typename TensorType::EigenValuesArrayType eigenvalues;
194 typename TensorType::EigenVectorsMatrixType eigenvectors;
195 if(tensor.GetTrace()!=0)
197 tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors);
199 vnl_vector_fixed<double,3> vec;
200 vec[0] = eigenvectors(2,0);
201 vec[1] = eigenvectors(2,1);
202 vec[2] = eigenvectors(2,2);
204 if (!m_NormalizeVectors)
205 vec *= eigenvalues[2];
207 if (eigenvalues[2]>m_MaxEigenvalue)
208 m_MaxEigenvalue = eigenvalues[2];
210 dir[0] = (TPDPixelType)vec[0];
211 dir[1] = (TPDPixelType)vec[1];
212 dir[2] = (TPDPixelType)vec[2];
215 m_NumDirectionsImage->SetPixel(index, 1);
223 std::cout << "One Thread finished extraction" << std::endl;
226 template< class TTensorPixelType,
228 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
230 ::PrintSelf(std::ostream& os, Indent indent) const
235 #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx