Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
ManualSegmentationEvaluation.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 
13 #include <mitkCoreObjectFactory.h>
14 #include "mitkImage.h"
15 #include <mitkLexicalCast.h>
16 #include <vnl/vnl_random.h>
17 #include "mitkCommandLineParser.h"
18 #include <mitkIOUtil.h>
19 #include <itksys/SystemTools.hxx>
20 #include <mitkITKImageImport.h>
21 #include <mitkImageCast.h>
22 #include <mitkProperties.h>
23 
24 // ITK
25 #include <itkImageRegionIterator.h>
26 
27 // MITK
28 #include <mitkIOUtil.h>
29 
30 // Classification
31 #include <mitkCLUtil.h>
33 
34 #include <QDir>
35 #include <QString>
36 #include <QStringList>
37 
38 #include <itkLabelOverlapMeasuresImageFilter.h>
39 #include <itkImageRegionIteratorWithIndex.h>
42 #include <itkNeighborhood.h>
46 
47 using namespace mitk;
48 
49 std::vector<mitk::Image::Pointer> m_FeatureImageVector;
50 
51 void ProcessFeatureImages(const mitk::Image::Pointer & raw_image, const mitk::Image::Pointer & brain_mask)
52 {
53  typedef itk::Image<double,3> DoubleImageType;
54  typedef itk::Image<short,3> ShortImageType;
55  typedef itk::ConstNeighborhoodIterator<DoubleImageType> NeighborhoodType; // Neighborhood iterator to access image
58 
59  m_FeatureImageVector.clear();
60 
61  // RAW
62  m_FeatureImageVector.push_back(raw_image);
63 
64  // GAUSS
65  mitk::Image::Pointer smoothed;
66  mitk::CLUtil::GaussianFilter(raw_image,smoothed,1);
67  m_FeatureImageVector.push_back(smoothed);
68 
69  // Calculate Probability maps (parameters used from literatur)
70  // CSF
72  mitk::CLUtil::ProbabilityMap(smoothed,13.9, 8.3,csf_prob);
73  m_FeatureImageVector.push_back(csf_prob);
74 
75  // Lesion
76 
78  mitk::CLUtil::ProbabilityMap(smoothed,59, 11.6,les_prob);
79  m_FeatureImageVector.push_back(les_prob);
80 
81  // Barin (GM/WM)
82 
84  mitk::CLUtil::ProbabilityMap(smoothed,32, 5.6,brain_prob);
85  m_FeatureImageVector.push_back(brain_prob);
86 
87  std::vector<unsigned int> FOS_sizes;
88  FOS_sizes.push_back(1);
89 
90  DoubleImageType::Pointer input;
91  ShortImageType::Pointer mask;
92  mitk::CastToItkImage(smoothed, input);
93  mitk::CastToItkImage(brain_mask, mask);
94 
95  for(unsigned int i = 0 ; i < FOS_sizes.size(); i++)
96  {
97  FOSFilerType::Pointer filter = FOSFilerType::New();
98  filter->SetNeighborhoodSize(FOS_sizes[i]);
99  filter->SetInput(input);
100  filter->SetMask(mask);
101 
102  filter->Update();
103  FOSFilerType::DataObjectPointerArray array = filter->GetOutputs();
104 
105  for( unsigned int i = 0; i < FunctorType::OutputCount; i++)
106  {
107  mitk::Image::Pointer featureimage;
108  mitk::CastToMitkImage(dynamic_cast<DoubleImageType *>(array[i].GetPointer()),featureimage);
109  m_FeatureImageVector.push_back(featureimage);
110  // AddImageAsDataNode(featureimage,FunctorType::GetFeatureName(i))->SetVisibility(show_nodes);
111  }
112  }
113 
114  {
116  filter->SetInput(input);
117  filter->SetImageMask(mask);
118  filter->SetSigma(3);
119  filter->Update();
120 
121  mitk::Image::Pointer o1,o2,o3;
122  mitk::CastToMitkImage(filter->GetOutput(0),o1);
123  mitk::CastToMitkImage(filter->GetOutput(1),o2);
124  mitk::CastToMitkImage(filter->GetOutput(2),o3);
125 
126  m_FeatureImageVector.push_back(o1);
127  m_FeatureImageVector.push_back(o2);
128  m_FeatureImageVector.push_back(o3);
129  }
130 
131  {
133  filter->SetInput(input);
134  filter->SetImageMask(mask);
135  filter->SetInnerScale(1.5);
136  filter->SetOuterScale(3);
137  filter->Update();
138 
139  mitk::Image::Pointer o1,o2,o3;
140  mitk::CastToMitkImage(filter->GetOutput(0),o1);
141  mitk::CastToMitkImage(filter->GetOutput(1),o2);
142  mitk::CastToMitkImage(filter->GetOutput(2),o3);
143 
144  m_FeatureImageVector.push_back(o1);
145  m_FeatureImageVector.push_back(o2);
146  m_FeatureImageVector.push_back(o3);
147  }
148 
149  {
151  filter->SetInput(input);
152  filter->SetImageMask(mask);
153  filter->Update();
154 
156  mitk::CastToMitkImage(filter->GetOutput(0),o1);
157 
158  m_FeatureImageVector.push_back(o1);
159  }
160 }
161 
162 std::vector<mitk::Point3D> PointSetToVector(const mitk::PointSet::Pointer & mps)
163 {
164  std::vector<mitk::Point3D> result;
165  for(int i = 0 ; i < mps->GetSize(); i++)
166  result.push_back(mps->GetPoint(i));
167  return result;
168 }
169 
170 int main(int argc, char* argv[])
171 {
172  mitkCommandLineParser parser;
173  parser.setArgumentPrefix("--", "-");
174 
175  // required params
176  parser.addArgument("inputdir", "i", mitkCommandLineParser::Directory, "Input Directory", "Contains input feature files.", us::Any(), false, false, false, mitkCommandLineParser::Input);
177  parser.addArgument("outputdir", "o", mitkCommandLineParser::Directory, "Output Directory", "Destination of output files.", us::Any(), false, false, false, mitkCommandLineParser::Output);
178  parser.addArgument("mitkprojectdata", "d", mitkCommandLineParser::File, "original class mask and raw image", "Orig. data.", us::Any(), false, false, false, mitkCommandLineParser::Input);
179  parser.addArgument("csfmps", "csf", mitkCommandLineParser::File, "CSF Pointset", ".", us::Any(), false, false, false, mitkCommandLineParser::Input);
180  parser.addArgument("lesmps", "les", mitkCommandLineParser::File, "LES Pointset", ".", us::Any(), false, false, false, mitkCommandLineParser::Input);
181  parser.addArgument("bramps", "bra", mitkCommandLineParser::File, "BRA Pointset", ".", us::Any(), false, false, false, mitkCommandLineParser::Input);
182  // parser.addArgument("points", "p", mitkCommandLineParser::Int, "Ensure that p points are selected", ".", us::Any(), false);
183 
184  // Miniapp Infos
185  parser.setCategory("Classification Tools");
186  parser.setTitle("Evaluationtool for Manual-Segmentation");
187  parser.setDescription("Uses Datacollection to calculate DICE scores for CSF LES BRA");
188  parser.setContributor("German Cancer Research Center (DKFZ)");
189 
190  // Params parsing
191  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
192  if (parsedArgs.size()==0)
193  return EXIT_FAILURE;
194 
195  std::string inputdir = us::any_cast<std::string>(parsedArgs["inputdir"]);
196  std::string outputdir = us::any_cast<std::string>(parsedArgs["outputdir"]);
197  std::string mitkprojectdata = us::any_cast<std::string>(parsedArgs["mitkprojectdata"]);
198  std::string csf_mps_name = us::any_cast<std::string>(parsedArgs["csfmps"]);
199  std::string les_mps_name = us::any_cast<std::string>(parsedArgs["lesmps"]);
200  std::string bra_mps_name = us::any_cast<std::string>(parsedArgs["bramps"]);
201 
202  mitk::Image::Pointer class_mask_sampled, raw_image, class_mask;
203  mitk::PointSet::Pointer CSF_mps, LES_mps, BRA_mps;
204 
205  // Load from mitk-project
206  auto so = mitk::IOUtil::Load(inputdir + "/" + mitkprojectdata);
207  std::map<unsigned int, unsigned int> map;
208  mitk::CLUtil::CountVoxel(dynamic_cast<mitk::Image *>(so[1].GetPointer()), map);
209  raw_image = map.size() <= 7 ? dynamic_cast<mitk::Image *>(so[0].GetPointer()) : dynamic_cast<mitk::Image *>(so[1].GetPointer());
210  class_mask = map.size() <= 7 ? dynamic_cast<mitk::Image *>(so[1].GetPointer()) : dynamic_cast<mitk::Image *>(so[0].GetPointer());
211 
212  CSF_mps = mitk::IOUtil::Load<mitk::PointSet>(inputdir + "/" + csf_mps_name);
213  LES_mps = mitk::IOUtil::Load<mitk::PointSet>(inputdir + "/" + les_mps_name);
214  BRA_mps = mitk::IOUtil::Load<mitk::PointSet>(inputdir + "/" + bra_mps_name);
215 
216  unsigned int num_points = CSF_mps->GetSize() + LES_mps->GetSize() + BRA_mps->GetSize();
217  MITK_INFO << "Found #" << num_points << " points over all classes.";
218 
219  ProcessFeatureImages(raw_image, class_mask);
220 
221  std::map<unsigned int, unsigned int> tmpMap;
222  tmpMap[0] = 0;
223  tmpMap[1] = 1;
224  tmpMap[2] = 1;
225  tmpMap[3] = 1;
226  tmpMap[4] = 2;
227  tmpMap[5] = 3;
228  tmpMap[6] = 3;
229  mitk::CLUtil::MergeLabels( class_mask, tmpMap);
230  class_mask_sampled = class_mask->Clone();
231  itk::Image<short,3>::Pointer itk_classmask_sampled;
232  mitk::CastToItkImage(class_mask_sampled,itk_classmask_sampled);
233  itk::ImageRegionIteratorWithIndex<itk::Image<short,3> >::IndexType index;
234  itk::ImageRegionIteratorWithIndex<itk::Image<short,3> > iit(itk_classmask_sampled,itk_classmask_sampled->GetLargestPossibleRegion());
235 
236  std::ofstream myfile;
237  myfile.open (inputdir + "/results_3.csv");
238 
239  Eigen::MatrixXd X_test;
240 
241  unsigned int count_test = 0;
242  mitk::CLUtil::CountVoxel(class_mask, count_test);
243  X_test = Eigen::MatrixXd(count_test, m_FeatureImageVector.size());
244 
245  unsigned int pos = 0;
246  for( const auto & image : m_FeatureImageVector)
247  {
248  X_test.col(pos) = mitk::CLUtil::Transform<double>(image,class_mask);
249  ++pos;
250  }
251 
252  unsigned int runs = 20;
253  for(unsigned int k = 0 ; k < runs; k++)
254  {
255  auto CSF_vec = PointSetToVector(CSF_mps);
256  auto LES_vec = PointSetToVector(LES_mps);
257  auto BRA_vec = PointSetToVector(BRA_mps);
258 
259  itk_classmask_sampled->FillBuffer(0);
260 
261  // initial draws
262  std::random_shuffle(CSF_vec.begin(), CSF_vec.end());
263  class_mask->GetGeometry()->WorldToIndex(CSF_vec.back(),index);
264  iit.SetIndex(index);
265  iit.Set(1);
266  CSF_vec.pop_back();
267 
268  std::random_shuffle(LES_vec.begin(), LES_vec.end());
269  class_mask->GetGeometry()->WorldToIndex(LES_vec.back(),index);
270  iit.SetIndex(index);
271  iit.Set(2);
272  LES_vec.pop_back();
273 
274  std::random_shuffle(BRA_vec.begin(), BRA_vec.end());
275  class_mask->GetGeometry()->WorldToIndex(BRA_vec.back(),index);
276  iit.SetIndex(index);
277  iit.Set(3);
278  BRA_vec.pop_back();
279 
280  std::stringstream ss;
281 
282  while(!(CSF_vec.empty() && LES_vec.empty() && BRA_vec.empty()))
283  {
284  mitk::CastToMitkImage(itk_classmask_sampled, class_mask_sampled);
285 
286  // Train forest
288  classifier->SetTreeCount(40);
289  classifier->SetSamplesPerTree(0.66);
290 
291  Eigen::MatrixXd X_train;
292 
293  unsigned int count_train = 0;
294  mitk::CLUtil::CountVoxel(class_mask_sampled, count_train);
295  X_train = Eigen::MatrixXd(count_train, m_FeatureImageVector.size() );
296 
297  unsigned int pos = 0;
298  for( const auto & image : m_FeatureImageVector)
299  {
300  X_train.col(pos) = mitk::CLUtil::Transform<double>(image,class_mask_sampled);
301  ++pos;
302  }
303 
304  Eigen::MatrixXi Y = mitk::CLUtil::Transform<int>(class_mask_sampled,class_mask_sampled);
305  classifier->Train(X_train,Y);
306 
307  Eigen::MatrixXi Y_test = classifier->Predict(X_test);
308 
309  mitk::Image::Pointer result_mask = mitk::CLUtil::Transform<int>(Y_test, class_mask);
310  itk::Image<short,3>::Pointer itk_result_mask, itk_class_mask;
311 
312  mitk::CastToItkImage(result_mask,itk_result_mask);
313  mitk::CastToItkImage(class_mask, itk_class_mask);
314 
315  itk::LabelOverlapMeasuresImageFilter<itk::Image<short,3> >::Pointer overlap_filter = itk::LabelOverlapMeasuresImageFilter<itk::Image<short,3> >::New();
316  overlap_filter->SetInput(0,itk_result_mask);
317  overlap_filter->SetInput(1,itk_class_mask);
318  overlap_filter->Update();
319 
320  MITK_INFO << "DICE (" << num_points - (CSF_vec.size() + LES_vec.size() + BRA_vec.size()) << "): " << overlap_filter->GetDiceCoefficient();
321  ss << overlap_filter->GetDiceCoefficient() <<",";
322 
323  // random class selection
324 
325  if(!CSF_vec.empty())
326  {
327  std::random_shuffle(CSF_vec.begin(), CSF_vec.end());
328  class_mask->GetGeometry()->WorldToIndex(CSF_vec.back(),index);
329  iit.SetIndex(index);
330  iit.Set(1);
331  CSF_vec.pop_back();
332  }
333 
334  if(!LES_vec.empty())
335  {
336  std::random_shuffle(LES_vec.begin(), LES_vec.end());
337  class_mask->GetGeometry()->WorldToIndex(LES_vec.back(),index);
338  iit.SetIndex(index);
339  iit.Set(2);
340  LES_vec.pop_back();
341  }
342 
343  if(!BRA_vec.empty())
344  {
345  std::random_shuffle(BRA_vec.begin(), BRA_vec.end());
346  class_mask->GetGeometry()->WorldToIndex(BRA_vec.back(),index);
347  iit.SetIndex(index);
348  iit.Set(3);
349  BRA_vec.pop_back();
350  }
351  }
352 
353  myfile << ss.str() << "\n";
354  myfile.flush();
355  }
356 
357  myfile.close();
358 
359  return EXIT_SUCCESS;
360 }
float k(1.0)
int main(int argc, char *argv[])
#define MITK_INFO
Definition: mitkLogMacros.h:18
std::vector< mitk::Image::Pointer > m_FeatureImageVector
std::vector< mitk::Point3D > PointSetToVector(const mitk::PointSet::Pointer &mps)
void setContributor(std::string contributor)
DataCollection - Class to facilitate loading/accessing structured data.
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
static void GaussianFilter(mitk::Image::Pointer image, mitk::Image::Pointer &smoothed, double sigma)
GaussianFilter.
Definition: mitkCLUtil.cpp:119
void addArgument(const std::string &longarg, const std::string &shortarg, Type type, const std::string &argLabel, const std::string &argHelp=std::string(), const us::Any &defaultValue=us::Any(), bool optional=true, bool ignoreRest=false, bool deprecated=false, mitkCommandLineParser::Channel channel=mitkCommandLineParser::Channel::None)
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
void ProcessFeatureImages(const mitk::Image::Pointer &raw_image, const mitk::Image::Pointer &brain_mask)
Image class for storing images.
Definition: mitkImage.h:72
Definition: usAny.h:163
static void ProbabilityMap(const mitk::Image::Pointer &sourceImage, double mean, double std_dev, mitk::Image::Pointer &resultImage)
ProbabilityMap.
Definition: mitkCLUtil.cpp:48
void setCategory(std::string category)
mitk::Image::Pointer image
static Pointer New()
MITKMATCHPOINTREGISTRATION_EXPORT ResultImageType::Pointer map(const InputImageType *input, const RegistrationType *registration, bool throwOnOutOfInputAreaError=false, const double &paddingValue=0, const ResultImageGeometryType *resultGeometry=nullptr, bool throwOnMappingError=true, const double &errorValue=0, mitk::ImageMappingInterpolator::Type interpolatorType=mitk::ImageMappingInterpolator::Linear)
void setArgumentPrefix(const std::string &longPrefix, const std::string &shortPrefix)
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
static void MergeLabels(mitk::Image::Pointer &img, const std::map< unsigned int, unsigned int > &map)
MergeLabels.
Definition: mitkCLUtil.cpp:83
mitk::Image::Pointer mask
void setTitle(std::string title)
static void CountVoxel(mitk::Image::Pointer image, std::map< unsigned int, unsigned int > &map)
CountVoxel.
Definition: mitkCLUtil.cpp:88
void setDescription(std::string description)
static DataStorage::SetOfObjects::Pointer Load(const std::string &path, DataStorage &storage, const ReaderOptionsFunctorBase *optionsCallback=nullptr)
Load a file into the given DataStorage.
Definition: mitkIOUtil.cpp:489