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