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
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.