X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractMediastinumFilter.txx;h=675718082192ea3663ada290675da3648b91a5e2;hb=595624eb4297e747630105d45017de69733caef2;hp=d3ba557b97d863cb7a486aab5c59c742a4f371b0;hpb=573d80d0f7a17607d2ee883c21c940c0ba020282;p=clitk.git diff --git a/segmentation/clitkExtractMediastinumFilter.txx b/segmentation/clitkExtractMediastinumFilter.txx index d3ba557..6757180 100644 --- a/segmentation/clitkExtractMediastinumFilter.txx +++ b/segmentation/clitkExtractMediastinumFilter.txx @@ -45,27 +45,16 @@ template clitk::ExtractMediastinumFilter:: ExtractMediastinumFilter(): - clitk::FilterBase(), - clitk::FilterWithAnatomicalFeatureDatabaseManagement(), - itk::ImageToImageFilter() + clitk::StructuresExtractionFilter() { this->SetNumberOfRequiredInputs(1); - SetBackgroundValuePatient(0); SetBackgroundValueLung(0); SetBackgroundValueBones(0); SetForegroundValueLeftLung(1); SetForegroundValueRightLung(2); - SetBackgroundValue(0); - SetForegroundValue(1); - - SetIntermediateSpacing(6); - SetFuzzyThreshold1(0.5); - SetFuzzyThreshold2(0.6); - SetFuzzyThreshold3(0.05); - - SetOutputMediastinumFilename("mediastinum.mhd"); - + SetDistanceMaxToAnteriorPartOfTheVertebralBody(10); + SetOutputMediastinumFilename("mediastinum.mha"); UseBonesOff(); } //-------------------------------------------------------------------- @@ -121,6 +110,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 @@ -133,7 +148,6 @@ GenerateInputRequestedRegion() // Superclass::GenerateInputRequestedRegion(); // DD("End GenerateInputRequestedRegion"); } - //-------------------------------------------------------------------- @@ -148,6 +162,7 @@ SetInput(const ImageType * image) //-------------------------------------------------------------------- + //-------------------------------------------------------------------- template void @@ -158,22 +173,21 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- // Get input pointers - clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK - LoadAFDB(); + this->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 = this->GetAFDB()->template GetImage ("Patient"); + lung = this->GetAFDB()->template GetImage ("Lungs"); if (GetUseBones()) { - bones = GetAFDB()->template GetImage ("Bones"); + bones = this->GetAFDB()->template GetImage ("Bones"); } - MaskImagePointer trachea = GetAFDB()->template GetImage ("Trachea"); - - clitk::PrintMemory(GetVerboseMemoryFlag(), "After read patient, lung"); - + trachea = this->GetAFDB()->template GetImage ("Trachea"); + + //this->VerboseImageSizeFlagOn(); + //-------------------------------------------------------------------- - // Step 1: Crop support (patient) to lung extend in RL - StartNewStep("Crop support like lungs along LR"); + // 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); @@ -181,11 +195,11 @@ 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"); + 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 @@ -195,8 +209,8 @@ GenerateOutputInformation() { } //-------------------------------------------------------------------- - // Step 3: patient minus lungs, minus bones, minus trachea - StartNewStep("Patient contours minus lungs, trachea [and bones]"); + // 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(); @@ -217,70 +231,60 @@ GenerateOutputInformation() { output = boolFilter->GetOutput(); // Auto crop to gain support area - output = clitk::AutoCrop(output, GetBackgroundValue()); + output = clitk::AutoCrop(output, this->GetBackgroundValue()); 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 // NOT NEED TO SEPARATE ? ) - StartNewStep("Left/Right limits with lungs"); + 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); - right_lung = clitk::ResizeImageLike(right_lung, output, GetBackgroundValue()); - left_lung = clitk::ResizeImageLike(left_lung, output, GetBackgroundValue()); - // 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->VerboseStepFlagOff(); - 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->Update(); - output = relPosFilter->GetOutput(); - //writeImage(right_lung, "step4-left.mhd"); - - relPosFilter = RelPosFilterType::New(); - relPosFilter->SetInput(output); - relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(right_lung); - //relPosFilter->SetInputObject(lung); - relPosFilter->AddOrientationType(RelPosFilterType::AtLeftTo); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); - relPosFilter->Update(); - output = relPosFilter->GetOutput(); + 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 5: AP limits from bones + // Step : AP limits from bones // Separate the bones in the ant-post middle - MaskImageConstPointer bones_ant; - MaskImageConstPointer bones_post; - MaskImagePointType middle_AntPost__position; + MaskImagePointer bones_ant; + MaskImagePointer bones_post; + MaskImagePointType middle_AntPost_position; + middle_AntPost_position.Fill(NumericTraits::Zero); if (GetUseBones()) { - StartNewStep("Ant/Post limits with bones"); - 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; + 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; - bones->TransformPhysicalPointToIndex(middle_AntPost__position, 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; @@ -291,8 +295,8 @@ GenerateOutputInformation() { 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; + // 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); @@ -302,7 +306,7 @@ GenerateOutputInformation() { roiFilter->ReleaseDataFlagOff(); roiFilter->Update(); bones_ant = roiFilter->GetOutput(); - // writeImage(bones_ant, "b_ant.mhd"); + // writeImage(bones_ant, "b_ant.mhd"); // POST part roiFilter = ROIFilterType::New(); index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1; @@ -314,154 +318,69 @@ GenerateOutputInformation() { 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->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(bones_post); - relPosFilter->AddOrientationType(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->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(bones_ant); - relPosFilter->AddOrientationType(RelPosFilterType::PostTo); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); - relPosFilter->Update(); - output = relPosFilter->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); } //-------------------------------------------------------------------- - // Step 6: Get CCL - StartNewStep("Keep main connected component"); - output = clitk::Labelize(output, GetBackgroundValue(), false, 500); - // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); - output = clitk::KeepLabels(output, GetBackgroundValue(), - GetForegroundValue(), 1, 1, 0); + // 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); - //-------------------------------------------------------------------- - // 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 ? - + // Generic RelativePosition processes + output = this->ApplyRelativePositionList("Mediastinum", output); - //TODO - // Confidence connected ? + //-------------------------------------------------------------------- + // 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 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: 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 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 post part from VertebralBody + this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody"); + RemovePostPartOfVertebralBody(); + this->template StopCurrentStep(output); //-------------------------------------------------------------------- - // Step 10: Slice by Slice CCL - StartNewStep("Slice by Slice keep only one component"); + // 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(); @@ -486,59 +405,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->StartNewStep("[Mediastinum] AutoCrop"); + output = clitk::AutoCrop(output, this->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()); @@ -556,10 +428,106 @@ GenerateData() { this->GraftOutput(output); // Store image filenames into AFDB - GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename()); - WriteAFDB(); + 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