Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkGIFVolumetricStatistics.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 
18 
19 // MITK
20 #include <mitkITKImageImport.h>
21 #include <mitkImageCast.h>
22 #include <mitkImageAccessByItk.h>
23 
24 // ITK
25 #include <itkLabelStatisticsImageFilter.h>
26 #include <itkNeighborhoodIterator.h>
27 
28 // VTK
29 #include <vtkSmartPointer.h>
30 #include <vtkImageMarchingCubes.h>
31 #include <vtkMassProperties.h>
32 
33 // STL
34 #include <sstream>
35 #include <vnl/vnl_math.h>
36 
37 template<typename TPixel, unsigned int VImageDimension>
38 void
39  CalculateVolumeStatistic(itk::Image<TPixel, VImageDimension>* itkImage, mitk::Image::Pointer mask, mitk::GIFVolumetricStatistics::FeatureListType & featureList)
40 {
41  typedef itk::Image<TPixel, VImageDimension> ImageType;
42  typedef itk::Image<int, VImageDimension> MaskType;
43  typedef itk::LabelStatisticsImageFilter<ImageType, MaskType> FilterType;
44 
45  typename MaskType::Pointer maskImage = MaskType::New();
46  mitk::CastToItkImage(mask, maskImage);
47 
48  typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New();
49  labelStatisticsImageFilter->SetInput( itkImage );
50  labelStatisticsImageFilter->SetLabelInput(maskImage);
51  labelStatisticsImageFilter->Update();
52 
53  double volume = labelStatisticsImageFilter->GetCount(1);
54  for (int i = 0; i < (int)(VImageDimension); ++i)
55  {
56  volume *= itkImage->GetSpacing()[i];
57  }
58 
59  featureList.push_back(std::make_pair("Volumetric Features Volume (pixel based)",volume));
60 }
61 
62 template<typename TPixel, unsigned int VImageDimension>
63 void
64  CalculateLargestDiameter(itk::Image<TPixel, VImageDimension>* mask, mitk::GIFVolumetricStatistics::FeatureListType & featureList)
65 {
66  typedef itk::Image<TPixel, VImageDimension> ImageType;
67  typedef typename ImageType::PointType PointType;
68  typename ImageType::SizeType radius;
69  for (int i=0; i < (int)VImageDimension; ++i)
70  radius[i] = 1;
71  itk::NeighborhoodIterator<ImageType> iterator(radius,mask, mask->GetRequestedRegion());
72  std::vector<PointType> borderPoints;
73  while(!iterator.IsAtEnd())
74  {
75  if (iterator.GetCenterPixel() == 0)
76  {
77  ++iterator;
78  continue;
79  }
80 
81  bool border = false;
82  for (int i = 0; i < (int)(iterator.Size()); ++i)
83  {
84  if (iterator.GetPixel(i) == 0)
85  {
86  border = true;
87  break;
88  }
89  }
90  if (border)
91  {
92  auto centerIndex = iterator.GetIndex();
93  PointType centerPoint;
94  mask->TransformIndexToPhysicalPoint(centerIndex, centerPoint );
95  borderPoints.push_back(centerPoint);
96  }
97  ++iterator;
98  }
99 
100  double longestDiameter = 0;
101  unsigned long numberOfBorderPoints = borderPoints.size();
102  for (int i = 0; i < (int)numberOfBorderPoints; ++i)
103  {
104  auto point = borderPoints[i];
105  for (int j = i; j < (int)numberOfBorderPoints; ++j)
106  {
107  double newDiameter=point.EuclideanDistanceTo(borderPoints[j]);
108  if (newDiameter > longestDiameter)
109  longestDiameter = newDiameter;
110  }
111  }
112 
113  featureList.push_back(std::make_pair("Volumetric Features Maximum 3D diameter",longestDiameter));
114 }
115 
117 {
118 }
119 
121 {
122  FeatureListType featureList;
123 
124  AccessByItk_2(image, CalculateVolumeStatistic, mask, featureList);
125  AccessByItk_1(mask, CalculateLargestDiameter, featureList);
126 
127  vtkSmartPointer<vtkImageMarchingCubes> mesher = vtkSmartPointer<vtkImageMarchingCubes>::New();
128  vtkSmartPointer<vtkMassProperties> stats = vtkSmartPointer<vtkMassProperties>::New();
129  mesher->SetInputData(mask->GetVtkImageData());
130  stats->SetInputConnection(mesher->GetOutputPort());
131  stats->Update();
132 
133  double pi = vnl_math::pi;
134 
135  double meshVolume = stats->GetVolume();
136  double meshSurf = stats->GetSurfaceArea();
137  double pixelVolume = featureList[0].second;
138 
139  double compactness1 = pixelVolume / ( std::sqrt(pi) * std::pow(meshSurf, 2.0/3.0));
140  double compactness2 = 36*pi*pixelVolume*pixelVolume/meshSurf/meshSurf/meshSurf;
141 
142  double sphericity=std::pow(pi,1/3.0) *std::pow(6*pixelVolume, 2.0/3.0) / meshSurf;
143  double surfaceToVolume = meshSurf / pixelVolume;
144  double sphericalDisproportion = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0);
145 
146  featureList.push_back(std::make_pair("Volumetric Features Volume (mesh based)",meshVolume));
147  featureList.push_back(std::make_pair("Volumetric Features Surface area",meshSurf));
148  featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio",surfaceToVolume));
149  featureList.push_back(std::make_pair("Volumetric Features Sphericity",sphericity));
150  featureList.push_back(std::make_pair("Volumetric Features Compactness 1",compactness1));
151  featureList.push_back(std::make_pair("Volumetric Features Compactness 2",compactness2));
152  featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion",sphericalDisproportion));
153 
154  return featureList;
155 }
156 
158 {
159  FeatureNameListType featureList;
160  featureList.push_back("Volumetric Features Compactness 1");
161  featureList.push_back("Volumetric Features Compactness 2");
162  featureList.push_back("Volumetric Features Sphericity");
163  featureList.push_back("Volumetric Features Surface to volume ratio");
164  featureList.push_back("Volumetric Features Surface area");
165  featureList.push_back("Volumetric Features Volume (mesh based)");
166  featureList.push_back("Volumetric Features Volume (pixel based)");
167  featureList.push_back("Volumetric Features Spherical disproportion");
168  featureList.push_back("Volumetric Features Maximum 3D diameter");
169  return featureList;
170 }
void CalculateVolumeStatistic(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, mitk::GIFVolumetricStatistics::FeatureListType &featureList)
mitk::Point3D PointType
itk::SmartPointer< Self > Pointer
itk::Image< unsigned char, 3 > ImageType
#define AccessByItk_1(mitkImage, itkImageTypeFunction, arg1)
std::vector< std::pair< std::string, double > > FeatureListType
void CalculateLargestDiameter(itk::Image< TPixel, VImageDimension > *mask, mitk::GIFVolumetricStatistics::FeatureListType &featureList)
virtual FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature)
Calculates the Cooccurence-Matrix based features for this class.
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
virtual FeatureNameListType GetFeatureNames()
Returns a list of the names of all features that are calculated from this class.
#define AccessByItk_2(mitkImage, itkImageTypeFunction, arg1, arg2)
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.