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