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
itkTractsToRgbaImageFilter.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 #include "itkTractsToRgbaImageFilter.h"
17 
18 // VTK
19 #include <vtkPolyLine.h>
20 #include <vtkCellArray.h>
21 #include <vtkCellData.h>
22 
23 // misc
24 #include <math.h>
25 #include <boost/progress.hpp>
26 
27 namespace itk{
28 
29  template< class OutputImageType >
30  TractsToRgbaImageFilter< OutputImageType >::TractsToRgbaImageFilter()
31  : m_UpsamplingFactor(1)
32  , m_InputImage(NULL)
33  , m_UseImageGeometry(false)
34  {
35 
36  }
37 
38  template< class OutputImageType >
39  TractsToRgbaImageFilter< OutputImageType >::~TractsToRgbaImageFilter()
40  {
41  }
42 
43  template< class OutputImageType >
44  itk::Point<float, 3> TractsToRgbaImageFilter< OutputImageType >::GetItkPoint(double point[3])
45  {
46  itk::Point<float, 3> itkPoint;
47  itkPoint[0] = point[0];
48  itkPoint[1] = point[1];
49  itkPoint[2] = point[2];
50  return itkPoint;
51  }
52 
53  template< class OutputImageType >
54  void TractsToRgbaImageFilter< OutputImageType >::GenerateData()
55  {
56  if(&typeid(OutPixelType) != &typeid(itk::RGBAPixel<unsigned char>))
57  return;
58 
59  // generate upsampled image
60  mitk::BaseGeometry::Pointer geometry = m_FiberBundle->GetGeometry();
61  typename OutputImageType::Pointer outImage = this->GetOutput();
62 
63  // calculate new image parameters
64  itk::Vector<double,3> newSpacing;
65  mitk::Point3D newOrigin;
66  itk::Matrix<double, 3, 3> newDirection;
67  ImageRegion<3> upsampledRegion;
68  if (m_UseImageGeometry && !m_InputImage.IsNull())
69  {
70  newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor;
71  upsampledRegion = m_InputImage->GetLargestPossibleRegion();
72  newOrigin = m_InputImage->GetOrigin();
73  typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize();
74  size[0] *= m_UpsamplingFactor;
75  size[1] *= m_UpsamplingFactor;
76  size[2] *= m_UpsamplingFactor;
77  upsampledRegion.SetSize(size);
78  newDirection = m_InputImage->GetDirection();
79  }
80  else
81  {
82  newSpacing = geometry->GetSpacing()/m_UpsamplingFactor;
83  newOrigin = geometry->GetOrigin();
84  mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds();
85  newOrigin[0] += bounds.GetElement(0);
86  newOrigin[1] += bounds.GetElement(2);
87  newOrigin[2] += bounds.GetElement(4);
88 
89  for (int i=0; i<3; i++)
90  for (int j=0; j<3; j++)
91  newDirection[j][i] = geometry->GetMatrixColumn(i)[j];
92  upsampledRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor);
93  upsampledRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor);
94  upsampledRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor);
95  }
96  typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize();
97 
98  // apply new image parameters
99  outImage->SetSpacing( newSpacing );
100  outImage->SetOrigin( newOrigin );
101  outImage->SetDirection( newDirection );
102  outImage->SetRegions( upsampledRegion );
103  outImage->Allocate();
104 
105  int w = upsampledSize[0];
106  int h = upsampledSize[1];
107  int d = upsampledSize[2];
108 
109  // set/initialize output
110  unsigned char* outImageBufferPointer = (unsigned char*)outImage->GetBufferPointer();
111  float* buffer = new float[w*h*d*4];
112  for (int i=0; i<w*h*d*4; i++)
113  buffer[i] = 0;
114 
115  // resample fiber bundle
116  float minSpacing = 1;
117  if(newSpacing[0]<newSpacing[1] && newSpacing[0]<newSpacing[2])
118  minSpacing = newSpacing[0];
119  else if (newSpacing[1] < newSpacing[2])
120  minSpacing = newSpacing[1];
121  else
122  minSpacing = newSpacing[2];
123 
124  m_FiberBundle = m_FiberBundle->GetDeepCopy();
125  m_FiberBundle->ResampleSpline(minSpacing);
126 
127  vtkSmartPointer<vtkPolyData> fiberPolyData = m_FiberBundle->GetFiberPolyData();
128  vtkSmartPointer<vtkCellArray> vLines = fiberPolyData->GetLines();
129  vLines->InitTraversal();
130 
131  int numFibers = m_FiberBundle->GetNumFibers();
132  boost::progress_display disp(numFibers);
133  for( int i=0; i<numFibers; i++ )
134  {
135  ++disp;
136  vtkIdType numPoints(0);
137  vtkIdType* points(NULL);
138  vLines->GetNextCell ( numPoints, points );
139 
140  // calc directions (which are used as weights)
141  std::list< itk::Point<float, 3> > rgbweights;
142  std::list<float> intensities;
143 
144  for( int j=0; j<numPoints-1; j++)
145  {
146 
147  itk::Point<float, 3> vertex = GetItkPoint(fiberPolyData->GetPoint(points[j]));
148  itk::Point<float, 3> vertexPost = GetItkPoint(fiberPolyData->GetPoint(points[j+1]));
149 
150  itk::Point<float, 3> dir;
151  dir[0] = fabs((vertexPost[0] - vertex[0]) * outImage->GetSpacing()[0]);
152  dir[1] = fabs((vertexPost[1] - vertex[1]) * outImage->GetSpacing()[1]);
153  dir[2] = fabs((vertexPost[2] - vertex[2]) * outImage->GetSpacing()[2]);
154 
155  rgbweights.push_back(dir);
156 
157  float intensity = sqrt(dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]);
158  intensities.push_back(intensity);
159 
160  // last point gets same as previous one
161  if(j==numPoints-2)
162  {
163  rgbweights.push_back(dir);
164  intensities.push_back(intensity);
165  }
166  }
167 
168  // fill output image
169  for( int j=0; j<numPoints; j++)
170  {
171  itk::Point<float, 3> vertex = GetItkPoint(fiberPolyData->GetPoint(points[j]));
172  itk::Index<3> index;
173  itk::ContinuousIndex<float, 3> contIndex;
174  outImage->TransformPhysicalPointToIndex(vertex, index);
175  outImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex);
176 
177  float frac_x = contIndex[0] - index[0];
178  float frac_y = contIndex[1] - index[1];
179  float frac_z = contIndex[2] - index[2];
180 
181  int px = index[0];
182  if (frac_x<0)
183  {
184  px -= 1;
185  frac_x += 1;
186  }
187 
188  int py = index[1];
189  if (frac_y<0)
190  {
191  py -= 1;
192  frac_y += 1;
193  }
194 
195  int pz = index[2];
196  if (frac_z<0)
197  {
198  pz -= 1;
199  frac_z += 1;
200  }
201 
202  // int coordinates inside image?
203  if (px < 0 || px >= w-1)
204  continue;
205  if (py < 0 || py >= h-1)
206  continue;
207  if (pz < 0 || pz >= d-1)
208  continue;
209 
210  float scale = 100 * pow((float)m_UpsamplingFactor,3);
211  itk::Point<float, 3> rgbweight = rgbweights.front();
212  rgbweights.pop_front();
213  float intweight = intensities.front();
214  intensities.pop_front();
215 
216  // add to r-channel in output image
217  buffer[0+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[0] * scale;
218  buffer[0+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * rgbweight[0] * scale;
219  buffer[0+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * rgbweight[0] * scale;
220  buffer[0+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * rgbweight[0] * scale;
221  buffer[0+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[0] * scale;
222  buffer[0+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * rgbweight[0] * scale;
223  buffer[0+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * rgbweight[0] * scale;
224  buffer[0+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * rgbweight[0] * scale;
225 
226  // add to g-channel in output image
227  buffer[1+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[1] * scale;
228  buffer[1+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * rgbweight[1] * scale;
229  buffer[1+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * rgbweight[1] * scale;
230  buffer[1+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * rgbweight[1] * scale;
231  buffer[1+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[1] * scale;
232  buffer[1+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * rgbweight[1] * scale;
233  buffer[1+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * rgbweight[1] * scale;
234  buffer[1+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * rgbweight[1] * scale;
235 
236  // add to b-channel in output image
237  buffer[2+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[2] * scale;
238  buffer[2+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * rgbweight[2] * scale;
239  buffer[2+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * rgbweight[2] * scale;
240  buffer[2+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * rgbweight[2] * scale;
241  buffer[2+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[2] * scale;
242  buffer[2+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * rgbweight[2] * scale;
243  buffer[2+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * rgbweight[2] * scale;
244  buffer[2+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * rgbweight[2] * scale;
245 
246  // add to a-channel in output image
247  buffer[3+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * intweight * scale;
248  buffer[3+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * intweight * scale;
249  buffer[3+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * intweight * scale;
250  buffer[3+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * intweight * scale;
251  buffer[3+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * intweight * scale;
252  buffer[3+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * intweight * scale;
253  buffer[3+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * intweight * scale;
254  buffer[3+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * intweight * scale;
255  }
256  }
257  float maxRgb = 0.000000001;
258  float maxInt = 0.000000001;
259  int numPix;
260 
261  numPix = w*h*d*4;
262  // calc maxima
263  for(int i=0; i<numPix; i++)
264  {
265  if((i-3)%4 != 0)
266  {
267  if(buffer[i] > maxRgb)
268  maxRgb = buffer[i];
269  }
270  else
271  {
272  if(buffer[i] > maxInt)
273  maxInt = buffer[i];
274  }
275  }
276 
277  // write output, normalized uchar 0..255
278  for(int i=0; i<numPix; i++)
279  {
280  if((i-3)%4 != 0)
281  outImageBufferPointer[i] = (unsigned char) (255.0 * buffer[i] / maxRgb);
282  else
283  outImageBufferPointer[i] = (unsigned char) (255.0 * buffer[i] / maxInt);
284  }
285  }
286 }
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
itk::TractsToRgbaImageFilter::OutPixelType
OutputImageType::PixelType OutPixelType
Definition: itkTractsToRgbaImageFilter.h:40
mitk::BaseGeometry::BoundsArrayType
BoundingBoxType::BoundsArrayType BoundsArrayType
Definition: mitkBaseGeometry.h:111
itk::TractsToRgbaImageFilter::GetItkPoint
itk::Point< float, 3 > GetItkPoint(double point[3])
Definition: itkTractsToRgbaImageFilter.cpp:44