Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkTractDensityImageFilter.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Coindex[1]right (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 #include "itkTractDensityImageFilter.h"
17 
18 // VTK
19 #include <vtkPolyLine.h>
20 #include <vtkCellArray.h>
21 #include <vtkCellData.h>
22 #include <vtkCell.h>
23 
24 // misc
25 #include <math.h>
26 #include <boost/progress.hpp>
27 
28 namespace itk{
29 
30 template< class OutputImageType >
31 TractDensityImageFilter< OutputImageType >::TractDensityImageFilter()
32  : m_InvertImage(false)
33  , m_FiberBundle(NULL)
34  , m_UpsamplingFactor(1)
35  , m_InputImage(NULL)
36  , m_BinaryOutput(false)
37  , m_UseImageGeometry(false)
38  , m_OutputAbsoluteValues(false)
39  , m_UseTrilinearInterpolation(false)
40  , m_DoFiberResampling(true)
41 {
42 
43 }
44 
45 template< class OutputImageType >
46 TractDensityImageFilter< OutputImageType >::~TractDensityImageFilter()
47 {
48 }
49 
50 template< class OutputImageType >
51 itk::Point<float, 3> TractDensityImageFilter< OutputImageType >::GetItkPoint(double point[3])
52 {
53  itk::Point<float, 3> itkPoint;
54  itkPoint[0] = point[0];
55  itkPoint[1] = point[1];
56  itkPoint[2] = point[2];
57  return itkPoint;
58 }
59 
60 template< class OutputImageType >
61 void TractDensityImageFilter< OutputImageType >::GenerateData()
62 {
63  // generate upsampled image
64  mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry();
65  typename OutputImageType::Pointer outImage = this->GetOutput();
66 
67  // calculate new image parameters
68  itk::Vector<double,3> newSpacing;
69  mitk::Point3D newOrigin;
70  itk::Matrix<double, 3, 3> newDirection;
71  ImageRegion<3> upsampledRegion;
72  if (m_UseImageGeometry && !m_InputImage.IsNull())
73  {
74  MITK_INFO << "TractDensityImageFilter: using image geometry";
75  newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor;
76  upsampledRegion = m_InputImage->GetLargestPossibleRegion();
77  newOrigin = m_InputImage->GetOrigin();
78  typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize();
79  size[0] *= m_UpsamplingFactor;
80  size[1] *= m_UpsamplingFactor;
81  size[2] *= m_UpsamplingFactor;
82  upsampledRegion.SetSize(size);
83  newDirection = m_InputImage->GetDirection();
84  }
85  else
86  {
87  MITK_INFO << "TractDensityImageFilter: using fiber bundle geometry";
88  newSpacing = geometry->GetSpacing()/m_UpsamplingFactor;
89  newOrigin = geometry->GetOrigin();
90  mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
91 
92  // we retrieve the origin from a vtk-polydata (corner-based) and hance have to translate it to an image geometry
93  // i.e. center-based
94  newOrigin[0] += bounds.GetElement(0) + 0.5 * newSpacing[0];
95  newOrigin[1] += bounds.GetElement(2) + 0.5 * newSpacing[1];
96  newOrigin[2] += bounds.GetElement(4) + 0.5 * newSpacing[2];
97 
98  for (int i=0; i<3; i++)
99  for (int j=0; j<3; j++)
100  newDirection[j][i] = geometry->GetMatrixColumn(i)[j];
101  upsampledRegion.SetSize(0, ceil( geometry->GetExtent(0)*m_UpsamplingFactor ) );
102  upsampledRegion.SetSize(1, ceil( geometry->GetExtent(1)*m_UpsamplingFactor ) );
103  upsampledRegion.SetSize(2, ceil( geometry->GetExtent(2)*m_UpsamplingFactor ) );
104  }
105  typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize();
106 
107  // apply new image parameters
108  outImage->SetSpacing( newSpacing );
109  outImage->SetOrigin( newOrigin );
110  outImage->SetDirection( newDirection );
111  outImage->SetLargestPossibleRegion( upsampledRegion );
112  outImage->SetBufferedRegion( upsampledRegion );
113  outImage->SetRequestedRegion( upsampledRegion );
114  outImage->Allocate();
115  outImage->FillBuffer(0.0);
116 
117  int w = upsampledSize[0];
118  int h = upsampledSize[1];
119  int d = upsampledSize[2];
120 
121  // set/initialize output
122  OutPixelType* outImageBufferPointer = (OutPixelType*)outImage->GetBufferPointer();
123 
124  // resample fiber bundle
125  float minSpacing = 1;
126  if(newSpacing[0]<newSpacing[1] && newSpacing[0]<newSpacing[2])
127  minSpacing = newSpacing[0];
128  else if (newSpacing[1] < newSpacing[2])
129  minSpacing = newSpacing[1];
130  else
131  minSpacing = newSpacing[2];
132 
133  MITK_INFO << "TractDensityImageFilter: resampling fibers to ensure sufficient voxel coverage";
134  if (m_DoFiberResampling)
135  {
136  m_FiberBundle = m_FiberBundle->GetDeepCopy();
137  m_FiberBundle->ResampleSpline(minSpacing/10);
138  }
139 
140  MITK_INFO << "TractDensityImageFilter: starting image generation";
141 
142  vtkSmartPointer<vtkPolyData> fiberPolyData = m_FiberBundle->GetFiberPolyData();
143  vtkSmartPointer<vtkCellArray> vLines = fiberPolyData->GetLines();
144  vLines->InitTraversal();
145  int numFibers = m_FiberBundle->GetNumFibers();
146  boost::progress_display disp(numFibers);
147  for( int i=0; i<numFibers; i++ )
148  {
149  ++disp;
150  vtkIdType numPoints(0);
151  vtkIdType* points(NULL);
152  vLines->GetNextCell ( numPoints, points );
153  float weight = m_FiberBundle->GetFiberWeight(i);
154 
155  // fill output image
156  for( int j=0; j<numPoints; j++)
157  {
158  itk::Point<float, 3> vertex = GetItkPoint(fiberPolyData->GetPoint(points[j]));
159  itk::Index<3> index;
160  itk::ContinuousIndex<float, 3> contIndex;
161  outImage->TransformPhysicalPointToIndex(vertex, index);
162  outImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
163 
164  if (!m_UseTrilinearInterpolation && outImage->GetLargestPossibleRegion().IsInside(index))
165  {
166  if (m_BinaryOutput)
167  outImage->SetPixel(index, 1);
168  else
169  outImage->SetPixel(index, outImage->GetPixel(index)+0.01*weight);
170  continue;
171  }
172 
173  float frac_x = contIndex[0] - index[0];
174  float frac_y = contIndex[1] - index[1];
175  float frac_z = contIndex[2] - index[2];
176 
177  if (frac_x<0)
178  {
179  index[0] -= 1;
180  frac_x += 1;
181  }
182  if (frac_y<0)
183  {
184  index[1] -= 1;
185  frac_y += 1;
186  }
187  if (frac_z<0)
188  {
189  index[2] -= 1;
190  frac_z += 1;
191  }
192 
193  frac_x = 1-frac_x;
194  frac_y = 1-frac_y;
195  frac_z = 1-frac_z;
196 
197  // int coordinates inside image?
198  if (index[0] < 0 || index[0] >= w-1)
199  continue;
200  if (index[1] < 0 || index[1] >= h-1)
201  continue;
202  if (index[2] < 0 || index[2] >= d-1)
203  continue;
204 
205  if (m_BinaryOutput)
206  {
207  outImageBufferPointer[( index[0] + w*(index[1] + h*index[2] ))] = 1;
208  outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2] ))] = 1;
209  outImageBufferPointer[( index[0] + w*(index[1] + h*index[2]+h))] = 1;
210  outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2]+h))] = 1;
211  outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2] ))] = 1;
212  outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2]+h))] = 1;
213  outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2] ))] = 1;
214  outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2]+h))] = 1;
215  }
216  else
217  {
218  outImageBufferPointer[( index[0] + w*(index[1] + h*index[2] ))] += ( frac_x)*( frac_y)*( frac_z);
219  outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2] ))] += ( frac_x)*(1-frac_y)*( frac_z);
220  outImageBufferPointer[( index[0] + w*(index[1] + h*index[2]+h))] += ( frac_x)*( frac_y)*(1-frac_z);
221  outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2]+h))] += ( frac_x)*(1-frac_y)*(1-frac_z);
222  outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2] ))] += (1-frac_x)*( frac_y)*( frac_z);
223  outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2]+h))] += (1-frac_x)*( frac_y)*(1-frac_z);
224  outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2] ))] += (1-frac_x)*(1-frac_y)*( frac_z);
225  outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2]+h))] += (1-frac_x)*(1-frac_y)*(1-frac_z);
226  }
227  }
228  }
229 
230  if (!m_OutputAbsoluteValues && !m_BinaryOutput)
231  {
232  MITK_INFO << "TractDensityImageFilter: max-normalizing output image";
233  OutPixelType max = 0;
234  for (int i=0; i<w*h*d; i++)
235  if (max < outImageBufferPointer[i])
236  max = outImageBufferPointer[i];
237  if (max>0)
238  for (int i=0; i<w*h*d; i++)
239  {
240  outImageBufferPointer[i] /= max;
241  }
242  }
243  if (m_InvertImage)
244  {
245  MITK_INFO << "TractDensityImageFilter: inverting image";
246  for (int i=0; i<w*h*d; i++)
247  outImageBufferPointer[i] = 1-outImageBufferPointer[i];
248  }
249  MITK_INFO << "TractDensityImageFilter: finished processing";
250 }
251 }
itk::TractDensityImageFilter::OutPixelType
OutputImageType::PixelType OutPixelType
Definition: itkTractDensityImageFilter.h:40
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
mitk::BaseGeometry::BoundsArrayType
BoundingBoxType::BoundsArrayType BoundsArrayType
Definition: mitkBaseGeometry.h:111
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
max
static T max(T x, T y)
Definition: svm.cpp:70
itk::TractDensityImageFilter::GetItkPoint
itk::Point< float, 3 > GetItkPoint(double point[3])
Definition: itkTractDensityImageFilter.cpp:51