Medical Imaging Interaction Toolkit  2018.4.99-b585543d
Medical Imaging Interaction Toolkit
NativeHeadCTSegmentation.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 "mitkCommandLineParser.h"
17 #include <mitkIOUtil.h>
18 #include <itksys/SystemTools.hxx>
19 #include <mitkITKImageImport.h>
20 #include <mitkImageCast.h>
21 #include <mitkProperties.h>
22 
23 // ITK
24 #include <itkImageRegionIterator.h>
25 
26 // MITK
27 #include <mitkIOUtil.h>
28 
29 // Classification
30 #include <mitkCLUtil.h>
32 #include <mitkAbstractFileReader.h>
33 
34 #include <QDir>
35 #include <QString>
36 #include <QStringList>
37 
38 #include <mitkImageAccessByItk.h>
39 
40 using namespace mitk;
41 
42 typedef unsigned int uint;
43 void ReadMitkProjectImageAndMask(std::string input_file, mitk::Image::Pointer & raw_image, mitk::Image::Pointer & class_mask, Image::Pointer &brain_mask);
44 std::map<unsigned int, double> VolumeUnderMaskByLabel(mitk::Image::Pointer mask);
45 
49 int main(int argc, char* argv[])
50 {
51  mitkCommandLineParser parser;
52  parser.setArgumentPrefix("--", "-");
53 
54  // required params
55  parser.addArgument("inputdir", "i", mitkCommandLineParser::Directory, "Input Directory", "Contains input feature files.", us::Any(), false, false, false, mitkCommandLineParser::Input);
56  parser.addArgument("outputdir", "o", mitkCommandLineParser::Directory, "Output Directory", "Destination of output files.", us::Any(), false, false, false, mitkCommandLineParser::Output);
57 
58  // optional params
59  parser.addArgument("select", "s", mitkCommandLineParser::String, "Item selection", "Using Regular expression, seperated by space e.g.: '*.nrrd *.vtk *test*'",std::string("*.mitk"),true);
60  // parser.addArgument("treecount", "tc", mitkCommandLineParser::Int, "Treecount", "Number of trees.",50,true);
61  // parser.addArgument("treedepth", "td", mitkCommandLineParser::Int, "Treedepth", "Maximal tree depth.",50,true);
62  // parser.addArgument("minsplitnodesize", "min", mitkCommandLineParser::Int, "Minimum split node size.", "Minimum split node size.",2,true);
63  // parser.addArgument("precision", "p", mitkCommandLineParser::Float, "Split precision.", "Precision.", mitk::eps,true);
64  // parser.addArgument("fraction", "f", mitkCommandLineParser::Float, "Fraction of samples per tree.", "Fraction of samples per tree.", 0.6f,true);
65  // parser.addArgument("replacment", "r", mitkCommandLineParser::Bool, "Sample with replacement.", "Sample with replacement.", true,true);
66 
67  // Miniapp Infos
68  parser.setCategory("Classification Tools");
69  parser.setTitle("Native Head CT Segmentation");
70  parser.setDescription("Using vigra random forest");
71  parser.setContributor("Jonas Cordes");
72 
73  // Params parsing
74  std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
75  if (parsedArgs.size()==0)
76  return EXIT_FAILURE;
77 
78  std::string inputdir = us::any_cast<std::string>(parsedArgs["inputdir"]);
79  std::string outputdir = us::any_cast<std::string>(parsedArgs["outputdir"]);
80 
81  // int treecount = parsedArgs.count("treecount") ? us::any_cast<int>(parsedArgs["treecount"]) : 50;
82  // int treedepth = parsedArgs.count("treedepth") ? us::any_cast<int>(parsedArgs["treedepth"]) : 50;
83  // int minsplitnodesize = parsedArgs.count("minsplitnodesize") ? us::any_cast<int>(parsedArgs["minsplitnodesize"]) : 2;
84  // float precision = parsedArgs.count("precision") ? us::any_cast<float>(parsedArgs["precision"]) : mitk::eps;
85  // float fraction = parsedArgs.count("fraction") ? us::any_cast<float>(parsedArgs["fraction"]) : 0.6;
86  // bool withreplacement = parsedArgs.count("replacment") ? us::any_cast<float>(parsedArgs["replacment"]) : true;
87  std::string filt_select = parsedArgs.count("select") ? us::any_cast<std::string>(parsedArgs["select"]) : "*.mitk";
88 
89  QString filter(filt_select.c_str());
90  QDir dir(inputdir.c_str());
91  auto strl = dir.entryList(filter.split(" "),QDir::Files);
92 
93  // create the one
95  uint n_samples = 45000;
96  uint n_samples_per_image = n_samples / strl.size();
97  uint n_features = 2;
98 
99  {// Training
100  Eigen::MatrixXd feature_matrix(n_samples, n_features);
101  Eigen::MatrixXi label_matrix(n_samples, 1);
102 
103  uint pos = 0;
104 
105  for(auto entry : strl)
106  {
107  mitk::Image::Pointer raw_image, class_mask, brain_mask;
108  ReadMitkProjectImageAndMask(inputdir + entry.toStdString(), raw_image, class_mask, brain_mask);
109 
110  mitk::Image::Pointer brain_mask_sampled;
111  AccessFixedDimensionByItk_2(brain_mask, mitk::CLUtil::itkSampleLabel, 3, brain_mask_sampled, n_samples_per_image);
112 
113  mitk::Image::Pointer csf_prob;
114  mitk::CLUtil::ProbabilityMap(raw_image,13.9, 8.3,csf_prob);
117  mitk::CLUtil::FillHoleGrayscale(csf_prob,csf_prob);
118 
119  feature_matrix.block(pos, 0, n_samples_per_image, 1) = mitk::CLUtil::Transform<double>(raw_image,brain_mask_sampled);
120  feature_matrix.block(pos, 1, n_samples_per_image, 1) = mitk::CLUtil::Transform<double>(csf_prob,brain_mask_sampled);
121 
122  label_matrix.block(pos, 0, n_samples_per_image, 1) = mitk::CLUtil::Transform<int>(class_mask, brain_mask_sampled);
123  pos += n_samples_per_image;
124  }
125  classifier->Train(feature_matrix, label_matrix);
126  classifier->PrintParameter();
127  }
128 
129  std::map<std::string, std::pair<double, double> > map_error;
130 
131  for(auto entry: strl)
132  {
133  mitk::Image::Pointer raw_image, class_mask, brain_mask;
134  ReadMitkProjectImageAndMask(inputdir + entry.toStdString(), raw_image, class_mask, brain_mask);
135 
136  auto map_true = VolumeUnderMaskByLabel(class_mask);
137 
138  mitk::Image::Pointer csf_prob;
139  mitk::CLUtil::ProbabilityMap(raw_image,13.9, 8.3,csf_prob);
142  mitk::CLUtil::FillHoleGrayscale(csf_prob,csf_prob);
143 
144  uint count = 0;
145  mitk::CLUtil::CountVoxel(brain_mask,count);
146 
147  Eigen::MatrixXd feature_matrix(count,n_features);
148  feature_matrix.block(0, 0, count, 1) = mitk::CLUtil::Transform<double>(raw_image,brain_mask);
149  feature_matrix.block(0, 1, count, 1) = mitk::CLUtil::Transform<double>(csf_prob,brain_mask);
150 
151  mitk::Image::Pointer result_mask = mitk::CLUtil::Transform<int>(classifier->Predict(feature_matrix),brain_mask);
152 
153  std::string name = itksys::SystemTools::GetFilenameWithoutExtension(entry.toStdString());
154  mitk::IOUtil::Save(result_mask,outputdir + name + ".nrrd");
155 
156  auto map_pred = VolumeUnderMaskByLabel(result_mask);
157 
158  map_error[entry.toStdString()] = std::make_pair(std::abs(map_true[1] - map_pred[1]) / map_true[1], map_true[2] != 0 ? std::abs(map_true[2] - map_pred[2]) / map_true[2]: 0);
159  }
160 
161  double mean_error_csf = 0;
162  double mean_error_les = 0;
163  double num_subjects = map_error.size();
164  for(auto entry: map_error)
165  {
166  MITK_INFO(entry.first.c_str()) << "CSF error: " << entry.second.first << "%\t LES error: " << entry.second.second << "%";
167  mean_error_csf += entry.second.first;
168  mean_error_les += entry.second.second;
169  }
170  MITK_INFO("Mean") << "CSF error: " << mean_error_csf/num_subjects << "%\t LES error: " << mean_error_les/num_subjects << "%";
171 
172  return EXIT_SUCCESS;
173 }
174 
175 void ReadMitkProjectImageAndMask(std::string input_file, mitk::Image::Pointer & raw_image, mitk::Image::Pointer & class_mask, mitk::Image::Pointer & brain_mask)
176 {
177  auto so = mitk::IOUtil::Load(input_file);
178 
179  std::map<uint, uint> map;
180 
181  mitk::CLUtil::CountVoxel(dynamic_cast<mitk::Image *>(so[1].GetPointer()), map);
182 
183  raw_image = map.size() <= 7 ? dynamic_cast<mitk::Image *>(so[0].GetPointer()) : dynamic_cast<mitk::Image *>(so[1].GetPointer());
184  class_mask = map.size() <= 7 ? dynamic_cast<mitk::Image *>(so[1].GetPointer()) : dynamic_cast<mitk::Image *>(so[0].GetPointer());
185 
186  std::map<uint, uint> merge_instructions;// = {{0,0},{1,1},{2,1},{3,1},{4,2},{5,3},{6,3}};
187  merge_instructions[0] = 0;
188  merge_instructions[1] = 1;
189  merge_instructions[2] = 1;
190  merge_instructions[3] = 1;
191  merge_instructions[4] = 2;
192  merge_instructions[5] = 3;
193  merge_instructions[6] = 3;
194  mitk::CLUtil::MergeLabels(class_mask, merge_instructions);
195 
196  brain_mask = class_mask->Clone();
197  //merge_instructions = {{0,0},{1,1},{2,1},{3,1},{4,1},{5,1},{6,1}};
198  merge_instructions[0] = 0;
199  merge_instructions[1] = 1;
200  merge_instructions[2] = 1;
201  merge_instructions[3] = 1;
202  merge_instructions[4] = 1;
203  merge_instructions[5] = 1;
204  merge_instructions[6] = 1;
205  mitk::CLUtil::MergeLabels(brain_mask, merge_instructions);
206 }
207 
208 std::map<unsigned int, double> VolumeUnderMaskByLabel(mitk::Image::Pointer mask)
209 {
210  itk::Image<unsigned int, 3>::Pointer img;
211  mitk::CastToItkImage(mask,img);
212 
213  std::map<unsigned int, double> volume_map;
214  itk::ImageRegionConstIterator<itk::Image<unsigned int, 3>> it(img,img->GetLargestPossibleRegion());
215  while(!it.IsAtEnd())
216  {
217  itk::Image<unsigned int, 3>::PixelType value = it.Value();
218  if(volume_map.find(value)== volume_map.end())
219  volume_map[value] = 0;
220  volume_map[value]++;
221  ++it;
222  }
223 
224  auto spacing = img->GetSpacing();
225  double volumeUnit = spacing[0] * spacing[1] * spacing[2];
226  volumeUnit *= 0.001;
227 
228  for(auto entry: volume_map)
229  entry.second *= volumeUnit;
230 
231  return volume_map;
232 }
#define MITK_INFO
Definition: mitkLogMacros.h:18
static void itkSampleLabel(TImageType1 *image, TImageType2 *output, double acceptrate, unsigned int label)
Definition: mitkCLUtil.h:355
void setContributor(std::string contributor)
DataCollection - Class to facilitate loading/accessing structured data.
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
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)
void ReadMitkProjectImageAndMask(std::string input_file, mitk::Image::Pointer &raw_image, mitk::Image::Pointer &class_mask, Image::Pointer &brain_mask)
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
int main(int argc, char *argv[])
static void DilateGrayscale(mitk::Image::Pointer &image, unsigned int radius, mitk::CLUtil::MorphologicalDimensions d, mitk::Image::Pointer &outimage)
DilateGrayscale.
Definition: mitkCLUtil.cpp:58
static void FillHoleGrayscale(mitk::Image::Pointer &image, mitk::Image::Pointer &outimage)
FillHoleGrayscale.
Definition: mitkCLUtil.cpp:63
static void ErodeGrayscale(mitk::Image::Pointer &image, unsigned int radius, mitk::CLUtil::MorphologicalDimensions d, mitk::Image::Pointer &outimage)
ErodeGrayscale.
Definition: mitkCLUtil.cpp:53
#define AccessFixedDimensionByItk_2(mitkImage, itkImageTypeFunction, dimension, arg1, arg2)
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)
unsigned int uint
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)
std::map< unsigned int, double > VolumeUnderMaskByLabel(mitk::Image::Pointer mask)
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 Save(const mitk::BaseData *data, const std::string &path, bool setPathProperty=false)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:774
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