Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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