1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
19 #ifndef CLITKEXTRACTMEDIASTINUMFILTER_TXX
20 #define CLITKEXTRACTMEDIASTINUMFILTER_TXX
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"
35 #include "itkStatisticsLabelMapFilter.h"
36 #include "itkLabelImageToStatisticsLabelMapFilter.h"
37 #include "itkRegionOfInterestImageFilter.h"
38 #include "itkBinaryThresholdImageFilter.h"
39 #include "itkScalarImageKmeansImageFilter.h"
42 #include "RelativePositionPropImageFilter.h"
44 //--------------------------------------------------------------------
45 template <class ImageType>
46 clitk::ExtractMediastinumFilter<ImageType>::
47 ExtractMediastinumFilter():
48 clitk::StructuresExtractionFilter<ImageType>()
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");
60 //--------------------------------------------------------------------
63 //--------------------------------------------------------------------
64 template <class ImageType>
66 clitk::ExtractMediastinumFilter<ImageType>::
67 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg)
69 this->SetNthInput(0, const_cast<MaskImageType *>(image));
70 m_BackgroundValuePatient = bg;
72 //--------------------------------------------------------------------
75 //--------------------------------------------------------------------
76 template <class ImageType>
78 clitk::ExtractMediastinumFilter<ImageType>::
79 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg,
80 MaskImagePixelType fgLeft, MaskImagePixelType fgRight)
82 this->SetNthInput(1, const_cast<MaskImageType *>(image));
83 m_BackgroundValueLung = bg;
84 m_ForegroundValueLeftLung = fgLeft;
85 m_ForegroundValueRightLung = fgRight;
87 //--------------------------------------------------------------------
90 //--------------------------------------------------------------------
91 template <class ImageType>
93 clitk::ExtractMediastinumFilter<ImageType>::
94 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg)
96 this->SetNthInput(2, const_cast<MaskImageType *>(image));
97 m_BackgroundValueBones = bg;
99 //--------------------------------------------------------------------
102 //--------------------------------------------------------------------
103 template <class ImageType>
105 clitk::ExtractMediastinumFilter<ImageType>::
106 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg)
108 this->SetNthInput(3, const_cast<MaskImageType *>(image));
109 m_BackgroundValueTrachea = bg;
111 //--------------------------------------------------------------------
113 //--------------------------------------------------------------------
114 template <class ImageType>
116 clitk::ExtractMediastinumFilter<ImageType>::
117 SetFuzzyThreshold(std::string tag, double value)
119 m_FuzzyThreshold[tag] = value;
121 //--------------------------------------------------------------------
124 //--------------------------------------------------------------------
125 template <class ImageType>
127 clitk::ExtractMediastinumFilter<ImageType>::
128 GetFuzzyThreshold(std::string tag)
130 if (m_FuzzyThreshold.find(tag) == m_FuzzyThreshold.end()) {
131 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThresholds.");
135 return m_FuzzyThreshold[tag];
137 //--------------------------------------------------------------------
140 //--------------------------------------------------------------------
141 template <class ImageType>
143 clitk::ExtractMediastinumFilter<ImageType>::
144 GenerateInputRequestedRegion()
146 // DD("GenerateInputRequestedRegion");
147 // Do not call default
148 // Superclass::GenerateInputRequestedRegion();
149 // DD("End GenerateInputRequestedRegion");
151 //--------------------------------------------------------------------
154 //--------------------------------------------------------------------
155 template <class ImageType>
157 clitk::ExtractMediastinumFilter<ImageType>::
158 SetInput(const ImageType * image)
160 this->SetNthInput(0, const_cast<ImageType *>(image));
162 //--------------------------------------------------------------------
166 //--------------------------------------------------------------------
167 template <class ImageType>
169 clitk::ExtractMediastinumFilter<ImageType>::
170 GenerateOutputInformation() {
171 // Do not call default
172 // Superclass::GenerateOutputInformation();
174 //--------------------------------------------------------------------
175 // Get input pointers
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");
182 bones = this->GetAFDB()->template GetImage <MaskImageType>("Bones");
184 trachea = this->GetAFDB()->template GetImage <MaskImageType>("Trachea");
186 //this->VerboseImageSizeFlagOn();
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);
199 //--------------------------------------------------------------------
200 // Step : Crop support (previous) to bones extend in AP
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);
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();
222 boolFilter->SetInput1(boolFilter->GetOutput());
223 boolFilter->SetInput2(bones);
224 boolFilter->SetOperationType(BoolFilterType::AndNot);
225 boolFilter->Update();
227 boolFilter->SetInput1(boolFilter->GetOutput());
228 boolFilter->SetInput2(trachea);
229 boolFilter->SetOperationType(BoolFilterType::AndNot);
230 boolFilter->Update();
231 output = boolFilter->GetOutput();
233 // Auto crop to gain support area
234 output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
235 this->template StopCurrentStep<MaskImageType>(output);
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");
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",
255 this->GetAFDB()->template SetImage<MaskImageType>("LeftLung", "seg/LeftLung.mha",
257 this->GetAFDB()->Write();
258 this->template StopCurrentStep<MaskImageType>(output);
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 middle_AntPost_position.Fill(NumericTraits<MaskImagePointType::ValueType>::Zero);
268 this->StartNewStep("[Mediastinum] Ant/Post limits with bones");
270 // To define ant and post part of the bones with a single horizontal line
271 MaskImageIndexType index_bones_middle;
273 // Method1: cut in the middle (not optimal)
275 middle_AntPost_position[0] = middle_AntPost_position[2] = 0;
276 middle_AntPost_position[1] = bones->GetOrigin()[1]+
277 (bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
278 DD(middle_AntPost_position);
279 bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
282 // Method2: Use VertebralBody, take most ant point
283 MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
284 FindExtremaPointInAGivenDirection<MaskImageType>(VertebralBody, this->GetBackgroundValue(),
285 1, true, middle_AntPost_position);
286 bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
288 // Split bone image first into two parts (ant and post), and crop
289 // lateraly to get vertebral
290 typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
291 // typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
292 typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
293 MaskImageRegionType region = bones->GetLargestPossibleRegion();
294 MaskImageSizeType size = region.GetSize();
295 MaskImageIndexType index = region.GetIndex();
297 // crop LR to keep 1/4 center part
298 // index[0] = size[0]/4+size[0]/8;
299 //size[0] = size[0]/4;
300 // crop AP to keep first (ant) part
301 size[1] = index_bones_middle[1]; //size[1]/2.0;
302 region.SetSize(size);
303 region.SetIndex(index);
304 roiFilter->SetInput(bones);
305 roiFilter->SetRegionOfInterest(region);
306 roiFilter->ReleaseDataFlagOff();
308 bones_ant = roiFilter->GetOutput();
309 // writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
311 roiFilter = ROIFilterType::New();
312 index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
313 size[1] = bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
314 region.SetIndex(index);
315 region.SetSize(size);
316 roiFilter->SetInput(bones);
317 roiFilter->SetRegionOfInterest(region);
318 roiFilter->ReleaseDataFlagOff();
320 bones_post = roiFilter->GetOutput();
321 // writeImage<MaskImageType>(bones_post, "b_post.mhd");
323 // Now, insert this image in the AFDB ==> (needed because used in the RelativePosition config file)
324 this->GetAFDB()->template SetImage<MaskImageType>("Bones_Post", "seg/Bones_Post.mha",
326 this->GetAFDB()->template SetImage<MaskImageType>("Bones_Ant", "seg/Bones_Ant.mha",
328 this->GetAFDB()->Write();
330 this->template StopCurrentStep<MaskImageType>(output);
333 //--------------------------------------------------------------------
334 // Remove VertebralBody part
335 this->StartNewStep("[Mediastinum] Remove VertebralBody");
336 MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
337 clitk::AndNot<MaskImageType>(output, VertebralBody, this->GetBackgroundValue());
338 this->template StopCurrentStep<MaskImageType>(output);
340 //--------------------------------------------------------------------
341 // Generic RelativePosition processes
342 output = this->ApplyRelativePositionList("Mediastinum", output);
345 //--------------------------------------------------------------------
346 // FIXME --> do not put this limits here !
348 // Step : SI limits It is better to do this limit *AFTER* the
349 // RelativePosition to avoid some issue due to superior boundaries.
350 this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag");
351 // load Cricoid, get centroid, cut above (or below), lower bound
352 MaskImagePointType p;
354 MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage <MaskImageType>("CricoidCartilag");
355 p[0] = p[1] = p[2] = 0.0; // to avoid warning
356 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(CricoidCartilag,
357 this->GetBackgroundValue(), 2, true, p);
358 } catch (clitk::ExceptionObject e) {
359 //DD("CricoidCartilag image not found, try CricoidCartilagZ");
360 this->GetAFDB()->GetPoint3D("CricoidCartilagPoint", p);
362 output = clitk::CropImageRemoveGreaterThan<MaskImageType>(output, 2, p[2], true, this->GetBackgroundValue());
363 this->template StopCurrentStep<MaskImageType>(output);
366 //--------------------------------------------------------------------
368 this->StartNewStep("[Mediastinum] Keep main connected component");
369 output = clitk::Labelize<MaskImageType>(output, this->GetBackgroundValue(), false, 500);
370 // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
371 output = clitk::KeepLabels<MaskImageType>(output, this->GetBackgroundValue(),
372 this->GetForegroundValue(), 1, 1, 0);
373 this->template StopCurrentStep<MaskImageType>(output);
375 //--------------------------------------------------------------------
376 // Step: Remove post part from VertebralBody
377 this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody");
378 RemovePostPartOfVertebralBody();
379 this->template StopCurrentStep<MaskImageType>(output);
381 //--------------------------------------------------------------------
382 // Step: Slice by Slice CCL
383 this->StartNewStep("[Mediastinum] Slice by Slice keep only one component");
384 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
385 // typename ExtractSliceFilterType::Pointer
386 ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
387 extractSliceFilter->SetInput(output);
388 extractSliceFilter->SetDirection(2);
389 extractSliceFilter->Update();
390 typedef typename ExtractSliceFilterType::SliceType SliceType;
391 std::vector<typename SliceType::Pointer> mSlices;
392 extractSliceFilter->GetOutputSlices(mSlices);
393 for(unsigned int i=0; i<mSlices.size(); i++) {
394 mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
395 mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
397 typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
398 typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
399 joinFilter->SetOrigin(output->GetOrigin()[2]);
400 joinFilter->SetSpacing(output->GetSpacing()[2]);
401 for(unsigned int i=0; i<mSlices.size(); i++) {
402 joinFilter->PushBackInput(mSlices[i]);
404 joinFilter->Update();
405 output = joinFilter->GetOutput();
406 this->template StopCurrentStep<MaskImageType>(output);
408 //--------------------------------------------------------------------
409 // Step 10 : AutoCrop
410 this->StartNewStep("[Mediastinum] AutoCrop");
411 output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
412 this->template StopCurrentStep<MaskImageType>(output);
414 // End, set the real size
415 this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
416 this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
417 this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
418 this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
420 //--------------------------------------------------------------------
423 //--------------------------------------------------------------------
424 template <class ImageType>
426 clitk::ExtractMediastinumFilter<ImageType>::
429 this->GraftOutput(output);
430 // Store image filenames into AFDB
431 this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());
434 //--------------------------------------------------------------------
437 //--------------------------------------------------------------------
438 template <class ImageType>
440 clitk::ExtractMediastinumFilter<ImageType>::
441 RemovePostPartOfVertebralBody()
445 Posteriorly, Station 8 abuts the descending aorta and the anterior
446 aspect of the vertebral body until an imaginary horizontal line
447 running 1 cm posterior to the anterior border of the vertebral
450 => We use this definition for all the mediastinum
452 Find most Ant point in VertebralBody. Consider the horizontal line
453 which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from
457 // Get VertebralBody mask image
458 MaskImagePointer VertebralBody =
459 this->GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
461 // Consider vertebral body slice by slice
462 std::vector<MaskSlicePointer> vertebralSlices;
463 clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vertebralSlices);
465 // For each slice, compute the most anterior point
466 std::map<int, MaskSlicePointType> vertebralAntPositionBySlice;
467 for(uint i=0; i<vertebralSlices.size(); i++) {
468 MaskSlicePointType p;
469 bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vertebralSlices[i],
470 this->GetBackgroundValue(),
473 vertebralAntPositionBySlice[i] = p;
476 // It should not happen ! But sometimes, a contour is missing or
477 // the VertebralBody is not delineated enough inferiorly ... in
478 // those cases, we consider the first found slice.
479 // std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
480 // [ Possible alternative -> consider previous limit ]
484 // Convert 2D points in slice into 3D points
485 std::vector<MaskImagePointType> vertebralAntPositions;
486 clitk::PointsUtils<MaskImageType>::Convert2DMapTo3DList(vertebralAntPositionBySlice,
488 vertebralAntPositions);
490 // DEBUG : write list of points
491 clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
492 "VertebralBodyMostAnteriorPoints.txt");
494 // Cut support posteriorly 1cm the most anterior point of the
495 // VertebralBody. Do nothing below and above the VertebralBody. To
496 // do that compute several region, slice by slice and fill.
497 MaskImageRegionType region;
498 MaskImageSizeType size;
499 MaskImageIndexType start;
500 size[2] = 1; // a single slice
501 start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
502 size[0] = output->GetLargestPossibleRegion().GetSize()[0];
503 for(uint i=0; i<vertebralAntPositions.size(); i++) {
504 typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
506 IteratorType(output, output->GetLargestPossibleRegion());
507 MaskImageIndexType index;
508 // Consider some cm posterior to most anterior positions (usually
510 vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody();
511 // Get index of this point
512 output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
513 // Compute region (a single slice)
515 start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
516 size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
517 region.SetSize(size);
518 region.SetIndex(start);
520 if (output->GetLargestPossibleRegion().IsInside(start)) {
521 itk::ImageRegionIterator<MaskImageType> it(output, region);
523 while (!it.IsAtEnd()) {
524 it.Set(this->GetBackgroundValue());
530 //--------------------------------------------------------------------
533 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX