ANIMA  4.0
animaGcStremMsLesionsSegmentation.cxx
Go to the documentation of this file.
1 #include <tclap/CmdLine.h>
2 
3 #include <itkTimeProbe.h>
6 
7 int main(int argc, const char** argv)
8 {
9  const unsigned int Dimension = 3;
10  typedef itk::Image <double,Dimension> InputImageTypeD;
11  typedef itk::Image <unsigned char,Dimension> InputImageTypeUC;
13 
14  // Parsing arguments
15  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
16 
17  // Setting up parameters
18 
19  // Input filenames
20  TCLAP::ValueArg<std::string> inputFileT1Arg("i","input-t1","T1 input image",false,"","t1 input image",cmd);
21  TCLAP::ValueArg<std::string> inputFileT2Arg("j","input-t2","T2 input image",false,"","t2 input image",cmd);
22  TCLAP::ValueArg<std::string> inputFileDPArg("k","input-dp","DP input image",false,"","dp input image",cmd);
23  TCLAP::ValueArg<std::string> inputFileFLAIRArg("l","input-flair","FLAIR input image",false,"","flair input image",cmd);
24  TCLAP::ValueArg<std::string> inputFileT1GdArg("","input-t1Gd","T1Gd input image",false,"","t1gd input image",cmd);
25 
26  TCLAP::ValueArg<std::string> atlasFileCSFArg("x","atlas-csf","CSF atlas image",false,"","csf atlas image",cmd);
27  TCLAP::ValueArg<std::string> atlasFileGMArg("y","atlas-gm","GM atlas image",false,"","gm atlas image",cmd);
28  TCLAP::ValueArg<std::string> atlasFileWMArg("z","atlas-wm","WM atlas image",false,"","wm atlas image",cmd);
29 
30  TCLAP::ValueArg<std::string> lesionPriorArg("","lesion-prior","Lesions prior image",false,"","lesions prior image",cmd);
31  TCLAP::ValueArg<double> lesionPriorProportionArg("","lesion-prior-rate","Lesions prior proportion (default: 0.23)",false,0.23,"lesions prior proportion",cmd);
32 
33  TCLAP::ValueArg<std::string> maskFileArg("m","mask","Brain mask",true,"","brain mask",cmd);
34 
35  TCLAP::ValueArg<std::string> sourcesMaskFileArg("","so","Binary source input image for graph cut",false,"","sources input image",cmd);
36  TCLAP::ValueArg<std::string> sinksMaskFileArg("","si","Binary sink input image for graph cut",false,"","sinks input image",cmd);
37 
38  // global
39  TCLAP::ValueArg<unsigned int> numThreadsArg("T","threads","Number of execution threads (default: 0 = all cores)",false,0,"number of threads",cmd);
40  TCLAP::SwitchArg verboseArg("v","verbose","verbose mode (default: false)",cmd,false);
41  TCLAP::ValueArg<double> tolArg("","tol","Filter tolerance (default: 0.0001)",false,0.0001,"filter tolerance",cmd);
42 
43  TCLAP::ValueArg<unsigned int> lesionsSelectionArg("a","lesion-select-method","Lesions selection method (0: strem, 1: gcem, 2: gcem + manual graph cut, 3: manual graph cut, default: 0)",false,0,"lesions selection",cmd);
44 
45  // Automatic segmentation parameters
46  TCLAP::ValueArg<unsigned int> initMethodArg("","ini","Initialisation method (0: Atlas, 1: Hierar DP, 2: Hierar FLAIR, default: 1)",false,1,"initialisation method",cmd);
47  TCLAP::ValueArg<int> emIterArg("","iter","Maximum number of iterations in EM algorithm (default: 100)",false,100,"maximum number of iterations",cmd);
48  TCLAP::ValueArg<double> minDistanceArg("","min-dist","Minimum distance in EM algorithm (default: 0.0001)",false,0.0001,"minimum distance",cmd);
49  TCLAP::ValueArg<double> rejRatioArg("","rej","Rejection ratio for EM algorithm (default: 0.2)",false,0.2,"rejection ratio",cmd);
50  TCLAP::ValueArg<int> emIter_concentrationArg("","iter-c","Maximum number of iterations between concentration step (default: 100)",false,100,"maximum number of iterations between concentration step ",cmd);
51  TCLAP::SwitchArg em_before_concentrationArg("","em-c","Use EM before concentration (default: false)",cmd,false);
52 
53  TCLAP::ValueArg<double> thStremCSFArg("","th-stremCSF","Strem threshold in CSF image (default: 0.5)",false,0.5,"strem threshold CSF",cmd);
54  TCLAP::ValueArg<double> thStremGMArg("","th-stremGM","Strem threshold in GM image (default: 0.5)",false,0.5,"strem threshold GM",cmd);
55  TCLAP::ValueArg<double> thStremWMArg("","th-stremWM","Strem threshold in WM image (default: 0.5)",false,0.5,"strem threshold WM",cmd);
56  TCLAP::ValueArg<double> fuzzyRuleMinArg("","min-fuzzy","Minimum factor in fuzzy rules (default: 1)",false,1,"Minimum in fuzzy rules ",cmd);
57  TCLAP::ValueArg<double> fuzzyRuleMaxArg("","max-fuzzy","Maximum factor in fuzzy rules (default: 2)",false,2,"Maximum in fuzzy rules",cmd);
58 
59  TCLAP::ValueArg<std::string> readSolutionFileArg("","read-sol","Filename to read automatic segmentation solution",false,"","read solution filename",cmd);
60  TCLAP::ValueArg<std::string> writeSolutionFileArg("","write-sol","Filename to write automatic segmentation solution",false,"","write solution filename",cmd);
61 
62  TCLAP::SwitchArg useT2Arg("","useT2","Use T2 image in automatic segmentation (default: false)",cmd,false);
63  TCLAP::SwitchArg useDPArg("","useDP","Use DP image in automatic segmentation (default: false)",cmd,false);
64  TCLAP::SwitchArg useFLAIRArg("","useFLAIR","Use FLAIR image in automatic segmentation(default: false)",cmd,false);
65 
66  // Graph cut parameters
67  TCLAP::SwitchArg notUseSpecGradArg("","no-usg","Do not use spectral gradient (default: false)",cmd,false);
68  TCLAP::ValueArg<double> multiVarSourcesArg("","mv","Coefficient to multiply the variance value of the source seeds (default: 1)",false,1,"sources multiply variance",cmd);
69  TCLAP::ValueArg<double> multiVarSinksArg("","ms","Coefficient to multiply the variance value of the sink seeds (default: 1)",false,1,"sinks multiply variance",cmd);
70  TCLAP::ValueArg<double> sigmaArg("","sigma","Sigma value (default: 0.6)",false,0.6,"sigma",cmd);
71  TCLAP::ValueArg<double> alphaArg("","alpha","Alpha value (default: 10)",false,10,"alpha",cmd);
72  TCLAP::ValueArg<std::string> matrixGradArg("","mat","Spectral gradient matrix file",false,"","spectral gradient matrix file",cmd);
73 
74  // Heuristic rules
75  TCLAP::ValueArg<double> minLesionSizeArg("","ml","Minimum lesion size in mm3 (default: 0)",false,0,"minimum lesion size",cmd);
76  TCLAP::SwitchArg removeBorderArg("","rb","Remove lesions that do not touch mask border (default: false)",cmd,false);
77  TCLAP::ValueArg<double> intensityT2Arg("","intT2","T2 intensity (default: 0)",false,0,"T2 intensity ",cmd);
78  TCLAP::ValueArg<double> intensityDPArg("","intDP","DP intensity (default: 0)",false,0,"DP intensity",cmd);
79  TCLAP::ValueArg<double> intensityFLAIRArg("","intFLAIR","FLAIR intensity (default: 0)",false,0,"FLAIR intensity",cmd);
80  TCLAP::ValueArg<double> thWMmapArg("","th-wm-map","Probability threshold to define white matter map (default: 0.2)",false,0.2,"threshold white matter map",cmd);
81  TCLAP::ValueArg<double> contourWMRatioArg("","WMratio","Minimum percentage of WM that must surround lesions",false,0,"WM ratio contour",cmd);
82 
83  // Output filenames
84  TCLAP::ValueArg<std::string> outputLesionFileArg("o","out-le","Output containing lesion segmentation",false,"","lesion segmentation output filename",cmd);
85  TCLAP::ValueArg<std::string> outputCSFFileArg("","out-csf","Output containing CSF region",false,"","csf output filename",cmd);
86  TCLAP::ValueArg<std::string> outputGMFileArg("","out-gm","Output containing GM region",false,"","gm output filename",cmd);
87  TCLAP::ValueArg<std::string> outputWMFileArg("","out-wm","Output containing WM region",false,"","wm output filename",cmd);
88  TCLAP::ValueArg<std::string> outputWholeFileArg("","out-wo","Output containing the whole brain segmentation",false,"","whole segmentation output filename",cmd);
89 
90  TCLAP::ValueArg<std::string> outputGCFileArg("","out-gc","Output containing graph cut segmentation",false,"","graph cut output filename",cmd);
91  TCLAP::ValueArg<std::string> outputStremFileArg("","out-strem","Output of automatic strem computation",false,"","strem output",cmd);
92  TCLAP::ValueArg<std::string> outputStremCSFFileArg("","out-strem-csf","Output of automatic csf strem computation",false,"","strem csf output",cmd);
93  TCLAP::ValueArg<std::string> outputStremGMFileArg("","out-strem-gm","Output of automatic grey matter strem computation",false,"","strem gm output",cmd);
94  TCLAP::ValueArg<std::string> outputStremWMFileArg("","out-strem-wm","Output of automatic white matter strem computation",false,"","strem wm output",cmd);
95 
96  TCLAP::ValueArg<std::string> outputObjectFileArg("","out-so", "Fuzzy object probability map (sources input for graph cut)",false,"","sources filename",cmd);
97 
98  TCLAP::ValueArg<std::string> outputMahaMaxiFileArg("","out-maha-maxi"," Mahalanobis maximum image output filename",false,"","mahalanobis maximum output filename",cmd);
99  TCLAP::ValueArg<std::string> outputMahaMiniFileArg("","out-maha-mini","Mahalanobis minimum image output filename",false,"","mahalanobis minimum output filename",cmd);
100  TCLAP::ValueArg<std::string> outputMahaCSFFileArg("","out-maha-csf","Mahalanobis CSF image output filename",false,"","mahalanobis csf output filename",cmd);
101  TCLAP::ValueArg<std::string> outputMahaGMFileArg("","out-maha-gm","Mahalanobis GM image output filename",false,"","mahalanobis gm output filename",cmd);
102  TCLAP::ValueArg<std::string> outputMahaWMFileArg("","out-maha-wm","Mahalanobis WM image output filename",false,"","mahalanobis wm output filename",cmd);
103 
104  TCLAP::ValueArg<std::string> outputHyperIntensity1FileArg("","out-hyper-im1","Hyper-intensity map in image 1 filename",false,"","output hyper-intensity filename 1",cmd);
105  TCLAP::ValueArg<std::string> outputHyperIntensity2FileArg("","out-hyper-im2","Hyper-intensity map in image 2 filename",false,"","output hyper-intensity filename 2",cmd);
106 
107  try
108  {
109  cmd.parse(argc,argv);
110  }
111  catch (TCLAP::ArgException& e)
112  {
113  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
114  return EXIT_FAILURE;
115  }
116 
117  FilterTypeSeg::Pointer segFilter = FilterTypeSeg::New();
118 
119  if (inputFileT1Arg.getValue() != "")
120  segFilter->SetInputImageT1( anima::readImage<InputImageTypeD>( inputFileT1Arg.getValue() ) );
121 
122  if(inputFileT2Arg.getValue() != "")
123  segFilter->SetInputImageT2( anima::readImage<InputImageTypeD>( inputFileT2Arg.getValue() ) );
124 
125  if (inputFileDPArg.getValue() != "")
126  segFilter->SetInputImageDP( anima::readImage<InputImageTypeD>( inputFileDPArg.getValue() ) );
127 
128  if (inputFileFLAIRArg.getValue() != "")
129  segFilter->SetInputImageFLAIR( anima::readImage<InputImageTypeD>( inputFileFLAIRArg.getValue() ) );
130 
131  if (inputFileT1GdArg.getValue() != "")
132  segFilter->SetInputImageT1Gd( anima::readImage<InputImageTypeD>( inputFileT1GdArg.getValue() ) );
133 
134  if( atlasFileCSFArg.getValue()!="" )
135  {
136  segFilter->SetInputCSFAtlas( anima::readImage<InputImageTypeD>( atlasFileCSFArg.getValue() ) );
137  }
138 
139  if( atlasFileGMArg.getValue()!="" )
140  {
141  segFilter->SetInputGMAtlas( anima::readImage<InputImageTypeD>( atlasFileGMArg.getValue() ) );
142  }
143 
144  if( atlasFileWMArg.getValue()!="" )
145  {
146  segFilter->SetInputWMAtlas( anima::readImage<InputImageTypeD>( atlasFileWMArg.getValue() ) );
147  }
148 
149  if (lesionPriorArg.getValue() != "")
150  {
151  segFilter->SetInputLesionPrior(anima::readImage<InputImageTypeD> (lesionPriorArg.getValue()));
152  segFilter->SetLesionPriorProportion(lesionPriorProportionArg.getValue());
153  }
154 
155  if( maskFileArg.getValue()!="" )
156  {
157  segFilter->SetMask( anima::readImage<InputImageTypeD>( maskFileArg.getValue() ) );
158  }
159 
160  if( sourcesMaskFileArg.getValue()!="" )
161  {
162  segFilter->SetSourcesMask( anima::readImage<InputImageTypeUC>( sourcesMaskFileArg.getValue()) );
163  }
164 
165  if( sinksMaskFileArg.getValue()!="" )
166  {
167  segFilter->SetSinksMask( anima::readImage<InputImageTypeUC>( sinksMaskFileArg.getValue()) );
168  }
169 
170  // Set parameters
171 
172  segFilter->SetLesionSegmentationType( (LesionSegmentationType) lesionsSelectionArg.getValue() );
173  segFilter->SetVerbose( verboseArg.getValue() );
174  segFilter->SetNumberOfWorkUnits( numThreadsArg.getValue() );
175  segFilter->SetTol( tolArg.getValue() );
176 
177  segFilter->SetUseT2( useT2Arg.getValue() );
178  segFilter->SetUseDP( useDPArg.getValue() );
179  segFilter->SetUseFLAIR( useFLAIRArg.getValue() );
180  segFilter->SetInitMethodType( (InitializationType) initMethodArg.getValue() );
181  segFilter->SetEmIter( emIterArg.getValue() );
182  segFilter->SetMinDistance( minDistanceArg.getValue() );
183  segFilter->SetRejRatio( rejRatioArg.getValue());
184  segFilter->SetEmIter_concentration( emIter_concentrationArg.getValue() );
185  segFilter->SetEM_before_concentration( em_before_concentrationArg.getValue() );
186  segFilter->SetMahalanobisThCSF( thStremCSFArg.getValue() );
187  segFilter->SetMahalanobisThGM( thStremGMArg.getValue() );
188  segFilter->SetMahalanobisThWM( thStremWMArg.getValue() );
189  segFilter->SetFuzzyRuleMin( fuzzyRuleMinArg.getValue() );
190  segFilter->SetFuzzyRuleMax( fuzzyRuleMaxArg.getValue() );
191  segFilter->SetSolutionReadFilename( readSolutionFileArg.getValue() );
192  segFilter->SetSolutionWriteFilename( writeSolutionFileArg.getValue() );
193 
194  segFilter->SetUseSpecGrad( !(notUseSpecGradArg.getValue()) );
195  segFilter->SetMultiVarSources( multiVarSourcesArg.getValue() );
196  segFilter->SetMultiVarSinks( multiVarSinksArg.getValue() );
197  segFilter->SetAlpha( alphaArg.getValue() );
198  segFilter->SetSigma( sigmaArg.getValue() );
199  segFilter->SetMatrixGradFilename( matrixGradArg.getValue() );
200 
201  segFilter->SetIntensityT2Factor( intensityT2Arg.getValue() );
202  segFilter->SetIntensityDPFactor( intensityDPArg.getValue() );
203  segFilter->SetIntensityFLAIRFactor( intensityFLAIRArg.getValue() );
204  segFilter->SetRemoveBorder( removeBorderArg.getValue() );
205  segFilter->SetMinLesionsSize( minLesionSizeArg.getValue() );
206  segFilter->SetThresoldWMmap( thWMmapArg.getValue());
207  segFilter->SetRatioContourWM( contourWMRatioArg.getValue());
208 
209  segFilter->SetOutputLesionFilename( outputLesionFileArg.getValue() );
210  segFilter->SetOutputCSFFilename( outputCSFFileArg.getValue() );
211  segFilter->SetOutputGMFilename( outputGMFileArg.getValue() );
212  segFilter->SetOutputWMFilename( outputWMFileArg.getValue() );
213  segFilter->SetOutputWholeFilename( outputWholeFileArg.getValue() );
214 
215  segFilter->SetOutputGCFilename( outputGCFileArg.getValue() );
216  segFilter->SetOutputStremFilename( outputStremFileArg.getValue() );
217  segFilter->SetOutputStremCSFFilename( outputStremCSFFileArg.getValue() );
218  segFilter->SetOutputStremGMFilename( outputStremGMFileArg.getValue() );
219  segFilter->SetOutputStremWMFilename( outputStremWMFileArg.getValue() );
220 
221  segFilter->SetOutputFuzzyObjectFilename( outputObjectFileArg.getValue() );
222  segFilter->SetOutputMahaMaximumFilename( outputMahaMaxiFileArg.getValue() );
223  segFilter->SetOutputMahaMinimumFilename( outputMahaMiniFileArg.getValue() );
224  segFilter->SetOutputMahaCSFFilename( outputMahaCSFFileArg.getValue() );
225  segFilter->SetOutputMahaGMFilename( outputMahaGMFileArg.getValue() );
226  segFilter->SetOutputMahaWMFilename( outputMahaWMFileArg.getValue() );
227 
228  segFilter->SetOutputIntensityImage1Filename( outputHyperIntensity1FileArg.getValue() );
229  segFilter->SetOutputIntensityImage2Filename( outputHyperIntensity2FileArg.getValue() );
230 
231  // Process
232  itk::TimeProbe timer;
233 
234  timer.Start();
235 
236  try
237  {
238  segFilter->Update();
239  segFilter->WriteOutputs();
240  }
241  catch (itk::ExceptionObject &e)
242  {
243  std::cerr << e << std::endl;
244  return(1);
245  }
246 
247  timer.Stop();
248 
249  std::cout << "Elapsed Time: " << timer.GetTotal() << timer.GetUnit() << std::endl;
250 
251  return EXIT_SUCCESS;
252 
253 }
int main(int argc, const char **argv)