Medical Imaging Interaction Toolkit  2018.4.99-c0f884b2
Medical Imaging Interaction Toolkit
mitkIntensityProfile.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 (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #include <itkLinearInterpolateImageFunction.h>
14 #include <itkNearestNeighborInterpolateImageFunction.h>
15 #include <itkPolyLineParametricPath.h>
16 #include <itkWindowedSincInterpolateImageFunction.h>
17 #include <mitkImageAccessByItk.h>
19 #include <mitkPixelTypeMultiplex.h>
21 #include "mitkIntensityProfile.h"
22 
23 using namespace mitk;
24 
25 template <class T>
26 static void ReadPixel(const PixelType&, Image::Pointer image, const itk::Index<3>& index, ScalarType* returnValue)
27 {
28  switch (image->GetDimension())
29  {
30  case 2:
31  {
32  ImagePixelReadAccessor<T, 2> readAccess(image, image->GetSliceData(0));
33  *returnValue = readAccess.GetPixelByIndex(reinterpret_cast<const itk::Index<2>&>(index));
34  break;
35  }
36 
37  case 3:
38  {
39  ImagePixelReadAccessor<T, 3> readAccess(image, image->GetVolumeData(0));
40  *returnValue = readAccess.GetPixelByIndex(index);
41  break;
42  }
43 
44  default:
45  *returnValue = 0;
46  break;
47  }
48 }
49 
50 static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path)
51 {
52  if (image->GetDimension() == 4)
53  {
54  mitkThrow() << "computation of intensity profiles not supported for 4D images";
55  }
56 
57  IntensityProfile::Pointer intensityProfile = IntensityProfile::New();
58  itk::PolyLineParametricPath<3>::InputType input = path->StartOfInput();
59  BaseGeometry* imageGeometry = image->GetGeometry();
60  const PixelType pixelType = image->GetPixelType();
61 
62  IntensityProfile::MeasurementVectorType measurementVector;
63  itk::PolyLineParametricPath<3>::OffsetType offset;
64  Point3D worldPoint;
65  itk::Index<3> index;
66 
67  do
68  {
69  imageGeometry->IndexToWorld(path->Evaluate(input), worldPoint);
70  imageGeometry->WorldToIndex(worldPoint, index);
71 
72  mitkPixelTypeMultiplex3(ReadPixel, pixelType, image, index, measurementVector.GetDataPointer());
73  intensityProfile->PushBack(measurementVector);
74 
75  offset = path->IncrementInput(input);
76  } while ((offset[0] | offset[1] | offset[2]) != 0);
77 
78  return intensityProfile;
79 }
80 
81 template <class TInputImage>
82 static typename itk::InterpolateImageFunction<TInputImage>::Pointer CreateInterpolateImageFunction(InterpolateImageFunction::Enum interpolator)
83 {
84  switch (interpolator)
85  {
87  return itk::NearestNeighborInterpolateImageFunction<TInputImage>::New().GetPointer();
88 
90  return itk::LinearInterpolateImageFunction<TInputImage>::New().GetPointer();
91 
93  return itk::WindowedSincInterpolateImageFunction<TInputImage, 3, itk::Function::BlackmanWindowFunction<3> >::New().GetPointer();
94 
96  return itk::WindowedSincInterpolateImageFunction<TInputImage, 4, itk::Function::BlackmanWindowFunction<4> >::New().GetPointer();
97 
99  return itk::WindowedSincInterpolateImageFunction<TInputImage, 5, itk::Function::BlackmanWindowFunction<5> >::New().GetPointer();
100 
102  return itk::WindowedSincInterpolateImageFunction<TInputImage, 3, itk::Function::CosineWindowFunction<3> >::New().GetPointer();
103 
105  return itk::WindowedSincInterpolateImageFunction<TInputImage, 4, itk::Function::CosineWindowFunction<4> >::New().GetPointer();
106 
108  return itk::WindowedSincInterpolateImageFunction<TInputImage, 5, itk::Function::CosineWindowFunction<5> >::New().GetPointer();
109 
111  return itk::WindowedSincInterpolateImageFunction<TInputImage, 3, itk::Function::HammingWindowFunction<3> >::New().GetPointer();
112 
114  return itk::WindowedSincInterpolateImageFunction<TInputImage, 4, itk::Function::HammingWindowFunction<4> >::New().GetPointer();
115 
117  return itk::WindowedSincInterpolateImageFunction<TInputImage, 5, itk::Function::HammingWindowFunction<5> >::New().GetPointer();
118 
120  return itk::WindowedSincInterpolateImageFunction<TInputImage, 3, itk::Function::LanczosWindowFunction<3> >::New().GetPointer();
121 
123  return itk::WindowedSincInterpolateImageFunction<TInputImage, 4, itk::Function::LanczosWindowFunction<4> >::New().GetPointer();
124 
126  return itk::WindowedSincInterpolateImageFunction<TInputImage, 5, itk::Function::LanczosWindowFunction<5> >::New().GetPointer();
127 
129  return itk::WindowedSincInterpolateImageFunction<TInputImage, 3, itk::Function::WelchWindowFunction<3> >::New().GetPointer();
130 
132  return itk::WindowedSincInterpolateImageFunction<TInputImage, 4, itk::Function::WelchWindowFunction<4> >::New().GetPointer();
133 
135  return itk::WindowedSincInterpolateImageFunction<TInputImage, 5, itk::Function::WelchWindowFunction<5> >::New().GetPointer();
136 
137  default:
138  return itk::NearestNeighborInterpolateImageFunction<TInputImage>::New().GetPointer();
139  }
140 }
141 
142 template <class TPixel, unsigned int VImageDimension>
143 static void ComputeIntensityProfile(itk::Image<TPixel, VImageDimension>* image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator, IntensityProfile::Pointer intensityProfile)
144 {
145  typename itk::InterpolateImageFunction<itk::Image<TPixel, VImageDimension> >::Pointer interpolateImageFunction = CreateInterpolateImageFunction<itk::Image<TPixel, VImageDimension> >(interpolator);
146  interpolateImageFunction->SetInputImage(image);
147 
148  const itk::PolyLineParametricPath<3>::InputType startOfInput = path->StartOfInput();
149  const itk::PolyLineParametricPath<3>::InputType delta = 1.0 / (numSamples - 1);
150 
151  IntensityProfile::MeasurementVectorType measurementVector;
152 
153  for (unsigned int i = 0; i < numSamples; ++i)
154  {
155  measurementVector[0] = interpolateImageFunction->EvaluateAtContinuousIndex(path->Evaluate(startOfInput + i * delta));
156  intensityProfile->PushBack(measurementVector);
157  }
158 }
159 
160 static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath<3>::Pointer path, unsigned int numSamples, InterpolateImageFunction::Enum interpolator)
161 {
162  IntensityProfile::Pointer intensityProfile = IntensityProfile::New();
163  AccessFixedDimensionByItk_n(image, ComputeIntensityProfile, 3, (path, numSamples, interpolator, intensityProfile));
164  return intensityProfile;
165 }
166 
167 class AddPolyLineElementToPath
168 {
169 public:
170  AddPolyLineElementToPath(const PlaneGeometry* planarFigureGeometry, const BaseGeometry* imageGeometry, itk::PolyLineParametricPath<3>::Pointer path)
171  : m_PlanarFigureGeometry(planarFigureGeometry),
172  m_ImageGeometry(imageGeometry),
173  m_Path(path)
174  {
175  }
176 
177  void operator()(const PlanarFigure::PolyLineElement& polyLineElement)
178  {
179  m_PlanarFigureGeometry->Map(polyLineElement, m_WorldPoint);
180  m_ImageGeometry->WorldToIndex(m_WorldPoint, m_ContinuousIndexPoint);
181  m_Vertex.CastFrom(m_ContinuousIndexPoint);
182  m_Path->AddVertex(m_Vertex);
183  }
184 
185 private:
186  const PlaneGeometry* m_PlanarFigureGeometry;
187  const BaseGeometry* m_ImageGeometry;
188  itk::PolyLineParametricPath<3>::Pointer m_Path;
189 
190  Point3D m_WorldPoint;
191  Point3D m_ContinuousIndexPoint;
192  itk::PolyLineParametricPath<3>::ContinuousIndexType m_Vertex;
193 };
194 
195 static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPlanarFigure(BaseGeometry* imageGeometry, PlanarFigure* planarFigure)
196 {
197  itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New();
198  const PlanarFigure::PolyLineType polyLine = planarFigure->GetPolyLine(0);
199 
200  std::for_each(polyLine.begin(), polyLine.end(),
201  AddPolyLineElementToPath(planarFigure->GetPlaneGeometry(), imageGeometry, path));
202 
203  return path;
204 }
205 
206 static void AddPointToPath(const BaseGeometry* imageGeometry, const Point3D& point, itk::PolyLineParametricPath<3>::Pointer path)
207 {
208  Point3D continuousIndexPoint;
209  imageGeometry->WorldToIndex(point, continuousIndexPoint);
210 
211  itk::PolyLineParametricPath<3>::ContinuousIndexType vertex;
212  vertex.CastFrom(continuousIndexPoint);
213 
214  path->AddVertex(vertex);
215 }
216 
217 static itk::PolyLineParametricPath<3>::Pointer CreatePathFromPoints(BaseGeometry* imageGeometry, const Point3D& startPoint, const Point3D& endPoint)
218 {
219  itk::PolyLineParametricPath<3>::Pointer path = itk::PolyLineParametricPath<3>::New();
220 
221  AddPointToPath(imageGeometry, startPoint, path);
222  AddPointToPath(imageGeometry, endPoint, path);
223 
224  return path;
225 }
226 
228 {
229  return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarFigure));
230 }
231 
232 IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, PlanarLine::Pointer planarLine, unsigned int numSamples, InterpolateImageFunction::Enum interpolator)
233 {
234  return ::ComputeIntensityProfile(image, CreatePathFromPlanarFigure(image->GetGeometry(), planarLine.GetPointer()), numSamples, interpolator);
235 }
236 
237 IntensityProfile::Pointer mitk::ComputeIntensityProfile(Image::Pointer image, const Point3D& startPoint, const Point3D& endPoint, unsigned int numSamples, InterpolateImageFunction::Enum interpolator)
238 {
239  return ::ComputeIntensityProfile(image, CreatePathFromPoints(image->GetGeometry(), startPoint, endPoint), numSamples, interpolator);
240 }
241 
242 IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMaximum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &max)
243 {
245  IntensityProfile::InstanceIdentifier maxIndex = 0;
246 
247  IntensityProfile::ConstIterator end = intensityProfile->End();
248  IntensityProfile::MeasurementType measurement;
249 
250  for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it)
251  {
252  measurement = it.GetMeasurementVector()[0];
253 
254  if (measurement > max)
255  {
256  max = measurement;
257  maxIndex = it.GetInstanceIdentifier();
258  }
259  }
260 
261  return maxIndex;
262 }
263 
264 IntensityProfile::InstanceIdentifier mitk::ComputeGlobalMinimum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &min)
265 {
267  IntensityProfile::InstanceIdentifier minIndex = 0;
268 
269  IntensityProfile::ConstIterator end = intensityProfile->End();
270  IntensityProfile::MeasurementType measurement;
271 
272  for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it)
273  {
274  measurement = it.GetMeasurementVector()[0];
275 
276  if (measurement < min)
277  {
278  min = measurement;
279  minIndex = it.GetInstanceIdentifier();
280  }
281  }
282 
283  return minIndex;
284 }
285 
286 IntensityProfile::InstanceIdentifier mitk::ComputeCenterOfMaximumArea(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::InstanceIdentifier radius)
287 {
288  //const IntensityProfile::MeasurementType min = intensityProfile->GetMeasurementVector(ComputeGlobalMinimum(intensityProfile))[0];
289  IntensityProfile::MeasurementType min;
290  ComputeGlobalMinimum(intensityProfile, min);
291  const IntensityProfile::InstanceIdentifier areaWidth = 1 + 2 * radius;
292 
293  IntensityProfile::MeasurementType maxArea = 0;
294 
295  for (IntensityProfile::InstanceIdentifier i = 0; i < areaWidth; ++i)
296  maxArea += intensityProfile->GetMeasurementVector(i)[0] - min;
297 
298  const IntensityProfile::InstanceIdentifier lastIndex = intensityProfile->Size() - areaWidth;
299  IntensityProfile::InstanceIdentifier centerOfMaxArea = radius;
300  IntensityProfile::MeasurementType area = maxArea;
301 
302  for (IntensityProfile::InstanceIdentifier i = 1; i <= lastIndex; ++i)
303  {
304  area += intensityProfile->GetMeasurementVector(i + areaWidth - 1)[0] - min;
305  area -= intensityProfile->GetMeasurementVector(i - 1)[0] - min;
306 
307  if (area > maxArea)
308  {
309  maxArea = area;
310  centerOfMaxArea = i + radius; // TODO: If multiple areas in the neighborhood have the same intensity chose the middle one instead of the first one.
311  }
312  }
313 
314  return centerOfMaxArea;
315 }
316 
317 std::vector<IntensityProfile::MeasurementType> mitk::CreateVectorFromIntensityProfile(IntensityProfile::ConstPointer intensityProfile)
318 {
319  std::vector<IntensityProfile::MeasurementType> result;
320  result.reserve(intensityProfile->Size());
321 
322  IntensityProfile::ConstIterator end = intensityProfile->End();
323 
324  for (IntensityProfile::ConstIterator it = intensityProfile->Begin(); it != end; ++it)
325  result.push_back(it.GetMeasurementVector()[0]);
326 
327  return result;
328 }
329 
330 IntensityProfile::Pointer mitk::CreateIntensityProfileFromVector(const std::vector<IntensityProfile::MeasurementType>& vector)
331 {
332  const IntensityProfile::InstanceIdentifier size = vector.size();
333 
334  IntensityProfile::Pointer result = IntensityProfile::New();
335  result->Resize(size);
336 
337  for (IntensityProfile::InstanceIdentifier i = 0; i < size; ++i)
338  result->SetMeasurement(i, 0, vector[i]);
339 
340  return result;
341 }
342 
343 void mitk::ComputeIntensityProfileStatistics(IntensityProfile::ConstPointer intensityProfile, ImageStatisticsContainer::ImageStatisticsObject& stats)
344 {
345  typedef std::vector<IntensityProfile::MeasurementType> StatsVecType;
346 
347  StatsVecType statsVec = mitk::CreateVectorFromIntensityProfile( intensityProfile );
348 
349  IntensityProfile::MeasurementType min;
350  IntensityProfile::MeasurementType max;
351  mitk::ComputeGlobalMinimum( intensityProfile, min );
352  mitk::ComputeGlobalMaximum( intensityProfile, max );
353  auto numSamples = static_cast<mitk::ImageStatisticsContainer::VoxelCountType>(statsVec.size());
354 
355  double mean = 0.0;
356  double rms = 0.0;
357  for ( StatsVecType::const_iterator it = statsVec.begin(); it != statsVec.end(); ++it )
358  {
359  double val = *it;
360  mean += val;
361  rms += val*val;
362  }
363  mean /= static_cast<double>(numSamples);
364  rms /= static_cast<double>(numSamples);
365 
366  double var = 0.0;
367  for ( StatsVecType::const_iterator it = statsVec.begin(); it != statsVec.end(); ++it )
368  {
369  double diff = *it - mean;
370  var += diff*diff;
371  }
372  var /= (static_cast<double>(numSamples) - 1 );
373 
374  rms = sqrt( rms );
375 
383 }
void IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const
Convert (continuous or discrete) index coordinates of a vector vec_units to world coordinates (in mm)...
Gives locked and index-based read access for a particular image part. The class provides several set-...
virtual const PlaneGeometry * GetPlaneGeometry() const
Returns (previously set) 2D geometry of this figure.
double ScalarType
MITKIMAGESTATISTICS_EXPORT IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, PlanarFigure::Pointer planarFigure)
Compute intensity profile of an image for each pixel along the first PolyLine of a given planar figur...
#define AccessFixedDimensionByItk_n(mitkImage, itkImageTypeFunction, dimension, va_tuple)
Access a mitk-image with known dimension by an itk-image with one or more parameters.
MITKIMAGESTATISTICS_EXPORT IntensityProfile::InstanceIdentifier ComputeCenterOfMaximumArea(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::InstanceIdentifier radius)
Compute center of maximum area under the curve of an intensity profile.
const TPixel & GetPixelByIndex(const itk::Index< VDimension > &idx) const
static IntensityProfile::Pointer ComputeIntensityProfile(Image::Pointer image, itk::PolyLineParametricPath< 3 >::Pointer path)
static itk::PolyLineParametricPath< 3 >::Pointer CreatePathFromPlanarFigure(BaseGeometry *imageGeometry, PlanarFigure *planarFigure)
DataCollection - Class to facilitate loading/accessing structured data.
MITKIMAGESTATISTICS_EXPORT IntensityProfile::InstanceIdentifier ComputeGlobalMaximum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &max)
Compute global maximum of an intensity profile.
void AddStatistic(const std::string &key, StatisticsVariantType value)
Adds a statistic to the statistics object.
const PolyLineType GetPolyLine(unsigned int index)
Returns the polyline representing the planar figure (for rendering, measurements, etc...
static void AddPointToPath(const BaseGeometry *imageGeometry, const Point3D &point, itk::PolyLineParametricPath< 3 >::Pointer path)
MITKIMAGESTATISTICS_EXPORT std::vector< IntensityProfile::MeasurementType > CreateVectorFromIntensityProfile(IntensityProfile::ConstPointer intensityProfile)
Convert an intensity profile to a standard library vector.
MITKIMAGESTATISTICS_EXPORT IntensityProfile::InstanceIdentifier ComputeGlobalMinimum(IntensityProfile::ConstPointer intensityProfile, IntensityProfile::MeasurementType &min)
Compute global minimum of an intensity profile.
static Vector3D offset
static itk::PolyLineParametricPath< 3 >::Pointer CreatePathFromPoints(BaseGeometry *imageGeometry, const Point3D &startPoint, const Point3D &endPoint)
#define mitkPixelTypeMultiplex3(function, ptype, param1, param2, param3)
#define mitkThrow()
MITKIMAGESTATISTICS_EXPORT IntensityProfile::Pointer CreateIntensityProfileFromVector(const std::vector< IntensityProfile::MeasurementType > &vector)
Convert a standard library vector to an intensity profile.
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
void WorldToIndex(const mitk::Point3D &pt_mm, mitk::Point3D &pt_units) const
Convert world coordinates (in mm) of a point to (continuous!) index coordinates.
static T min(T x, T y)
Definition: svm.cpp:53
static const std::string STANDARDDEVIATION()
Base-class for geometric planar (2D) figures, such as lines, circles, rectangles, polygons...
MITKIMAGESTATISTICS_EXPORT void ComputeIntensityProfileStatistics(IntensityProfile::ConstPointer intensityProfile, ImageStatisticsContainer::ImageStatisticsObject &stats)
Compute statistics of an intensity profile.
std::vector< PolyLineElement > PolyLineType
Describes a two-dimensional, rectangular plane.
static void ReadPixel(const PixelType &, Image::Pointer image, const itk::Index< 3 > &index, ScalarType *returnValue)
Container class for storing the computed image statistics.
static itk::InterpolateImageFunction< TInputImage >::Pointer CreateInterpolateImageFunction(InterpolateImageFunction::Enum interpolator)
BaseGeometry Describes the geometry of a data object.
Class for defining the data type of pixels.
Definition: mitkPixelType.h:51