Medical Imaging Interaction Toolkit  2018.4.99-9a29ffc6
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 (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 #ifndef mitkCLVoxeFeatures_cpp
13 #define mitkCLVoxeFeatures_cpp
14 
15 #include "time.h"
16 #include <sstream>
17 #include <fstream>
18 
19 #include <mitkIOUtil.h>
20 #include <mitkImageAccessByItk.h>
21 #include <mitkImageCast.h>
22 #include "mitkCommandLineParser.h"
23 
24 #include "itkDiscreteGaussianImageFilter.h"
25 #include <itkLaplacianRecursiveGaussianImageFilter.h>
26 #include "itkHessianRecursiveGaussianImageFilter.h"
27 #include "itkUnaryFunctorImageFilter.h"
28 #include "vnl/algo/vnl_symmetric_eigensystem.h"
30 #include <itkSubtractImageFilter.h>
32 
33 static std::vector<double> splitDouble(std::string str, char delimiter) {
34  std::vector<double> internal;
35  std::stringstream ss(str); // Turn the string into a stream.
36  std::string tok;
37  double val;
38  while (std::getline(ss, tok, delimiter)) {
39  std::stringstream s2(tok);
40  s2 >> val;
41  internal.push_back(val);
42  }
43 
44  return internal;
45 }
46 
47 namespace Functor
48 {
49  template <class TInput, class TOutput>
50  class MatrixFirstEigenvalue
51  {
52  public:
53  MatrixFirstEigenvalue() {}
54  virtual ~MatrixFirstEigenvalue() {}
55 
56  int order;
57 
58  inline TOutput operator ()(const TInput& input)
59  {
60  double a,b,c;
61  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)
62  return 0;
63  vnl_symmetric_eigensystem_compute_eigenvals(input[0], input[1],input[2],input[3],input[4],input[5],a,b,c);
64  switch (order)
65  {
66  case 0: return a;
67  case 1: return b;
68  case 2: return c;
69  default: return a;
70  }
71  }
72  bool operator !=(const MatrixFirstEigenvalue) const
73  {
74  return false;
75  }
76  bool operator ==(const MatrixFirstEigenvalue& other) const
77  {
78  return !(*this != other);
79  }
80  };
81 }
82 
83 template<typename TPixel, unsigned int VImageDimension>
84 void
85  GaussianFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, mitk::Image::Pointer &output)
86 {
87  typedef itk::Image<TPixel, VImageDimension> ImageType;
88  typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > GaussFilterType;
89 
90  typename GaussFilterType::Pointer gaussianFilter = GaussFilterType::New();
91  gaussianFilter->SetInput( itkImage );
92  gaussianFilter->SetVariance(variance);
93  gaussianFilter->Update();
94  mitk::CastToMitkImage(gaussianFilter->GetOutput(), output);
95 }
96 
97 template<typename TPixel, unsigned int VImageDimension>
98 void
99  DifferenceOfGaussFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, mitk::Image::Pointer &output)
100 {
101  typedef itk::Image<TPixel, VImageDimension> ImageType;
102  typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > GaussFilterType;
103  typedef itk::SubtractImageFilter<ImageType, ImageType, ImageType> SubFilterType;
104 
105  typename GaussFilterType::Pointer gaussianFilter1 = GaussFilterType::New();
106  gaussianFilter1->SetInput( itkImage );
107  gaussianFilter1->SetVariance(variance);
108  gaussianFilter1->Update();
109  typename GaussFilterType::Pointer gaussianFilter2 = GaussFilterType::New();
110  gaussianFilter2->SetInput( itkImage );
111  gaussianFilter2->SetVariance(variance*0.66*0.66);
112  gaussianFilter2->Update();
113  typename SubFilterType::Pointer subFilter = SubFilterType::New();
114  subFilter->SetInput1(gaussianFilter1->GetOutput());
115  subFilter->SetInput2(gaussianFilter2->GetOutput());
116  subFilter->Update();
117  mitk::CastToMitkImage(subFilter->GetOutput(), output);
118 }
119 
120 template<typename TPixel, unsigned int VImageDimension>
121 void
122  LaplacianOfGaussianFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, mitk::Image::Pointer &output)
123 {
124  typedef itk::Image<TPixel, VImageDimension> ImageType;
125  typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > GaussFilterType;
126  typedef itk::LaplacianRecursiveGaussianImageFilter<ImageType, ImageType> LaplacianFilter;
127 
128  typename GaussFilterType::Pointer gaussianFilter = GaussFilterType::New();
129  gaussianFilter->SetInput( itkImage );
130  gaussianFilter->SetVariance(variance);
131  gaussianFilter->Update();
132  typename LaplacianFilter::Pointer laplaceFilter = LaplacianFilter::New();
133  laplaceFilter->SetInput(gaussianFilter->GetOutput());
134  laplaceFilter->Update();
135  mitk::CastToMitkImage(laplaceFilter->GetOutput(), output);
136 }
137 
138 template<typename TPixel, unsigned int VImageDimension>
139 void
140  HessianOfGaussianFilter(itk::Image<TPixel, VImageDimension>* itkImage, double variance, std::vector<mitk::Image::Pointer> &out)
141 {
142  typedef itk::Image<TPixel, VImageDimension> ImageType;
143  typedef itk::Image<double, VImageDimension> FloatImageType;
144  typedef itk::HessianRecursiveGaussianImageFilter <ImageType> HessianFilterType;
145  typedef typename HessianFilterType::OutputImageType VectorImageType;
146  typedef Functor::MatrixFirstEigenvalue<typename VectorImageType::PixelType, double> DeterminantFunctorType;
147  typedef itk::UnaryFunctorImageFilter<VectorImageType, FloatImageType, DeterminantFunctorType> DetFilterType;
148 
149  typename HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
150  hessianFilter->SetInput(itkImage);
151  hessianFilter->SetSigma(std::sqrt(variance));
152  for (unsigned int i = 0; i < VImageDimension; ++i)
153  {
154  typename DetFilterType::Pointer detFilter = DetFilterType::New();
155  detFilter->SetInput(hessianFilter->GetOutput());
156  detFilter->GetFunctor().order = i;
157  detFilter->Update();
158  mitk::CastToMitkImage(detFilter->GetOutput(), out[i]);
159  }
160 }
161 
162 template<typename TPixel, unsigned int VImageDimension>
163 void
164 LocalHistograms2(itk::Image<TPixel, VImageDimension>* itkImage, std::vector<mitk::Image::Pointer> &out, std::vector<double> params)
165 {
166  typedef itk::Image<TPixel, VImageDimension> ImageType;
167  typedef itk::MultiHistogramFilter <ImageType, ImageType> MultiHistogramType;
168 
169  double minimum = params[0];
170  double maximum = params[1];
171  int bins = std::round(params[2]);
172 
173  double offset = minimum;
174  double delta = (maximum - minimum) / bins;
175 
176  typename MultiHistogramType::Pointer filter = MultiHistogramType::New();
177  filter->SetInput(itkImage);
178  filter->SetOffset(offset);
179  filter->SetDelta(delta);
180  filter->SetBins(bins);
181  filter->Update();
182  for (int i = 0; i < bins; ++i)
183  {
185  mitk::CastToMitkImage(filter->GetOutput(i), img);
186  out.push_back(img);
187  }
188 }
189 
190 template<typename TPixel, unsigned int VImageDimension>
191 void
192  LocalHistograms(itk::Image<TPixel, VImageDimension>* itkImage, std::vector<mitk::Image::Pointer> &out, double offset, double delta)
193 {
194  typedef itk::Image<TPixel, VImageDimension> ImageType;
195  typedef itk::MultiHistogramFilter <ImageType,ImageType> MultiHistogramType;
196 
197  typename MultiHistogramType::Pointer filter = MultiHistogramType::New();
198  filter->SetInput(itkImage);
199  filter->SetOffset(offset);
200  filter->SetDelta(delta);
201  filter->Update();
202  for (int i = 0; i < 11; ++i)
203  {
205  mitk::CastToMitkImage(filter->GetOutput(i), img);
206  out.push_back(img);
207  }
208 }
209 
210 template<typename TPixel, unsigned int VImageDimension>
211 void
212 localStatistic(itk::Image<TPixel, VImageDimension>* itkImage, std::vector<mitk::Image::Pointer> &out, int size)
213 {
214  typedef itk::Image<TPixel, VImageDimension> ImageType;
215  typedef itk::LocalStatisticFilter <ImageType, ImageType> MultiHistogramType;
216 
217  typename MultiHistogramType::Pointer filter = MultiHistogramType::New();
218  filter->SetInput(itkImage);
219  filter->SetSize(size);
220  filter->Update();
221  for (int i = 0; i < 5; ++i)
222  {
224  mitk::CastToMitkImage(filter->GetOutput(i), img);
225  out.push_back(img);
226  }
227 }
228 
229 
230 int main(int argc, char* argv[])
231 {
232  mitkCommandLineParser parser;
233  parser.setArgumentPrefix("--", "-");
234  // required params
235  parser.addArgument("image", "i", mitkCommandLineParser::Image, "Input Image", "Path to the input VTK polydata", us::Any(), false, false, false, mitkCommandLineParser::Input);
236  parser.addArgument("output", "o", mitkCommandLineParser::File, "Output text file", "Target file. The output statistic is appended to this file.", us::Any(), false, false, false, mitkCommandLineParser::Output);
237  parser.addArgument("extension", "e", mitkCommandLineParser::File, "Extension", "File extension. Default is .nii.gz", us::Any(), true, false, false, mitkCommandLineParser::Output);
238 
239  parser.addArgument("gaussian","g",mitkCommandLineParser::String, "Gaussian Filtering of the input images", "Gaussian Filter. Followed by the used variances seperated by ';' ",us::Any());
240  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());
241  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());
242  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());
243  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());
244  parser.addArgument("local-histogram2", "lh2", mitkCommandLineParser::String, "Local Histograms", "Calculate the local histogram based feature. Specify Minimum;Maximum;Bins, for exampel -3;3;6 ", us::Any());
245  parser.addArgument("local-statistic", "ls", mitkCommandLineParser::String, "Local Histograms", "Calculate the local histogram based feature. Specify Offset and Delta, for exampel -3;0.6 ", us::Any());
246  // Miniapp Infos
247  parser.setCategory("Classification Tools");
248  parser.setTitle("Global Image Feature calculator");
249  parser.setDescription("Calculates different global statistics for a given segmentation / image combination");
250  parser.setContributor("German Cancer Research Center (DKFZ)");
251 
252  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
253 
254  if (parsedArgs.size()==0)
255  {
256  return EXIT_FAILURE;
257  }
258  if ( parsedArgs.count("help") || parsedArgs.count("h"))
259  {
260  return EXIT_SUCCESS;
261  }
262 
263  mitk::Image::Pointer image = mitk::IOUtil::Load<mitk::Image>(parsedArgs["image"].ToString());
264  std::string filename=parsedArgs["output"].ToString();
265 
266  std::string extension = ".nii.gz";
267  if (parsedArgs.count("extension"))
268  {
269  extension = parsedArgs["extension"].ToString();
270  }
271 
273  // CAlculate Local Histogram
275  MITK_INFO << "Check for Local Histogram...";
276  if (parsedArgs.count("local-histogram"))
277  {
278  std::vector<mitk::Image::Pointer> outs;
279  auto ranges = splitDouble(parsedArgs["local-histogram"].ToString(), ';');
280  if (ranges.size() < 2)
281  {
282  MITK_INFO << "Missing Delta and Offset for Local Histogram";
283  }
284  else
285  {
286  AccessByItk_3(image, LocalHistograms, outs, ranges[0], ranges[1]);
287  for (std::size_t i = 0; i < outs.size(); ++i)
288  {
289  std::string name = filename + "-lh" + us::any_value_to_string<int>(i)+extension;
290  mitk::IOUtil::Save(outs[i], name);
291  }
292  }
293  }
294 
296  // CAlculate Local Histogram 2
298  MITK_INFO << "Check for Local Histogram...";
299  if (parsedArgs.count("local-histogram2"))
300  {
301  std::vector<mitk::Image::Pointer> outs;
302  auto ranges = splitDouble(parsedArgs["local-histogram2"].ToString(), ';');
303  if (ranges.size() < 3)
304  {
305  MITK_INFO << "Missing Delta and Offset for Local Histogram";
306  }
307  else
308  {
309  AccessByItk_2(image, LocalHistograms2, outs, ranges);
310  for (std::size_t i = 0; i < outs.size(); ++i)
311  {
312  std::string name = filename + "-lh2" + us::any_value_to_string<int>(i)+extension;
313  mitk::IOUtil::Save(outs[i], name);
314  }
315  }
316  }
317 
319  // CAlculate Local Statistic
321  MITK_INFO << "Check for Local Histogram...";
322  if (parsedArgs.count("local-statistic"))
323  {
324  std::vector<mitk::Image::Pointer> outs;
325  auto ranges = splitDouble(parsedArgs["local-statistic"].ToString(), ';');
326  if (ranges.size() < 1)
327  {
328  MITK_INFO << "Missing Rage";
329  }
330  else
331  {
332  for (std::size_t j = 0; j < ranges.size(); ++j)
333  {
334  AccessByItk_2(image, localStatistic, outs, ranges[j]);
335  for (std::size_t i = 0; i < outs.size(); ++i)
336  {
337  std::string name = filename + "-lstat" + us::any_value_to_string<int>(ranges[j])+ "_" +us::any_value_to_string<int>(i)+extension;
338  mitk::IOUtil::Save(outs[i], name);
339  }
340  outs.clear();
341  }
342  }
343  }
344 
345 
347  // CAlculate Gaussian Features
349  MITK_INFO << "Check for Gaussian...";
350  if (parsedArgs.count("gaussian"))
351  {
352  MITK_INFO << "Calculate Gaussian... " << parsedArgs["gaussian"].ToString();
353  auto ranges = splitDouble(parsedArgs["gaussian"].ToString(),';');
354 
355  for (std::size_t i = 0; i < ranges.size(); ++i)
356  {
357  MITK_INFO << "Gaussian with sigma: " << ranges[i];
358  mitk::Image::Pointer output;
359  AccessByItk_2(image, GaussianFilter, ranges[i], output);
360  MITK_INFO << "Write output:";
361  std::string name = filename + "-gaussian-" + us::any_value_to_string(ranges[i]) + extension;
362  mitk::IOUtil::Save(output, name);
363  }
364  }
365 
367  // CAlculate Difference of Gaussian Features
369  MITK_INFO << "Check for DoG...";
370  if (parsedArgs.count("difference-of-gaussian"))
371  {
372  MITK_INFO << "Calculate Difference of Gaussian... " << parsedArgs["difference-of-gaussian"].ToString();
373  auto ranges = splitDouble(parsedArgs["difference-of-gaussian"].ToString(),';');
374 
375  for (std::size_t i = 0; i < ranges.size(); ++i)
376  {
377  mitk::Image::Pointer output;
378  AccessByItk_2(image, DifferenceOfGaussFilter, ranges[i], output);
379  std::string name = filename + "-dog-" + us::any_value_to_string(ranges[i]) + extension;
380  mitk::IOUtil::Save(output, name);
381  }
382  }
383 
384  MITK_INFO << "Check for LoG...";
386  // CAlculate Laplacian Of Gauss Features
388  if (parsedArgs.count("laplace-of-gauss"))
389  {
390  MITK_INFO << "Calculate LoG... " << parsedArgs["laplace-of-gauss"].ToString();
391  auto ranges = splitDouble(parsedArgs["laplace-of-gauss"].ToString(),';');
392 
393  for (std::size_t i = 0; i < ranges.size(); ++i)
394  {
395  mitk::Image::Pointer output;
396  AccessByItk_2(image, LaplacianOfGaussianFilter, ranges[i], output);
397  std::string name = filename + "-log-" + us::any_value_to_string(ranges[i]) + extension;
398  mitk::IOUtil::Save(output, name);
399  }
400  }
401 
402  MITK_INFO << "Check for HoG...";
404  // CAlculate Hessian Of Gauss Features
406  if (parsedArgs.count("hessian-of-gauss"))
407  {
408  MITK_INFO << "Calculate HoG... " << parsedArgs["hessian-of-gauss"].ToString();
409  auto ranges = splitDouble(parsedArgs["hessian-of-gauss"].ToString(),';');
410 
411  for (std::size_t i = 0; i < ranges.size(); ++i)
412  {
413  std::vector<mitk::Image::Pointer> outs;
414  outs.push_back(mitk::Image::New());
415  outs.push_back(mitk::Image::New());
416  outs.push_back(mitk::Image::New());
417  AccessByItk_2(image, HessianOfGaussianFilter, ranges[i], outs);
418  std::string name = filename + "-hog0-" + us::any_value_to_string(ranges[i]) + extension;
419  mitk::IOUtil::Save(outs[0], name);
420  name = filename + "-hog1-" + us::any_value_to_string(ranges[i]) + extension;
421  mitk::IOUtil::Save(outs[1], name);
422  name = filename + "-hog2-" + us::any_value_to_string(ranges[i]) + extension;
423  mitk::IOUtil::Save(outs[2], name);
424  }
425  }
426 
427  return 0;
428 }
429 
430 #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)
MITKCORE_EXPORT bool operator!=(const InteractionEvent &a, const InteractionEvent &b)
int main(int argc, char *argv[])
#define MITK_INFO
Definition: mitkLogMacros.h:18
void LocalHistograms2(itk::Image< TPixel, VImageDimension > *itkImage, std::vector< mitk::Image::Pointer > &out, std::vector< double > params)
itk::Image< unsigned char, 3 > ImageType
void setContributor(std::string contributor)
itk::Image< double, 3 > FloatImageType
Definition: CLBrainMask.cpp:31
void HessianOfGaussianFilter(itk::Image< TPixel, VImageDimension > *itkImage, double variance, std::vector< mitk::Image::Pointer > &out)
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)
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
void DifferenceOfGaussFilter(itk::Image< TPixel, VImageDimension > *itkImage, double variance, mitk::Image::Pointer &output)
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)
mitk::Image::Pointer image
void localStatistic(itk::Image< TPixel, VImageDimension > *itkImage, std::vector< mitk::Image::Pointer > &out, int size)
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:74
static void Save(const mitk::BaseData *data, const std::string &path, bool setPathProperty=false)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:774
#define AccessByItk_2(mitkImage, itkImageTypeFunction, arg1, arg2)
void setTitle(std::string title)
void setDescription(std::string description)