X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractMediastinumFilter.txx;h=b76e617946593c7faaca4924beb61f6d7a43e796;hb=3b5ab5fac2a4b3d4c752722687a04a0dd08e80b0;hp=b5535cf23af74e19f32275fca7ba5cc89e9fe8b1;hpb=e008d74b0ecdc4ca2eaae8c429901a78f9ef5c31;p=clitk.git diff --git a/segmentation/clitkExtractMediastinumFilter.txx b/segmentation/clitkExtractMediastinumFilter.txx index b5535cf..b76e617 100644 --- a/segmentation/clitkExtractMediastinumFilter.txx +++ b/segmentation/clitkExtractMediastinumFilter.txx @@ -23,10 +23,14 @@ #include "clitkCommon.h" #include "clitkExtractMediastinumFilter.h" #include "clitkAddRelativePositionConstraintToLabelImageFilter.h" -#include "clitkSegmentationFunctions.h" +#include "clitkSegmentationUtils.h" +#include "clitkExtractAirwaysTreeInfoFilter.h" +#include "clitkCropLikeImageFilter.h" -// itk +// std #include + +// itk #include "itkStatisticsLabelMapFilter.h" #include "itkLabelImageToStatisticsLabelMapFilter.h" #include "itkRegionOfInterestImageFilter.h" @@ -36,13 +40,14 @@ #include "RelativePositionPropImageFilter.h" //-------------------------------------------------------------------- -template -clitk::ExtractMediastinumFilter:: +template +clitk::ExtractMediastinumFilter:: ExtractMediastinumFilter(): clitk::FilterBase(), - itk::ImageToImageFilter() + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), + itk::ImageToImageFilter() { - this->SetNumberOfRequiredInputs(3); + this->SetNumberOfRequiredInputs(4); SetBackgroundValuePatient(0); SetBackgroundValueLung(0); @@ -53,30 +58,33 @@ ExtractMediastinumFilter(): SetForegroundValue(1); SetIntermediateSpacing(6); - SetFuzzyThreshold1(0.6); - SetFuzzyThreshold2(0.7); + SetFuzzyThreshold1(0.4); + SetFuzzyThreshold2(0.6); + SetFuzzyThreshold3(0.2); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -SetInputPatientLabelImage(const TImageType * image, ImagePixelType bg) { - this->SetNthInput(0, const_cast(image)); +clitk::ExtractMediastinumFilter:: +SetInputPatientLabelImage(const ImageType * image, ImagePixelType bg) +{ + this->SetNthInput(0, const_cast(image)); m_BackgroundValuePatient = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -SetInputLungLabelImage(const TImageType * image, ImagePixelType bg, - ImagePixelType fgLeft, ImagePixelType fgRight) { - this->SetNthInput(1, const_cast(image)); +clitk::ExtractMediastinumFilter:: +SetInputLungLabelImage(const ImageType * image, ImagePixelType bg, + ImagePixelType fgLeft, ImagePixelType fgRight) +{ + this->SetNthInput(1, const_cast(image)); m_BackgroundValueLung = bg; m_ForegroundValueLeftLung = fgLeft; m_ForegroundValueRightLung = fgRight; @@ -85,21 +93,34 @@ SetInputLungLabelImage(const TImageType * image, ImagePixelType bg, //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -SetInputBonesLabelImage(const TImageType * image, ImagePixelType bg) { - this->SetNthInput(2, const_cast(image)); +clitk::ExtractMediastinumFilter:: +SetInputBonesLabelImage(const ImageType * image, ImagePixelType bg) +{ + this->SetNthInput(2, const_cast(image)); m_BackgroundValueBones = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template +void +clitk::ExtractMediastinumFilter:: +SetInputTracheaLabelImage(const ImageType * image, ImagePixelType bg) +{ + this->SetNthInput(3, const_cast(image)); + m_BackgroundValueTrachea = bg; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template template void -clitk::ExtractMediastinumFilter:: +clitk::ExtractMediastinumFilter:: SetArgsInfo(ArgsInfoType mArgsInfo) { SetVerboseOption_GGO(mArgsInfo); @@ -109,7 +130,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetBackgroundValuePatient_GGO(mArgsInfo); SetBackgroundValueLung_GGO(mArgsInfo); - SetBackgroundValueBones_GGO(mArgsInfo); + SetBackgroundValueTrachea_GGO(mArgsInfo); SetForegroundValueLeftLung_GGO(mArgsInfo); SetForegroundValueRightLung_GGO(mArgsInfo); @@ -117,16 +138,19 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetIntermediateSpacing_GGO(mArgsInfo); SetFuzzyThreshold1_GGO(mArgsInfo); SetFuzzyThreshold2_GGO(mArgsInfo); + SetFuzzyThreshold3_GGO(mArgsInfo); + + SetAFDBFilename_GGO(mArgsInfo); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: +clitk::ExtractMediastinumFilter:: GenerateOutputInformation() { - ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); ImagePointer outputImage = this->GetOutput(0); outputImage->SetRegions(outputImage->GetLargestPossibleRegion()); } @@ -134,43 +158,69 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -GenerateInputRequestedRegion() { +clitk::ExtractMediastinumFilter:: +GenerateInputRequestedRegion() +{ // Call default Superclass::GenerateInputRequestedRegion(); - // Get input pointers - ImagePointer patient = dynamic_cast(itk::ProcessObject::GetInput(0)); - ImagePointer lung = dynamic_cast(itk::ProcessObject::GetInput(1)); - ImagePointer bones = dynamic_cast(itk::ProcessObject::GetInput(2)); + + // Get input pointers + LoadAFDB(); + ImagePointer patient = GetAFDB()->template GetImage ("patient"); + ImagePointer lung = GetAFDB()->template GetImage ("lungs"); + ImagePointer bones = GetAFDB()->template GetImage ("bones"); + ImagePointer trachea = GetAFDB()->template GetImage ("trachea"); patient->SetRequestedRegion(patient->GetLargestPossibleRegion()); lung->SetRequestedRegion(lung->GetLargestPossibleRegion()); bones->SetRequestedRegion(bones->GetLargestPossibleRegion()); + trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion()); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -GenerateData() { +clitk::ExtractMediastinumFilter:: +GenerateData() +{ // Get input pointers - ImageConstPointer patient = dynamic_cast(itk::ProcessObject::GetInput(0)); - ImageConstPointer lung = dynamic_cast(itk::ProcessObject::GetInput(1)); - ImageConstPointer bones = dynamic_cast(itk::ProcessObject::GetInput(2)); + ImagePointer patient = GetAFDB()->template GetImage ("patient"); + ImagePointer lung = GetAFDB()->template GetImage ("lungs"); + ImagePointer bones = GetAFDB()->template GetImage ("bones"); + ImagePointer trachea = GetAFDB()->template GetImage ("trachea"); // Get output pointer ImagePointer output; + // Step 0: Crop support (patient) to lung extend in RL + StartNewStep("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 0: Crop support (previous) to bones extend in AP + StartNewStep("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 1: patient minus lungs, minus bones StartNewStep("Patient contours minus lungs and minus bones"); typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); boolFilter->InPlaceOn(); - boolFilter->SetInput1(patient); + boolFilter->SetInput1(output); boolFilter->SetInput2(lung); boolFilter->SetOperationType(BoolFilterType::AndNot); boolFilter->Update(); @@ -180,79 +230,158 @@ GenerateData() { boolFilter->Update(); output = boolFilter->GetOutput(); + // Auto crop to gain support area output = clitk::AutoCrop(output, GetBackgroundValue()); - ////autoCropFilter->GetOutput(); typedef clitk::AutoCropFilter CropFilterType; - //typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - //cropFilter->SetInput(output); - //cropFilter->Update(); - //output = cropFilter->GetOutput(); - - this->template StopCurrentStep(output); + this->template StopCurrentStep(output); // Step 2: LR limits from lung (need separate lung ?) StartNewStep("Left/Right limits with lungs"); - typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + + /* + // WE DO NOT NEED THE FOLLOWING ? + // 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). + ImagePointer right_lung = clitk::SetBackground(lung, lung, 2, 0); + ImagePointer left_lung = clitk::SetBackground(lung, lung, 1, 0); + writeImage(right_lung, "right.mhd"); + writeImage(left_lung, "left.mhd"); + */ + + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); relPosFilter->SetInput(output); + //relPosFilter->SetInputObject(left_lung); relPosFilter->SetInputObject(lung); - relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); + relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); // warning left lung is at right ;) relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); relPosFilter->Update(); output = relPosFilter->GetOutput(); + // writeImage(right_lung, "step4-left.mhd"); relPosFilter->SetInput(output); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); + //relPosFilter->SetInputObject(right_lung); relPosFilter->SetInputObject(lung); relPosFilter->SetOrientationType(RelPosFilterType::RightTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); relPosFilter->Update(); output = relPosFilter->GetOutput(); - this->template StopCurrentStep(output); + this->template StopCurrentStep(output); // Step 3: AP limits from bones StartNewStep("Ant/Post limits with bones"); + ImageConstPointer bones_ant; + ImageConstPointer bones_post; + + // Find ant-post coordinate of trachea (assume the carena position is a + // good estimation of the ant-post position of the trachea) + ImagePointType carina; + LoadAFDB(); + GetAFDB()->GetPoint3D("Carina", carina); + DD(carina); + ImageIndexType index_trachea; + bones->TransformPhysicalPointToIndex(carina, index_trachea); + DD(index_trachea); + + // Split bone image first into two parts (ant and post) + typedef itk::RegionOfInterestImageFilter ROIFilterType; + // typedef itk::ExtractImageFilter ROIFilterType; + typename ROIFilterType::Pointer roiFilter = ROIFilterType::New(); + ImageRegionType region = bones->GetLargestPossibleRegion(); + ImageSizeType size = region.GetSize(); + DD(size); + size[1] = index_trachea[1]; //size[1]/2.0; + DD(size); + region.SetSize(size); + roiFilter->SetInput(bones); + // roiFilter->SetExtractionRegion(region); + roiFilter->SetRegionOfInterest(region); + roiFilter->ReleaseDataFlagOff(); + roiFilter->Update(); + bones_ant = roiFilter->GetOutput(); + writeImage(bones_ant, "b_ant.mhd"); + + // roiFilter->ResetPipeline();// = ROIFilterType::New(); + roiFilter = ROIFilterType::New(); + ImageIndexType index = region.GetIndex(); + index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1; + size[1] = bones->GetLargestPossibleRegion().GetSize()[1] - size[1]; + DD(size); + region.SetIndex(index); + region.SetSize(size); + roiFilter->SetInput(bones); + // roiFilter->SetExtractionRegion(region); + roiFilter->SetRegionOfInterest(region); + roiFilter->ReleaseDataFlagOff(); + roiFilter->Update(); + bones_post = roiFilter->GetOutput(); + writeImage(bones_post, "b_post.mhd"); + + // Go ! relPosFilter->SetCurrentStepNumber(0); relPosFilter->ResetPipeline();// = RelPosFilterType::New(); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); relPosFilter->SetInput(output); - relPosFilter->SetInputObject(bones); + relPosFilter->SetInputObject(bones_post); relPosFilter->SetOrientationType(RelPosFilterType::AntTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); relPosFilter->Update(); + output = relPosFilter->GetOutput(); + writeImage(output, "post.mhd"); relPosFilter->SetInput(relPosFilter->GetOutput()); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); - relPosFilter->SetInputObject(bones); + relPosFilter->SetInput(output); + relPosFilter->SetInputObject(bones_ant); relPosFilter->SetOrientationType(RelPosFilterType::PostTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); relPosFilter->Update(); output = relPosFilter->GetOutput(); - this->template StopCurrentStep(output); + this->template StopCurrentStep(output); // Get CCL - output = clitk::Labelize(output, GetBackgroundValue(), true, 100); - // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); - output = clitk::KeepLabels(output, GetBackgroundValue(), - GetForegroundValue(), 1, 1, 0); + StartNewStep("Keep main connected component"); + output = clitk::Labelize(output, GetBackgroundValue(), true, 500); + // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); + output = clitk::KeepLabels(output, GetBackgroundValue(), + GetForegroundValue(), 1, 1, 0); + this->template StopCurrentStep(output); + + // Step : Lower limits from lung (need separate lung ?) + StartNewStep("Lower limits with lungs"); + relPosFilter = RelPosFilterType::New(); + relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); + relPosFilter->VerboseStepOff(); + relPosFilter->WriteStepOff(); + relPosFilter->SetInput(output); + DD(output->GetLargestPossibleRegion().GetIndex()); + // relPosFilter->SetInputObject(left_lung); + relPosFilter->SetInputObject(lung); + relPosFilter->SetOrientationType(RelPosFilterType::SupTo); + relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3()); + relPosFilter->Update(); + output = relPosFilter->GetOutput(); + DD(output->GetLargestPossibleRegion()); output = clitk::AutoCrop(output, GetBackgroundValue()); - // cropFilter = CropFilterType::New(); - //cropFilter->SetInput(output); - //cropFilter->Update(); - //output = cropFilter->GetOutput(); + // roiFilter = ROIFilterType::New(); + //roiFilter->SetInput(output); + //roiFilter->Update(); + //output = roiFilter->GetOutput(); // Final Step -> set output this->SetNthOutput(0, output);