Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkDiffusionTensorPrincipalDirectionImageFilter.txx
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 __itkDiffusionTensorPrincipalDirectionImageFilter_txx
18 #define __itkDiffusionTensorPrincipalDirectionImageFilter_txx
19 
20 #include <time.h>
21 #include <stdio.h>
22 #include <stdlib.h>
23 
24 #include "itkDiffusionTensorPrincipalDirectionImageFilter.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionConstIteratorWithIndex.h"
27 #include "itkImageRegionIterator.h"
28 #include "itkArray.h"
29 #include "vnl/vnl_vector.h"
30 #include <boost/progress.hpp>
31 
32 #include <vtkSmartPointer.h>
33 #include <vtkPolyData.h>
34 #include <vtkCellArray.h>
35 #include <vtkPoints.h>
36 #include <vtkPolyLine.h>
37 
38 #define _USE_MATH_DEFINES
39 #include <math.h>
40 
41 namespace itk {
42 
43 //#define QBALL_RECON_PI M_PI
44 
45 template< class TTensorPixelType, class TPDPixelType>
46 DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
47 TPDPixelType>
48 ::DiffusionTensorPrincipalDirectionImageFilter()
49  : m_NormalizeVectors(true)
50  , m_MaxEigenvalue(0.0)
51 {
52  this->SetNumberOfRequiredInputs( 1 );
53 }
54 
55 template< class TTensorPixelType,
56  class TPDPixelType>
57 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
58 TPDPixelType>
59 ::BeforeThreadedGenerateData()
60 {
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();
66 
67  if (m_MaskImage.IsNull())
68  {
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);
76  }
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);
84 
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);
94 }
95 
96 template< class TTensorPixelType,
97  class TPDPixelType>
98 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
99 TPDPixelType>
100 ::AfterThreadedGenerateData()
101 {
102  vtkSmartPointer<vtkCellArray> m_VtkCellArray = vtkSmartPointer<vtkCellArray>::New();
103  vtkSmartPointer<vtkPoints> m_VtkPoints = vtkSmartPointer<vtkPoints>::New();
104 
105  typename OutputImageType::Pointer directionImage = static_cast< OutputImageType* >( this->ProcessObject::GetPrimaryOutput() );
106  ImageRegionConstIterator< OutputImageType > it(directionImage, directionImage->GetLargestPossibleRegion() );
107 
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];
114 
115  while( !it.IsAtEnd() )
116  {
117  typename OutputImageType::IndexType index = it.GetIndex();
118  if (m_MaskImage->GetPixel(index)==0)
119  {
120  ++it;
121  continue;
122  }
123 
124  itk::Vector< float, 3 > pixel = directionImage->GetPixel(index);
125  DirectionType dir; dir[0] = pixel[0]; dir[1] = pixel[1]; dir[2] = pixel[2];
126 
127  if (!m_NormalizeVectors && m_MaxEigenvalue>0)
128  dir /= m_MaxEigenvalue;
129 
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 );
137 
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);
151 
152  ++it;
153  }
154 
155  vtkSmartPointer<vtkPolyData> directionsPolyData = vtkSmartPointer<vtkPolyData>::New();
156  directionsPolyData->SetPoints(m_VtkPoints);
157  directionsPolyData->SetLines(m_VtkCellArray);
158  m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData);
159 }
160 
161 template< class TTensorPixelType,
162  class TPDPixelType>
163 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
164 TPDPixelType>
165 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType )
166 {
167 
168  typedef itk::DiffusionTensor3D<TTensorPixelType> TensorType;
169  typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
170  typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
171 
172  typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
173 
174  ImageRegionIterator< OutputImageType > outIt(outputImage, outputRegionForThread);
175  InputIteratorType inIt(inputImagePointer, outputRegionForThread );
176  ImageRegionIterator< ItkUcharImgType > nIt(m_NumDirectionsImage, outputRegionForThread );
177 
178  while( !inIt.IsAtEnd() )
179  {
180  typename InputImageType::IndexType index = inIt.GetIndex();
181  if (m_MaskImage->GetPixel(index)==0)
182  {
183  ++inIt;
184  ++nIt;
185  ++outIt;
186  continue;
187  }
188 
189  typename InputImageType::PixelType b = inIt.Get();
190  TensorType tensor = b.GetDataPointer();
191 
192  typename OutputImageType::PixelType dir;
193  typename TensorType::EigenValuesArrayType eigenvalues;
194  typename TensorType::EigenVectorsMatrixType eigenvectors;
195  if(tensor.GetTrace()!=0)
196  {
197  tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors);
198 
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);
203 
204  if (!m_NormalizeVectors)
205  vec *= eigenvalues[2];
206 
207  if (eigenvalues[2]>m_MaxEigenvalue)
208  m_MaxEigenvalue = eigenvalues[2];
209 
210  dir[0] = (TPDPixelType)vec[0];
211  dir[1] = (TPDPixelType)vec[1];
212  dir[2] = (TPDPixelType)vec[2];
213 
214  outIt.Set( dir );
215  m_NumDirectionsImage->SetPixel(index, 1);
216  }
217 
218  ++outIt;
219  ++inIt;
220  ++nIt;
221  }
222 
223  std::cout << "One Thread finished extraction" << std::endl;
224 }
225 
226 template< class TTensorPixelType,
227  class TPDPixelType>
228 void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType,
229 TPDPixelType>
230 ::PrintSelf(std::ostream& os, Indent indent) const
231 {
232 }
233 
234 }
235 #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx