/*========================================================================= Program: vv http://www.creatis.insa-lyon.fr/rio/vv Authors belong to: - University of LYON http://www.universite-lyon.fr/ - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the copyright notices for more information. It is distributed under dual licence - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html ===========================================================================**/ #ifndef CLITKEXTRACTMEDIASTINUMFILTER_TXX #define CLITKEXTRACTMEDIASTINUMFILTER_TXX // clitk #include "clitkCommon.h" #include "clitkExtractMediastinumFilter.h" #include "clitkAddRelativePositionConstraintToLabelImageFilter.h" #include "clitkSliceBySliceRelativePositionFilter.h" #include "clitkSegmentationUtils.h" #include "clitkExtractAirwaysTreeInfoFilter.h" #include "clitkCropLikeImageFilter.h" // std #include // itk #include "itkStatisticsLabelMapFilter.h" #include "itkLabelImageToStatisticsLabelMapFilter.h" #include "itkRegionOfInterestImageFilter.h" #include "itkBinaryThresholdImageFilter.h" #include "itkScalarImageKmeansImageFilter.h" // itk ENST #include "RelativePositionPropImageFilter.h" //-------------------------------------------------------------------- template clitk::ExtractMediastinumFilter:: ExtractMediastinumFilter(): clitk::StructuresExtractionFilter() { this->SetNumberOfRequiredInputs(1); SetBackgroundValuePatient(0); SetBackgroundValueLung(0); SetBackgroundValueBones(0); SetForegroundValueLeftLung(1); SetForegroundValueRightLung(2); SetDistanceMaxToAnteriorPartOfTheVertebralBody(10); SetOutputMediastinumFilename("mediastinum.mha"); UseBonesOff(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg) { this->SetNthInput(0, const_cast(image)); m_BackgroundValuePatient = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg, MaskImagePixelType fgLeft, MaskImagePixelType fgRight) { this->SetNthInput(1, const_cast(image)); m_BackgroundValueLung = bg; m_ForegroundValueLeftLung = fgLeft; m_ForegroundValueRightLung = fgRight; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg) { this->SetNthInput(2, const_cast(image)); m_BackgroundValueBones = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg) { this->SetNthInput(3, const_cast(image)); m_BackgroundValueTrachea = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: SetFuzzyThreshold(std::string tag, double value) { m_FuzzyThreshold[tag] = value; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template double clitk::ExtractMediastinumFilter:: GetFuzzyThreshold(std::string tag) { if (m_FuzzyThreshold.find(tag) == m_FuzzyThreshold.end()) { clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThresholds."); return 0.0; } return m_FuzzyThreshold[tag]; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: GenerateInputRequestedRegion() { // DD("GenerateInputRequestedRegion"); // Do not call default // Superclass::GenerateInputRequestedRegion(); // DD("End GenerateInputRequestedRegion"); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: SetInput(const ImageType * image) { this->SetNthInput(0, const_cast(image)); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: GenerateOutputInformation() { // Do not call default // Superclass::GenerateOutputInformation(); //-------------------------------------------------------------------- // Get input pointers this->LoadAFDB(); ImageConstPointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); MaskImagePointer patient, lung, bones, trachea; patient = this->GetAFDB()->template GetImage ("Patient"); lung = this->GetAFDB()->template GetImage ("Lungs"); if (GetUseBones()) { bones = this->GetAFDB()->template GetImage ("Bones"); } trachea = this->GetAFDB()->template GetImage ("Trachea"); //this->VerboseImageSizeFlagOn(); //-------------------------------------------------------------------- // Step : Crop support (patient) to lung extend in RL this->StartNewStep("[Mediastinum] Crop support like lungs along LR"); typedef clitk::CropLikeImageFilter CropFilterType; typename CropFilterType::Pointer cropFilter = CropFilterType::New(); cropFilter->SetInput(patient); cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe cropFilter->Update(); output = cropFilter->GetOutput(); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step : Crop support (previous) to bones extend in AP if (GetUseBones()) { this->StartNewStep("[Mediastinum] Crop support like bones along AP"); cropFilter = CropFilterType::New(); cropFilter->SetInput(output); cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe cropFilter->Update(); output = cropFilter->GetOutput(); this->template StopCurrentStep(output); } //-------------------------------------------------------------------- // Step : patient minus lungs, minus bones, minus trachea this->StartNewStep("[Mediastinum] Patient contours minus lungs, trachea [and bones]"); typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); boolFilter->InPlaceOn(); boolFilter->SetInput1(output); boolFilter->SetInput2(lung); boolFilter->SetOperationType(BoolFilterType::AndNot); boolFilter->Update(); if (GetUseBones()) { boolFilter->SetInput1(boolFilter->GetOutput()); boolFilter->SetInput2(bones); boolFilter->SetOperationType(BoolFilterType::AndNot); boolFilter->Update(); } boolFilter->SetInput1(boolFilter->GetOutput()); boolFilter->SetInput2(trachea); boolFilter->SetOperationType(BoolFilterType::AndNot); boolFilter->Update(); output = boolFilter->GetOutput(); // Auto crop to gain support area output = clitk::AutoCrop(output, this->GetBackgroundValue()); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step : LR limits from lung (need separate lung ?) // Get separate lung images to get only the right and left lung // (because RelativePositionPropImageFilter only consider fg=1); // (label must be '1' because right is greater than left). (WE DO // NOT NEED TO SEPARATE ? ) this->StartNewStep("[Mediastinum] Left/Right limits with lungs"); // The following cannot be "inplace" because mask is the same than input ... MaskImagePointer right_lung = clitk::SetBackground(lung, lung, 2, 0, false); MaskImagePointer left_lung = clitk::SetBackground(lung, lung, 1, 0, false); left_lung = clitk::SetBackground(left_lung, left_lung, 2, 1, false); right_lung = clitk::ResizeImageLike(right_lung, output, this->GetBackgroundValue()); left_lung = clitk::ResizeImageLike(left_lung, output, this->GetBackgroundValue()); this->GetAFDB()->template SetImage("RightLung", "seg/RightLung.mha", right_lung, true); this->GetAFDB()->template SetImage("LeftLung", "seg/LeftLung.mha", left_lung, true); this->GetAFDB()->Write(); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step : AP limits from bones // Separate the bones in the ant-post middle MaskImagePointer bones_ant; MaskImagePointer bones_post; MaskImagePointType middle_AntPost_position; middle_AntPost_position.Fill(NumericTraits::Zero); if (GetUseBones()) { this->StartNewStep("[Mediastinum] Ant/Post limits with bones"); // To define ant and post part of the bones with a single horizontal line MaskImageIndexType index_bones_middle; // Method1: cut in the middle (not optimal) /* middle_AntPost_position[0] = middle_AntPost_position[2] = 0; middle_AntPost_position[1] = bones->GetOrigin()[1]+ (bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0; DD(middle_AntPost_position); bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle); */ // Method2: Use VertebralBody, take most ant point MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage("VertebralBody"); FindExtremaPointInAGivenDirection(VertebralBody, this->GetBackgroundValue(), 1, true, middle_AntPost_position); bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle); // Split bone image first into two parts (ant and post), and crop // lateraly to get vertebral typedef itk::RegionOfInterestImageFilter ROIFilterType; // typedef itk::ExtractImageFilter ROIFilterType; typename ROIFilterType::Pointer roiFilter = ROIFilterType::New(); MaskImageRegionType region = bones->GetLargestPossibleRegion(); MaskImageSizeType size = region.GetSize(); MaskImageIndexType index = region.GetIndex(); // ANT part // crop LR to keep 1/4 center part // index[0] = size[0]/4+size[0]/8; //size[0] = size[0]/4; // crop AP to keep first (ant) part size[1] = index_bones_middle[1]; //size[1]/2.0; region.SetSize(size); region.SetIndex(index); roiFilter->SetInput(bones); roiFilter->SetRegionOfInterest(region); roiFilter->ReleaseDataFlagOff(); roiFilter->Update(); bones_ant = roiFilter->GetOutput(); // writeImage(bones_ant, "b_ant.mhd"); // POST part roiFilter = ROIFilterType::New(); index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1; size[1] = bones->GetLargestPossibleRegion().GetSize()[1] - size[1]; region.SetIndex(index); region.SetSize(size); roiFilter->SetInput(bones); roiFilter->SetRegionOfInterest(region); roiFilter->ReleaseDataFlagOff(); roiFilter->Update(); bones_post = roiFilter->GetOutput(); // writeImage(bones_post, "b_post.mhd"); // Now, insert this image in the AFDB ==> (needed because used in the RelativePosition config file) this->GetAFDB()->template SetImage("Bones_Post", "seg/Bones_Post.mha", bones_post, true); this->GetAFDB()->template SetImage("Bones_Ant", "seg/Bones_Ant.mha", bones_ant, true); this->GetAFDB()->Write(); this->template StopCurrentStep(output); } //-------------------------------------------------------------------- // Remove VertebralBody part this->StartNewStep("[Mediastinum] Remove VertebralBody"); MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage("VertebralBody"); clitk::AndNot(output, VertebralBody, this->GetBackgroundValue()); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Generic RelativePosition processes output = this->ApplyRelativePositionList("Mediastinum", output); //-------------------------------------------------------------------- // FIXME --> do not put this limits here ! /* // Step : SI limits It is better to do this limit *AFTER* the // RelativePosition to avoid some issue due to superior boundaries. this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag"); // load Cricoid, get centroid, cut above (or below), lower bound MaskImagePointType p; try { MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage ("CricoidCartilag"); p[0] = p[1] = p[2] = 0.0; // to avoid warning clitk::FindExtremaPointInAGivenDirection(CricoidCartilag, this->GetBackgroundValue(), 2, true, p); } catch (clitk::ExceptionObject e) { //DD("CricoidCartilag image not found, try CricoidCartilagZ"); this->GetAFDB()->GetPoint3D("CricoidCartilagPoint", p); } output = clitk::CropImageRemoveGreaterThan(output, 2, p[2], true, this->GetBackgroundValue()); this->template StopCurrentStep(output); */ //-------------------------------------------------------------------- // Step: Get CCL this->StartNewStep("[Mediastinum] Keep main connected component"); output = clitk::Labelize(output, this->GetBackgroundValue(), false, 500); // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); output = clitk::KeepLabels(output, this->GetBackgroundValue(), this->GetForegroundValue(), 1, 1, 0); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step: Remove post part from VertebralBody this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody"); RemovePostPartOfVertebralBody(); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step: Slice by Slice CCL this->StartNewStep("[Mediastinum] Slice by Slice keep only one component"); typedef clitk::ExtractSliceFilter ExtractSliceFilterType; // typename ExtractSliceFilterType::Pointer ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); extractSliceFilter->SetInput(output); extractSliceFilter->SetDirection(2); extractSliceFilter->Update(); typedef typename ExtractSliceFilterType::SliceType SliceType; std::vector mSlices; extractSliceFilter->GetOutputSlices(mSlices); for(unsigned int i=0; i(mSlices[i], 0, true, 100); mSlices[i] = KeepLabels(mSlices[i], 0, 1, 1, 1, true); } typedef itk::JoinSeriesImageFilter JoinSeriesFilterType; typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New(); joinFilter->SetOrigin(output->GetOrigin()[2]); joinFilter->SetSpacing(output->GetSpacing()[2]); for(unsigned int i=0; iPushBackInput(mSlices[i]); } joinFilter->Update(); output = joinFilter->GetOutput(); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step 10 : AutoCrop this->StartNewStep("[Mediastinum] AutoCrop"); output = clitk::AutoCrop(output, this->GetBackgroundValue()); this->template StopCurrentStep(output); // End, set the real size this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion()); this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion()); this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion()); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: GenerateData() { this->GraftOutput(output); // Store image filenames into AFDB this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename()); this->WriteAFDB(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template void clitk::ExtractMediastinumFilter:: RemovePostPartOfVertebralBody() { /* Posteriorly, Station 8 abuts the descending aorta and the anterior aspect of the vertebral body until an imaginary horizontal line running 1 cm posterior to the anterior border of the vertebral body (Fig. 3C). => We use this definition for all the mediastinum Find most Ant point in VertebralBody. Consider the horizontal line which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from the most ant point. */ // Get VertebralBody mask image MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage ("VertebralBody"); // Consider vertebral body slice by slice std::vector vertebralSlices; clitk::ExtractSlices(VertebralBody, 2, vertebralSlices); // For each slice, compute the most anterior point std::map vertebralAntPositionBySlice; for(uint i=0; i(vertebralSlices[i], this->GetBackgroundValue(), 1, true, p); if (found) { vertebralAntPositionBySlice[i] = p; } else { // It should not happen ! But sometimes, a contour is missing or // the VertebralBody is not delineated enough inferiorly ... in // those cases, we consider the first found slice. // std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl; // [ Possible alternative -> consider previous limit ] } } // Convert 2D points in slice into 3D points std::vector vertebralAntPositions; clitk::PointsUtils::Convert2DMapTo3DList(vertebralAntPositionBySlice, VertebralBody, vertebralAntPositions); // DEBUG : write list of points clitk::WriteListOfLandmarks(vertebralAntPositions, "VertebralBodyMostAnteriorPoints.txt"); // Cut support posteriorly 1cm the most anterior point of the // VertebralBody. Do nothing below and above the VertebralBody. To // do that compute several region, slice by slice and fill. MaskImageRegionType region; MaskImageSizeType size; MaskImageIndexType start; size[2] = 1; // a single slice start[0] = output->GetLargestPossibleRegion().GetIndex()[0]; size[0] = output->GetLargestPossibleRegion().GetSize()[0]; for(uint i=0; i IteratorType; IteratorType iter = IteratorType(output, output->GetLargestPossibleRegion()); MaskImageIndexType index; // Consider some cm posterior to most anterior positions (usually // 1 cm). vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody(); // Get index of this point output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index); // Compute region (a single slice) start[2] = index[2]; start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1]; size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1]; region.SetSize(size); region.SetIndex(start); // Fill region if (output->GetLargestPossibleRegion().IsInside(start)) { itk::ImageRegionIterator it(output, region); it.GoToBegin(); while (!it.IsAtEnd()) { it.Set(this->GetBackgroundValue()); ++it; } } } } //-------------------------------------------------------------------- #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX