Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
mitkPartialVolumeAnalysisClusteringCalculator.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 
17 
18 #include "mitkPartialVolumeAnalysisClusteringCalculator.h"
19 
20 #include "mitkImageAccessByItk.h"
21 #include "mitkImageToItk.h"
22 #include "itkScalarImageToHistogramGenerator.h"
23 #include "itkListSample.h"
24 
25 #define _USE_MATH_DEFINES
26 #include <math.h>
27 
28 #define PVA_PI M_PI
29 
30 namespace mitk
31 {
32 
33  PartialVolumeAnalysisClusteringCalculator::PartialVolumeAnalysisClusteringCalculator()
34  : m_MaxIt(100), m_StepsNumIntegration(100)
35  {
36  }
37 
38 
39  PartialVolumeAnalysisClusteringCalculator::~PartialVolumeAnalysisClusteringCalculator()
40  {
41  }
42 
43  PartialVolumeAnalysisClusteringCalculator::ParamsType*
44  PartialVolumeAnalysisClusteringCalculator::InitialGuess(HistType h) const
45  {
46  int s = h.xVals.size();
47  int s30 = 0.3 * s;
48  int s70 = 0.7 * s;
49  auto params = new ParamsType();
50  params->means[0] = h.xVals(s30);
51  params->means[1] = h.xVals(s70);
52  params->sigmas[0] = (h.xVals(s-1) - h.xVals(0)) / 5.0;
53  params->sigmas[1] = (h.xVals(s-1) - h.xVals(0)) / 5.0;
54  params->ps[0] = 0.4;
55  params->ps[1] = 0.4;
56 
57  return params;
58  }
59 
60  PartialVolumeAnalysisClusteringCalculator::ParamsType*
61  PartialVolumeAnalysisClusteringCalculator::Cluster(HistType h, ParamsType *initialGuess) const
62  {
63  auto params = new ParamsType();
64  params->Initialize(initialGuess);
65 
66  double ll = 9999999999999999999.9;
67  double oll;
68  double sml = (h.xVals[1] - h.xVals[0]) / 1000;
69 
70  int arraysz = h.xVals.size();
71 
72  for (unsigned int it = 0; it < m_MaxIt; it++)
73  {
74  // wie sehen basisfunktionen zu aktuellen params aus?
75  ClusterResultType curves = CalculateCurves(*params,h.xVals);
76 
77  // summe der basisfunktionen
78  for(int i=0; i<arraysz; i++)
79  {
80  if(curves.combiVals(i) == 0)
81  {
82  curves.combiVals(i) = 0.00000000000000000001;
83  }
84  }
85 
86  // alte log likelihood
87  oll = ll;
88  ll = 0;
89 
90  // neue likelihood
91  for(int i=0; i<arraysz; i++)
92  {
93  ll -= h.hVals(i) * log(curves.combiVals(i));
94  }
95 
96  if( it>2 && oll-ll < arraysz/2*1e-15 )
97  {
98  break;
99  }
100 
101  for(int j=0; j<2; j++)
102  {
103  VecType array, p(arraysz);
104  array = curves.vals[j];
105 
106  for(int i=0; i<arraysz; i++)
107  {
108  p(i) = h.hVals(i)*array(i)/curves.combiVals(i);
109  }
110 
111  params->ps[j] = 0;
112  params->means[j] = 0;
113  for(int i=0; i<arraysz; i++)
114  {
115  params->ps[j] += p(i);
116  params->means[j] += h.xVals(i)*p(i);
117  }
118  params->means[j] /= params->ps[j];
119 
120  VecType vr = h.xVals;
121  for(int i=0; i<arraysz; i++)
122  {
123  vr(i) -= params->means[j];
124  }
125 
126  params->sigmas[j] = 0;
127  for(int i=0; i<arraysz; i++)
128  {
129  params->sigmas[j] += vr(i)*vr(i)*p(i);
130  }
131  params->sigmas[j] /= params->ps[j];
132  params->sigmas[j] += sml;
133 
134  }
135 
136  double p3 = 0;
137  for(int i=0; i<arraysz; i++)
138  {
139  p3 += h.hVals(i)*curves.mixedVals[0](i)/curves.combiVals(i);
140  }
141  p3 += 1e-3;
142 
143  double sum = 0;
144  for(int j=0; j<2; j++)
145  {
146  params->ps[j] = params->ps[j] + 1e-3;
147  sum += params->ps[j];
148  }
149  sum += p3;
150 
151  for(int j=0; j<2; j++)
152  {
153  params->ps[j] = params->ps[j] / sum;
154  }
155  p3 /= sum;
156 
157  }
158 
159  return params;
160 
161  }
162 
163  void PartialVolumeAnalysisClusteringCalculator::Normalize(ParamsType params, ClusterResultType* curves) const
164  {
165  double sum1=0, sum2=0, sum3=0;
166  int arraysz = curves->vals[0].size();
167 
168  for(int i=0; i<arraysz; i++)
169  {
170  sum1 += curves->vals[0](i);
171  sum2 += curves->vals[1](i);
172  sum3 += curves->mixedVals[0](i);
173  }
174 
175  sum1 /= params.ps[0];
176  sum2 /= params.ps[1];
177  sum3 /= 1.0 -params.ps[0] -params.ps[1];
178 
179  for(int i=0; i<arraysz; i++)
180  {
181  curves->vals[0](i) /= sum1;
182  curves->vals[1](i) /= sum2;
183  curves->mixedVals[0](i) /= sum3;
184  }
185 
186  for(int i=0; i<arraysz; i++)
187  {
188  curves->combiVals(i) = curves->mixedVals[0](i) + curves->vals[0](i) + curves->vals[1](i);
189  }
190  }
191 
192  PartialVolumeAnalysisClusteringCalculator::ClusterResultType
193  PartialVolumeAnalysisClusteringCalculator::CalculateCurves(ParamsType params, VecType xVals) const
194  {
195 
196  int arraysz = xVals.size();
197  ClusterResultType result(arraysz);
198 
199  for( int j=0; j<2; j++)
200  {
201  for(int i=0; i<arraysz; i++)
202  {
203  double d = xVals(i)-params.means[j];
204  double amp = params.ps[j]/sqrt(2*PVA_PI*params.sigmas[j]);
205 
206  result.vals[j](i) = amp*exp(-0.5 * (d*d)/params.sigmas[j]);
207  }
208  }
209 
210  for(int i=0; i<arraysz; i++)
211  {
212  result.mixedVals[0](i) = 0;
213  }
214 
215  for(double t=0; t<=1; t = t + 1.0/(m_StepsNumIntegration-1.0))
216  {
217  for(int i=0; i<arraysz; i++)
218  {
219  double d = xVals(i)-(t*params.means[0]+(1-t)*params.means[1]);
220  double v = t*params.sigmas[0]+(1-t)*params.sigmas[1];
221  double p = 1.0 - params.ps[0] - params.ps[1];
222  double amp = (1.0/m_StepsNumIntegration) * p / sqrt(2.0*PVA_PI*v);
223 
224  result.mixedVals[0](i) = result.mixedVals[0](i) + amp*exp(-0.5 * (d*d)/v);
225 
226  // MITK_INFO << "d=" << d << "v=" <<v << "p=" << p << "amp=" << amp << "result=" << amp*exp(-0.5 * (d*d)/v) << std::endl;
227  }
228  // MITK_INFO << "t=" << t << std::endl;
229  // params.Print();
230  // result.Print();
231  }
232 
233  for(int i=0; i<arraysz; i++)
234  {
235  result.combiVals(i) = result.mixedVals[0](i) + result.vals[0](i) + result.vals[1](i);
236  }
237 
238  return result;
239  }
240 
241  PartialVolumeAnalysisClusteringCalculator::HelperStructPerformRGBClusteringRetval*
242  PartialVolumeAnalysisClusteringCalculator::PerformRGBClustering(mitk::Image::ConstPointer image, const MitkHistType *histogram) const
243  {
244 
245  auto rgbChannels = new HelperStructRGBChannels();
246 
247  HelperStructPerformClusteringRetval *resultr = PerformClustering(image, histogram, 2);
248  rgbChannels->r = resultr->clusteredImage;
249 
250  HelperStructPerformClusteringRetval *resultg = PerformClustering(image, histogram, 0, resultr);
251  rgbChannels->g = resultg->clusteredImage;
252 
253  HelperStructPerformClusteringRetval *resultb = PerformClustering(image, histogram, 1, resultr);
254  rgbChannels->b = resultb->clusteredImage;
255 
256  mitk::Image::Pointer outImage = mitk::Image::New();
257 
258  switch(rgbChannels->r->GetDimension())
259  {
260  case 2:
261  InternalGenerateRGB<2>(rgbChannels, outImage);
262  break;
263  case 3:
264  InternalGenerateRGB<3>(rgbChannels, outImage);
265  break;
266  case 4:
267  InternalGenerateRGB<4>(rgbChannels, outImage);
268  break;
269  default:
270  InternalGenerateRGB<3>(rgbChannels, outImage);
271  }
272 
273  auto retval
274  = new HelperStructPerformRGBClusteringRetval();
275  retval->rgbChannels = rgbChannels;
276  retval->rgb = outImage;
277  retval->params = resultr->params;
278  retval->result = resultr->result;
279  retval->hist = resultr->hist;
280 
281  delete resultr;
282  delete resultg;
283 
284  return retval;
285 
286  }
287 
288  template < unsigned int VImageDimension >
289  void PartialVolumeAnalysisClusteringCalculator::InternalGenerateRGB( HelperStructRGBChannels *rgbin, mitk::Image::Pointer retval ) const
290  {
291  typedef itk::Image< float, VImageDimension > ProbImageType;
292  typedef itk::Image< typename itk::RGBAPixel<unsigned char>, VImageDimension > RGBImageType;
293 
294  typedef mitk::ImageToItk<ProbImageType> CastFilterType;
295  typename CastFilterType::Pointer castFilter = CastFilterType::New();
296  castFilter->SetInput( rgbin->r );
297  castFilter->Update();
298 
299  typename ProbImageType::Pointer r = castFilter->GetOutput();
300 
301  castFilter = CastFilterType::New();
302  castFilter->SetInput( rgbin->g );
303  castFilter->Update();
304  typename ProbImageType::Pointer g = castFilter->GetOutput();
305 
306  typename RGBImageType::Pointer rgb = RGBImageType::New();
307  rgb->SetSpacing( g->GetSpacing() ); // Set the image spacing
308  rgb->SetOrigin( g->GetOrigin() ); // Set the image origin
309  rgb->SetDirection( g->GetDirection() ); // Set the image direction
310  rgb->SetRegions( g->GetLargestPossibleRegion() );
311  rgb->Allocate();
312 
313  itk::ImageRegionConstIterator<ProbImageType>
314  itr(r, r->GetLargestPossibleRegion());
315  itk::ImageRegionConstIterator<ProbImageType>
316  itg(g, g->GetLargestPossibleRegion());
317 
318  itk::ImageRegionIterator<RGBImageType>
319  itrgb(rgb, rgb->GetLargestPossibleRegion());
320 
321  itr.GoToBegin();
322  itg.GoToBegin();
323 
324  float maxr = 0;
325  float maxg = 0;
326 
327  while( !itr.IsAtEnd() )
328  {
329  typename ProbImageType::PixelType pr = itr.Get();
330  typename ProbImageType::PixelType pg = itg.Get();
331 
332  if(pr > maxr)
333  {
334  maxr = pr;
335  }
336 
337  if(pg > maxg)
338  {
339  maxg = pg;
340  }
341 
342  ++itr;
343  ++itg;
344  }
345 
346  itr.GoToBegin();
347  itg.GoToBegin();
348  itrgb.GoToBegin();
349 
350  while( !itr.IsAtEnd() )
351  {
352  typename ProbImageType::PixelType pr = itr.Get();
353  typename ProbImageType::PixelType pg = itg.Get();
354 
355  typename RGBImageType::PixelType prgb;
356 
357  float valr = (pr/maxr)*255.0f;
358  float valg = (pg/maxg)*255.0f;
359  float alpha = valr>valg ? valr : valg;
360  prgb.Set(valr, valg, 0.0f, alpha);
361 
362  itrgb.Set(prgb);
363 
364  ++itr;
365  ++itg;
366  ++itrgb;
367  }
368 
369  retval->InitializeByItk(rgb.GetPointer());
370  retval->SetVolume(rgb->GetBufferPointer());
371 
372  }
373 
374 
375  PartialVolumeAnalysisClusteringCalculator::HelperStructPerformClusteringRetval*
376  PartialVolumeAnalysisClusteringCalculator::PerformClustering(mitk::Image::ConstPointer image, const MitkHistType *histogram, int classIdent, HelperStructPerformClusteringRetval* precResult) const
377  {
378 
379  auto retval =
380  new HelperStructPerformClusteringRetval();
381 
382  if(precResult == nullptr)
383  {
384  retval->hist = new HistType();
385  retval->hist->InitByMitkHistogram(histogram);
386 
387  ParamsType params;
388  params.Initialize( Cluster(*(retval->hist)) );
389  ClusterResultType result = CalculateCurves(params,retval->hist->xVals);
390  Normalize(params, &result);
391 
392  retval->params = new ParamsType();
393  retval->params->Initialize(&params);
394  retval->result = new ClusterResultType(10);
395  retval->result->Initialize(&result);
396  }
397  else
398  {
399  retval->params = new ParamsType();
400  retval->params->Initialize(precResult->params);
401  retval->result = new ClusterResultType(10);
402  retval->result->Initialize(precResult->result);
403  }
404 
405  VecType totalProbs = retval->result->combiVals;
406  VecType pvProbs = retval->result->mixedVals[0];
407  VecType fiberProbs;
408  VecType nonFiberProbs;
409  VecType interestingProbs;
410  double p_fiber;
411  double p_nonFiber;
412  double p_interesting;
413 // if(retval->params->means[0]<retval->params->means[1])
414 // {
415  fiberProbs = retval->result->vals[1];
416  nonFiberProbs = retval->result->vals[0];
417  p_fiber = retval->params->ps[1];
418  p_nonFiber = retval->params->ps[0];
419 // }
420 // else
421 // {
422 // fiberProbs = retval->result->vals[0];
423 // nonFiberProbs = retval->result->vals[1];
424 // p_fiber = retval->params->ps[0];
425 // p_nonFiber = retval->params->ps[1];
426 // }
427 
428  switch(classIdent)
429  {
430  case 0:
431  interestingProbs = nonFiberProbs;
432  p_interesting = p_nonFiber;
433  break;
434  case 1:
435  interestingProbs = pvProbs;
436  p_interesting = 1 - p_fiber - p_nonFiber;
437  break;
438  case 2:
439  default:
440  interestingProbs = fiberProbs;
441  p_interesting = p_fiber;
442  break;
443  }
444 
445  double sum = histogram->GetTotalFrequency();
446 
447  // initialize two histograms for class and total probabilities
448  MitkHistType::MeasurementVectorType min(1);
449  MitkHistType::MeasurementVectorType max(1);
450  min.Fill(histogram->GetDimensionMins(0)[0]);
451  max.Fill(histogram->GetDimensionMaxs(0)[histogram->GetDimensionMaxs(0).size()-1]);
452 
453  MitkHistType::Pointer interestingHist = MitkHistType::New();
454  interestingHist->SetMeasurementVectorSize(1);
455  interestingHist->Initialize(histogram->GetSize(),min,max);
456  MitkHistType::Iterator newIt = interestingHist->Begin();
457  MitkHistType::Iterator newEnd = interestingHist->End();
458 
459  MitkHistType::Pointer totalHist = MitkHistType::New();
460  totalHist->SetMeasurementVectorSize(1);
461  totalHist->Initialize(histogram->GetSize(),min,max);
462  MitkHistType::Iterator totalIt = totalHist->Begin();
463 
464  int i=0;
465  while (newIt != newEnd)
466  {
467  newIt.SetFrequency(interestingProbs(i)*sum);
468  totalIt.SetFrequency(totalProbs(i)*sum);
469  ++newIt;
470  ++totalIt;
471  ++i;
472  }
473 
474  mitk::Image::Pointer outImage1 = mitk::Image::New();
475  mitk::Image::Pointer outImage2 = mitk::Image::New();
476 
477  HelperStructClusteringResults clusterResults;
478  clusterResults.interestingHist = interestingHist;
479  clusterResults.totalHist = totalHist;
480  clusterResults.p_interesting = p_interesting;
481 
482  AccessFixedDimensionByItk_3(
483  image.GetPointer(),
484  InternalGenerateProbabilityImage,
485  3,
486  clusterResults,
487  outImage1, outImage2);
488 
489  retval->clusteredImage = outImage1;
490  retval->displayImage = outImage2;
491 
492  return retval;
493 
494  }
495 
496  template < typename TPixel, unsigned int VImageDimension >
497  void PartialVolumeAnalysisClusteringCalculator::InternalGenerateProbabilityImage(
498  const itk::Image< TPixel, VImageDimension > *image,
499  const HelperStructClusteringResults clusterResults,
500  mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2 ) const
501  {
502 
503  typedef itk::Image< TPixel, VImageDimension > ImageType;
504  typedef itk::Image< itk::RGBAPixel<unsigned char>, VImageDimension > DisplayImageType;
505  typedef itk::Image< float, VImageDimension > ProbImageType;
506 
507  typename ProbImageType::Pointer probimage = ProbImageType::New();
508  probimage->SetSpacing( image->GetSpacing() ); // Set the image spacing
509  probimage->SetOrigin( image->GetOrigin() ); // Set the image origin
510  probimage->SetDirection( image->GetDirection() ); // Set the image direction
511  probimage->SetRegions( image->GetLargestPossibleRegion() );
512  probimage->Allocate();
513  probimage->FillBuffer(0);
514 
515  typename DisplayImageType::Pointer displayimage = DisplayImageType::New();
516  displayimage->SetSpacing( image->GetSpacing() ); // Set the image spacing
517  displayimage->SetOrigin( image->GetOrigin() ); // Set the image origin
518  displayimage->SetDirection( image->GetDirection() ); // Set the image direction
519  displayimage->SetRegions( image->GetLargestPossibleRegion() );
520  displayimage->Allocate();
521 
522  typename DisplayImageType::PixelType rgba;
523  rgba.Set(0.0f, 0.0f, 0.0f, 0.0f);
524  displayimage->FillBuffer(rgba);
525 
526  itk::ImageRegionConstIterator<ImageType>
527  itimage(image, image->GetLargestPossibleRegion());
528 
529  itk::ImageRegionIterator<ProbImageType>
530  itprob(probimage, probimage->GetLargestPossibleRegion());
531 
532  itk::ImageRegionIterator<DisplayImageType>
533  itdisp(displayimage, displayimage->GetLargestPossibleRegion());
534 
535  itimage.GoToBegin();
536  itprob.GoToBegin();
537 
538  MitkHistType::IndexType index(1);
539  float maxp = 0;
540  while( !itimage.IsAtEnd() )
541  {
542  if(itimage.Get())
543  {
544  MitkHistType::MeasurementVectorType meas(1);
545  meas.Fill(itimage.Get());
546  double aposteriori = 0;
547  bool success = clusterResults.interestingHist->GetIndex(meas, index );
548  if(success)
549  {
550  double aprioriProb = clusterResults.interestingHist->GetFrequency(index);
551  double intensityProb = clusterResults.totalHist->GetFrequency(index);
552  double p_interesting = clusterResults.p_interesting;
553  aposteriori = p_interesting * aprioriProb / intensityProb;
554  }
555  else
556  {
557  MITK_ERROR << "index not found in histogram";
558  }
559 
560  if(aposteriori > 0.0000000000000001)
561  {
562  itprob.Set( aposteriori );
563  maxp = aposteriori > maxp ? aposteriori : maxp;
564  }
565  else
566  {
567  itprob.Set(0.0f);
568  }
569  }
570 
571  ++itimage;
572  ++itprob;
573  }
574 
575  itprob.GoToBegin();
576  itdisp.GoToBegin();
577 
578  while( !itprob.IsAtEnd() )
579  {
580  if(itprob.Get())
581  {
582  typename DisplayImageType::PixelType rgba;
583  rgba.Set(255.0f, 0.0f, 0.0f, 255.0f*(itprob.Get()/maxp));
584  itdisp.Set( rgba );
585  }
586  ++itprob;
587  ++itdisp;
588  }
589 
590  outImage1->InitializeByItk(probimage.GetPointer());
591  outImage1->SetVolume(probimage->GetBufferPointer());
592 
593  outImage2->InitializeByItk(displayimage.GetPointer());
594  outImage2->SetVolume(displayimage->GetBufferPointer());
595 
596  }
597 
598 
599  double* PartialVolumeAnalysisClusteringCalculator::PerformQuantification(
600  mitk::Image::ConstPointer image, mitk::Image::Pointer clusteredImage, mitk::Image::Pointer mask) const
601  {
602 
603  auto retval = new double[2];
604 
605  AccessFixedDimensionByItk_3(
606  image.GetPointer(),
607  InternalQuantify,
608  3,
609  clusteredImage.GetPointer(),
610  retval, mask );
611 
612  return retval;
613 
614  }
615 
616  template < typename TPixel, unsigned int VImageDimension >
617  void PartialVolumeAnalysisClusteringCalculator::InternalQuantify(
618  const itk::Image< TPixel, VImageDimension > *image,
619  mitk::Image::Pointer clusteredImage, double* retval, mitk::Image::Pointer mask ) const
620  {
621  typedef itk::Image< TPixel, VImageDimension > ImageType;
622  typedef itk::Image< float, VImageDimension > ProbImageType;
623  typedef itk::Image< unsigned char, VImageDimension > MaskImageType;
624 
625  typedef mitk::ImageToItk<ProbImageType> CastFilterType;
626  typename CastFilterType::Pointer castFilter = CastFilterType::New();
627  castFilter->SetInput( clusteredImage );
628  castFilter->Update();
629  typename ProbImageType::Pointer clusterImage = castFilter->GetOutput();
630 
631  typename MaskImageType::Pointer itkmask = nullptr;
632  if(mask.IsNotNull())
633  {
634  typedef mitk::ImageToItk<MaskImageType> CastFilterType2;
635  typename CastFilterType2::Pointer castFilter2 = CastFilterType2::New();
636  castFilter2->SetInput( mask );
637  castFilter2->Update();
638  itkmask = castFilter2->GetOutput();
639  }
640  else
641  {
642  itkmask = MaskImageType::New();
643  itkmask->SetSpacing( clusterImage->GetSpacing() ); // Set the image spacing
644  itkmask->SetOrigin( clusterImage->GetOrigin() ); // Set the image origin
645  itkmask->SetDirection( clusterImage->GetDirection() ); // Set the image direction
646  itkmask->SetRegions( clusterImage->GetLargestPossibleRegion() );
647  itkmask->Allocate();
648  itkmask->FillBuffer(1);
649  }
650 
651  itk::ImageRegionConstIterator<ImageType>
652  itimage(image, image->GetLargestPossibleRegion());
653 
654  itk::ImageRegionConstIterator<ProbImageType>
655  itprob(clusterImage, clusterImage->GetLargestPossibleRegion());
656 
657  itk::ImageRegionConstIterator<MaskImageType>
658  itmask(itkmask, itkmask->GetLargestPossibleRegion());
659 
660  itimage.GoToBegin();
661  itprob.GoToBegin();
662  itmask.GoToBegin();
663 
664  double totalProb = 0;
665  double measurement = 0;
666  double error = 0;
667 
668  while( !itimage.IsAtEnd() && !itprob.IsAtEnd() && !itmask.IsAtEnd() )
669  {
670  double valImag = itimage.Get();
671  double valProb = itprob.Get();
672  double valMask = itmask.Get();
673 
674  typename ProbImageType::PixelType prop = valProb * valMask;
675 
676  totalProb += prop;
677  measurement += valImag * prop;
678  error += valImag * valImag * prop;
679 
680  ++itimage;
681  ++itprob;
682  ++itmask;
683  }
684 
685  measurement = measurement / totalProb;
686  error = error / totalProb;
687  retval[0] = measurement;
688  retval[1] = sqrt( error - measurement*measurement );
689 
690  }
691 
692 
693  mitk::Image::Pointer PartialVolumeAnalysisClusteringCalculator::CaculateAngularErrorImage(
694  mitk::Image::Pointer comp1, mitk::Image::Pointer comp2, mitk::Image::Pointer probImg) const
695  {
696 
697  // cast input images to itk
698  typedef itk::Image<float, 3> ImageType;
699  typedef mitk::ImageToItk<ImageType> CastType;
700  CastType::Pointer caster = CastType::New();
701  caster->SetInput(comp1);
702  caster->Update();
703  ImageType::Pointer comp1Image = caster->GetOutput();
704 
705  caster = CastType::New();
706  caster->SetInput(comp2);
707  caster->Update();
708  ImageType::Pointer comp2Image = caster->GetOutput();
709 
710  caster = CastType::New();
711  caster->SetInput(probImg);
712  caster->Update();
713  ImageType::Pointer probImage = caster->GetOutput();
714 
715  // figure out maximum probability for fiber class
716  float maxProb = 0;
717  itk::ImageRegionConstIterator<ImageType>
718  itprob(probImage, probImage->GetLargestPossibleRegion());
719  itprob.GoToBegin();
720  while( !itprob.IsAtEnd() )
721  {
722  maxProb = itprob.Get() > maxProb ? itprob.Get() : maxProb;
723  ++itprob;
724  }
725 
726  // generate a list sample of angles at positions
727  // where the fiber-prob is higher than .2*maxprob
728  typedef float MeasurementType;
729  const unsigned int MeasurementVectorLength = 2;
730  typedef itk::Vector< MeasurementType , MeasurementVectorLength >
731  MeasurementVectorType;
732  typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType;
733  ListSampleType::Pointer listSample = ListSampleType::New();
734  listSample->SetMeasurementVectorSize( MeasurementVectorLength );
735 
736  itk::ImageRegionIterator<ImageType>
737  it1(comp1Image, comp1Image->GetLargestPossibleRegion());
738  itk::ImageRegionIterator<ImageType>
739  it2(comp2Image, comp2Image->GetLargestPossibleRegion());
740 
741  it1.GoToBegin();
742  it2.GoToBegin();
743  itprob.GoToBegin();
744 
745  while( !itprob.IsAtEnd() )
746  {
747  if(itprob.Get() > 0.2 * maxProb)
748  {
749  MeasurementVectorType mv;
750  mv[0] = ( MeasurementType ) it1.Get();
751  mv[1] = ( MeasurementType ) it2.Get();
752  listSample->PushBack(mv);
753  }
754  ++it1;
755  ++it2;
756  ++itprob;
757  }
758 
759  // generate a histogram from the list sample
760  typedef float HistogramMeasurementType;
761  typedef itk::Statistics::Histogram< HistogramMeasurementType, itk::Statistics::DenseFrequencyContainer2 > HistogramType;
762  typedef itk::Statistics::SampleToHistogramFilter< ListSampleType, HistogramType > GeneratorType;
763  GeneratorType::Pointer generator = GeneratorType::New();
764 
765  GeneratorType::HistogramType::SizeType size(2);
766  size.Fill(30);
767  generator->SetHistogramSize( size );
768 
769  generator->SetInput( listSample );
770  generator->SetMarginalScale( 10.0 );
771 
772  generator->Update();
773 
774  // look for frequency mode in the histogram
775  GeneratorType::HistogramType::ConstPointer histogram = generator->GetOutput();
776  GeneratorType::HistogramType::ConstIterator iter = histogram->Begin();
777  float maxFreq = 0;
778  MeasurementVectorType maxValue;
779  maxValue.Fill(0);
780  while ( iter != histogram->End() )
781  {
782  if(iter.GetFrequency() > maxFreq)
783  {
784  maxFreq = iter.GetFrequency();
785  maxValue[0] = iter.GetMeasurementVector()[0];
786  maxValue[1] = iter.GetMeasurementVector()[1];
787  }
788  ++iter;
789  }
790 
791  // generate return image that contains the angular
792  // error of the voxels to the histogram max measurement
793  ImageType::Pointer returnImage = ImageType::New();
794  returnImage->SetSpacing( comp1Image->GetSpacing() ); // Set the image spacing
795  returnImage->SetOrigin( comp1Image->GetOrigin() ); // Set the image origin
796  returnImage->SetDirection( comp1Image->GetDirection() ); // Set the image direction
797  returnImage->SetRegions( comp1Image->GetLargestPossibleRegion() );
798  returnImage->Allocate();
799 
800  itk::ImageRegionConstIterator<ImageType>
801  cit1(comp1Image, comp1Image->GetLargestPossibleRegion());
802 
803  itk::ImageRegionConstIterator<ImageType>
804  cit2(comp2Image, comp2Image->GetLargestPossibleRegion());
805 
806  itk::ImageRegionIterator<ImageType>
807  itout(returnImage, returnImage->GetLargestPossibleRegion());
808 
809  cit1.GoToBegin();
810  cit2.GoToBegin();
811  itout.GoToBegin();
812 
813  vnl_vector<float> v(3);
814  v[0] = cos( maxValue[0] ) * sin( maxValue[1] );
815  v[1] = sin( maxValue[0] ) * sin( maxValue[1] );
816  v[2] = cos( maxValue[1] );
817 // MITK_INFO << "max vector: " << v;
818  while( !cit1.IsAtEnd() )
819  {
820  vnl_vector<float> v1(3);
821  v1[0] = cos( cit1.Get() ) * sin( cit2.Get() );
822  v1[1] = sin( cit1.Get() ) * sin( cit2.Get() );
823  v1[2] = cos( cit2.Get() );
824 
825  itout.Set(fabs(angle(v,v1)));
826 // MITK_INFO << "ang_error " << v1 << ": " << fabs(angle(v,v1));
827 
828  ++cit1;
829  ++cit2;
830  ++itout;
831  }
832 
833  mitk::Image::Pointer retval = mitk::Image::New();
834  retval->InitializeByItk(returnImage.GetPointer());
835  retval->SetVolume(returnImage->GetBufferPointer());
836  return retval;
837  }
838 
839  PartialVolumeAnalysisClusteringCalculator::HelperStructPerformRGBClusteringRetval*
840  PartialVolumeAnalysisClusteringCalculator::PerformRGBQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2) const
841  {
842 
843  auto rgbChannels = new HelperStructRGBChannels();
844 
845  HelperStructPerformClusteringRetval *resultr = PerformQuantiles(image, histogram, p2, 999999999.0 );
846  rgbChannels->r = resultr->clusteredImage;
847 
848  HelperStructPerformClusteringRetval *resultg = PerformQuantiles(image, histogram, -999999999.0, p1 );
849  rgbChannels->g = resultg->clusteredImage;
850 
851  HelperStructPerformClusteringRetval *resultb = PerformQuantiles(image, histogram, p1, p2 );
852  rgbChannels->b = resultb->clusteredImage;
853 
854  mitk::Image::Pointer outImage = mitk::Image::New();
855 
856  switch(rgbChannels->r->GetDimension())
857  {
858  case 2:
859  InternalGenerateRGB<2>(rgbChannels, outImage);
860  break;
861  case 3:
862  InternalGenerateRGB<3>(rgbChannels, outImage);
863  break;
864  case 4:
865  InternalGenerateRGB<4>(rgbChannels, outImage);
866  break;
867  default:
868  InternalGenerateRGB<3>(rgbChannels, outImage);
869  }
870 
871  auto retval
872  = new HelperStructPerformRGBClusteringRetval();
873  retval->rgbChannels = rgbChannels;
874  retval->rgb = outImage;
875  retval->params = resultr->params;
876  retval->result = resultr->result;
877  retval->hist = resultr->hist;
878 
879  delete resultr;
880  delete resultg;
881 
882  return retval;
883 
884  }
885 
886  PartialVolumeAnalysisClusteringCalculator::HelperStructPerformClusteringRetval*
887  PartialVolumeAnalysisClusteringCalculator::PerformQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2 ) const
888  {
889 
890  auto retval =
891  new HelperStructPerformClusteringRetval();
892 
893  retval->hist = new HistType();
894  retval->hist->InitByMitkHistogram(histogram);
895 
896  auto q = new double[2];
897  q[0] = histogram->Quantile(0, p1);
898  q[1] = histogram->Quantile(0, p2);
899 
900  mitk::Image::Pointer outImage1 = mitk::Image::New();
901  mitk::Image::Pointer outImage2 = mitk::Image::New();
902 
903  AccessFixedDimensionByItk_3(
904  image.GetPointer(),
905  InternalGenerateQuantileImage,
906  3, q,
907  outImage1, outImage2);
908 
909  retval->clusteredImage = outImage1;
910  retval->displayImage = outImage2;
911 
912  delete[] q;
913  return retval;
914 
915  }
916 
917  template < typename TPixel, unsigned int VImageDimension >
918  void PartialVolumeAnalysisClusteringCalculator::InternalGenerateQuantileImage(
919  const itk::Image< TPixel, VImageDimension > *image, double* q,
920  mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2 ) const
921  {
922 
923  typedef itk::Image< TPixel, VImageDimension > ImageType;
924  typedef itk::Image< itk::RGBAPixel<unsigned char>, VImageDimension > DisplayImageType;
925  typedef itk::Image< float, VImageDimension > ProbImageType;
926 
927  typename ProbImageType::Pointer probimage = ProbImageType::New();
928  probimage->SetSpacing( image->GetSpacing() ); // Set the image spacing
929  probimage->SetOrigin( image->GetOrigin() ); // Set the image origin
930  probimage->SetDirection( image->GetDirection() ); // Set the image direction
931  probimage->SetRegions( image->GetLargestPossibleRegion() );
932  probimage->Allocate();
933  probimage->FillBuffer(0);
934 
935  typename DisplayImageType::Pointer displayimage = DisplayImageType::New();
936  displayimage->SetSpacing( image->GetSpacing() ); // Set the image spacing
937  displayimage->SetOrigin( image->GetOrigin() ); // Set the image origin
938  displayimage->SetDirection( image->GetDirection() ); // Set the image direction
939  displayimage->SetRegions( image->GetLargestPossibleRegion() );
940  displayimage->Allocate();
941 
942  typename DisplayImageType::PixelType rgba;
943  rgba.Set(0.0f, 0.0f, 0.0f, 0.0f);
944  displayimage->FillBuffer(rgba);
945 
946  itk::ImageRegionConstIterator<ImageType>
947  itimage(image, image->GetLargestPossibleRegion());
948 
949  itk::ImageRegionIterator<ProbImageType>
950  itprob(probimage, probimage->GetLargestPossibleRegion());
951 
952  itk::ImageRegionIterator<DisplayImageType>
953  itdisp(displayimage, displayimage->GetLargestPossibleRegion());
954 
955  itimage.GoToBegin();
956  itprob.GoToBegin();
957 
958  while( !itimage.IsAtEnd() )
959  {
960  if(itimage.Get() > q[0] && itimage.Get() < q[1])
961  {
962  itprob.Set(1.0f);
963  }
964 
965  ++itimage;
966  ++itprob;
967  }
968 
969  itprob.GoToBegin();
970  itdisp.GoToBegin();
971 
972  while( !itprob.IsAtEnd() )
973  {
974  if(itprob.Get())
975  {
976  typename DisplayImageType::PixelType rgba;
977  rgba.Set(255.0f, 0.0f, 0.0f, 255.0f);
978  itdisp.Set( rgba );
979  }
980  ++itprob;
981  ++itdisp;
982  }
983 
984  outImage1->InitializeByItk(probimage.GetPointer());
985  outImage1->SetVolume(probimage->GetBufferPointer());
986 
987  outImage2->InitializeByItk(displayimage.GetPointer());
988  outImage2->SetVolume(displayimage->GetBufferPointer());
989 
990  }
991 
992 }
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
mitk::PartialVolumeAnalysisClusteringCalculator::Normalize
void Normalize(ParamsType params, ClusterResultType *curves) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:163
mitk::PartialVolumeAnalysisClusteringCalculator::PerformQuantiles
HelperStructPerformClusteringRetval * PerformQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:887
MITK_ERROR
#define MITK_ERROR
Definition: mitkLogMacros.h:24
mitk::PartialVolumeAnalysisClusteringCalculator::PerformClustering
HelperStructPerformClusteringRetval * PerformClustering(mitk::Image::ConstPointer image, const MitkHistType *histogram, int classIdent, HelperStructPerformClusteringRetval *precResult=nullptr) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:376
mitk
DataCollection - Class to facilitate loading/accessing structured data.
Definition: GeometryOverview.dox:1
AccessFixedDimensionByItk_3
#define AccessFixedDimensionByItk_3(mitkImage, itkImageTypeFunction, dimension, arg1, arg2, arg3)
Definition: mitkImageAccessByItk.h:592
mitk::ConstPointer
itk::SmartPointer< const Self > ConstPointer
Definition: mitkRenderingManager.h:389
mitk::ImageType
map::core::discrete::Elements< 3 >::InternalImageType ImageType
Definition: mitkMultiModalAffine_default.cpp:26
mitk::PartialVolumeAnalysisClusteringCalculator::InternalQuantify
void InternalQuantify(const itk::Image< TPixel, VImageDimension > *image, mitk::Image::Pointer clusteredImage, double *retval, mitk::Image::Pointer mask) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:617
mitk::PartialVolumeAnalysisClusteringCalculator::PerformRGBQuantiles
HelperStructPerformRGBClusteringRetval * PerformRGBQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:840
mitk::PartialVolumeAnalysisClusteringCalculator::PerformQuantification
double * PerformQuantification(mitk::Image::ConstPointer image, mitk::Image::Pointer clusteredImage, mitk::Image::Pointer mask=nullptr) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:599
mitk::PartialVolumeAnalysisClusteringCalculator::CaculateAngularErrorImage
mitk::Image::Pointer CaculateAngularErrorImage(mitk::Image::Pointer comp1, mitk::Image::Pointer comp2, mitk::Image::Pointer probImg) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:693
MaskImageType
itk::Image< unsigned char, 3 > MaskImageType
Definition: CLBrainMask.cpp:36
max
static T max(T x, T y)
Definition: svm.cpp:70
mitk::Image::New
static Pointer New()
min
static T min(T x, T y)
Definition: svm.cpp:67
mitk::PartialVolumeAnalysisClusteringCalculator::PerformRGBClustering
HelperStructPerformRGBClusteringRetval * PerformRGBClustering(mitk::Image::ConstPointer image, const MitkHistType *histogram) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:242
mitk::PartialVolumeAnalysisClusteringCalculator::Cluster
ParamsType * Cluster(const HistType h, ParamsType *initialGuess) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:61
mitk::PartialVolumeAnalysisClusteringCalculator::InternalGenerateRGB
void InternalGenerateRGB(HelperStructRGBChannels *rgb, mitk::Image::Pointer retval) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:289
PixelType
unsigned short PixelType
Definition: TumorProgressionMapping.cpp:55
mitk::PartialVolumeAnalysisClusteringCalculator::InternalGenerateQuantileImage
void InternalGenerateQuantileImage(const itk::Image< TPixel, VImageDimension > *image, double *q, mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:918
mitk::PartialVolumeAnalysisClusteringCalculator::InternalGenerateProbabilityImage
void InternalGenerateProbabilityImage(const itk::Image< TPixel, VImageDimension > *image, const HelperStructClusteringResults clusterResults, mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:497
mitk::PartialVolumeAnalysisClusteringCalculator::CalculateCurves
ClusterResultType CalculateCurves(ParamsType params, VecType xVals) const
Definition: mitkPartialVolumeAnalysisClusteringCalculator.cpp:193
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.