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;
267 this->StartNewStep("[Mediastinum] Ant/Post limits with bones");
269 // To define ant and post part of the bones with a single horizontal line
270 MaskImageIndexType index_bones_middle;
272 // Method1: cut in the middle (not optimal)
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);
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);
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();
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();
307 bones_ant = roiFilter->GetOutput();
308 // writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
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();
319 bones_post = roiFilter->GetOutput();
320 // writeImage<MaskImageType>(bones_post, "b_post.mhd");
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",
325 this->GetAFDB()->template SetImage<MaskImageType>("Bones_Ant", "seg/Bones_Ant.mha",
327 this->GetAFDB()->Write();
329 this->template StopCurrentStep<MaskImageType>(output);
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);
339 //--------------------------------------------------------------------
340 // Generic RelativePosition processes
341 output = this->ApplyRelativePositionList("Mediastinum", output);
344 //--------------------------------------------------------------------
345 // FIXME --> do not put this limits here !
347 // Step : SI limits It is better to do this limit *AFTER* the
348 // RelativePosition to avoid some issue due to superior boundaries.
349 this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag");
350 // load Cricoid, get centroid, cut above (or below), lower bound
351 MaskImagePointType p;
353 MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage <MaskImageType>("CricoidCartilag");
354 p[0] = p[1] = p[2] = 0.0; // to avoid warning
355 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(CricoidCartilag,
356 this->GetBackgroundValue(), 2, true, p);
357 } catch (clitk::ExceptionObject e) {
358 //DD("CricoidCartilag image not found, try CricoidCartilagZ");
359 this->GetAFDB()->GetPoint3D("CricoidCartilagPoint", p);
361 output = clitk::CropImageRemoveGreaterThan<MaskImageType>(output, 2, p[2], true, this->GetBackgroundValue());
362 this->template StopCurrentStep<MaskImageType>(output);
365 //--------------------------------------------------------------------
367 this->StartNewStep("[Mediastinum] Keep main connected component");
368 output = clitk::Labelize<MaskImageType>(output, this->GetBackgroundValue(), false, 500);
369 // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
370 output = clitk::KeepLabels<MaskImageType>(output, this->GetBackgroundValue(),
371 this->GetForegroundValue(), 1, 1, 0);
372 this->template StopCurrentStep<MaskImageType>(output);
374 //--------------------------------------------------------------------
375 // Step: Remove post part from VertebralBody
376 this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody");
377 RemovePostPartOfVertebralBody();
378 this->template StopCurrentStep<MaskImageType>(output);
380 //--------------------------------------------------------------------
381 // Step: Slice by Slice CCL
382 this->StartNewStep("[Mediastinum] Slice by Slice keep only one component");
383 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
384 // typename ExtractSliceFilterType::Pointer
385 ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
386 extractSliceFilter->SetInput(output);
387 extractSliceFilter->SetDirection(2);
388 extractSliceFilter->Update();
389 typedef typename ExtractSliceFilterType::SliceType SliceType;
390 std::vector<typename SliceType::Pointer> mSlices;
391 extractSliceFilter->GetOutputSlices(mSlices);
392 for(unsigned int i=0; i<mSlices.size(); i++) {
393 mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
394 mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
396 typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
397 typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
398 joinFilter->SetOrigin(output->GetOrigin()[2]);
399 joinFilter->SetSpacing(output->GetSpacing()[2]);
400 for(unsigned int i=0; i<mSlices.size(); i++) {
401 joinFilter->PushBackInput(mSlices[i]);
403 joinFilter->Update();
404 output = joinFilter->GetOutput();
405 this->template StopCurrentStep<MaskImageType>(output);
407 //--------------------------------------------------------------------
408 // Step 10 : AutoCrop
409 this->StartNewStep("[Mediastinum] AutoCrop");
410 output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
411 this->template StopCurrentStep<MaskImageType>(output);
413 // End, set the real size
414 this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
415 this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
416 this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
417 this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
419 //--------------------------------------------------------------------
422 //--------------------------------------------------------------------
423 template <class ImageType>
425 clitk::ExtractMediastinumFilter<ImageType>::
428 this->GraftOutput(output);
429 // Store image filenames into AFDB
430 this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());
433 //--------------------------------------------------------------------
436 //--------------------------------------------------------------------
437 template <class ImageType>
439 clitk::ExtractMediastinumFilter<ImageType>::
440 RemovePostPartOfVertebralBody()
444 Posteriorly, Station 8 abuts the descending aorta and the anterior
445 aspect of the vertebral body until an imaginary horizontal line
446 running 1 cm posterior to the anterior border of the vertebral
449 => We use this definition for all the mediastinum
451 Find most Ant point in VertebralBody. Consider the horizontal line
452 which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from
456 // Get VertebralBody mask image
457 MaskImagePointer VertebralBody =
458 this->GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
460 // Consider vertebral body slice by slice
461 std::vector<MaskSlicePointer> vertebralSlices;
462 clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vertebralSlices);
464 // For each slice, compute the most anterior point
465 std::map<int, MaskSlicePointType> vertebralAntPositionBySlice;
466 for(uint i=0; i<vertebralSlices.size(); i++) {
467 MaskSlicePointType p;
468 bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vertebralSlices[i],
469 this->GetBackgroundValue(),
472 vertebralAntPositionBySlice[i] = p;
475 // It should not happen ! But sometimes, a contour is missing or
476 // the VertebralBody is not delineated enough inferiorly ... in
477 // those cases, we consider the first found slice.
478 // std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
479 // [ Possible alternative -> consider previous limit ]
483 // Convert 2D points in slice into 3D points
484 std::vector<MaskImagePointType> vertebralAntPositions;
485 clitk::PointsUtils<MaskImageType>::Convert2DMapTo3DList(vertebralAntPositionBySlice,
487 vertebralAntPositions);
489 // DEBUG : write list of points
490 clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
491 "VertebralBodyMostAnteriorPoints.txt");
493 // Cut support posteriorly 1cm the most anterior point of the
494 // VertebralBody. Do nothing below and above the VertebralBody. To
495 // do that compute several region, slice by slice and fill.
496 MaskImageRegionType region;
497 MaskImageSizeType size;
498 MaskImageIndexType start;
499 size[2] = 1; // a single slice
500 start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
501 size[0] = output->GetLargestPossibleRegion().GetSize()[0];
502 for(uint i=0; i<vertebralAntPositions.size(); i++) {
503 typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
505 IteratorType(output, output->GetLargestPossibleRegion());
506 MaskImageIndexType index;
507 // Consider some cm posterior to most anterior positions (usually
509 vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody();
510 // Get index of this point
511 output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
512 // Compute region (a single slice)
514 start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
515 size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
516 region.SetSize(size);
517 region.SetIndex(start);
519 if (output->GetLargestPossibleRegion().IsInside(start)) {
520 itk::ImageRegionIterator<MaskImageType> it(output, region);
522 while (!it.IsAtEnd()) {
523 it.Set(this->GetBackgroundValue());
529 //--------------------------------------------------------------------
532 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX