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
mitkTumorInvasionClassification.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 #include "mitkTumorInvasionClassification.h"
18 
19 #include "mitkCLUtil.h"
20 #include "mitkCollectionDilatation.h"
21 #include "mitkCollectionGrayOpening.h"
22 #include "mitkDataCollectionImageIterator.h"
23 #include "mitkDataCollectionUtilities.h"
24 #include "mitkIOUtil.h"
25 #include <fstream>
26 // To initialize random number generator
27 #include <time.h>
28 
29 static void EnsureDataImageInCollection(mitk::DataCollection::Pointer collection,
30  std::string origin,
31  std::string target)
32 {
33  typedef itk::Image<double, 3> FeatureImage;
34 
35  if (collection->HasElement(origin))
36  {
37  mitk::Image::Pointer originImage = dynamic_cast<mitk::Image *>(collection->GetMitkImage(origin).GetPointer());
38  FeatureImage::Pointer itkOriginImage = FeatureImage::New();
39  mitk::CastToItkImage(originImage, itkOriginImage);
40 
41  if (!collection->HasElement(target) && itkOriginImage.IsNotNull())
42  {
43  FeatureImage::Pointer image = FeatureImage::New();
44  image->SetRegions(itkOriginImage->GetLargestPossibleRegion());
45  image->SetSpacing(itkOriginImage->GetSpacing());
46  image->SetOrigin(itkOriginImage->GetOrigin());
47  image->SetDirection(itkOriginImage->GetDirection());
48  image->Allocate();
49 
50  collection->AddData(dynamic_cast<itk::DataObject *>(image.GetPointer()), target, "");
51  }
52  }
53  for (std::size_t i = 0; i < collection->Size(); ++i)
54  {
55  mitk::DataCollection *newCol = dynamic_cast<mitk::DataCollection *>(collection->GetData(i).GetPointer());
56  if (newCol != 0)
57  {
58  EnsureDataImageInCollection(newCol, origin, target);
59  }
60  }
61 }
62 
63 mitk::TumorInvasionClassification::TumorInvasionClassification()
64  : m_TargetID("TARGET"),
65  m_TumorID("GTV"),
66  m_MaskID("BRAINMASK"),
67  m_ResultID("RESULT"),
68  m_Randomize(false),
69  m_WeightSamples(false)
70 {
71 }
72 
73 void mitk::TumorInvasionClassification::SelectTrainingSamples(mitk::DataCollection *collection, unsigned int mode)
74 {
75  srand(time(NULL));
76  MITK_INFO << "LearnProgressionFeatures: Selecting training voxels.";
77  switch (mode)
78  {
79  case 0:
80  {
81  MITK_INFO << " Selection Mode " << mode << " use all tumor voxels, healthy: 50% vicinity / 50% far away";
82 
83  CollectionDilation::DilateBinaryByName(collection, m_TargetID, 2, 0, "EXCLUDE");
84  CollectionDilation::ErodeBinaryByName(collection, m_TargetID, 1, 0, "ERODE");
85  CollectionDilation::DilateBinaryByName(collection, m_TargetID, m_TargetDil2D, m_TargetDil3D, "TRAIN");
86  DataCollectionImageIterator<unsigned char, 3> gtvIter(collection, m_TumorID);
87  DataCollectionImageIterator<unsigned char, 3> brainMaskIter(collection, m_MaskID);
88  DataCollectionImageIterator<unsigned char, 3> targetIter(collection, m_TargetID);
89  DataCollectionImageIterator<unsigned char, 3> targetDil(collection, m_TargetID + "TRAIN");
90 
91  // Count Healthy/ Tumor voxels
92  // i.o. to decide how many of each are taken
93  unsigned int totalTumor = 0;
94  unsigned int totalHealthyClose = 0;
95  unsigned int totalHealthy = 0;
96 
97  while (!brainMaskIter.IsAtEnd())
98  {
99  if (brainMaskIter.GetVoxel() != 0)
100  {
101  if (targetIter.GetVoxel() == 1 && gtvIter.GetVoxel() == 0)
102  ++totalTumor;
103 
104  if (targetIter.GetVoxel() == 0 && targetDil.GetVoxel() == 1 && gtvIter.GetVoxel() == 0)
105  ++totalHealthyClose;
106 
107  if (targetIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 0 && targetDil.GetVoxel() == 0)
108  ++totalHealthy; // healthy but not close
109  }
110  ++brainMaskIter;
111  ++targetIter;
112  ++targetDil;
113  ++gtvIter;
114  }
115  brainMaskIter.ToBegin();
116  targetIter.ToBegin();
117  targetDil.ToBegin();
118  gtvIter.ToBegin();
119  // Total of healthy samples that is to be collected
120  unsigned int targetHealthy = totalTumor * m_ClassRatio;
121  // Determines which portion of the healthy samples is drawn from the immediate vicinity of the newly grown tumor
122  ScalarType ratioClose = .5;
123 
124  // Compute probabilities thresholds for choosing a close healthy voxel / any healthy voxel
125  ScalarType thHealthyClose = std::min(1.0, (targetHealthy * ratioClose) / totalHealthyClose);
126  ScalarType thHealthyAny = std::min(1.0, (targetHealthy * (1.0 - ratioClose)) / totalHealthy);
127  // Some stats
128  {
129  MITK_INFO << "Total # Tumor Voxels " << totalTumor;
130  MITK_INFO << "Total # Healthy Voxels" << totalHealthy;
131  MITK_INFO << "Target Ratio " << m_ClassRatio;
132  MITK_INFO << "Ratio for healthy close: " << thHealthyClose;
133  MITK_INFO << "Ratio for healthy any other: " << thHealthyAny;
134  }
135 
136  // DEBUG count occurances to compare
137  unsigned int tumor = 0;
138  unsigned int healthy = 0;
139  while (!brainMaskIter.IsAtEnd())
140  {
141  if (brainMaskIter.GetVoxel() != 0)
142  {
143  if (targetIter.GetVoxel() == 1)
144  { // choose tumor voxels for training
145  if (gtvIter.GetVoxel() == 0) // tumor always
146  {
147  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
148  ++tumor;
149  }
150  else
151  targetIter.SetVoxel(0);
152  }
153  else
154  { // choose healty tissue voxels for training
155  ScalarType rndVal = (float)rand() / (float)(RAND_MAX); //(0..1)
156  if (gtvIter.GetVoxel() == 0 && ((targetDil.GetVoxel() == 1 && rndVal <= thHealthyClose)))
157  {
158  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
159  ++healthy;
160  }
161  else if (((targetDil.GetVoxel() == 0 && rndVal <= thHealthyAny)))
162  {
163  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
164  ++healthy;
165  }
166  else
167  targetIter.SetVoxel(0);
168  }
169  }
170  else
171  targetIter.SetVoxel(0);
172 
173  ++brainMaskIter;
174  ++targetIter;
175  ++gtvIter;
176  ++targetDil;
177  }
178  MITK_INFO << "Training with Samples #Tumor " << tumor << " / healthy # " << healthy;
179  }
180  break;
181  case 2:
182  {
183  MITK_INFO << " Selection Mode " << mode << " Weighted with ratio";
184 
185  EnsureDataImageInCollection(collection, m_TumorID, "WEIGHTS");
186 
187  CollectionDilation::DilateBinaryByName(collection, m_TargetID, 1, 0, "EXCLUDE");
188  CollectionDilation::ErodeBinaryByName(collection, m_TargetID, 1, 0, "ERODE");
189  CollectionDilation::DilateBinaryByName(collection, m_TargetID + "EXCLUDE", m_TargetDil2D, m_TargetDil3D, "TRAIN");
190  DataCollectionImageIterator<unsigned char, 3> gtvIter(collection, m_TumorID);
191  DataCollectionImageIterator<unsigned char, 3> brainMaskIter(collection, m_MaskID);
192  DataCollectionImageIterator<unsigned char, 3> targetIter(collection, m_TargetID);
193  DataCollectionImageIterator<unsigned char, 3> targetDil(collection, m_TargetID + "EXCLUDETRAIN");
194  DataCollectionImageIterator<unsigned char, 3> excludeTumorIter(collection, m_TargetID + "ERODE");
195  DataCollectionImageIterator<unsigned char, 3> excludeHealthyIter(collection, m_TargetID + "EXCLUDE");
196 
197  DataCollectionImageIterator<double, 3> weightsIter(collection, "WEIGHTS");
198 
199  // Count Healthy/ Tumor voxels
200  // i.o. to decide how many of each are taken
201  double totalTumor = 0;
202  double totalHealthyClose = 0;
203  double totalHealthyNonClose = 0;
204  double totalHealthy = 0;
205 
206  while (!brainMaskIter.IsAtEnd())
207  {
208  if (brainMaskIter.GetVoxel() != 0)
209  {
210  if (targetIter.GetVoxel() == 1 && gtvIter.GetVoxel() == 0 && excludeTumorIter.GetVoxel() == 1)
211  ++totalTumor;
212 
213  if (excludeHealthyIter.GetVoxel() == 0 && targetDil.GetVoxel() == 1 && gtvIter.GetVoxel() == 0)
214  ++totalHealthyClose;
215 
216  if (excludeHealthyIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 0 && targetDil.GetVoxel() == 0)
217  ++totalHealthyNonClose; // healthy but not close
218 
219  if (excludeHealthyIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 0 && targetIter.GetVoxel() == 0)
220  ++totalHealthy; // healthy
221  }
222  ++brainMaskIter;
223  ++targetIter;
224  ++targetDil;
225  ++gtvIter;
226  ++excludeHealthyIter;
227  ++excludeTumorIter;
228  }
229  brainMaskIter.ToBegin();
230  targetIter.ToBegin();
231  targetDil.ToBegin();
232  gtvIter.ToBegin();
233  excludeHealthyIter.ToBegin();
234  excludeTumorIter.ToBegin();
235 
236  // Compute probabilities thresholds for choosing a close healthy voxel / any healthy voxel
237  ScalarType wHealthyClose = 10000.0 / totalHealthyClose;
238  ScalarType wHealthyAny = 5000.0 / totalHealthyNonClose;
239  ScalarType wTumor = (m_ClassRatio * (10000.0 + 5000.0)) / totalTumor;
240 
241  // DEBUG count occurances to compare
242  double potentialClose = 0;
243  double selectedClose = 0;
244  double tumor = 0;
245  double healthy = 0;
246  while (!brainMaskIter.IsAtEnd())
247  {
248  weightsIter.SetVoxel(0);
249  if (brainMaskIter.GetVoxel() != 0)
250  {
251  if (targetIter.GetVoxel() == 1)
252  { // choose tumor voxels for training
253  if (gtvIter.GetVoxel() == 0 && excludeTumorIter.GetVoxel() == 1) // tumor always
254  {
255  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
256  weightsIter.SetVoxel(wTumor);
257  ++tumor;
258  }
259  else if (gtvIter.GetVoxel() == 0)
260  {
261  weightsIter.SetVoxel(0); //.1);
262  targetIter.SetVoxel(0); // 2);
263  }
264  }
265  else
266  { // choose healty tissue voxels for training
267  if (gtvIter.GetVoxel() == 0 && ((excludeHealthyIter.GetVoxel() == 0 && targetDil.GetVoxel() == 1)))
268  {
269  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
270  weightsIter.SetVoxel(wHealthyClose);
271  healthy += wHealthyClose;
272  }
273  else if (((targetDil.GetVoxel() == 0 && excludeHealthyIter.GetVoxel() == 0)))
274  {
275  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
276  weightsIter.SetVoxel(wHealthyAny);
277  healthy += wHealthyAny;
278  }
279  else if ((gtvIter.GetVoxel() == 0 && excludeHealthyIter.GetVoxel() == 1))
280  {
281  targetIter.SetVoxel(0); // brainMaskIter.GetVoxel()+targetIter.GetVoxel());
282  weightsIter.SetVoxel(0); //.1);
283  }
284  else
285  {
286  targetIter.SetVoxel(0);
287  weightsIter.SetVoxel(0);
288  }
289  }
290  }
291  else
292  {
293  targetIter.SetVoxel(0);
294  weightsIter.SetVoxel(0);
295  }
296 
297  if (gtvIter.GetVoxel() != 0)
298  targetIter.SetVoxel(0);
299 
300  ++brainMaskIter;
301  ++targetIter;
302  ++gtvIter;
303  ++targetDil;
304  ++excludeHealthyIter;
305  ++excludeTumorIter;
306  ++weightsIter;
307  }
308  MITK_INFO << "Training with Samples #Tumor " << tumor << " / healthy # " << healthy;
309  MITK_INFO << "Potential Close" << potentialClose;
310  MITK_INFO << "Selected Close" << selectedClose;
311  }
312  break;
313  default:
314  {
315  MITK_INFO << " Selection Mode " << mode
316  << " Exclude voxels in border regions, healthy: 50% vicinity / 50% far away";
317 
318  // weights
319  ScalarType tumorWeight = 1;
320  // ScalarType unsureRegion = .25;
321  ScalarType healthyTissue = 1;
322 
323  EnsureDataImageInCollection(collection, m_TumorID, "WEIGHTS");
324 
325  CollectionDilation::DilateBinaryByName(collection, m_TargetID, 1, 0, "EXCLUDE");
326  CollectionDilation::ErodeBinaryByName(collection, m_TargetID, 1, 0, "ERODE");
327  CollectionDilation::DilateBinaryByName(collection, m_TargetID + "EXCLUDE", m_TargetDil2D, m_TargetDil3D, "TRAIN");
328  DataCollectionImageIterator<unsigned char, 3> gtvIter(collection, m_TumorID);
329  DataCollectionImageIterator<unsigned char, 3> brainMaskIter(collection, m_MaskID);
330  DataCollectionImageIterator<unsigned char, 3> targetIter(collection, m_TargetID);
331  DataCollectionImageIterator<unsigned char, 3> targetDil(collection, m_TargetID + "EXCLUDETRAIN");
332  DataCollectionImageIterator<unsigned char, 3> excludeTumorIter(collection, m_TargetID + "ERODE");
333  DataCollectionImageIterator<unsigned char, 3> excludeHealthyIter(collection, m_TargetID + "EXCLUDE");
334 
335  DataCollectionImageIterator<double, 3> weightsIter(collection, "WEIGHTS");
336 
337  // Count Healthy/ Tumor voxels
338  // i.o. to decide how many of each are taken
339  double totalTumor = 0;
340  double totalHealthyClose = 0;
341  double totalHealthyNonClose = 0;
342  double totalHealthy = 0;
343 
344  while (!brainMaskIter.IsAtEnd())
345  {
346  if (brainMaskIter.GetVoxel() != 0)
347  {
348  if (targetIter.GetVoxel() == 1 && gtvIter.GetVoxel() == 0 && excludeTumorIter.GetVoxel() == 1)
349  ++totalTumor;
350 
351  if (excludeHealthyIter.GetVoxel() == 0 && targetDil.GetVoxel() == 1 && gtvIter.GetVoxel() == 0)
352  ++totalHealthyClose;
353 
354  if (excludeHealthyIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 0 && targetDil.GetVoxel() == 0)
355  ++totalHealthyNonClose; // healthy but not close
356 
357  if (excludeHealthyIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 0 && targetIter.GetVoxel() == 0)
358  ++totalHealthy; // healthy
359  }
360  ++brainMaskIter;
361  ++targetIter;
362  ++targetDil;
363  ++gtvIter;
364  ++excludeHealthyIter;
365  ++excludeTumorIter;
366  }
367  brainMaskIter.ToBegin();
368  targetIter.ToBegin();
369  targetDil.ToBegin();
370  gtvIter.ToBegin();
371  excludeHealthyIter.ToBegin();
372  excludeTumorIter.ToBegin();
373 
374  // Total of healthy samples that is to be collected
375  unsigned int targetHealthy = (tumorWeight / healthyTissue) * totalTumor * m_ClassRatio;
376  // Determines which portion of the healthy samples is drawn from the immediate vicinity of the newly grown tumor
377  ScalarType ratioClose = .5;
378 
379  // Compute probabilities thresholds for choosing a close healthy voxel / any healthy voxel
380  ScalarType thHealthyClose = std::min(1.0, ((double)targetHealthy * ratioClose) / totalHealthyClose);
381  ScalarType thHealthyAny = std::min(1.0, ((double)targetHealthy * (1.0 - ratioClose)) / totalHealthyNonClose);
382  // Some stats
383  {
384  MITK_INFO << "Total Tumor " << totalTumor;
385  MITK_INFO << "Total healthy " << totalHealthyNonClose;
386  MITK_INFO << "Total healthy close " << totalHealthyClose;
387  MITK_INFO << "Total healthy non-close " << totalHealthyNonClose;
388  MITK_INFO << "Target Ratio " << m_ClassRatio;
389 
390  MITK_INFO << "Target Healthy " << targetHealthy;
391  MITK_INFO << "Probabilty close " << thHealthyClose;
392  MITK_INFO << "Probabilty any other " << thHealthyAny;
393  }
394 
395  // DEBUG count occurances to compare
396  double potentialClose = 0;
397  double selectedClose = 0;
398  double tumor = 0;
399  double healthy = 0;
400  while (!brainMaskIter.IsAtEnd())
401  {
402  weightsIter.SetVoxel(0);
403  if (brainMaskIter.GetVoxel() != 0)
404  {
405  if (targetIter.GetVoxel() == 1)
406  { // choose tumor voxels for training
407  if (gtvIter.GetVoxel() == 0 && excludeTumorIter.GetVoxel() == 1) // tumor always
408  {
409  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
410  weightsIter.SetVoxel(tumorWeight);
411  ++tumor;
412  }
413  else if (gtvIter.GetVoxel() == 0)
414  {
415  weightsIter.SetVoxel(0); //.1);
416  targetIter.SetVoxel(0); // 2);
417  }
418  }
419  else
420  { // choose healty tissue voxels for training
421  if (gtvIter.GetVoxel() == 0 && ((excludeHealthyIter.GetVoxel() == 0 && targetDil.GetVoxel() == 1)))
422  {
423  if (!m_WeightSamples && ((float)rand() / (float)(RAND_MAX) < thHealthyClose))
424  {
425  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
426  weightsIter.SetVoxel(1);
427  ++healthy;
428  }
429  else
430  {
431  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
432  weightsIter.SetVoxel(thHealthyClose);
433  healthy += thHealthyClose;
434  }
435  }
436  else if (((targetDil.GetVoxel() == 0 && excludeHealthyIter.GetVoxel() == 0)))
437  {
438  if (!m_WeightSamples && ((float)rand() / (float)(RAND_MAX) < thHealthyAny))
439  {
440  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
441  weightsIter.SetVoxel(1);
442  ++healthy;
443  }
444  else
445  {
446  targetIter.SetVoxel(brainMaskIter.GetVoxel() + targetIter.GetVoxel());
447  weightsIter.SetVoxel(thHealthyAny);
448  healthy += thHealthyAny;
449  }
450  }
451  else if ((gtvIter.GetVoxel() == 0 && excludeHealthyIter.GetVoxel() == 1))
452  {
453  targetIter.SetVoxel(0); // brainMaskIter.GetVoxel()+targetIter.GetVoxel());
454  weightsIter.SetVoxel(0); //.1);
455  }
456  else
457  {
458  targetIter.SetVoxel(0);
459  weightsIter.SetVoxel(0);
460  }
461 
462  if (gtvIter.GetVoxel() == 0 && ((excludeHealthyIter.GetVoxel() == 0 && targetDil.GetVoxel() == 1)))
463  {
464  potentialClose++;
465 
466  if ((float)rand() / (float)(RAND_MAX) < thHealthyClose)
467  selectedClose++;
468  }
469  }
470  }
471  else
472  {
473  targetIter.SetVoxel(0);
474  weightsIter.SetVoxel(0);
475  }
476 
477  if (gtvIter.GetVoxel() != 0)
478  targetIter.SetVoxel(0);
479 
480  ++brainMaskIter;
481  ++targetIter;
482  ++gtvIter;
483  ++targetDil;
484  ++excludeHealthyIter;
485  ++excludeTumorIter;
486  ++weightsIter;
487  }
488  MITK_INFO << "Training with Samples #Tumor " << tumor << " / healthy # " << healthy;
489  MITK_INFO << "Potential Close" << potentialClose;
490  MITK_INFO << "Selected Close" << selectedClose;
491  }
492  break;
493  }
494 }
495 
496 void mitk::TumorInvasionClassification::PrepareResponseSamples(mitk::DataCollection *collection)
497 {
498  srand(time(NULL));
499  MITK_INFO << "PrepareResponseSamples: Selecting training voxels.";
500 
501  EnsureDataImageInCollection(collection, m_TumorID, "WEIGHTS");
502 
503  CollectionDilation::DilateBinaryByName(collection, m_TargetID, 10, 4, "TRAIN");
504  DataCollectionImageIterator<unsigned char, 3> gtvIter(collection, m_TumorID);
505  DataCollectionImageIterator<unsigned char, 3> brainMaskIter(collection, m_MaskID);
506  DataCollectionImageIterator<unsigned char, 3> targetIter(collection, m_TargetID);
507  DataCollectionImageIterator<unsigned char, 3> dilIter(collection, m_TargetID + "TRAIN");
508 
509  // Count Non-Involved/Responsive/Rezidiv Voxels
510 
511  double nonInvolved = 0;
512  double responsive = 0;
513  double rezidiv = 0;
514 
515  while (!brainMaskIter.IsAtEnd())
516  {
517  if (brainMaskIter.GetVoxel() != 0 && dilIter.GetVoxel() != 0)
518  {
519  if (targetIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 0)
520  {
521  targetIter.SetVoxel(1);
522  ++nonInvolved;
523  }
524 
525  if (targetIter.GetVoxel() == 0 && gtvIter.GetVoxel() == 1)
526  {
527  targetIter.SetVoxel(3);
528  ++responsive;
529  }
530 
531  if (targetIter.GetVoxel() == 1 && gtvIter.GetVoxel() == 0)
532  {
533  targetIter.SetVoxel(2);
534  ++rezidiv;
535  }
536  }
537  else
538  targetIter.SetVoxel(0);
539 
540  ++brainMaskIter;
541  ++targetIter;
542  ++gtvIter;
543  ++dilIter;
544  }
545  brainMaskIter.ToBegin();
546  targetIter.ToBegin();
547  gtvIter.ToBegin();
548  dilIter.ToBegin();
549 
550  // weight for groups
551 
552  double wNonInvolved = 10000.0 / nonInvolved;
553  double wResponsive = 10000.0 / responsive;
554  double wRezidiv = 10000.0 / rezidiv;
555 
556  std::cout << "Weights are " << wNonInvolved << "/ " << wResponsive << " / " << wRezidiv << std::endl;
557 
558  // Assign weight for each voxel
559  DataCollectionImageIterator<double, 3> weightsIter(collection, "WEIGHTS");
560 
561  while (!brainMaskIter.IsAtEnd())
562  {
563  if (brainMaskIter.GetVoxel() != 0 && dilIter.GetVoxel() != 0)
564  {
565  if (targetIter.GetVoxel() == 1)
566  weightsIter.SetVoxel(wNonInvolved);
567 
568  if (targetIter.GetVoxel() == 2)
569  weightsIter.SetVoxel(wRezidiv);
570 
571  if (targetIter.GetVoxel() == 3)
572  weightsIter.SetVoxel(wResponsive);
573  }
574 
575  ++dilIter;
576  ++brainMaskIter;
577  ++targetIter;
578  ++weightsIter;
579  }
580 }
581 
582 void mitk::TumorInvasionClassification::LearnProgressionFeatures(mitk::DataCollection *collection,
583  std::vector<std::string> modalitiesList,
584  size_t forestSize,
585  size_t treeDepth)
586 {
587  auto trainDataX = mitk::DCUtilities::DC3dDToMatrixXd(collection, modalitiesList, m_TargetID);
588  auto trainDataY = mitk::DCUtilities::DC3dDToMatrixXi(collection, m_TargetID, m_TargetID);
589 
590  if (treeDepth != 0)
591  m_Forest.SetMaximumTreeDepth(treeDepth);
592 
593  m_Forest.SetTreeCount(forestSize);
594 
595  if (m_WeightSamples)
596  {
597  auto trainDataW = mitk::DCUtilities::DC3dDToMatrixXd(collection, "WEIGHTS", m_TargetID);
598  m_Forest.SetPointWiseWeight(trainDataW);
599  m_Forest.UsePointWiseWeight(true);
600  }
601 
602  MITK_INFO << "Start Training:";
603  try
604  {
605  m_Forest.Train(trainDataX, trainDataY);
606  }
607  catch (std::exception &e)
608  {
609  MITK_INFO << "Exception while training forest: " << e.what();
610  }
611  std::cout << "Training finished" << std::endl;
612 }
613 
614 void mitk::TumorInvasionClassification::PredictInvasion(mitk::DataCollection *collection,
615  std::vector<std::string> modalitiesList)
616 {
617  if (collection != NULL)
618  {
619  MITK_INFO << "Evaluating Forest";
620 
621  auto testDataX = mitk::DCUtilities::DC3dDToMatrixXd(collection, modalitiesList, m_MaskID);
622  auto testDataNewY = m_Forest.Predict(testDataX);
623  mitk::DCUtilities::MatrixToDC3d(testDataNewY, collection, m_ResultID, m_MaskID);
624 
625  Eigen::MatrixXd Probs = m_Forest.GetPointWiseProbabilities();
626 
627  Eigen::MatrixXd prob0 = Probs.col(0);
628  Eigen::MatrixXd prob1 = Probs.col(1);
629 
630  mitk::DCUtilities::MatrixToDC3d(prob0, collection, "prob0", m_MaskID);
631  mitk::DCUtilities::MatrixToDC3d(prob1, collection, "prob1", m_MaskID);
632  if (Probs.cols() >= 3)
633  {
634  Eigen::MatrixXd prob2 = Probs.col(2);
635  mitk::DCUtilities::MatrixToDC3d(prob2, collection, "prob2", m_MaskID);
636  }
637  }
638  else
639  MITK_ERROR << "TumorInvasionClassification::PredictInvasion - provided collection is NULL.";
640 }
641 
642 void mitk::TumorInvasionClassification::SanitizeResults(mitk::DataCollection *collection, std::string resultID)
643 {
644  CollectionGrayOpening::PerformGrayOpening(collection, resultID, "OPEN");
645 }
646 
647 void mitk::TumorInvasionClassification::DescriptionToLogFile(mitk::DataCollection *collection, std::string outputFile)
648 {
649  std::ofstream fileStream;
650  fileStream.open(outputFile.c_str(), std::ios_base::app);
651 
652  fileStream << "Configuration:" << std::endl;
653  for (size_t patient = 0; patient < collection->Size(); ++patient)
654  {
655  DataCollection *dataPatient = dynamic_cast<DataCollection *>(collection->GetData(patient).GetPointer());
656  fileStream << "Name : " << dataPatient->GetName() << " with time steps:" << std::endl;
657  for (size_t timeStep = 0; timeStep < dataPatient->Size(); ++timeStep)
658  {
659  fileStream << dataPatient->IndexToName(timeStep) << "-" << std::endl;
660  DataCollection *dataTS = dynamic_cast<DataCollection *>(dataPatient->GetData(timeStep).GetPointer());
661  for (size_t ts = 0; ts < dataTS->Size(); ++ts)
662  fileStream << dataTS->IndexToName(ts) << "-" << std::endl;
663  }
664  }
665 }
666 
667 void mitk::TumorInvasionClassification::SaveRandomForest(std::string filename)
668 {
669  mitk::IOUtil::Save(&m_Forest, filename);
670 }
671 
672 void mitk::TumorInvasionClassification::LoadRandomForest(std::string /*filename*/)
673 {
674  // TBD
675  MITK_ERROR << "not yet implemented ...";
676 }
mitk::CollectionDilation::ErodeBinaryByName
static void ErodeBinaryByName(DataCollection *dataCollection, std::string name, unsigned int xy=5, unsigned int z=0, std::string suffix="ERODE")
ErodeBinaryByName - Erode all occurances of a modality within a mitk::DataCollection.
Definition: mitkCollectionDilatation.cpp:76
mitk::IOUtil::Save
static void Save(const mitk::BaseData *data, const std::string &path)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:824
FeatureImage
itk::Image< mitk::ScalarType, 3 > FeatureImage
Definition: HistogramAdaption.cpp:26
mitk::Pointer
itk::SmartPointer< Self > Pointer
Definition: mitkRenderingManager.h:389
MITK_INFO
#define MITK_INFO
Definition: mitkLogMacros.h:22
mitk::TumorInvasionClassification::SaveRandomForest
void SaveRandomForest(std::string filename)
SaveRandomForest - Saves a trained random forest.
Definition: mitkTumorInvasionClassification.cpp:667
MITK_ERROR
#define MITK_ERROR
Definition: mitkLogMacros.h:24
mitk::ScalarType
double ScalarType
Definition: mitkNumericConstants.h:24
mitk::DataCollection::Size
vcl_size_t Size() const
Size - number of data items in collection.
Definition: mitkDataCollection.cpp:116
mitk::DataCollection::GetName
std::string GetName() const
Definition: mitkDataCollection.cpp:61
mitk::TumorInvasionClassification::DescriptionToLogFile
void DescriptionToLogFile(DataCollection *collection, std::string outputFile)
DescriptionToLogFile - Append to log file: patient ids and time steps contained in this data collecti...
Definition: mitkTumorInvasionClassification.cpp:647
mitk::TumorInvasionClassification::SelectTrainingSamples
void SelectTrainingSamples(DataCollection *collection, unsigned int mode=0)
SelectTrainingSamples.
Definition: mitkTumorInvasionClassification.cpp:73
mitk::TumorInvasionClassification::PrepareResponseSamples
void PrepareResponseSamples(DataCollection *collection)
PrepareResponseSamples.
Definition: mitkTumorInvasionClassification.cpp:496
mitk::TumorInvasionClassification::PredictInvasion
void PredictInvasion(DataCollection *collection, std::vector< std::string > modalitiesList)
PredictGrowth - Classify voxels into remaining healthy / turning into tumor.
Definition: mitkTumorInvasionClassification.cpp:614
mitk::DataCollection::IndexToName
std::string IndexToName(vcl_size_t index) const
IndexToName - Get name from index.
Definition: mitkDataCollection.cpp:90
mitk::CollectionGrayOpening::PerformGrayOpening
static void PerformGrayOpening(mitk::DataCollection *dataCollection, std::string name, std::string suffix)
PerformGrayOpening - Opening operation on a specific modality type wihtin the DataCollection. Creates a new item.
Definition: mitkCollectionGrayOpening.cpp:26
mitk::DCUtilities::MatrixToDC3d
static void MatrixToDC3d(const Eigen::MatrixXd &matrix, mitk::DataCollection::Pointer dc, const std::vector< std::string > &names, std::string mask)
Definition: mitkDataCollectionUtilities.cpp:125
filename
static const std::string filename
Definition: mitkTinyXMLTest.cpp:28
mitk::DCUtilities::DC3dDToMatrixXi
static Eigen::MatrixXi DC3dDToMatrixXi(mitk::DataCollection::Pointer dc, std::string name, std::string mask)
Definition: mitkDataCollectionUtilities.cpp:80
mitk::Image
Image class for storing images.
Definition: mitkImage.h:76
mitk::DCUtilities::DC3dDToMatrixXd
static Eigen::MatrixXd DC3dDToMatrixXd(mitk::DataCollection::Pointer dc, std::string names, std::string mask)
Definition: mitkDataCollectionUtilities.cpp:36
min
static T min(T x, T y)
Definition: svm.cpp:67
mitk::TumorInvasionClassification::LearnProgressionFeatures
void LearnProgressionFeatures(DataCollection *collection, std::vector< std::string > modalitiesList, vcl_size_t forestSize=300, vcl_size_t treeDepth=10)
LearnProgressionFeatures.
Definition: mitkTumorInvasionClassification.cpp:582
mitk::CastToItkImage
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
EnsureDataImageInCollection
static void EnsureDataImageInCollection(mitk::DataCollection::Pointer collection, std::string origin, std::string target)
Definition: mitkTumorInvasionClassification.cpp:29
mitk::DataCollection::GetData
itk::DataObject::Pointer GetData(vcl_size_t index)
GetData Get original data by index.
Definition: mitkDataCollection.cpp:141
mitk::TumorInvasionClassification::SanitizeResults
void SanitizeResults(DataCollection *collection, std::string resultID="RESULT")
SanitizeResults - Performs an Opening Operation on tha data to remove isolated misclassifications.
Definition: mitkTumorInvasionClassification.cpp:642
mitk::TumorInvasionClassification::LoadRandomForest
void LoadRandomForest(std::string)
LoadRandomForest - loads a random forest to be used for classification.
Definition: mitkTumorInvasionClassification.cpp:672
mitk::CollectionDilation::DilateBinaryByName
static void DilateBinaryByName(DataCollection *dataCollection, std::string name, unsigned int xy=5, unsigned int z=0, std::string suffix="DILATE")
DilateBinaryByName - Dilate all occurances of a modality within a mitk::DataCollection.
Definition: mitkCollectionDilatation.cpp:27
mitk::New
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.