Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkLRDensityEstimation.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
14 
17 
18 #include <itkImageRegionIterator.h>
19 
21 
22 #include <mitkIOUtil.h>
23 #include "itkImageFileWriter.h"
24 
25 static void EnsureDataImageInCollection(mitk::DataCollection::Pointer collection, std::string origin, std::string target);
26 //static std::size_t CountElementsInCollection(mitk::DataCollection::Pointer collection, std::string maskName);
27 
28 void
30 {
31  m_Collection = data;
32 }
33 
36 {
37  return m_Collection;
38 }
39 
40 void
42 {
43  m_TestMask = name;
44 }
45 
46 std::string
48 {
49  return m_TestMask;
50 }
51 
52 void
54 {
55  m_TrainMask = name;
56 }
57 
58 std::string
60 {
61  return m_TrainMask;
62 }
63 
64 void
66 {
67  m_WeightName = name;
68 }
69 
70 std::string
72 {
73  return m_WeightName;
74 }
75 
76 void
78 {
79  typedef itk::Image<unsigned char, 3> MaskImageType;
80  typedef itk::Image<double , 3> FeatureImageType;
81 
83  EnsureDataImageInCollection(m_Collection, m_TrainMask, m_WeightName);
84 
86  DataCollectionSingleImageIterator<unsigned char, 3> train(m_Collection, m_TrainMask);
88  DataCollectionSingleImageIterator<double, 3> weight (m_Collection, m_WeightName);
89  std::vector<DataCollectionSingleImageIterator<double, 3> > featureList;
90  for (unsigned int i = 0; i < m_Modalities.size(); ++i)
91  {
92  DataCollectionSingleImageIterator<double, 3> iter(m_Collection, m_Modalities[i]);
93  featureList.push_back(iter);
94  }
95 
96  // Do for each image...
97  // char buchstabe = 'A';
98  while (!train.IsAtEnd())
99  {
100  itk::ImageRegionIterator<MaskImageType> trainIter(train.GetImage(), train.GetImage()->GetLargestPossibleRegion() );
101  itk::ImageRegionIterator<MaskImageType> testIter(test.GetImage(), test.GetImage()->GetLargestPossibleRegion() );
102 
103  int trainElements = 0;
104  int testElements = 0;
105  while (!trainIter.IsAtEnd())
106  {
107  if (trainIter.Get() > 0)
108  {
109  ++trainElements;
110  }
111  if (testIter.Get() > 0)
112  {
113  ++testElements;
114  }
115  ++trainIter;
116  ++testIter;
117  }
118  trainIter.GoToBegin();
119  testIter.GoToBegin();
120 
121  std::vector<itk::ImageRegionIterator<FeatureImageType> > featureIter;
122  for (unsigned int i = 0; i < featureList.size(); ++i)
123  {
124  itk::ImageRegionIterator<FeatureImageType> iter(featureList[i].GetImage(), featureList[i].GetImage()->GetLargestPossibleRegion() );
125  featureIter.push_back(iter);
126  }
127 
128  vnl_vector<double> label(trainElements + testElements);
129  vnl_matrix<double> feature(trainElements + testElements, featureList.size() );
130  vnl_matrix<double> trainFeature(trainElements, featureList.size() );
131 
132  int index = 0;
133  int trainIndex = 0;
134  while (!trainIter.IsAtEnd())
135  {
136  if (trainIter.Get() > 0)
137  {
138  label(index) = 1;
139  for (unsigned int i = 0; i < featureIter.size(); ++i)
140  {
141  feature(index, i) = featureIter[i].Get();
142  trainFeature(trainIndex, i) = featureIter[i].Get();
143  }
144  ++index;
145  ++trainIndex;
146  }
147  if (testIter.Get() > 0)
148  {
149  label(index) = 0;
150  for (unsigned int i = 0; i < featureIter.size(); ++i)
151  {
152  feature(index, i) = featureIter[i].Get();
153  }
154  ++index;
155  }
156 
157  ++trainIter;
158  ++testIter;
159  for (unsigned int i = 0; i < featureIter.size(); ++i)
160  {
161  ++(featureIter[i]);
162  }
163  }
164 
165  mitk::GeneralizedLinearModel glm(feature, label);
166  vnl_vector<double> weightVector = glm.ExpMu(trainFeature);
167 
168  itk::ImageRegionIterator<FeatureImageType> weightIter(weight.GetImage(), weight.GetImage()->GetLargestPossibleRegion() );
169 
170  trainIter.GoToBegin();
171  index = 0;
172  //double ratio = (trainElements * 1.0) / (testElements * 1.0);
173  //ratio = 1.0 / weightVector.one_norm()*weightVector.size();
174  while(!trainIter.IsAtEnd())
175  {
176  if (trainIter.Get() > 0)
177  {
178  weightIter.Set(weightVector(index) );
179  ++index;
180  }
181  ++trainIter;
182  ++weightIter;
183  }
184 
185  //std::stringstream s;
186  // s<<"d:/tmp/img/datei"<<buchstabe<<".nrrd";
187  //std::string blub =s.str();
188  //MITK_INFO << blub;
189  //FeatureImageType::Pointer image = weight.GetImage();
190  //typedef itk::ImageFileWriter< FeatureImageType > WriterType;
191  //WriterType::Pointer writer = WriterType::New();
192  //writer->SetFileName(blub);
193  //writer->SetInput(image);
194  //writer->Update();
195  //++buchstabe;
196 
197  ++train;
198  ++test;
199  ++weight;
200  for (unsigned int i = 0; i < featureList.size(); ++i)
201  {
202  ++(featureList[i]);
203  }
204  }
205 }
206 
207 void
208  mitk::LRDensityEstimation::SetModalities(std::vector<std::string> modalities)
209 {
210  m_Modalities = modalities;
211 }
212 
213 std::vector<std::string>
215 {
216  return m_Modalities;
217 }
218 
219 //static std::size_t CountElementsInCollection(mitk::DataCollection::Pointer collection, std::string maskName)
220 //{
221 // std::size_t count = 0;
222 // mitk::DataCollectionImageIterator<unsigned char,3> iter(collection, maskName);
223 // while (!iter.IsAtEnd())
224 // {
225 // if (iter.GetVoxel() > 0)
226 // {
227 // ++count;
228 // }
229 // ++iter;
230 // }
231 // return count;
232 //}
233 
234 static void EnsureDataImageInCollection(mitk::DataCollection::Pointer collection, std::string origin, std::string target)
235 {
236  typedef itk::Image<double, 3> FeatureImage;
237  typedef itk::Image<unsigned char, 3> LabelImage;
238 
240  while (!iter.IsAtEnd())
241  {
242  ++iter;
243  }
244 
245  if (collection->HasElement(origin))
246  {
247  LabelImage::Pointer originImage = dynamic_cast<LabelImage*>(collection->GetData(origin).GetPointer());
248  if (!collection->HasElement(target) && originImage.IsNotNull())
249  {
250  MITK_INFO << "New image necessary";
251  FeatureImage::Pointer image = FeatureImage::New();
252  image->SetRegions(originImage->GetLargestPossibleRegion());
253  image->SetSpacing(originImage->GetSpacing());
254  image->SetOrigin(originImage->GetOrigin());
255  image->SetDirection(originImage->GetDirection());
256  image->Allocate();
257 
258  collection->AddData(dynamic_cast<itk::DataObject*>(image.GetPointer()),target,"");
259  }
260  }
261  for (std::size_t i = 0; i < collection->Size();++i)
262  {
263  mitk::DataCollection* newCol = dynamic_cast<mitk::DataCollection*>(collection->GetData(i).GetPointer());
264  if (newCol != nullptr)
265  {
266  EnsureDataImageInCollection(newCol, origin, target);
267  }
268  }
269 }
270 
271 void
273 {
275  EnsureDataImageInCollection(train, m_TrainMask, m_WeightName);
276 
278  DataCollectionImageIterator<unsigned char, 3> trainIter(train, m_TrainMask);
279  DataCollectionImageIterator<unsigned char, 3> testIter (test, m_TestMask);
280  DataCollectionImageIterator<double, 3> weightIter (train, m_WeightName);
281  std::vector<DataCollectionImageIterator<double, 3> > trainFeatureList;
282  std::vector<DataCollectionImageIterator<double, 3> > testFeatureList;
283  for (unsigned int i = 0; i < m_Modalities.size(); ++i)
284  {
285  DataCollectionImageIterator<double, 3> iter(train, m_Modalities[i]);
286  trainFeatureList.push_back(iter);
287  DataCollectionImageIterator<double, 3> iter2(train, m_Modalities[i]);
288  testFeatureList.push_back(iter2);
289  }
290 
291  int trainElements = 0;
292  int testElements = 0;
293  while (!trainIter.IsAtEnd())
294  {
295  if (trainIter.GetVoxel() > 0)
296  ++trainElements;
297  ++trainIter;
298  }
299  while (!testIter.IsAtEnd() )
300  {
301  if (testIter.GetVoxel() > 0)
302  ++testElements;
303  ++testIter;
304  }
305  trainIter.ToBegin();
306  testIter.ToBegin();
307 
308  vnl_vector<double> label(trainElements + testElements);
309  vnl_matrix<double> feature(trainElements + testElements, trainFeatureList.size() );
310  vnl_matrix<double> trainFeature(trainElements, trainFeatureList.size() );
311 
312  int trainIndex = 0;
313  int testIndex = 0;
314  while (!trainIter.IsAtEnd())
315  {
316  if (trainIter.GetVoxel() > 0)
317  {
318  label(trainIndex) = 1;
319  for (unsigned int i = 0; i < trainFeatureList.size(); ++i)
320  {
321  feature(trainIndex, i) = trainFeatureList[i].GetVoxel();
322  trainFeature(trainIndex, i) = trainFeatureList[i].GetVoxel();
323  }
324 
325  ++trainIndex;
326  }
327 
328  ++trainIter;
329  for (unsigned int i = 0; i < trainFeatureList.size(); ++i)
330  {
331  ++(trainFeatureList[i]);
332  }
333  }
334 
335  while (!testIter.IsAtEnd())
336  {
337  if (testIter.GetVoxel() > 0)
338  {
339  label(trainIndex) = 1;
340  for (unsigned int i = 0; i < testFeatureList.size(); ++i)
341  {
342  feature(trainIndex+testIndex, i) = testFeatureList[i].GetVoxel();
343  }
344  ++testIndex;
345  }
346 
347  ++testIter;
348  for (unsigned int i = 0; i < testFeatureList.size(); ++i)
349  {
350  ++(testFeatureList[i]);
351  }
352  }
353 
354  mitk::GeneralizedLinearModel glm(feature, label);
355  vnl_vector<double> weightVector = glm.ExpMu(trainFeature);
356 
357  trainIter.ToBegin();
358  int index = 0;
359  while (!trainIter.IsAtEnd())
360  {
361  if (trainIter.GetVoxel() > 0)
362  {
363  weightIter.SetVoxel(weightVector(index));
364  ++index;
365  }
366  ++trainIter;
367  ++weightIter;
368  }
369 }
370 
372 {
373  typedef itk::Image<unsigned char, 3> MaskImageType;
374  typedef itk::Image<double , 3> FeatureImageType;
375 
377  EnsureDataImageInCollection(m_Collection, m_TrainMask, m_WeightName);
378 
380  DataCollectionSingleImageIterator<unsigned char, 3> train(m_Collection, m_TrainMask);
382  DataCollectionSingleImageIterator<double, 3> weight (m_Collection, m_WeightName);
383  std::vector<DataCollectionSingleImageIterator<double, 3> > featureList;
384  for (unsigned int i = 0; i < m_Modalities.size(); ++i)
385  {
386  DataCollectionSingleImageIterator<double, 3> iter(m_Collection, m_Modalities[i]);
387  featureList.push_back(iter);
388  }
389 
390  while (!train.IsAtEnd())
391  {
392  itk::ImageRegionIterator<MaskImageType> trainIter(train.GetImage(), train.GetImage()->GetLargestPossibleRegion() );
393  itk::ImageRegionIterator<MaskImageType> testIter(test.GetImage(), test.GetImage()->GetLargestPossibleRegion() );
394 
395  int trainElements = 0;
396  int testElements = 0;
397  while (!trainIter.IsAtEnd())
398  {
399  if (trainIter.Get() > 0)
400  {
401  ++trainElements;
402  }
403  if (testIter.Get() > 0)
404  {
405  ++testElements;
406  }
407  ++trainIter;
408  ++testIter;
409  }
410  trainIter.GoToBegin();
411  testIter.GoToBegin();
412 
413  std::vector<itk::ImageRegionIterator<FeatureImageType> > featureIter;
414  for (unsigned int i = 0; i < featureList.size(); ++i)
415  {
416  itk::ImageRegionIterator<FeatureImageType> iter(featureList[i].GetImage(), featureList[i].GetImage()->GetLargestPossibleRegion() );
417  featureIter.push_back(iter);
418  }
419 
420  vnl_vector<double> label(trainElements + testElements);
421  vnl_matrix<double> feature(trainElements + testElements, featureList.size() );
422  vnl_matrix<double> trainFeature(trainElements, featureList.size() );
423 
424  int index = 0;
425  int trainIndex = 0;
426  while (!trainIter.IsAtEnd())
427  {
428  if (trainIter.Get() > 0)
429  {
430  label(index) = 1;
431  for (unsigned int i = 0; i < featureIter.size(); ++i)
432  {
433  feature(index, i) = featureIter[i].Get();
434  trainFeature(trainIndex, i) = featureIter[i].Get();
435  }
436  ++index;
437  ++trainIndex;
438  }
439  if (testIter.Get() > 0)
440  {
441  label(index) = 0;
442  for (unsigned int i = 0; i < featureIter.size(); ++i)
443  {
444  feature(index, i) = featureIter[i].Get();
445  }
446  ++index;
447  }
448 
449  ++trainIter;
450  ++testIter;
451  for (unsigned int i = 0; i < featureIter.size(); ++i)
452  {
453  ++(featureIter[i]);
454  }
455  }
456 
457  mitk::GeneralizedLinearModel glm(feature, label);
458  vnl_vector<double> weightVector = glm.Predict(trainFeature);
459 
460  itk::ImageRegionIterator<FeatureImageType> weightIter(weight.GetImage(), weight.GetImage()->GetLargestPossibleRegion() );
461 
462  trainIter.GoToBegin();
463  index = 0;
464 
465  while(!trainIter.IsAtEnd())
466  {
467  if (trainIter.Get() > 0)
468  {
469  weightIter.Set(weightVector(index) );
470  ++index;
471  }
472  ++trainIter;
473  ++weightIter;
474  }
475 
476  ++train;
477  ++test;
478  ++weight;
479  for (unsigned int i = 0; i < featureList.size(); ++i)
480  {
481  ++(featureList[i]);
482  }
483  }
484 }
void SetWeightName(std::string name)
#define MITK_INFO
Definition: mitkLogMacros.h:18
Follow Up Storage - Class to facilitate loading/accessing structured follow-up data.
Definition: testcase.h:28
double Predict(const vnl_vector< double > &c)
Predicts the value corresponding to the given vector.
void SetTrainMask(std::string name)
DataCollection::Pointer GetCollection()
Generalized Linear Model that allows linear models for non-gaussian data.
void SetTestMask(std::string name)
static void EnsureDataImageInCollection(mitk::DataCollection::Pointer collection, std::string origin, std::string target)
std::vector< std::string > GetModalities()
itk::Image< unsigned char, 3 > MaskImageType
Definition: CLBrainMask.cpp:32
mitk::Image::Pointer image
void SetModalities(std::vector< std::string > modalities)
void SetCollection(DataCollection::Pointer data)
vnl_vector< double > ExpMu(const vnl_matrix< double > &x)
Estimation of the exponential factor for a given function.
itk::SmartPointer< Self > Pointer
Definition: mitkBaseData.h:41
void WeightsForAll(mitk::DataCollection::Pointer train, mitk::DataCollection::Pointer test)