Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkEvaluateTractogramDirectionsFilter.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 __itkEvaluateTractogramDirectionsFilter_cpp
18 #define __itkEvaluateTractogramDirectionsFilter_cpp
19 
20 #include "itkEvaluateTractogramDirectionsFilter.h"
21 #include <itkImageRegionConstIterator.h>
22 #include <itkImageRegionIterator.h>
23 #include <itkImageDuplicator.h>
24 #include <boost/progress.hpp>
25 
26 #define _USE_MATH_DEFINES
27 #include <math.h>
28 
29 namespace itk {
30 
31 template< class PixelType >
32 EvaluateTractogramDirectionsFilter< PixelType >
33 ::EvaluateTractogramDirectionsFilter():
34  m_ReferenceImageSet(NULL),
35  m_IgnoreMissingDirections(false),
36  m_Eps(0.0001),
37  m_UseInterpolation(false)
38 {
39  this->SetNumberOfOutputs(1);
40 }
41 
42 template< class PixelType >
43 itk::Vector<PixelType, 3> EvaluateTractogramDirectionsFilter< PixelType >::GetItkVector(double point[3])
44 {
45  itk::Vector<PixelType, 3> itkVector;
46  itkVector[0] = point[0];
47  itkVector[1] = point[1];
48  itkVector[2] = point[2];
49  return itkVector;
50 }
51 
52 template< class PixelType >
53 vnl_vector_fixed<PixelType, 3> EvaluateTractogramDirectionsFilter< PixelType >::GetVnlVector(double point[3])
54 {
55  vnl_vector_fixed<PixelType, 3> vnlVector;
56  vnlVector[0] = point[0];
57  vnlVector[1] = point[1];
58  vnlVector[2] = point[2];
59  return vnlVector;
60 }
61 
62 template< class PixelType >
63 vnl_vector_fixed<PixelType, 3> EvaluateTractogramDirectionsFilter< PixelType >::GetVnlVector(Vector<PixelType,3>& vector)
64 {
65  vnl_vector_fixed<PixelType, 3> vnlVector;
66  vnlVector[0] = vector[0];
67  vnlVector[1] = vector[1];
68  vnlVector[2] = vector[2];
69  return vnlVector;
70 }
71 
72 template< class PixelType >
73 itk::Point<PixelType, 3> EvaluateTractogramDirectionsFilter< PixelType >::GetItkPoint(double point[3])
74 {
75  itk::Point<PixelType, 3> itkPoint;
76  itkPoint[0] = point[0];
77  itkPoint[1] = point[1];
78  itkPoint[2] = point[2];
79  return itkPoint;
80 }
81 
82 template< class PixelType >
83 void EvaluateTractogramDirectionsFilter< PixelType >::GenerateData()
84 {
85  if (m_Tractogram.IsNull() || m_ReferenceImageSet.IsNull())
86  return;
87 
88  if (m_UseInterpolation)
89  MITK_INFO << "Using trilinear interpolation";
90  else
91  MITK_INFO << "Using nearest neighbor interpolation";
92 
93  if (m_IgnoreMissingDirections)
94  MITK_INFO << "Ignoring missing directions";
95  else
96  MITK_INFO << "Penalizing missing directions";
97 
98  // angular error image
99  typename OutputImageType::Pointer outputImage = OutputImageType::New();
100  outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() );
101  outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() );
102  outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() );
103  outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() );
104  outputImage->Allocate();
105  outputImage->FillBuffer(0.0);
106 
107  DoubleImageType::Pointer counterImage = DoubleImageType::New();
108  counterImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() );
109  counterImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() );
110  counterImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() );
111  counterImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() );
112  counterImage->Allocate();
113  counterImage->FillBuffer(0.0);
114 
115  std::vector< DoubleImageType::Pointer > directionFound;
116  for (int i=0; i<m_ReferenceImageSet->Size(); i++)
117  {
118  DoubleImageType::Pointer check = DoubleImageType::New();
119  check->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() );
120  check->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() );
121  check->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() );
122  check->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() );
123  check->Allocate();
124  check->FillBuffer(90);
125  directionFound.push_back(check);
126  }
127 
128  if (m_MaskImage.IsNull())
129  {
130  m_MaskImage = UCharImageType::New();
131  m_MaskImage->SetOrigin( outputImage->GetOrigin() );
132  m_MaskImage->SetRegions( outputImage->GetLargestPossibleRegion() );
133  m_MaskImage->SetSpacing( outputImage->GetSpacing() );
134  m_MaskImage->SetDirection( outputImage->GetDirection() );
135  m_MaskImage->Allocate();
136  m_MaskImage->FillBuffer(1);
137  }
138 
139  m_MeanAngularError = 0.0;
140  m_MedianAngularError = 0;
141  m_MaxAngularError = 0.0;
142  m_MinAngularError = itk::NumericTraits<float>::max();
143  m_VarAngularError = 0.0;
144  m_AngularErrorVector.clear();
145 
146  float minSpacing = 1;
147  if(outputImage->GetSpacing()[0]<outputImage->GetSpacing()[1] && outputImage->GetSpacing()[0]<outputImage->GetSpacing()[2])
148  minSpacing = outputImage->GetSpacing()[0];
149  else if (outputImage->GetSpacing()[1] < outputImage->GetSpacing()[2])
150  minSpacing = outputImage->GetSpacing()[1];
151  else
152  minSpacing = outputImage->GetSpacing()[2];
153  FiberBundleType::Pointer fiberBundle = m_Tractogram->GetDeepCopy();
154  fiberBundle->ResampleFibers(minSpacing/10);
155 
156  MITK_INFO << "Evaluating tractogram";
157  vtkSmartPointer<vtkPolyData> fiberPolyData = fiberBundle->GetFiberPolyData();
158  boost::progress_display disp( fiberBundle->GetNumFibers() );
159  for( int i=0; i<fiberBundle->GetNumFibers(); i++ )
160  {
161  vtkCell* cell = fiberPolyData->GetCell(i);
162  int numPoints = cell->GetNumberOfPoints();
163  vtkPoints* points = cell->GetPoints();
164 
165  if (numPoints<2)
166  continue;
167 
168  for( int j=0; j<numPoints; j++)
169  {
170  double* temp = points->GetPoint(j);
171  itk::Point<PixelType, 3> vertex = GetItkPoint(temp);
172  itk::Vector<PixelType> v = GetItkVector(temp);
173 
174  itk::Vector<PixelType, 3> dir(3);
175  if (j<numPoints-1)
176  dir = GetItkVector(points->GetPoint(j+1))-v;
177  else
178  dir = v-GetItkVector(points->GetPoint(j-1));
179 
180  vnl_vector_fixed< PixelType, 3 > fiberDir = GetVnlVector(dir);
181  fiberDir.normalize();
182 
183  itk::Index<3> idx;
184  itk::ContinuousIndex<float, 3> contIndex;
185  m_MaskImage->TransformPhysicalPointToIndex(vertex, idx);
186  m_MaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
187 
188  if (!m_UseInterpolation) // use nearest neighbour interpolation
189  {
190  if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)<=0)
191  continue;
192 
193  double angle = 90;
194  int usedIndex = -1;
195  for (int k=0; k<m_ReferenceImageSet->Size(); k++) // and each test direction
196  {
197  vnl_vector_fixed< PixelType, 3 > refDir = m_ReferenceImageSet->GetElement(k)->GetPixel(idx).GetVnlVector();
198  if (refDir.magnitude()>m_Eps) // normalize if not null
199  refDir.normalize();
200  else
201  continue;
202 
203  // calculate angle between directions
204  double tempAngle = acos(fabs(dot_product(refDir, fiberDir)))*180.0/M_PI;
205  directionFound.at(k)->SetPixel(idx, tempAngle);
206 
207  if (tempAngle < angle)
208  {
209  angle = tempAngle;
210  usedIndex = k;
211  }
212  }
213  if (usedIndex>=0)
214  directionFound.at(usedIndex)->SetPixel(idx, -1);
215  else if (m_IgnoreMissingDirections)
216  angle = 0;
217 
218  counterImage->SetPixel(idx, counterImage->GetPixel(idx)+1);
219  outputImage->SetPixel(idx, outputImage->GetPixel(idx)+angle);
220  continue;
221  }
222 
223  double frac_x = contIndex[0] - idx[0];
224  double frac_y = contIndex[1] - idx[1];
225  double frac_z = contIndex[2] - idx[2];
226  if (frac_x<0)
227  {
228  idx[0] -= 1;
229  frac_x += 1;
230  }
231  if (frac_y<0)
232  {
233  idx[1] -= 1;
234  frac_y += 1;
235  }
236  if (frac_z<0)
237  {
238  idx[2] -= 1;
239  frac_z += 1;
240  }
241 
242  // use trilinear interpolation
243  itk::Index<3> newIdx;
244  for (int x=0; x<2; x++)
245  {
246  frac_x = 1-frac_x;
247  for (int y=0; y<2; y++)
248  {
249  frac_y = 1-frac_y;
250  for (int z=0; z<2; z++)
251  {
252  frac_z = 1-frac_z;
253 
254  newIdx[0] = idx[0]+x;
255  newIdx[1] = idx[1]+y;
256  newIdx[2] = idx[2]+z;
257 
258  double frac = frac_x*frac_y*frac_z;
259 
260  // is position valid?
261  if (!m_MaskImage->GetLargestPossibleRegion().IsInside(newIdx) || m_MaskImage->GetPixel(newIdx)<=0)
262  continue;
263 
264  double angle = 90;
265  int usedIndex = -1;
266  for (int k=0; k<m_ReferenceImageSet->Size(); k++) // and each test direction
267  {
268  vnl_vector_fixed< PixelType, 3 > refDir = m_ReferenceImageSet->GetElement(k)->GetPixel(newIdx).GetVnlVector();
269  if (refDir.magnitude()>m_Eps) // normalize if not null
270  refDir.normalize();
271  else
272  continue;
273 
274  // calculate angle between directions
275  double tempAngle = acos(fabs(dot_product(refDir, fiberDir)))*180.0/M_PI;
276  directionFound.at(k)->SetPixel(newIdx, tempAngle);
277 
278  if (tempAngle < angle)
279  {
280  angle = tempAngle;
281  usedIndex = k;
282  }
283  }
284  if (usedIndex>=0)
285  directionFound.at(usedIndex)->SetPixel(newIdx, -1);
286  else if (m_IgnoreMissingDirections)
287  angle = 0;
288 
289  counterImage->SetPixel(newIdx, counterImage->GetPixel(newIdx)+1);
290  outputImage->SetPixel(newIdx, outputImage->GetPixel(newIdx)+frac*angle);
291  }
292  }
293  }
294  }
295  ++disp;
296  }
297 
298  ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion());
299  ImageRegionIterator< UCharImageType > mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion());
300  ImageRegionIterator< DoubleImageType > cit(counterImage, counterImage->GetLargestPossibleRegion());
301 
302 
303  if (!m_IgnoreMissingDirections)
304  MITK_INFO << "Creatings statistics and accounting for missing directions";
305  else
306  MITK_INFO << "Creatings statistics";
307  boost::progress_display disp2(outputImage->GetLargestPossibleRegion().GetNumberOfPixels());
308  while( !oit.IsAtEnd() )
309  {
310  if ( mit.Get()>0)
311  {
312  if ( cit.Get()>m_Eps )
313  oit.Set(oit.Get()/cit.Get());
314 
315  if (!m_IgnoreMissingDirections)
316  {
317  int missingCount = 0;
318 
319  if ( cit.Get()>m_Eps )
320  missingCount++;
321 
322  for (int i=0; i<directionFound.size(); i++)
323  {
324  if ( directionFound.at(i)->GetPixel(oit.GetIndex())>0 && m_ReferenceImageSet->GetElement(i)->GetPixel(oit.GetIndex()).GetVnlVector().magnitude()>m_Eps )
325  {
326  oit.Set(oit.Get()+directionFound.at(i)->GetPixel(oit.GetIndex()));
327  missingCount++;
328  }
329  }
330  if (missingCount>1)
331  oit.Set(oit.Get()/missingCount);
332  }
333 
334  if (oit.Get()>m_MaxAngularError)
335  m_MaxAngularError = oit.Get();
336  if (oit.Get()<m_MinAngularError)
337  m_MinAngularError = oit.Get();
338  }
339 
340  ++oit;
341  ++cit;
342  ++mit;
343  ++disp2;
344  }
345 
346  this->SetNthOutput(0, outputImage);
347 }
348 
349 }
350 
351 #endif // __itkEvaluateTractogramDirectionsFilter_cpp
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
itk::EvaluateTractogramDirectionsFilter::GetItkVector
itk::Vector< PixelType, 3 > GetItkVector(double point[3])
Definition: itkEvaluateTractogramDirectionsFilter.cpp:43
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
itk::EvaluateTractogramDirectionsFilter::GetVnlVector
vnl_vector_fixed< PixelType, 3 > GetVnlVector(double point[3])
Definition: itkEvaluateTractogramDirectionsFilter.cpp:53
max
static T max(T x, T y)
Definition: svm.cpp:70
itk::EvaluateTractogramDirectionsFilter::GetItkPoint
itk::Point< PixelType, 3 > GetItkPoint(double point[3])
Definition: itkEvaluateTractogramDirectionsFilter.cpp:73
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.