Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkGIFVolumetricDensityStatistics.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 
14 
15 // MITK
16 #include <mitkITKImageImport.h>
17 #include <mitkImageCast.h>
18 #include <mitkImageAccessByItk.h>
19 #include <mitkPixelTypeMultiplex.h>
21 
22 // ITK
23 #include <itkLabelStatisticsImageFilter.h>
24 #include <itkNeighborhoodIterator.h>
25 #include <itkImageRegionConstIteratorWithIndex.h>
26 #include <itkLabelGeometryImageFilter.h>
27 
28 // VTK
29 #include <vtkSmartPointer.h>
30 #include <vtkImageMarchingCubes.h>
31 #include <vtkMassProperties.h>
32 #include <vtkDelaunay3D.h>
33 #include <vtkGeometryFilter.h>
34 #include <vtkDoubleArray.h>
35 #include <vtkPCAStatistics.h>
36 #include <vtkTable.h>
37 
38 // STL
39 #include <limits>
40 #include <vnl/vnl_math.h>
41 
42 // Eigen
43 #include <Eigen/Dense>
44 
45 struct GIFVolumetricDensityStatisticsParameters
46 {
47  double volume;
48  std::string prefix;
49 };
50 
51 template<typename TPixel, unsigned int VImageDimension>
52 void
53 CalculateVolumeDensityStatistic(itk::Image<TPixel, VImageDimension>* itkImage, mitk::Image::Pointer mask, GIFVolumetricDensityStatisticsParameters params, mitk::GIFVolumetricDensityStatistics::FeatureListType & featureList)
54 {
55  typedef itk::Image<TPixel, VImageDimension> ImageType;
56  typedef itk::Image<unsigned short, VImageDimension> MaskType;
57 
58  double volume = params.volume;
59  std::string prefix = params.prefix;
60 
61  typename MaskType::Pointer maskImage = MaskType::New();
62  mitk::CastToItkImage(mask, maskImage);
63 
64  itk::ImageRegionConstIteratorWithIndex<ImageType> imgA(itkImage, itkImage->GetLargestPossibleRegion());
65  itk::ImageRegionConstIteratorWithIndex<ImageType> imgB(itkImage, itkImage->GetLargestPossibleRegion());
66  itk::ImageRegionConstIteratorWithIndex<MaskType> maskA(maskImage, maskImage->GetLargestPossibleRegion());
67  itk::ImageRegionConstIteratorWithIndex<MaskType> maskB(maskImage, maskImage->GetLargestPossibleRegion());
68 
69  double moranA = 0;
70  double moranB = 0;
71  double geary = 0;
72  double Nv = 0;
73  double w_ij = 0;
74  double mean = 0;
75 
76  typename ImageType::PointType pointA;
77  typename ImageType::PointType pointB;
78 
79  while (!imgA.IsAtEnd())
80  {
81  if (maskA.Get() > 0)
82  {
83  Nv += 1;
84  mean += imgA.Get();
85  }
86  ++imgA;
87  ++maskA;
88  }
89  mean /= Nv;
90  imgA.GoToBegin();
91  maskA.GoToBegin();
92 
93  while (!imgA.IsAtEnd())
94  {
95  if (maskA.Get() > 0)
96  {
97  imgB.GoToBegin();
98  maskB.GoToBegin();
99  while (!imgB.IsAtEnd())
100  {
101  if ((imgA.GetIndex() == imgB.GetIndex()) ||
102  (maskB.Get() < 1))
103  {
104  ++imgB;
105  ++maskB;
106  continue;
107  }
108  itkImage->TransformIndexToPhysicalPoint(maskA.GetIndex(), pointA);
109  itkImage->TransformIndexToPhysicalPoint(maskB.GetIndex(), pointB);
110 
111  double w = 1 / pointA.EuclideanDistanceTo(pointB);
112  moranA += w*(imgA.Get() - mean)* (imgB.Get() - mean);
113  geary += w * (imgA.Get() - imgB.Get()) * (imgA.Get() - imgB.Get());
114 
115  w_ij += w;
116 
117  ++imgB;
118  ++maskB;
119  }
120  moranB += (imgA.Get() - mean)* (imgA.Get() - mean);
121  }
122  ++imgA;
123  ++maskA;
124  }
125 
126  MITK_INFO << "Volume: " << volume;
127  MITK_INFO << " Mean: " << mean;
128  featureList.push_back(std::make_pair(prefix + "Volume integrated intensity", volume* mean));
129  featureList.push_back(std::make_pair(prefix + "Volume Moran's I index", Nv / w_ij * moranA / moranB));
130  featureList.push_back(std::make_pair(prefix + "Volume Geary's C measure", ( Nv -1 ) / 2 / w_ij * geary/ moranB));
131 }
132 
133 void calculateMOBB(vtkPointSet *pointset, double &volume, double &surface)
134 {
136 
137  for (int cellID = 0; cellID < pointset->GetNumberOfCells(); ++cellID)
138  {
139  auto cell = pointset->GetCell(cellID);
140 
141  for (int edgeID = 0; edgeID < 3; ++edgeID)
142  {
143  auto edge = cell->GetEdge(edgeID);
144 
145  double pA[3], pB[3];
146  double pAA[3], pBB[3];
147 
148  vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
149  transform->PostMultiply();
150  pointset->GetPoint(edge->GetPointId(0), pA);
151  pointset->GetPoint(edge->GetPointId(1), pB);
152 
153  double angleZ = std::atan2((- pA[2] + pB[2]) ,(pA[1] - pB[1]));
154  angleZ *= 180 / vnl_math::pi;
155  if (pA[2] == pB[2])
156  angleZ = 0;
157 
158  transform->RotateX(angleZ);
159  transform->TransformPoint(pA, pAA);
160  transform->TransformPoint(pB, pBB);
161 
162  double angleY = std::atan2((pAA[1] -pBB[1]) ,-(pAA[0] - pBB[0]));
163  angleY *= 180 / vnl_math::pi;
164  if (pAA[1] == pBB[1])
165  angleY = 0;
166  transform->RotateZ(angleY);
167 
168  double p0[3];
169  pointset->GetPoint(edge->GetPointId(0), p0);
170 
171  double curMinX = std::numeric_limits<double>::max();
172  double curMaxX = std::numeric_limits<double>::lowest();
173  double curMinY = std::numeric_limits<double>::max();
174  double curMaxY = std::numeric_limits<double>::lowest();
175  double curMinZ = std::numeric_limits<double>::max();
176  double curMaxZ = std::numeric_limits<double>::lowest();
177  for (int pointID = 0; pointID < pointset->GetNumberOfPoints(); ++pointID)
178  {
179  double p[3];
180  double p2[3];
181  pointset->GetPoint(pointID, p);
182  p[0] -= p0[0]; p[1] -= p0[1]; p[2] -= p0[2];
183  transform->TransformPoint(p, p2);
184 
185  curMinX = std::min<double>(p2[0], curMinX);
186  curMaxX = std::max<double>(p2[0], curMaxX);
187  curMinY = std::min<double>(p2[1], curMinY);
188  curMaxY = std::max<double>(p2[1], curMaxY);
189  curMinZ = std::min<double>(p2[2], curMinZ);
190  curMaxZ = std::max<double>(p2[2], curMaxZ);
191  }
192 
193  if ((curMaxX - curMinX)*(curMaxY - curMinY)*(curMaxZ - curMinZ) < volume)
194  {
195  volume = (curMaxX - curMinX)*(curMaxY - curMinY)*(curMaxZ - curMinZ);
196  surface = (curMaxX - curMinX)*(curMaxX - curMinX) + (curMaxY - curMinY)*(curMaxY - curMinY) + (curMaxZ - curMinZ)*(curMaxZ - curMinZ);
197  surface *= 2;
198  }
199 
200 
201  }
202  }
203 }
204 
205 void calculateMEE(vtkPointSet *pointset, double &vol, double &surf, double tolerance=0.0001)
206 {
207  // Inspired by https://github.com/smdabdoub/ProkaryMetrics/blob/master/calc/fitting.py
208 
209  int numberOfPoints = pointset->GetNumberOfPoints();
210  int dimension = 3;
211  Eigen::MatrixXd points(3, numberOfPoints);
212  Eigen::MatrixXd Q(3+1, numberOfPoints);
213  double p[3];
214 
215  std::cout << "Initialize Q " << std::endl;
216  for (int i = 0; i < numberOfPoints; ++i)
217  {
218  pointset->GetPoint(i, p);
219  points(0, i) = p[0];
220  points(1, i) = p[1];
221  points(2, i) = p[2];
222  Q(0, i) = p[0];
223  Q(1, i) = p[1];
224  Q(2, i) = p[2];
225  Q(3, i) = 1.0;
226  }
227 
228  int count = 1;
229  double error = 1;
230  Eigen::VectorXd u_vector(numberOfPoints);
231  u_vector.fill(1.0 / numberOfPoints);
232  Eigen::DiagonalMatrix<double, Eigen::Dynamic> u = u_vector.asDiagonal();
233  Eigen::VectorXd ones(dimension + 1);
234  ones.fill(1);
235  Eigen::MatrixXd Ones = ones.asDiagonal();
236 
237  // Khachiyan Algorithm
238  while (error > tolerance)
239  {
240  auto Qt = Q.transpose();
241  Eigen::MatrixXd X = Q*u*Qt;
242  Eigen::FullPivHouseholderQR<Eigen::MatrixXd> qr(X);
243  Eigen::MatrixXd Xi = qr.solve(Ones);
244 
245  Eigen::MatrixXd M = Qt * Xi * Q;
246 
247  double maximumValue = M(0, 0);
248  int maximumPosition = 0;
249  for (int i = 0; i < numberOfPoints; ++i)
250  {
251  if (maximumValue < M(i, i))
252  {
253  maximumValue = M(i, i);
254  maximumPosition = i;
255  }
256  }
257  double stepsize = (maximumValue - dimension - 1) / ((dimension + 1) * (maximumValue - 1));
258  Eigen::DiagonalMatrix<double, Eigen::Dynamic> new_u = (1.0 - stepsize) * u;
259  new_u.diagonal()[maximumPosition] = (new_u.diagonal())(maximumPosition) + stepsize;
260  ++count;
261  error = (new_u.diagonal() - u.diagonal()).norm();
262  u.diagonal() = new_u.diagonal();
263  }
264 
265  // U = u
266 
267  Eigen::MatrixXd Ai = points * u * points.transpose() - points * u *(points * u).transpose();
268  Eigen::FullPivHouseholderQR<Eigen::MatrixXd> qr(Ai);
269  Eigen::VectorXd ones2(dimension);
270  ones2.fill(1);
271  Eigen::MatrixXd Ones2 = ones2.asDiagonal();
272  Eigen::MatrixXd A = qr.solve(Ones2)*1.0/dimension;
273 
274  Eigen::JacobiSVD<Eigen::MatrixXd> svd(A);
275  double c = 1 / sqrt(svd.singularValues()[0]);
276  double b = 1 / sqrt(svd.singularValues()[1]);
277  double a = 1 / sqrt(svd.singularValues()[2]);
278  double V = 4 * vnl_math::pi*a*b*c / 3;
279 
280  double ad_mvee= 0;
281  double alpha = std::sqrt(1 - b*b / a / a);
282  double beta = std::sqrt(1 - c*c / a / a);
283  for (int i = 0; i < 20; ++i)
284  {
285  ad_mvee += 4 * vnl_math::pi*a*b*(alpha*alpha + beta*beta) / (2 * alpha*beta) * (std::pow(alpha*beta, i)) / (1 - 4 * i*i);
286  }
287  vol = V;
288  surf = ad_mvee;
289 }
290 
292 {
293  FeatureListType featureList;
294  if (image->GetDimension() < 3)
295  {
296  return featureList;
297  }
298 
299  std::string prefix = FeatureDescriptionPrefix();
300 
301  vtkSmartPointer<vtkImageMarchingCubes> mesher = vtkSmartPointer<vtkImageMarchingCubes>::New();
302  vtkSmartPointer<vtkMassProperties> stats = vtkSmartPointer<vtkMassProperties>::New();
303  vtkSmartPointer<vtkMassProperties> stats2 = vtkSmartPointer<vtkMassProperties>::New();
304  mesher->SetInputData(mask->GetVtkImageData());
305  mesher->SetValue(0, 0.5);
306  stats->SetInputConnection(mesher->GetOutputPort());
307  stats->Update();
308 
309  vtkSmartPointer<vtkDelaunay3D> delaunay =
310  vtkSmartPointer< vtkDelaunay3D >::New();
311  delaunay->SetInputConnection(mesher->GetOutputPort());
312  delaunay->SetAlpha(0);
313  delaunay->Update();
314  vtkSmartPointer<vtkGeometryFilter> geometryFilter =
315  vtkSmartPointer<vtkGeometryFilter>::New();
316  geometryFilter->SetInputConnection(delaunay->GetOutputPort());
317  geometryFilter->Update();
318  stats2->SetInputConnection(geometryFilter->GetOutputPort());
319  stats2->Update();
320 
321  double vol_mvee;
322  double surf_mvee;
323  calculateMEE(mesher->GetOutput(), vol_mvee, surf_mvee);
324 
325  double vol_mobb;
326  double surf_mobb;
327  calculateMOBB(geometryFilter->GetOutput(), vol_mobb, surf_mobb);
328 
329  double pi = vnl_math::pi;
330 
331  double meshVolume = stats->GetVolume();
332  double meshSurf = stats->GetSurfaceArea();
333 
334  GIFVolumetricDensityStatisticsParameters params;
335  params.volume = meshVolume;
336  params.prefix = prefix;
337  AccessByItk_3(image, CalculateVolumeDensityStatistic, mask, params, featureList);
338 
339  //Calculate center of mass shift
340  int xx = mask->GetDimensions()[0];
341  int yy = mask->GetDimensions()[1];
342  int zz = mask->GetDimensions()[2];
343 
344  double xd = mask->GetGeometry()->GetSpacing()[0];
345  double yd = mask->GetGeometry()->GetSpacing()[1];
346  double zd = mask->GetGeometry()->GetSpacing()[2];
347 
348  int minimumX=xx;
349  int maximumX=0;
350  int minimumY=yy;
351  int maximumY=0;
352  int minimumZ=zz;
353  int maximumZ=0;
354 
355  vtkSmartPointer<vtkDoubleArray> dataset1Arr = vtkSmartPointer<vtkDoubleArray>::New();
356  vtkSmartPointer<vtkDoubleArray> dataset2Arr = vtkSmartPointer<vtkDoubleArray>::New();
357  vtkSmartPointer<vtkDoubleArray> dataset3Arr = vtkSmartPointer<vtkDoubleArray>::New();
358  dataset1Arr->SetNumberOfComponents(1);
359  dataset2Arr->SetNumberOfComponents(1);
360  dataset3Arr->SetNumberOfComponents(1);
361  dataset1Arr->SetName("M1");
362  dataset2Arr->SetName("M2");
363  dataset3Arr->SetName("M3");
364 
365  vtkSmartPointer<vtkDoubleArray> dataset1ArrU = vtkSmartPointer<vtkDoubleArray>::New();
366  vtkSmartPointer<vtkDoubleArray> dataset2ArrU = vtkSmartPointer<vtkDoubleArray>::New();
367  vtkSmartPointer<vtkDoubleArray> dataset3ArrU = vtkSmartPointer<vtkDoubleArray>::New();
368  dataset1ArrU->SetNumberOfComponents(1);
369  dataset2ArrU->SetNumberOfComponents(1);
370  dataset3ArrU->SetNumberOfComponents(1);
371  dataset1ArrU->SetName("M1");
372  dataset2ArrU->SetName("M2");
373  dataset3ArrU->SetName("M3");
374 
375  vtkSmartPointer<vtkPoints> points =
376  vtkSmartPointer< vtkPoints >::New();
377 
378  for (int x = 0; x < xx; x++)
379  {
380  for (int y = 0; y < yy; y++)
381  {
382  for (int z = 0; z < zz; z++)
383  {
384  itk::Image<int,3>::IndexType index;
385 
386  index[0] = x;
387  index[1] = y;
388  index[2] = z;
389 
390  mitk::ScalarType pxImage;
391  mitk::ScalarType pxMask;
392 
395  image->GetChannelDescriptor().GetPixelType(),
396  image,
397  image->GetVolumeData(),
398  index,
399  pxImage,
400  0);
401 
404  mask->GetChannelDescriptor().GetPixelType(),
405  mask,
406  mask->GetVolumeData(),
407  index,
408  pxMask,
409  0);
410 
411  //Check if voxel is contained in segmentation
412  if (pxMask > 0)
413  {
414  minimumX = std::min<int>(x, minimumX);
415  minimumY = std::min<int>(y, minimumY);
416  minimumZ = std::min<int>(z, minimumZ);
417  maximumX = std::max<int>(x, maximumX);
418  maximumY = std::max<int>(y, maximumY);
419  maximumZ = std::max<int>(z, maximumZ);
420  points->InsertNextPoint(x*xd, y*yd, z*zd);
421 
422  if (pxImage == pxImage)
423  {
424  dataset1Arr->InsertNextValue(x*xd);
425  dataset2Arr->InsertNextValue(y*yd);
426  dataset3Arr->InsertNextValue(z*zd);
427  }
428  }
429  }
430  }
431  }
432 
433  vtkSmartPointer<vtkTable> datasetTable = vtkSmartPointer<vtkTable>::New();
434  datasetTable->AddColumn(dataset1Arr);
435  datasetTable->AddColumn(dataset2Arr);
436  datasetTable->AddColumn(dataset3Arr);
437 
438  vtkSmartPointer<vtkPCAStatistics> pcaStatistics = vtkSmartPointer<vtkPCAStatistics>::New();
439  pcaStatistics->SetInputData(vtkStatisticsAlgorithm::INPUT_DATA, datasetTable);
440  pcaStatistics->SetColumnStatus("M1", 1);
441  pcaStatistics->SetColumnStatus("M2", 1);
442  pcaStatistics->SetColumnStatus("M3", 1);
443  pcaStatistics->RequestSelectedColumns();
444  pcaStatistics->SetDeriveOption(true);
445  pcaStatistics->Update();
446 
447  vtkSmartPointer<vtkDoubleArray> eigenvalues = vtkSmartPointer<vtkDoubleArray>::New();
448  pcaStatistics->GetEigenvalues(eigenvalues);
449 
450  std::vector<double> eigen_val(3);
451  eigen_val[2] = eigenvalues->GetValue(0);
452  eigen_val[1] = eigenvalues->GetValue(1);
453  eigen_val[0] = eigenvalues->GetValue(2);
454 
455  double major = 2*sqrt(eigen_val[2]);
456  double minor = 2*sqrt(eigen_val[1]);
457  double least = 2*sqrt(eigen_val[0]);
458 
459  double alpha = std::sqrt(1 - minor*minor / major / major);
460  double beta = std::sqrt(1 - least*least / major / major);
461 
462  double a = (maximumX - minimumX+1) * xd;
463  double b = (maximumY - minimumY+1) * yd;
464  double c = (maximumZ - minimumZ+1) * zd;
465 
466  double vd_aabb = meshVolume / (a*b*c);
467  double ad_aabb = meshSurf / (2 * a*b + 2 * a*c + 2 * b*c);
468 
469  double vd_aee = 3 * meshVolume / (4.0*pi*major*minor*least);
470  double ad_aee = 0;
471  for (int i = 0; i < 20; ++i)
472  {
473  ad_aee += 4 * pi*major*minor*(alpha*alpha + beta*beta) / (2 * alpha*beta) * (std::pow(alpha*beta, i)) / (1 - 4 * i*i);
474  }
475  ad_aee = meshSurf / ad_aee;
476 
477  double vd_ch = meshVolume / stats2->GetVolume();
478  double ad_ch = meshSurf / stats2->GetSurfaceArea();
479 
480  featureList.push_back(std::make_pair(prefix + "Volume density axis-aligned bounding box", vd_aabb));
481  featureList.push_back(std::make_pair(prefix + "Surface density axis-aligned bounding box", ad_aabb));
482  featureList.push_back(std::make_pair(prefix + "Volume density oriented minimum bounding box", meshVolume / vol_mobb));
483  featureList.push_back(std::make_pair(prefix + "Surface density oriented minimum bounding box", meshSurf / surf_mobb));
484  featureList.push_back(std::make_pair(prefix + "Volume density approx. enclosing ellipsoid", vd_aee));
485  featureList.push_back(std::make_pair(prefix + "Surface density approx. enclosing ellipsoid", ad_aee));
486  featureList.push_back(std::make_pair(prefix + "Volume density approx. minimum volume enclosing ellipsoid", meshVolume / vol_mvee));
487  featureList.push_back(std::make_pair(prefix + "Surface density approx. minimum volume enclosing ellipsoid", meshSurf / surf_mvee));
488  featureList.push_back(std::make_pair(prefix + "Volume density convex hull", vd_ch));
489  featureList.push_back(std::make_pair(prefix + "Surface density convex hull", ad_ch));
490 
491  return featureList;
492 }
493 
495 {
496  SetLongName("volume-density");
497  SetShortName("volden");
498  SetFeatureClassName("Morphological Density");
499 }
500 
502 {
503  FeatureNameListType featureList;
504  return featureList;
505 }
506 
507 
509 {
510  std::string name = GetOptionPrefix();
511 
512  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Volume-Density Statistic", "calculates volume density based features", us::Any());
513 }
514 
515 void
517 {
518  auto parsedArgs = GetParameter();
519  if (parsedArgs.count(GetLongName()))
520  {
521  MITK_INFO << "Start calculating volumetric density features ....";
522  auto localResults = this->CalculateFeatures(feature, mask);
523  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
524  MITK_INFO << "Finished calculating volumetric density features....";
525  }
526 }
527 
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
mitk::ScalarType FastSinglePixelAccess(mitk::PixelType, mitk::Image::Pointer im, ImageDataItem *item, itk::Index< 3 > idx, mitk::ScalarType &val, int component=0)
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
void calculateMEE(vtkPointSet *pointset, double &vol, double &surf, double tolerance=0.0001)
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
double ScalarType
void CalculateVolumeDensityStatistic(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, GIFVolumetricDensityStatisticsParameters params, mitk::GIFVolumetricDensityStatistics::FeatureListType &featureList)
void CalculateFeaturesUsingParameters(const Image::Pointer &feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) override
Calculates the feature of this abstact interface. Does not necessarily considers the parameter settin...
void addArgument(const std::string &longarg, const std::string &shortarg, Type type, const std::string &argLabel, const std::string &argHelp=std::string(), const us::Any &defaultValue=us::Any(), bool optional=true, bool ignoreRest=false, bool deprecated=false, mitkCommandLineParser::Channel channel=mitkCommandLineParser::Channel::None)
void calculateMOBB(vtkPointSet *pointset, double &volume, double &surface)
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for this class.
std::vector< std::pair< std::string, double > > FeatureListType
Definition: usAny.h:163
void AddArguments(mitkCommandLineParser &parser) override
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
mitk::Image::Pointer mask
#define mitkPixelTypeMultiplex5(function, ptype, param1, param2, param3, param4, param5)
void CalculateFeatures(mitk::CoocurenceMatrixHolder &holder, mitk::CoocurenceMatrixFeatures &results)