Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
CLVoxelFeatures.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 #ifndef mitkCLVoxeFeatures_cpp
17 #define mitkCLVoxeFeatures_cpp
18 
19 #include "time.h"
20 #include <sstream>
21 #include <fstream>
22 
23 #include <mitkIOUtil.h>
24 #include <mitkImageAccessByItk.h>
25 #include <mitkImageCast.h>
26 #include "mitkCommandLineParser.h"
27 
30 
31 #include "itkDiscreteGaussianImageFilter.h"
32 #include <itkLaplacianRecursiveGaussianImageFilter.h>
33 #include "itkHessianRecursiveGaussianImageFilter.h"
34 #include "itkUnaryFunctorImageFilter.h"
36 #include "vnl/algo/vnl_symmetric_eigensystem.h"
37 #include <itkSubtractImageFilter.h>
38 
39 static std::vector<double> splitDouble(std::string str, char delimiter) {
40  std::vector<double> internal;
41  std::stringstream ss(str); // Turn the string into a stream.
42  std::string tok;
43  double val;
44  while (std::getline(ss, tok, delimiter)) {
45  std::stringstream s2(tok);
46  s2 >> val;
47  internal.push_back(val);
48  }
49 
50  return internal;
51 }
52 
53 namespace Functor
54 {
55  template <class TInput, class TOutput>
56  class MatrixFirstEigenvalue
57  {
58  public:
59  MatrixFirstEigenvalue() {}
60  virtual ~MatrixFirstEigenvalue() {}
61 
62  int order;
63 
64  inline TOutput operator ()(const TInput& input)
65  {
66  double a,b,c;
67  if (input[0] < 0.01 && input[1] < 0.01 &&input[2] < 0.01 &&input[3] < 0.01 &&input[4] < 0.01 &&input[5] < 0.01)
68  return 0;
69  vnl_symmetric_eigensystem_compute_eigenvals(input[0], input[1],input[2],input[3],input[4],input[5],a,b,c);
70  switch (order)
71  {
72  case 0: return a;
73  case 1: return b;
74  case 2: return c;
75  default: return a;
76  }
77  }
78  bool operator !=(const MatrixFirstEigenvalue) const
79  {
80  return false;
81  }
82  bool operator ==(const MatrixFirstEigenvalue& other) const
83  {
84  return !(*this != other);
85  }
86  };
87 }
88 
89 template<typename TPixel, unsigned int VImageDimension>
90 void
91  GaussianFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, mitk::Image::Pointer &output)
92 {
93  typedef itk::Image<TPixel, VImageDimension> ImageType;
94  typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > GaussFilterType;
95 
96  typename GaussFilterType::Pointer gaussianFilter = GaussFilterType::New();
97  gaussianFilter->SetInput( itkImage );
98  gaussianFilter->SetVariance(variance);
99  gaussianFilter->Update();
100  mitk::CastToMitkImage(gaussianFilter->GetOutput(), output);
101 }
102 
103 template<typename TPixel, unsigned int VImageDimension>
104 void
105  DifferenceOfGaussFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, mitk::Image::Pointer &output)
106 {
107  typedef itk::Image<TPixel, VImageDimension> ImageType;
108  typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > GaussFilterType;
109  typedef itk::SubtractImageFilter<ImageType, ImageType, ImageType> SubFilterType;
110 
111  typename GaussFilterType::Pointer gaussianFilter1 = GaussFilterType::New();
112  gaussianFilter1->SetInput( itkImage );
113  gaussianFilter1->SetVariance(variance);
114  gaussianFilter1->Update();
115  typename GaussFilterType::Pointer gaussianFilter2 = GaussFilterType::New();
116  gaussianFilter2->SetInput( itkImage );
117  gaussianFilter2->SetVariance(variance*0.66*0.66);
118  gaussianFilter2->Update();
119  typename SubFilterType::Pointer subFilter = SubFilterType::New();
120  subFilter->SetInput1(gaussianFilter1->GetOutput());
121  subFilter->SetInput2(gaussianFilter2->GetOutput());
122  subFilter->Update();
123  mitk::CastToMitkImage(subFilter->GetOutput(), output);
124 }
125 
126 template<typename TPixel, unsigned int VImageDimension>
127 void
128  LaplacianOfGaussianFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, mitk::Image::Pointer &output)
129 {
130  typedef itk::Image<TPixel, VImageDimension> ImageType;
131  typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > GaussFilterType;
132  typedef itk::LaplacianRecursiveGaussianImageFilter<ImageType, ImageType> LaplacianFilter;
133 
134  typename GaussFilterType::Pointer gaussianFilter = GaussFilterType::New();
135  gaussianFilter->SetInput( itkImage );
136  gaussianFilter->SetVariance(variance);
137  gaussianFilter->Update();
138  typename LaplacianFilter::Pointer laplaceFilter = LaplacianFilter::New();
139  laplaceFilter->SetInput(gaussianFilter->GetOutput());
140  laplaceFilter->Update();
141  mitk::CastToMitkImage(laplaceFilter->GetOutput(), output);
142 }
143 
144 template<typename TPixel, unsigned int VImageDimension>
145 void
146  HessianOfGaussianFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, std::vector<mitk::Image::Pointer> &out)
147 {
148  typedef itk::Image<TPixel, VImageDimension> ImageType;
149  typedef itk::Image<double, VImageDimension> FloatImageType;
150  typedef itk::HessianRecursiveGaussianImageFilter <ImageType> HessianFilterType;
151  typedef typename HessianFilterType::OutputImageType VectorImageType;
152  typedef Functor::MatrixFirstEigenvalue<typename VectorImageType::PixelType, double> DeterminantFunctorType;
153  typedef itk::UnaryFunctorImageFilter<VectorImageType, FloatImageType, DeterminantFunctorType> DetFilterType;
154 
155  typename HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
156  hessianFilter->SetInput(itkImage);
157  hessianFilter->SetSigma(std::sqrt(variance));
158  for (int i = 0; i < VImageDimension; ++i)
159  {
160  typename DetFilterType::Pointer detFilter = DetFilterType::New();
161  detFilter->SetInput(hessianFilter->GetOutput());
162  detFilter->GetFunctor().order = i;
163  detFilter->Update();
164  mitk::CastToMitkImage(detFilter->GetOutput(), out[i]);
165  }
166 }
167 
168 template<typename TPixel, unsigned int VImageDimension>
169 void
170  LocalHistograms(itk::Image<TPixel, VImageDimension>* itkImage, std::vector<mitk::Image::Pointer> &out, double offset, double delta)
171 {
172  typedef itk::Image<TPixel, VImageDimension> ImageType;
173  typedef itk::Image<double, VImageDimension> FloatImageType;
174  typedef itk::MultiHistogramFilter <ImageType,ImageType> MultiHistogramType;
175 
177  filter->SetInput(itkImage);
178  filter->SetOffset(offset);
179  filter->SetDelta(delta);
180  filter->Update();
181  for (int i = 0; i < VImageDimension; ++i)
182  {
184  mitk::CastToMitkImage(filter->GetOutput(), img);
185  out.push_back(img);
186  }
187 }
188 
189 int main(int argc, char* argv[])
190 {
191  mitkCommandLineParser parser;
192  parser.setArgumentPrefix("--", "-");
193  // required params
194  parser.addArgument("image", "i", mitkCommandLineParser::InputImage, "Input Image", "Path to the input VTK polydata", us::Any(), false);
195  parser.addArgument("output", "o", mitkCommandLineParser::OutputFile, "Output text file", "Target file. The output statistic is appended to this file.", us::Any(), false);
196 
197  parser.addArgument("gaussian","g",mitkCommandLineParser::String, "Gaussian Filtering of the input images", "Gaussian Filter. Followed by the used variances seperated by ';' ",us::Any());
198  parser.addArgument("difference-of-gaussian","dog",mitkCommandLineParser::String, "Difference of Gaussian Filtering of the input images", "Difference of Gaussian Filter. Followed by the used variances seperated by ';' ",us::Any());
199  parser.addArgument("laplace-of-gauss","log",mitkCommandLineParser::String, "Laplacian of Gaussian Filtering", "Laplacian of Gaussian Filter. Followed by the used variances seperated by ';' ",us::Any());
200  parser.addArgument("hessian-of-gauss","hog",mitkCommandLineParser::String, "Hessian of Gaussian Filtering", "Hessian of Gaussian Filter. Followed by the used variances seperated by ';' ",us::Any());
201  parser.addArgument("local-histogram", "lh", mitkCommandLineParser::String, "Local Histograms", "Calculate the local histogram based feature. Specify Offset and Delta, for exampel -3;0.6 ", us::Any());
202  // Miniapp Infos
203  parser.setCategory("Classification Tools");
204  parser.setTitle("Global Image Feature calculator");
205  parser.setDescription("Calculates different global statistics for a given segmentation / image combination");
206  parser.setContributor("MBI");
207 
208  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
209 
210  if (parsedArgs.size()==0)
211  {
212  return EXIT_FAILURE;
213  }
214  if ( parsedArgs.count("help") || parsedArgs.count("h"))
215  {
216  return EXIT_SUCCESS;
217  }
218  bool useCooc = parsedArgs.count("cooccurence");
219 
220  mitk::Image::Pointer image = mitk::IOUtil::LoadImage(parsedArgs["image"].ToString());
221  std::string filename=parsedArgs["output"].ToString();
222 
224  // CAlculate Gaussian Features
226  MITK_INFO << "Check for Local Histogram...";
227  if (parsedArgs.count("local-histogram"))
228  {
229  std::vector<mitk::Image::Pointer> outs;
230  auto ranges = splitDouble(parsedArgs["local-histogram"].ToString(), ';');
231  if (ranges.size() < 2)
232  {
233  MITK_INFO << "Missing Delta and Offset for Local Histogram";
234  }
235  else
236  {
237  AccessByItk_3(image, LocalHistograms, outs, ranges[0], ranges[1]);
238  for (int i = 0; i < outs.size(); ++i)
239  {
240  std::string name = filename + "-lh" + us::any_value_to_string<int>(i)+".nii.gz";
241  mitk::IOUtil::SaveImage(outs[i], name);
242  }
243  }
244  }
245 
247  // CAlculate Gaussian Features
249  MITK_INFO << "Check for Gaussian...";
250  if (parsedArgs.count("gaussian"))
251  {
252  MITK_INFO << "Calculate Gaussian... " << parsedArgs["gaussian"].ToString();
253  auto ranges = splitDouble(parsedArgs["gaussian"].ToString(),';');
254 
255  for (int i = 0; i < ranges.size(); ++i)
256  {
257  mitk::Image::Pointer output;
258  AccessByItk_2(image, GaussianFilter, ranges[i], output);
259  std::string name = filename + "-gaussian-" + us::any_value_to_string(ranges[i])+".nii.gz";
260  mitk::IOUtil::SaveImage(output, name);
261  }
262  }
263 
265  // CAlculate Difference of Gaussian Features
267  MITK_INFO << "Check for DoG...";
268  if (parsedArgs.count("difference-of-gaussian"))
269  {
270  MITK_INFO << "Calculate Difference of Gaussian... " << parsedArgs["difference-of-gaussian"].ToString();
271  auto ranges = splitDouble(parsedArgs["difference-of-gaussian"].ToString(),';');
272 
273  for (int i = 0; i < ranges.size(); ++i)
274  {
275  mitk::Image::Pointer output;
276  AccessByItk_2(image, DifferenceOfGaussFilter, ranges[i], output);
277  std::string name = filename + "-dog-" + us::any_value_to_string(ranges[i])+".nii.gz";
278  mitk::IOUtil::SaveImage(output, name);
279  }
280  }
281 
282  MITK_INFO << "Check for LoG...";
284  // CAlculate Laplacian Of Gauss Features
286  if (parsedArgs.count("laplace-of-gauss"))
287  {
288  MITK_INFO << "Calculate LoG... " << parsedArgs["laplace-of-gauss"].ToString();
289  auto ranges = splitDouble(parsedArgs["laplace-of-gauss"].ToString(),';');
290 
291  for (int i = 0; i < ranges.size(); ++i)
292  {
293  mitk::Image::Pointer output;
294  AccessByItk_2(image, LaplacianOfGaussianFilter, ranges[i], output);
295  std::string name = filename + "-log-" + us::any_value_to_string(ranges[i])+".nii.gz";
296  mitk::IOUtil::SaveImage(output, name);
297  }
298  }
299 
300  MITK_INFO << "Check for HoG...";
302  // CAlculate Hessian Of Gauss Features
304  if (parsedArgs.count("hessian-of-gauss"))
305  {
306  MITK_INFO << "Calculate HoG... " << parsedArgs["hessian-of-gauss"].ToString();
307  auto ranges = splitDouble(parsedArgs["hessian-of-gauss"].ToString(),';');
308 
309  for (int i = 0; i < ranges.size(); ++i)
310  {
311  std::vector<mitk::Image::Pointer> outs;
312  outs.push_back(mitk::Image::New());
313  outs.push_back(mitk::Image::New());
314  outs.push_back(mitk::Image::New());
315  AccessByItk_2(image, HessianOfGaussianFilter, ranges[i], outs);
316  std::string name = filename + "-hog0-" + us::any_value_to_string(ranges[i])+".nii.gz";
317  mitk::IOUtil::SaveImage(outs[0], name);
318  name = filename + "-hog1-" + us::any_value_to_string(ranges[i])+".nii.gz";
319  mitk::IOUtil::SaveImage(outs[1], name);
320  name = filename + "-hog2-" + us::any_value_to_string(ranges[i])+".nii.gz";
321  mitk::IOUtil::SaveImage(outs[2], name);
322  }
323  }
324 
325  return 0;
326 }
327 
328 #endif
US_Core_EXPORT std::string any_value_to_string(const Any &any)
Definition: usAny.cpp:26
void LaplacianOfGaussianFilter(itk::Image< TPixel, VImageDimension > *itkImage, double variance, mitk::Image::Pointer &output)
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
itk::SmartPointer< Self > Pointer
MITKCORE_EXPORT bool operator!=(const InteractionEvent &a, const InteractionEvent &b)
int main(int argc, char *argv[])
#define MITK_INFO
Definition: mitkLogMacros.h:22
itk::Image< unsigned char, 3 > ImageType
void setContributor(std::string contributor)
itk::Image< double, 3 > FloatImageType
Definition: CLBrainMask.cpp:35
void HessianOfGaussianFilter(itk::Image< TPixel, VImageDimension > *itkImage, double variance, std::vector< mitk::Image::Pointer > &out)
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
MITKCORE_EXPORT bool operator==(const InteractionEvent &a, const InteractionEvent &b)
void GaussianFilter(itk::Image< TPixel, VImageDimension > *itkImage, double variance, mitk::Image::Pointer &output)
static Vector3D offset
static bool SaveImage(mitk::Image::Pointer image, const std::string &path)
SaveImage Convenience method to save an arbitrary mitkImage.
Definition: mitkIOUtil.cpp:870
itk::VectorImage< float, 3 > VectorImageType
void DifferenceOfGaussFilter(itk::Image< TPixel, VImageDimension > *itkImage, double variance, mitk::Image::Pointer &output)
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)
static const std::string filename
Definition: usAny.h:163
void LocalHistograms(itk::Image< TPixel, VImageDimension > *itkImage, std::vector< mitk::Image::Pointer > &out, double offset, double delta)
static std::vector< double > splitDouble(std::string str, char delimiter)
void setCategory(std::string category)
static Pointer New()
void setArgumentPrefix(const std::string &longPrefix, const std::string &shortPrefix)
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:78
#define AccessByItk_2(mitkImage, itkImageTypeFunction, arg1, arg2)
void setTitle(std::string title)
void setDescription(std::string description)
static mitk::Image::Pointer LoadImage(const std::string &path)
LoadImage Convenience method to load an arbitrary mitkImage.
Definition: mitkIOUtil.cpp:597
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.