]> Creatis software - clitk.git/blob - segmentation/clitkExtractMediastinumFilter.txx
Read CricoidCartilag image or Point position (simpler)
[clitk.git] / segmentation / clitkExtractMediastinumFilter.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://www.centreleonberard.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17   ===========================================================================**/
18
19 #ifndef CLITKEXTRACTMEDIASTINUMFILTER_TXX
20 #define CLITKEXTRACTMEDIASTINUMFILTER_TXX
21
22 // clitk
23 #include "clitkCommon.h"
24 #include "clitkExtractMediastinumFilter.h"
25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
26 #include "clitkSliceBySliceRelativePositionFilter.h"
27 #include "clitkSegmentationUtils.h"
28 #include "clitkExtractAirwaysTreeInfoFilter.h"
29 #include "clitkCropLikeImageFilter.h"
30
31 // std
32 #include <deque>
33
34 // itk
35 #include "itkStatisticsLabelMapFilter.h"
36 #include "itkLabelImageToStatisticsLabelMapFilter.h"
37 #include "itkRegionOfInterestImageFilter.h"
38 #include "itkBinaryThresholdImageFilter.h"
39 #include "itkScalarImageKmeansImageFilter.h"
40
41 // itk ENST
42 #include "RelativePositionPropImageFilter.h"
43
44 //--------------------------------------------------------------------
45 template <class ImageType>
46 clitk::ExtractMediastinumFilter<ImageType>::
47 ExtractMediastinumFilter():
48   clitk::StructuresExtractionFilter<ImageType>()
49 {
50   this->SetNumberOfRequiredInputs(1);
51   SetBackgroundValuePatient(0);
52   SetBackgroundValueLung(0);
53   SetBackgroundValueBones(0);
54   SetForegroundValueLeftLung(1);
55   SetForegroundValueRightLung(2);
56   SetDistanceMaxToAnteriorPartOfTheVertebralBody(10);  
57   SetOutputMediastinumFilename("mediastinum.mha");  
58   UseBonesOff();
59 }
60 //--------------------------------------------------------------------
61
62
63 //--------------------------------------------------------------------
64 template <class ImageType>
65 void 
66 clitk::ExtractMediastinumFilter<ImageType>::
67 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg) 
68 {
69   this->SetNthInput(0, const_cast<MaskImageType *>(image));
70   m_BackgroundValuePatient = bg;
71 }
72 //--------------------------------------------------------------------
73
74
75 //--------------------------------------------------------------------
76 template <class ImageType>
77 void 
78 clitk::ExtractMediastinumFilter<ImageType>::
79 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg, 
80                        MaskImagePixelType fgLeft, MaskImagePixelType fgRight) 
81 {
82   this->SetNthInput(1, const_cast<MaskImageType *>(image));
83   m_BackgroundValueLung = bg;
84   m_ForegroundValueLeftLung = fgLeft;
85   m_ForegroundValueRightLung = fgRight;
86 }
87 //--------------------------------------------------------------------
88
89
90 //--------------------------------------------------------------------
91 template <class ImageType>
92 void 
93 clitk::ExtractMediastinumFilter<ImageType>::
94 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg) 
95 {
96   this->SetNthInput(2, const_cast<MaskImageType *>(image));
97   m_BackgroundValueBones = bg;
98 }
99 //--------------------------------------------------------------------
100
101
102 //--------------------------------------------------------------------
103 template <class ImageType>
104 void 
105 clitk::ExtractMediastinumFilter<ImageType>::
106 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg) 
107 {
108   this->SetNthInput(3, const_cast<MaskImageType *>(image));
109   m_BackgroundValueTrachea = bg;
110 }
111 //--------------------------------------------------------------------
112
113 //--------------------------------------------------------------------
114 template <class ImageType>
115 void 
116 clitk::ExtractMediastinumFilter<ImageType>::
117 SetFuzzyThreshold(std::string tag, double value)
118 {
119   m_FuzzyThreshold[tag] = value;
120 }
121 //--------------------------------------------------------------------
122
123
124 //--------------------------------------------------------------------
125 template <class ImageType>
126 double 
127 clitk::ExtractMediastinumFilter<ImageType>::
128 GetFuzzyThreshold(std::string tag)
129 {
130   if (m_FuzzyThreshold.find(tag) == m_FuzzyThreshold.end()) {
131     clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThresholds.");
132     return 0.0;
133   }
134   
135   return m_FuzzyThreshold[tag]; 
136 }
137 //--------------------------------------------------------------------
138
139
140 //--------------------------------------------------------------------
141 template <class ImageType>
142 void 
143 clitk::ExtractMediastinumFilter<ImageType>::
144 GenerateInputRequestedRegion() 
145 {
146   // DD("GenerateInputRequestedRegion");
147   // Do not call default
148   // Superclass::GenerateInputRequestedRegion();  
149   // DD("End GenerateInputRequestedRegion");
150 }
151 //--------------------------------------------------------------------
152
153
154 //--------------------------------------------------------------------
155 template <class ImageType>
156 void 
157 clitk::ExtractMediastinumFilter<ImageType>::
158 SetInput(const ImageType * image) 
159 {
160   this->SetNthInput(0, const_cast<ImageType *>(image));
161 }
162 //--------------------------------------------------------------------
163
164
165
166 //--------------------------------------------------------------------
167 template <class ImageType>
168 void 
169 clitk::ExtractMediastinumFilter<ImageType>::
170 GenerateOutputInformation() { 
171   // Do not call default
172   // Superclass::GenerateOutputInformation();
173
174   //--------------------------------------------------------------------
175   // Get input pointers
176   this->LoadAFDB();
177   ImageConstPointer input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
178   MaskImagePointer patient, lung, bones, trachea;
179   patient = this->GetAFDB()->template GetImage <MaskImageType>("Patient");
180   lung = this->GetAFDB()->template GetImage <MaskImageType>("Lungs");
181   if (GetUseBones()) {
182     bones = this->GetAFDB()->template GetImage <MaskImageType>("Bones");  
183   }
184   trachea = this->GetAFDB()->template GetImage <MaskImageType>("Trachea");  
185
186   //this->VerboseImageSizeFlagOn();
187
188   //--------------------------------------------------------------------
189   // Step : Crop support (patient) to lung extend in RL
190   this->StartNewStep("[Mediastinum] Crop support like lungs along LR");
191   typedef clitk::CropLikeImageFilter<MaskImageType> CropFilterType;
192   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
193   cropFilter->SetInput(patient);
194   cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
195   cropFilter->Update();
196   output = cropFilter->GetOutput();
197   this->template StopCurrentStep<MaskImageType>(output);
198
199   //--------------------------------------------------------------------
200   // Step : Crop support (previous) to bones extend in AP
201   if (GetUseBones()) {
202     this->StartNewStep("[Mediastinum] Crop support like bones along AP");
203     cropFilter = CropFilterType::New();
204     cropFilter->SetInput(output);
205     cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
206     cropFilter->Update();
207     output = cropFilter->GetOutput();
208     this->template StopCurrentStep<MaskImageType>(output);
209   }
210
211   //--------------------------------------------------------------------
212   // Step : patient minus lungs, minus bones, minus trachea
213   this->StartNewStep("[Mediastinum] Patient contours minus lungs, trachea [and bones]");
214   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
215   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
216   boolFilter->InPlaceOn();
217   boolFilter->SetInput1(output);
218   boolFilter->SetInput2(lung);    
219   boolFilter->SetOperationType(BoolFilterType::AndNot);
220   boolFilter->Update();    
221   if (GetUseBones()) {
222     boolFilter->SetInput1(boolFilter->GetOutput());
223     boolFilter->SetInput2(bones);
224     boolFilter->SetOperationType(BoolFilterType::AndNot);
225     boolFilter->Update(); 
226   }
227   boolFilter->SetInput1(boolFilter->GetOutput());
228   boolFilter->SetInput2(trachea);
229   boolFilter->SetOperationType(BoolFilterType::AndNot);
230   boolFilter->Update(); 
231   output = boolFilter->GetOutput();
232
233   // Auto crop to gain support area
234   output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue()); 
235   this->template StopCurrentStep<MaskImageType>(output);
236
237   //--------------------------------------------------------------------
238   // Step : LR limits from lung (need separate lung ?)
239   // Get separate lung images to get only the right and left lung
240   // (because RelativePositionPropImageFilter only consider fg=1);
241   // (label must be '1' because right is greater than left).  (WE DO
242   // NOT NEED TO SEPARATE ? )
243   this->StartNewStep("[Mediastinum] Left/Right limits with lungs");
244   
245   // The following cannot be "inplace" because mask is the same than input ...
246   MaskImagePointer right_lung = 
247     clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 2, 0, false);
248   MaskImagePointer left_lung = 
249     clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 1, 0, false);
250   left_lung = clitk::SetBackground<MaskImageType, MaskImageType>(left_lung, left_lung, 2, 1, false);
251   right_lung = clitk::ResizeImageLike<MaskImageType>(right_lung, output, this->GetBackgroundValue());
252   left_lung = clitk::ResizeImageLike<MaskImageType>(left_lung, output, this->GetBackgroundValue());
253   this->GetAFDB()->template SetImage<MaskImageType>("RightLung", "seg/RightLung.mha", 
254                                                     right_lung, true);
255   this->GetAFDB()->template SetImage<MaskImageType>("LeftLung", "seg/LeftLung.mha", 
256                                                     left_lung, true);
257   this->GetAFDB()->Write();
258   this->template StopCurrentStep<MaskImageType>(output);
259
260   //--------------------------------------------------------------------
261   // Step : AP limits from bones
262   // Separate the bones in the ant-post middle
263   MaskImagePointer bones_ant;
264   MaskImagePointer bones_post;
265   MaskImagePointType middle_AntPost_position;
266   if (GetUseBones()) { 
267     this->StartNewStep("[Mediastinum] Ant/Post limits with bones");
268
269     // To define ant and post part of the bones with a single horizontal line 
270     MaskImageIndexType index_bones_middle;
271
272     // Method1: cut in the middle (not optimal)
273     /*
274     middle_AntPost_position[0] = middle_AntPost_position[2] = 0;
275     middle_AntPost_position[1] = bones->GetOrigin()[1]+
276     (bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
277     DD(middle_AntPost_position);
278     bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
279     */
280   
281     // Method2: Use VertebralBody, take most ant point
282     MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
283     FindExtremaPointInAGivenDirection<MaskImageType>(VertebralBody, this->GetBackgroundValue(), 
284                                                      1, true, middle_AntPost_position);
285     bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);  
286
287     // Split bone image first into two parts (ant and post), and crop
288     // lateraly to get vertebral 
289     typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
290     //  typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
291     typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
292     MaskImageRegionType region = bones->GetLargestPossibleRegion();
293     MaskImageSizeType size = region.GetSize();
294     MaskImageIndexType index = region.GetIndex();
295     // ANT part
296     // crop LR to keep 1/4 center part
297     //    index[0] = size[0]/4+size[0]/8;
298     //size[0] = size[0]/4; 
299     // crop AP to keep first (ant) part
300     size[1] =  index_bones_middle[1]; //size[1]/2.0;
301     region.SetSize(size);
302     region.SetIndex(index);
303     roiFilter->SetInput(bones);
304     roiFilter->SetRegionOfInterest(region);
305     roiFilter->ReleaseDataFlagOff();
306     roiFilter->Update();
307     bones_ant = roiFilter->GetOutput();
308     //      writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
309     // POST part
310     roiFilter = ROIFilterType::New();  
311     index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
312     size[1] =  bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
313     region.SetIndex(index);
314     region.SetSize(size);
315     roiFilter->SetInput(bones);
316     roiFilter->SetRegionOfInterest(region);
317     roiFilter->ReleaseDataFlagOff();
318     roiFilter->Update();
319     bones_post = roiFilter->GetOutput();
320     // writeImage<MaskImageType>(bones_post, "b_post.mhd");
321
322     // Now, insert this image in the AFDB ==> (needed because used in the RelativePosition config file)
323     this->GetAFDB()->template SetImage<MaskImageType>("Bones_Post", "seg/Bones_Post.mha", 
324                                                  bones_post, true);
325     this->GetAFDB()->template SetImage<MaskImageType>("Bones_Ant", "seg/Bones_Ant.mha", 
326                                                  bones_ant, true);
327     this->GetAFDB()->Write();
328
329     this->template StopCurrentStep<MaskImageType>(output);
330   }
331
332   //--------------------------------------------------------------------
333   // Remove VertebralBody part
334   this->StartNewStep("[Mediastinum] Remove VertebralBody");  
335   MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
336   clitk::AndNot<MaskImageType>(output, VertebralBody, this->GetBackgroundValue());
337   this->template StopCurrentStep<MaskImageType>(output);
338
339   //--------------------------------------------------------------------
340   // Generic RelativePosition processes
341   output = this->ApplyRelativePositionList("Mediastinum", output);
342
343   //--------------------------------------------------------------------
344   // Step : SI limits It is better to do this limit *AFTER* the
345   // RelativePosition to avoid some issue due to superior boundaries.
346   this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag");
347   // load Cricoid, get centroid, cut above (or below), lower bound
348   MaskImagePointType p;
349   try {
350     MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage <MaskImageType>("CricoidCartilag");
351     p[0] = p[1] = p[2] =  0.0; // to avoid warning
352     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(CricoidCartilag, 
353                                                             this->GetBackgroundValue(), 2, true, p);
354   } catch (clitk::ExceptionObject e) {
355     //DD("CricoidCartilag image not found, try CricoidCartilagZ");
356     this->GetAFDB()->GetPoint3D("CricoidCartilagPoint", p);
357   }
358   output = clitk::CropImageRemoveGreaterThan<MaskImageType>(output, 2, p[2], true, this->GetBackgroundValue());
359   this->template StopCurrentStep<MaskImageType>(output);
360
361   //--------------------------------------------------------------------
362   // Step: Get CCL
363   this->StartNewStep("[Mediastinum] Keep main connected component");
364   output = clitk::Labelize<MaskImageType>(output, this->GetBackgroundValue(), false, 500);
365   // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
366   output = clitk::KeepLabels<MaskImageType>(output, this->GetBackgroundValue(), 
367                                             this->GetForegroundValue(), 1, 1, 0);
368   this->template StopCurrentStep<MaskImageType>(output);
369
370   //--------------------------------------------------------------------
371   // Step: Remove post part from VertebralBody
372   this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody");
373   RemovePostPartOfVertebralBody();
374   this->template StopCurrentStep<MaskImageType>(output);
375
376   //--------------------------------------------------------------------
377   // Step: Slice by Slice CCL
378   this->StartNewStep("[Mediastinum] Slice by Slice keep only one component");
379   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
380   //  typename ExtractSliceFilterType::Pointer 
381   ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
382   extractSliceFilter->SetInput(output);
383   extractSliceFilter->SetDirection(2);
384   extractSliceFilter->Update();
385   typedef typename ExtractSliceFilterType::SliceType SliceType;
386   std::vector<typename SliceType::Pointer> mSlices;
387   extractSliceFilter->GetOutputSlices(mSlices);
388   for(unsigned int i=0; i<mSlices.size(); i++) {
389     mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
390     mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
391   }
392   typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
393   typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
394   joinFilter->SetOrigin(output->GetOrigin()[2]);
395   joinFilter->SetSpacing(output->GetSpacing()[2]);
396   for(unsigned int i=0; i<mSlices.size(); i++) {
397     joinFilter->PushBackInput(mSlices[i]);
398   }
399   joinFilter->Update();
400   output = joinFilter->GetOutput();
401   this->template StopCurrentStep<MaskImageType>(output);
402
403   //--------------------------------------------------------------------
404   // Step 10 : AutoCrop
405   this->StartNewStep("[Mediastinum] AutoCrop");
406   output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue()); 
407   this->template StopCurrentStep<MaskImageType>(output);
408
409   // End, set the real size
410   this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
411   this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
412   this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
413   this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
414 }
415 //--------------------------------------------------------------------
416
417
418 //--------------------------------------------------------------------
419 template <class ImageType>
420 void 
421 clitk::ExtractMediastinumFilter<ImageType>::
422 GenerateData() 
423 {
424   this->GraftOutput(output);
425   // Store image filenames into AFDB 
426   this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());  
427   this->WriteAFDB();
428 }
429 //--------------------------------------------------------------------
430
431   
432 //--------------------------------------------------------------------
433 template <class ImageType>
434 void 
435 clitk::ExtractMediastinumFilter<ImageType>::
436 RemovePostPartOfVertebralBody()
437 {
438   
439   /*
440     Posteriorly, Station 8 abuts the descending aorta and the anterior
441     aspect of the vertebral body until an imaginary horizontal line
442     running 1 cm posterior to the anterior border of the vertebral
443     body (Fig. 3C).
444     
445     => We use this definition for all the mediastinum
446
447    Find most Ant point in VertebralBody. Consider the horizontal line
448    which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from
449    the most ant point.
450   */
451
452   // Get VertebralBody mask image
453   MaskImagePointer VertebralBody = 
454     this->GetAFDB()->template GetImage <MaskImageType>("VertebralBody");  
455
456   // Consider vertebral body slice by slice
457   std::vector<MaskSlicePointer> vertebralSlices;
458   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vertebralSlices);
459
460   // For each slice, compute the most anterior point
461   std::map<int, MaskSlicePointType> vertebralAntPositionBySlice;
462   for(uint i=0; i<vertebralSlices.size(); i++) {
463     MaskSlicePointType p;
464     bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vertebralSlices[i], 
465                                                                      this->GetBackgroundValue(), 
466                                                                      1, true, p);
467     if (found) {
468       vertebralAntPositionBySlice[i] = p;
469     }
470     else { 
471       // It should not happen ! But sometimes, a contour is missing or
472       // the VertebralBody is not delineated enough inferiorly ... in
473       // those cases, we consider the first found slice.
474       //        std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
475       // [ Possible alternative -> consider previous limit ]
476     }
477   }
478
479   // Convert 2D points in slice into 3D points
480   std::vector<MaskImagePointType> vertebralAntPositions;
481   clitk::PointsUtils<MaskImageType>::Convert2DMapTo3DList(vertebralAntPositionBySlice, 
482                                                           VertebralBody, 
483                                                           vertebralAntPositions);
484
485   // DEBUG : write list of points
486   clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions, 
487                                              "VertebralBodyMostAnteriorPoints.txt");
488
489   // Cut support posteriorly 1cm the most anterior point of the
490   // VertebralBody. Do nothing below and above the VertebralBody. To
491   // do that compute several region, slice by slice and fill. 
492   MaskImageRegionType region;
493   MaskImageSizeType size;
494   MaskImageIndexType start;
495   size[2] = 1; // a single slice
496   start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
497   size[0] = output->GetLargestPossibleRegion().GetSize()[0];
498   for(uint i=0; i<vertebralAntPositions.size(); i++) {
499     typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
500     IteratorType iter = 
501       IteratorType(output, output->GetLargestPossibleRegion());
502     MaskImageIndexType index;
503     // Consider some cm posterior to most anterior positions (usually
504     // 1 cm).
505     vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody();
506     // Get index of this point
507     output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
508     // Compute region (a single slice)
509     start[2] = index[2];
510     start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
511     size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
512     region.SetSize(size);
513     region.SetIndex(start);
514     // Fill region
515     if (output->GetLargestPossibleRegion().IsInside(start))  {
516       itk::ImageRegionIterator<MaskImageType> it(output, region);
517       it.GoToBegin();
518       while (!it.IsAtEnd()) {
519         it.Set(this->GetBackgroundValue());
520         ++it;
521       }
522     }
523   }  
524 }
525 //--------------------------------------------------------------------
526
527
528 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX