X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractMediastinumFilter.txx;h=67bc305d2ae0bf73f5808727fa3750d22c19deff;hb=3c86758765bc9bcba20d439424bcf97091b5af6f;hp=f9cd006d27f3d20b743ac360e17ebb3e5b1e6463;hpb=a339bdc482ea9752ec53195bc9a47e8b05dba582;p=clitk.git diff --git a/segmentation/clitkExtractMediastinumFilter.txx b/segmentation/clitkExtractMediastinumFilter.txx index f9cd006..67bc305 100644 --- a/segmentation/clitkExtractMediastinumFilter.txx +++ b/segmentation/clitkExtractMediastinumFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.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 @@ -14,7 +14,7 @@ - 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 @@ -50,7 +50,6 @@ ExtractMediastinumFilter(): itk::ImageToImageFilter() { this->SetNumberOfRequiredInputs(1); - SetBackgroundValuePatient(0); SetBackgroundValueLung(0); SetBackgroundValueBones(0); @@ -58,14 +57,12 @@ ExtractMediastinumFilter(): SetForegroundValueRightLung(2); SetBackgroundValue(0); SetForegroundValue(1); - SetIntermediateSpacing(6); - SetFuzzyThreshold1(0.5); - SetFuzzyThreshold2(0.6); - SetFuzzyThreshold3(0.05); - - SetOutputMediastinumFilename("mediastinum.mhd"); - + SetFuzzyThreshold("LR_lungs", 0.3); + SetFuzzyThreshold("bones", 0.6); + SetFuzzyThreshold("inf_lungs", 0.05); + SetDistanceMaxToAnteriorPartOfTheVertebralBody(10); + SetOutputMediastinumFilename("mediastinum.mhd"); UseBonesOff(); } //-------------------------------------------------------------------- @@ -121,6 +118,32 @@ SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType 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 @@ -158,21 +181,18 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- // Get input pointers - clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK LoadAFDB(); ImageConstPointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); - MaskImagePointer patient = GetAFDB()->template GetImage ("Patient"); - MaskImagePointer lung = GetAFDB()->template GetImage ("Lungs"); - MaskImagePointer bones; + MaskImagePointer patient, lung, bones, trachea; + patient = GetAFDB()->template GetImage ("Patient"); + lung = GetAFDB()->template GetImage ("Lungs"); if (GetUseBones()) { bones = GetAFDB()->template GetImage ("Bones"); } - MaskImagePointer trachea = GetAFDB()->template GetImage ("Trachea"); - - clitk::PrintMemory(GetVerboseMemoryFlag(), "After read patient, lung"); + trachea = GetAFDB()->template GetImage ("Trachea"); //-------------------------------------------------------------------- - // Step 1: Crop support (patient) to lung extend in RL + // Step : 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(); @@ -181,9 +201,9 @@ GenerateOutputInformation() { cropFilter->Update(); output = cropFilter->GetOutput(); this->template StopCurrentStep(output); - + //-------------------------------------------------------------------- - // Step 2: Crop support (previous) to bones extend in AP + // Step : Crop support (previous) to bones extend in AP if (GetUseBones()) { StartNewStep("Crop support like bones along AP"); cropFilter = CropFilterType::New(); @@ -195,7 +215,7 @@ GenerateOutputInformation() { } //-------------------------------------------------------------------- - // Step 3: patient minus lungs, minus bones, minus trachea + // Step : patient minus lungs, minus bones, minus trachea StartNewStep("Patient contours minus lungs, trachea [and bones]"); typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); @@ -221,7 +241,7 @@ GenerateOutputInformation() { this->template StopCurrentStep(output); //-------------------------------------------------------------------- - // Step 4: LR limits from lung (need separate lung ?) + // 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 @@ -245,13 +265,11 @@ GenerateOutputInformation() { relPosFilter->WriteStepFlagOff(); relPosFilter->SetInput(output); relPosFilter->SetInputObject(left_lung); - // relPosFilter->SetInputObject(lung); relPosFilter->AddOrientationType(RelPosFilterType::AtRightTo); // warning left lung is at right ;) relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("LR_lungs")); relPosFilter->Update(); output = relPosFilter->GetOutput(); - //writeImage(right_lung, "step4-left.mhd"); relPosFilter = RelPosFilterType::New(); relPosFilter->SetInput(output); @@ -260,16 +278,26 @@ GenerateOutputInformation() { relPosFilter->WriteStepFlagOff(); relPosFilter->SetInput(output); relPosFilter->SetInputObject(right_lung); - //relPosFilter->SetInputObject(lung); relPosFilter->AddOrientationType(RelPosFilterType::AtLeftTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("LR_lungs")); relPosFilter->Update(); output = relPosFilter->GetOutput(); this->template StopCurrentStep(output); //-------------------------------------------------------------------- - // Step 5: AP limits from bones + // Step : superior limits + StartNewStep("Keep inferior to CricoidCartilag"); + // load Cricoid, get centroid, cut above (or below), lower bound + MaskImagePointer CricoidCartilag = GetAFDB()->template GetImage ("CricoidCartilag"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(CricoidCartilag, GetBackgroundValue(), 2, true, p); + output = clitk::CropImageRemoveGreaterThan(output, 2, p[2], true, GetBackgroundValue()); + this->template StopCurrentStep(output); + + //-------------------------------------------------------------------- + // Step : AP limits from bones // Separate the bones in the ant-post middle MaskImageConstPointer bones_ant; MaskImageConstPointer bones_post; @@ -326,7 +354,7 @@ GenerateOutputInformation() { relPosFilter->SetInputObject(bones_post); relPosFilter->AddOrientationType(RelPosFilterType::AntTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("bones")); relPosFilter->Update(); output = relPosFilter->GetOutput(); // writeImage(output, "post.mhd"); @@ -339,14 +367,14 @@ GenerateOutputInformation() { relPosFilter->SetInputObject(bones_ant); relPosFilter->AddOrientationType(RelPosFilterType::PostTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("bones")); relPosFilter->Update(); output = relPosFilter->GetOutput(); this->template StopCurrentStep(output); } //-------------------------------------------------------------------- - // Step 6: Get CCL + // Step: Get CCL StartNewStep("Keep main connected component"); output = clitk::Labelize(output, GetBackgroundValue(), false, 500); // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); @@ -354,113 +382,23 @@ GenerateOutputInformation() { GetForegroundValue(), 1, 1, 0); this->template StopCurrentStep(output); - - //-------------------------------------------------------------------- - // Step 8: Trial segmentation KMeans - if (0) { - StartNewStep("K means"); - // Take input, crop like current mask - typedef CropLikeImageFilter CropLikeFilterType; - typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New(); - cropLikeFilter->SetInput(input); - cropLikeFilter->SetCropLikeImage(output); - cropLikeFilter->Update(); - ImagePointer working_input = cropLikeFilter->GetOutput(); - writeImage(working_input, "crop-input.mhd"); - // Set bG at -1000 - working_input = clitk::SetBackground(working_input, output, GetBackgroundValue(), -1000, true); - writeImage(working_input, "crop-input2.mhd"); - // Kmeans - typedef itk::ScalarImageKmeansImageFilter KMeansFilterType; - typename KMeansFilterType::Pointer kmeansFilter = KMeansFilterType::New(); - kmeansFilter->SetInput(working_input); - // const unsigned int numberOfInitialClasses = 3; - // const unsigned int useNonContiguousLabels = 0; - kmeansFilter->AddClassWithInitialMean(-1000); - kmeansFilter->AddClassWithInitialMean(30); - kmeansFilter->AddClassWithInitialMean(-40); // ==> I want this one - DD("Go!"); - kmeansFilter->Update(); - DD("End"); - typename KMeansFilterType::ParametersType estimatedMeans = kmeansFilter->GetFinalMeans(); - const unsigned int numberOfClasses = estimatedMeans.Size(); - for ( unsigned int i = 0 ; i < numberOfClasses ; ++i ) { - std::cout << "cluster[" << i << "] "; - std::cout << " estimated mean : " << estimatedMeans[i] << std::endl; - } - MaskImageType::Pointer kmeans = kmeansFilter->GetOutput(); - kmeans = clitk::SetBackground(kmeans, kmeans, - 1, GetBackgroundValue(), true); - writeImage(kmeans, "kmeans.mhd"); - // Get final results, and remove from current mask - boolFilter = BoolFilterType::New(); - boolFilter->InPlaceOn(); - boolFilter->SetInput1(output); - boolFilter->SetInput2(kmeans); - boolFilter->SetOperationType(BoolFilterType::And); - boolFilter->Update(); - output = boolFilter->GetOutput(); - writeImage(output, "out-kmean.mhd"); - this->template StopCurrentStep(output); - - // TODO -> FillMASK ? - // comment speed ? mask ? 2 class ? - - - //TODO - // Confidence connected ? - - } - //-------------------------------------------------------------------- - // Step 8: Lower limits from lung (need separate lung ?) - if (0) { - // StartNewStep("Trial : minus segmented struct"); - // MaskImagePointer heart = GetAFDB()->template GetImage ("heart"); - // boolFilter = BoolFilterType::New(); - // boolFilter->InPlaceOn(); - // boolFilter->SetInput1(output); - // boolFilter->SetInput2(heart); - // boolFilter->SetOperationType(BoolFilterType::AndNot); - // boolFilter->Update(); - // output = boolFilter->GetOutput(); // not needed because InPlace - - // Not below the heart - // relPosFilter = RelPosFilterType::New(); - // relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - // relPosFilter->VerboseStepFlagOff(); - // relPosFilter->WriteStepFlagOff(); - // relPosFilter->SetInput(output); - // relPosFilter->SetInputObject(heart); - // relPosFilter->SetOrientationType(RelPosFilterType::SupTo); - // relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - // relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3()); - // relPosFilter->Update(); - // output = relPosFilter->GetOutput(); - } + // Step: Remove post part from VertebralBody + StartNewStep("Remove post part according to VertebralBody"); + RemovePostPartOfVertebralBody(); + this->template StopCurrentStep(output); //-------------------------------------------------------------------- - // Step 8: Lower limits from lung (need separate lung ?) - if (0) { - StartNewStep("Lower limits with lungs"); - // TODO BOFFF ???? - relPosFilter = RelPosFilterType::New(); - relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - // relPosFilter->SetInputObject(left_lung); - relPosFilter->SetInputObject(lung); - relPosFilter->AddOrientationType(RelPosFilterType::SupTo); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3()); - relPosFilter->Update(); - output = relPosFilter->GetOutput(); - this->template StopCurrentStep(output); - } + // Step: Remove ant part according to Sternum + StartNewStep("Remove ant part according to Sternum"); + MaskImagePointer Sternum = GetAFDB()->template GetImage ("Sternum"); + output = clitk::SliceBySliceRelativePosition(output, Sternum, 2, + GetFuzzyThreshold("ant_sternum"), + "PostTo", false, 3, true, false); + this->template StopCurrentStep(output); //-------------------------------------------------------------------- - // Step 10: Slice by Slice CCL + // Step: Slice by Slice CCL StartNewStep("Slice by Slice keep only one component"); typedef clitk::ExtractSliceFilter ExtractSliceFilterType; // typename ExtractSliceFilterType::Pointer @@ -486,59 +424,12 @@ GenerateOutputInformation() { output = joinFilter->GetOutput(); this->template StopCurrentStep(output); - //-------------------------------------------------------------------- - // Step 9: Binarize to remove too high HU - // --> warning CCL slice by slice must be done before - if (0) { - StartNewStep("Remove hypersignal (bones and injected part"); - // Crop initial ct like current support - typedef CropLikeImageFilter CropLikeFilterType; - typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New(); - cropLikeFilter->SetInput(input); - cropLikeFilter->SetCropLikeImage(output); - cropLikeFilter->Update(); - ImagePointer working_input = cropLikeFilter->GetOutput(); - // writeImage(working_input, "crop-ct.mhd"); - // Binarize - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; - typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New(); - binarizeFilter->SetInput(working_input); - binarizeFilter->SetLowerThreshold(GetLowerThreshold()); - binarizeFilter->SetUpperThreshold(GetUpperThreshold()); - binarizeFilter->SetInsideValue(this->GetBackgroundValue()); // opposite - binarizeFilter->SetOutsideValue(this->GetForegroundValue()); // opposite - binarizeFilter->Update(); - MaskImagePointer working_bin = binarizeFilter->GetOutput(); - // writeImage(working_bin, "bin.mhd"); - // Remove from support - boolFilter = BoolFilterType::New(); - boolFilter->InPlaceOn(); - boolFilter->SetInput1(output); - boolFilter->SetInput2(working_bin); - boolFilter->SetOperationType(BoolFilterType::AndNot); - boolFilter->Update(); - output = boolFilter->GetOutput(); - StopCurrentStep(output); - } - //-------------------------------------------------------------------- // Step 10 : AutoCrop StartNewStep("AutoCrop"); output = clitk::AutoCrop(output, GetBackgroundValue()); this->template StopCurrentStep(output); - // Bones ? pb with RAM ? FillHoles ? - - // how to do with post part ? spine /lung ? - // POST the spine (should be separated from the rest) - /// DO THAT ---->> - // histo Y on the whole bones_post (3D) -> result is the Y center on the spine (?) - // by slice on bones_post - // find the most ant point in the center - // from this point go to post until out of bones. - // - - // End, set the real size this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion()); @@ -560,6 +451,102 @@ GenerateData() 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 = + 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], + 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(GetBackgroundValue()); + ++it; + } + } + } +} +//-------------------------------------------------------------------- + #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX