Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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.