X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLymphStation_8.txx;h=fc39203ebf7ee646fbc2b169adfad864d9aabd92;hb=431fea9a7c825943b07fc4b4f3423355a5fee501;hp=469ff2f836f38a053739151f642c79dc461c8163;hpb=9de203e39ab0d5c9779bc0827abdb8b5d8a58975;p=clitk.git diff --git a/segmentation/clitkExtractLymphStation_8.txx b/segmentation/clitkExtractLymphStation_8.txx index 469ff2f..fc39203 100644 --- a/segmentation/clitkExtractLymphStation_8.txx +++ b/segmentation/clitkExtractLymphStation_8.txx @@ -2,6 +2,23 @@ #include #include +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_SetDefaultValues() +{ + SetDistanceMaxToAnteriorPartOfTheSpine(10); + MaskImagePointType p; + p[0] = 15; p[1] = 2; p[2] = 1; + SetEsophagusDiltationForAnt(p); + p[0] = 5; p[1] = 10; p[2] = 1; + SetEsophagusDiltationForRight(p); + SetFuzzyThreshold("8", "Esophagus", 0.5); + SetInjectedThresholdForS8(150); +} +//-------------------------------------------------------------------- + //-------------------------------------------------------------------- template void @@ -29,39 +46,25 @@ ExtractStation_8_SI_Limits() // Get Carina Z position MaskImagePointer Carina = GetAFDB()->template GetImage("Carina"); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Carina"); + std::vector centroids; clitk::ComputeCentroids(Carina, GetBackgroundValue(), centroids); - DD(centroids.size()); - DD(centroids[0]); // BG - DD(centroids[1]); m_CarinaZ = centroids[1][2]; - DD(m_CarinaZ); + // DD(m_CarinaZ); // add one slice to include carina ? m_CarinaZ += m_Mediastinum->GetSpacing()[2]; - DD(m_CarinaZ); - DD(Carina->GetReferenceCount()); // We dont need Carina structure from now Carina->Delete(); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after delete Carina"); + GetAFDB()->SetDouble("CarinaZ", m_CarinaZ); - // Get left lower lobe bronchus (L), take the upper border - // Get right bronchus intermedius (RML), take the lower border - - /* - double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) + - m_Mediastinum->GetSpacing()[2]; // add one slice to include carina - // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2); - */ - // Found most inferior part of the lung MaskImagePointer Lungs = GetAFDB()->template GetImage("Lungs"); // It should be already croped, so I took the origin and add 10mm above m_DiaphragmInferiorLimit = Lungs->GetOrigin()[2]+10; - DD(m_DiaphragmInferiorLimit); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Lung"); - Lungs->Delete(); // we don't need it, release memory - clitk::PrintMemory(GetVerboseMemoryFlag(), "after release Lung"); + // Lungs->Delete(); // we don't need it, release memory -> it we want to release, also free in AFDB + clitk::PrintMemory(GetVerboseMemoryFlag(), "after reading lungs"); + GetAFDB()->template ReleaseImage("Lungs"); + clitk::PrintMemory(GetVerboseMemoryFlag(), "after release lungs"); /* Crop support : Superior limit = carina @@ -84,11 +87,16 @@ ExtractStation_8_SI_Limits() boolFilter->SetInput1(m_Working_Support); boolFilter->SetInput2(VertebralBody); boolFilter->SetOperationType(BoolFilterType::AndNot); - boolFilter->Update(); - boolFilter->SetInput1(boolFilter->GetOutput()); - boolFilter->SetInput2(Aorta); - boolFilter->SetOperationType(BoolFilterType::AndNot); - boolFilter->Update(); + boolFilter->Update(); + /* + + DO NOT REMOVE AORTA YET (LATER) + + boolFilter->SetInput1(boolFilter->GetOutput()); + boolFilter->SetInput2(Aorta); + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + */ m_Working_Support = boolFilter->GetOutput(); // Done. @@ -102,7 +110,7 @@ ExtractStation_8_SI_Limits() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_8_AP_Limits() +ExtractStation_8_Post_Limits() { /* Station 8: paraeosphageal nodes @@ -155,13 +163,13 @@ ExtractStation_8_AP_Limits() // Convert 2D points in slice into 3D points std::vector vertebralAntPositions; - clitk::PointsUtils::Convert2DTo3DList(vertebralAntPositionBySlice, - VertebralBody, - vertebralAntPositions); + clitk::PointsUtils::Convert2DMapTo3DList(vertebralAntPositionBySlice, + VertebralBody, + vertebralAntPositions); // DEBUG : write list of points clitk::WriteListOfLandmarks(vertebralAntPositions, - "vertebralMostAntPositions.txt"); + "S8-vertebralMostAntPositions-points.txt"); // Cut support posteriorly 1cm the most anterior point of the // VertebralBody. Do nothing below and above the VertebralBody. To @@ -201,7 +209,16 @@ ExtractStation_8_AP_Limits() StopCurrentStep(m_Working_Support); m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_Ant_Sup_Limits() +{ //-------------------------------------------------------------------- StartNewStep("[Station8] Ant limits with S7 above Carina"); /* @@ -222,24 +239,17 @@ ExtractStation_8_AP_Limits() // Ant limit from carina (start) to end of S7 = originRMLB // Get Trachea MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Trachea"); - // Crop above Carina - //double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) + - // Trachea->GetSpacing()[2]; // add one slice to include carina; - MaskImagePointer m_Working_Trachea = - clitk::CropImageAbove(Trachea, 2, m_CarinaZ, true, // AutoCrop + clitk::CropImageRemoveGreaterThan(Trachea, 2, m_CarinaZ, true, // AutoCrop GetBackgroundValue()); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after crop"); // Seprate into two main bronchi - MaskImagePointer LeftBronchus; MaskImagePointer RightBronchus; + MaskImagePointer LeftBronchus; // Labelize and consider the two first (main) labels m_Working_Trachea = Labelize(m_Working_Trachea, 0, true, 1); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after labelize"); // Carina position must at the first slice that separate the two // main bronchus (not superiorly). We thus first check that the @@ -268,67 +278,56 @@ ExtractStation_8_AP_Limits() ImagePointType C1 = c[1]; ImagePointType C2 = c[2]; - ImagePixelType leftLabel; - ImagePixelType rightLabel; - if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; } - else { leftLabel = 2; rightLabel = 1; } + ImagePixelType rightLabel; + ImagePixelType leftLabel; + if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; } + else { rightLabel = 2; leftLabel = 1; } // Select LeftLabel (set one label to Backgroundvalue) - LeftBronchus = + RightBronchus = + clitk::Binarize(m_Working_Trachea, rightLabel, rightLabel, + GetBackgroundValue(), GetForegroundValue()); + /* SetBackground(m_Working_Trachea, m_Working_Trachea, - rightLabel, GetBackgroundValue(), false); - RightBronchus = + leftLabel, GetBackgroundValue(), false); + */ + LeftBronchus = clitk::Binarize(m_Working_Trachea, leftLabel, leftLabel, + GetBackgroundValue(), GetForegroundValue()); + /* SetBackground(m_Working_Trachea, m_Working_Trachea, - leftLabel, GetBackgroundValue(), false); + rightLabel, GetBackgroundValue(), false); + */ + // Crop images - LeftBronchus = clitk::AutoCrop(LeftBronchus, GetBackgroundValue()); RightBronchus = clitk::AutoCrop(RightBronchus, GetBackgroundValue()); + LeftBronchus = clitk::AutoCrop(LeftBronchus, GetBackgroundValue()); // Insert int AFDB if need after - GetAFDB()->template SetImage ("LeftBronchus", "seg/leftBronchus.mhd", - LeftBronchus, true); GetAFDB()->template SetImage ("RightBronchus", "seg/rightBronchus.mhd", RightBronchus, true); - + GetAFDB()->template SetImage ("LeftBronchus", "seg/leftBronchus.mhd", + LeftBronchus, true); // Now crop below OriginOfRightMiddleLobeBronchusZ // It is not done before to keep entire bronchi. - /* - double OriginOfRightMiddleLobeBronchusZ = - GetAFDB()->GetPoint3D("OriginOfRightMiddleLobeBronchus", 2)+ - LeftBronchus->GetSpacing()[2]; - // ^--> Add one slice because the origin is the first slice without S7 - DD(OriginOfRightMiddleLobeBronchusZ); - DD(OriginOfRightMiddleLobeBronchusZ-LeftBronchus->GetSpacing()[2]); - */ - MaskImagePointer OriginOfRightMiddleLobeBronchus = GetAFDB()->template GetImage("OriginOfRightMiddleLobeBronchus"); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after read OriginOfRightMiddleLobeBronchus"); std::vector centroids; clitk::ComputeCentroids(OriginOfRightMiddleLobeBronchus, GetBackgroundValue(), centroids); - DD(centroids.size()); - DD(centroids[0]); // BG - DD(centroids[1]); m_OriginOfRightMiddleLobeBronchusZ = centroids[1][2]; - DD(m_OriginOfRightMiddleLobeBronchusZ); // add one slice to include carina ? - m_OriginOfRightMiddleLobeBronchusZ += LeftBronchus->GetSpacing()[2]; - DD(m_OriginOfRightMiddleLobeBronchusZ); - DD(OriginOfRightMiddleLobeBronchus->GetReferenceCount()); + m_OriginOfRightMiddleLobeBronchusZ += RightBronchus->GetSpacing()[2]; // We dont need Carina structure from now OriginOfRightMiddleLobeBronchus->Delete(); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after delete OriginOfRightMiddleLobeBronchus"); - - LeftBronchus = - clitk::CropImageBelow(LeftBronchus, 2, + RightBronchus = + clitk::CropImageRemoveLowerThan(RightBronchus, 2, m_OriginOfRightMiddleLobeBronchusZ, true, // AutoCrop GetBackgroundValue()); - RightBronchus = - clitk::CropImageBelow(RightBronchus, 2, + LeftBronchus = + clitk::CropImageRemoveLowerThan(LeftBronchus, 2, m_OriginOfRightMiddleLobeBronchusZ, true, // AutoCrop GetBackgroundValue()); @@ -336,87 +335,42 @@ ExtractStation_8_AP_Limits() // Search for points that are the most left/post/ant and most // right/post/ant of the left and right bronchus // 15 = not 15 mm more distance than the middle point. - FindExtremaPointsInBronchus(LeftBronchus, 0, 10, m_RightMostInLeftBronchus, - m_AntMostInLeftBronchus, m_PostMostInLeftBronchus); - FindExtremaPointsInBronchus(RightBronchus, 1, 10, m_LeftMostInRightBronchus, + FindExtremaPointsInBronchus(RightBronchus, 0, 10, m_LeftMostInRightBronchus, m_AntMostInRightBronchus, m_PostMostInRightBronchus); + FindExtremaPointsInBronchus(LeftBronchus, 1, 10, m_RightMostInLeftBronchus, + m_AntMostInLeftBronchus, m_PostMostInLeftBronchus); + // DEBUG : write the list of points ListOfPointsType v; + v.reserve(m_LeftMostInRightBronchus.size()+m_AntMostInRightBronchus.size()+ + m_PostMostInRightBronchus.size()); + v.insert(v.end(), m_LeftMostInRightBronchus.begin(), m_LeftMostInRightBronchus.end() ); + v.insert(v.end(), m_AntMostInRightBronchus.begin(), m_AntMostInRightBronchus.end() ); + v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() ); + clitk::WriteListOfLandmarks(v, "S8-RightBronchus-points.txt"); + + v.clear(); v.reserve(m_RightMostInLeftBronchus.size()+m_AntMostInLeftBronchus.size()+ m_PostMostInLeftBronchus.size()); v.insert(v.end(), m_RightMostInLeftBronchus.begin(), m_RightMostInLeftBronchus.end() ); v.insert(v.end(), m_AntMostInLeftBronchus.begin(), m_AntMostInLeftBronchus.end() ); v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() ); - clitk::WriteListOfLandmarks(v, "LeftBronchusPoints.txt"); + clitk::WriteListOfLandmarks(v, "S8-LeftBronchus-points.txt"); v.clear(); - v.reserve(m_LeftMostInRightBronchus.size()+m_AntMostInRightBronchus.size()+ - m_PostMostInRightBronchus.size()); - v.insert(v.end(), m_LeftMostInRightBronchus.begin(), m_LeftMostInRightBronchus.end() ); - v.insert(v.end(), m_AntMostInRightBronchus.begin(), m_AntMostInRightBronchus.end() ); + v.reserve(m_PostMostInLeftBronchus.size()+m_PostMostInRightBronchus.size()); + v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() ); v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() ); - clitk::WriteListOfLandmarks(v, "RightBronchusPoints.txt"); + clitk::WriteListOfLandmarks(v, "S8-RightLeftBronchus-points.txt"); // Now uses these points to limit, slice by slice - // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870 - /* - Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute: - (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax) - This will equal zero if the point C is on the line formed by - points A and B, and will have a different sign depending on the - side. Which side this is depends on the orientation of your (x,y) - coordinates, but you can plug test values for A,B and C into this - formula to determine whether negative values are to the left or to - the right. - => to accelerate, start with formula, when change sign -> stop and fill - */ - typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; - SliceIteratorType siter = SliceIteratorType(m_Working_Support, - m_Working_Support->GetLargestPossibleRegion()); - siter.SetFirstDirection(0); - siter.SetSecondDirection(1); - siter.GoToBegin(); - int i=0; - MaskImageType::PointType A; - MaskImageType::PointType B; - MaskImageType::PointType C; - while (!siter.IsAtEnd()) { - // Check that the current slice correspond to the current point - m_Working_Support->TransformIndexToPhysicalPoint(siter.GetIndex(), C); - if (C[2] != m_PostMostInLeftBronchus[i][2]) { - // m_Working_Support start from GOjunction while list of point - // start at OriginOfRightMiddleLobeBronchusZ, so we must skip some slices. - } - else { - // Define A,B,C points - A = m_PostMostInLeftBronchus[i]; - B = m_PostMostInRightBronchus[i]; - C = A; - C[1] += 10; // I know I must keep this point - double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); - bool isPositive = s<0; - while (!siter.IsAtEndOfSlice()) { - while (!siter.IsAtEndOfLine()) { - // Very slow, I know ... but image should be very small - m_Working_Support->TransformIndexToPhysicalPoint(siter.GetIndex(), C); - double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); - if (s == 0) siter.Set(2); - if (isPositive) { - if (s > 0) siter.Set(GetBackgroundValue()); - } - else { - if (s < 0) siter.Set(GetBackgroundValue()); - } - ++siter; - } - siter.NextLine(); - } - ++i; - } - siter.NextSlice(); - } + // line is mainly horizontal, so mainDirection=1 + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + m_PostMostInRightBronchus, + m_PostMostInLeftBronchus, + GetBackgroundValue(), 1, 10); // Keep main 3D CCL : m_Working_Support = Labelize(m_Working_Support, 0, false, 10); @@ -431,8 +385,17 @@ ExtractStation_8_AP_Limits() StopCurrentStep(m_Working_Support); // m_ListOfStations["8"] = m_Working_Support; +} + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_Ant_Inf_Limits() +{ + //-------------------------------------------------------------------- - StartNewStep("[Station8] Ant limits arround esophagus below Carina"); + StartNewStep("[Station8] Ant part: not post to Esophagus"); /* Consider Esophagus, dilate it and remove ant part. It remains part on L & R, than can be partly removed by cutting what remains at @@ -440,259 +403,739 @@ ExtractStation_8_AP_Limits() */ // Get Esophagus - DD("Esophagus"); - MaskImagePointer Esophagus = GetAFDB()->template GetImage("Esophagus"); + m_Esophagus = GetAFDB()->template GetImage("Esophagus"); clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Esophagus"); + // In images from the original article, Atlas – UM, the oesophagus + //was included in nodal stations 3p and 8. Having said that, in the + //description for station 8, it indicates that “the delineation of + //station 8 is limited to the soft tissues surrounding the + //oesophagus”. In the recent article, The IASLC Lung Cancer Staging + //Project (J Thorac Oncol 4:5, 568-77), the images drawn by + //Dr. Aletta Frasier exclude this structure. From an oncological + //prospective, the oesophagus should be excluded from these nodal + //stations. + + /* NOT YET !! DO IT LATER + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(m_Esophagus); + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + + */ + // Crop Esophagus : keep only below the OriginOfRightMiddleLobeBronchusZ - Esophagus = - clitk::CropImageAbove(Esophagus, 2, + m_Esophagus = + clitk::CropImageRemoveGreaterThan(m_Esophagus, 2, m_OriginOfRightMiddleLobeBronchusZ, true, // AutoCrop GetBackgroundValue()); // Dilate to keep only not Anterior positions MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt(); - Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus); - Esophagus = clitk::Dilate(Esophagus, - radiusInMM, - GetBackgroundValue(), - GetForegroundValue(), true); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after dilate Esophagus"); - writeImage(Esophagus, "enlarged-eso.mhd"); - // Remove Anterior part according to this dilatated esophagus + // m_Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(m_Esophagus); + + m_Esophagus = clitk::Dilate(m_Esophagus, + radiusInMM, + GetBackgroundValue(), + GetForegroundValue(), true); + + // Remove Anterior part according to this dilatated esophagus. Note: + // because we crop Esophagus with ORML, the support will also be + // croped in the same way. Here it is a desired feature. If we dont + // want, use SetIgnoreEmptySliceObject(true) + + // In the new IASCL definition, it is not clear if sup limits is + // around carina or On the right, it is “the lower border of the + // bronchus intermedius”, indicated on the image set as a point + // (“lower border of the bronchus intermedius”) + typedef clitk::SliceBySliceRelativePositionFilter RelPosFilterType; typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); relPosFilter->VerboseStepFlagOff(); + relPosFilter->WriteStepFlagOff(); relPosFilter->SetInput(m_Working_Support); - relPosFilter->SetInputObject(Esophagus); - relPosFilter->AddOrientationTypeString("P"); - relPosFilter->InverseOrientationFlagOff(); + relPosFilter->SetInputObject(m_Esophagus); + relPosFilter->AddOrientationTypeString("PostTo"); + // relPosFilter->InverseOrientationFlagOff(); relPosFilter->SetDirection(2); // Z axis relPosFilter->UniqueConnectedComponentBySliceOff(); relPosFilter->SetIntermediateSpacing(3); - relPosFilter->ResampleBeforeRelativePositionFilterOn(); - relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8()); - relPosFilter->RemoveObjectFlagOff(); + relPosFilter->IntermediateSpacingFlagOn(); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("8", "Esophagus")); + relPosFilter->RemoveObjectFlagOff(); // Do not exclude here because it is dilated + relPosFilter->CombineWithOrFlagOff(); // NO ! + relPosFilter->IgnoreEmptySliceObjectFlagOn(); relPosFilter->Update(); m_Working_Support = relPosFilter->GetOutput(); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after SbS Rel pos Post"); // AutoCrop (OR SbS ?) m_Working_Support = clitk::AutoCrop(m_Working_Support, GetBackgroundValue()); - writeImage(m_Working_Support, "step1.4.1.mhd"); + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_Ant_Injected_Limits() +{ + + //-------------------------------------------------------------------- + StartNewStep("[Station8] Ant part (remove high density, injected part)"); + + // Consider initial image, crop to current support + ImagePointer working_input = + clitk::ResizeImageLike(m_Input, + m_Working_Support, + (short)GetBackgroundValue()); + + // Threshold + typedef itk::BinaryThresholdImageFilter BinarizeFilterType; + typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); + binarizeFilter->SetInput(working_input); + binarizeFilter->SetLowerThreshold(GetInjectedThresholdForS8()); + binarizeFilter->SetInsideValue(GetForegroundValue()); + binarizeFilter->SetOutsideValue(GetBackgroundValue()); + binarizeFilter->Update(); + MaskImagePointer injected = binarizeFilter->GetOutput(); + + // Combine with current support + clitk::AndNot(m_Working_Support, injected, GetBackgroundValue()); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_LR_1_Limits() +{ + //-------------------------------------------------------------------- + StartNewStep("[Station8] Left and Right (from Carina to PulmonaryTrunk): Right to LeftPulmonaryArtery"); + + /* + We remove LeftPulmonaryArtery structure and what is at Left to + this structure. + */ + MaskImagePointer LeftPulmonaryArtery = GetAFDB()->template GetImage("LeftPulmonaryArtery"); + + // Relative Position : not at Left + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->VerboseStepFlagOff(); + relPosFilter->WriteStepFlagOff(); + relPosFilter->SetBackgroundValue(GetBackgroundValue()); + relPosFilter->SetInput(m_Working_Support); + relPosFilter->SetInputObject(LeftPulmonaryArtery); + relPosFilter->RemoveObjectFlagOn(); // remove the object too + relPosFilter->AddOrientationTypeString("L"); + relPosFilter->InverseOrientationFlagOn(); // Not at Left + relPosFilter->SetIntermediateSpacing(3); + relPosFilter->IntermediateSpacingFlagOn(); + relPosFilter->SetFuzzyThreshold(0.7); + relPosFilter->AutoCropFlagOn(); + relPosFilter->Update(); + m_Working_Support = relPosFilter->GetOutput(); + + // Release LeftPulmonaryArtery + GetAFDB()->template ReleaseImage("LeftPulmonaryArtery"); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + - clitk::PrintMemory(GetVerboseMemoryFlag(), "after autocrop"); +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_LR_2_Limits() +{ + //-------------------------------------------------------------------- + StartNewStep("[Station8] Left and Right (from PulmTrunk to OriginMiddleLobeBronchus) Right to line from Aorta to PulmonaryTrunk"); + + /* + We consider a line from Left part of Aorta to left part of + PulmonaryTrunk and remove what is at Left. + */ + + // Load the structures + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + MaskImagePointer PulmonaryTrunk = GetAFDB()->template GetImage("PulmonaryTrunk"); + + // Resize like the PT and define the slices + MaskImagePointType min, max; + clitk::GetMinMaxPointPosition(PulmonaryTrunk, min, max); + Aorta = clitk::CropImageAlongOneAxis(Aorta, 2, min[2], max[2], false, GetBackgroundValue()); + std::vector slices_aorta; + clitk::ExtractSlices(Aorta, 2, slices_aorta); + std::vector slices_PT; + clitk::ExtractSlices(PulmonaryTrunk, 2, slices_PT); + + // Find the points at left + std::vector p_A; + std::vector p_B; + MaskImagePointType pA; + MaskImagePointType pB; + for(uint i=0; i(slices_aorta[i], GetBackgroundValue(), 0, false, ps); + clitk::PointsUtils::Convert2DTo3D(ps, Aorta, i, pA); + + if (found) { + // In PT : generally 2 CCL, we keep the most at left + found = + clitk::FindExtremaPointInAGivenDirection(slices_PT[i], GetBackgroundValue(), 0, false, ps); + clitk::PointsUtils::Convert2DTo3D(ps, PulmonaryTrunk, i, pB); + } + + if (found) { + p_A.push_back(pA); + p_B.push_back(pB); + } + } + clitk::WriteListOfLandmarks(p_A, "S8-Aorta-Left-points.txt"); + clitk::WriteListOfLandmarks(p_B, "S8-PT-Left-points.txt"); + + // Remove part at Left + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + p_A, p_B, + GetBackgroundValue(), + 0, -10); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_Single_CCL_Limits() +{ + //-------------------------------------------------------------------- + StartNewStep("[Station8 or 3P] Slice by slice, keep a single CCL (the closest to VertebralBody)"); + + // Consider slices + std::vector slices; + std::vector slices_vb; + clitk::ExtractSlices(m_Working_Support, 2, slices); + MaskImagePointer VertebralBody = + GetAFDB()->template GetImage ("VertebralBody"); + VertebralBody = clitk::ResizeImageLike(VertebralBody, m_Working_Support, GetBackgroundValue()); + clitk::ExtractSlices(VertebralBody, 2, slices_vb); + + for(uint i=0; i(slices[i], 0, true, 100); + // Compute centroids coordinate + std::vector centroids; + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), centroids); + clitk::ComputeCentroids(slices_vb[i], GetBackgroundValue(), c); + if ((c.size() > 1) && (centroids.size() > 1)) { + // keep the one which is the closer to vertebralbody centroid + double distance=1000000; + int index=0; + for(uint j=1; j(slices[i], GetBackgroundValue(), + GetForegroundValue(), index, index, true); + } + } + m_Working_Support = clitk::JoinSlices(slices, m_Working_Support, 2); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_LR_Limits_old2() +{ + + //-------------------------------------------------------------------- + StartNewStep("[Station8] Left and Right limits arround esophagus (below Carina)"); // Estract slices for current support for slice by slice processing std::vector slices; clitk::ExtractSlices(m_Working_Support, 2, slices); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after support slices"); + // Dilate the Esophagus to consider a margins around + MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt(); + m_Esophagus = clitk::Dilate(m_Esophagus, + radiusInMM, + GetBackgroundValue(), + GetForegroundValue(), true); + + // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select + // 'better' extrema points (not too post). + MaskImagePointer Lungs = GetAFDB()->template GetImage("Lungs"); + clitk::AndNot(m_Esophagus, Lungs, GetBackgroundValue()); + GetAFDB()->template ReleaseImage("Lungs"); + // Estract slices of Esophagus (resize like support before to have the same set of slices) - MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike(Esophagus, m_Working_Support, GetBackgroundValue()); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after Eso resize"); + MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike(m_Esophagus, m_Working_Support, GetBackgroundValue()); + std::vector eso_slices; clitk::ExtractSlices(EsophagusForSlice, 2, eso_slices); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after Eso slices"); // Estract slices of Vertebral (resize like support before to have the same set of slices) - VertebralBody = GetAFDB()->template GetImage("VertebralBody"); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after Read VertebralBody"); + MaskImagePointer VertebralBody = GetAFDB()->template GetImage("VertebralBody"); VertebralBody = clitk::ResizeImageLike(VertebralBody, m_Working_Support, GetBackgroundValue()); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after VertebralBody Resize"); std::vector vert_slices; clitk::ExtractSlices(VertebralBody, 2, vert_slices); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after VertebralBody slices"); + + // Estract slices of Aorta (resize like support before to have the same set of slices) + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + Aorta = clitk::ResizeImageLike(Aorta, m_Working_Support, GetBackgroundValue()); + std::vector aorta_slices; + clitk::ExtractSlices(Aorta, 2, aorta_slices); // Extract slices of Mediastinum (resize like support before to have the same set of slices) m_Mediastinum = GetAFDB()->template GetImage("Mediastinum"); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Mediastinum"); m_Mediastinum = clitk::ResizeImageLike(m_Mediastinum, m_Working_Support, GetBackgroundValue()); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after resize Mediastinum"); std::vector mediast_slices; clitk::ExtractSlices(m_Mediastinum, 2, mediast_slices); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after Mediastinum slices"); - - writeImage(EsophagusForSlice, "slices_eso.mhd"); - writeImage(VertebralBody, "slices_vert.mhd"); - writeImage(m_Mediastinum, "slices_medias.mhd"); - writeImage(m_Working_Support, "slices_support.mhd"); - - // Find common slices between Eso and m_Working_Support - // int s=0; - // MaskImageIndexType z_Eso = Esophagus->GetLargestPossibleRegion().GetIndex(); - // MaskImagePointType p_Eso; - // Esophagus->TransformIndexToPhysicalPoint(z_Eso, p_Eso); - // MaskImageIndexType z_Support; - // z_Support = m_Working_Support->GetLargestPossibleRegion().GetIndex(); - // MaskImagePointType p_Support; - // m_Working_Support->TransformIndexToPhysicalPoint(z_Support, p_Support); - // while (p_Eso[2] < p_Support[2]) { - // z_Eso[2] ++; - // Esophagus->TransformIndexToPhysicalPoint(z_Eso, p_Eso); - // } - // s = z_Eso[2] - Esophagus->GetLargestPossibleRegion().GetIndex()[2]; - // DD(s); - - // Find common slices between m_Working_Support and Mediastinum - // int sm=0; - // MaskImageIndexType z_Mediast = m_Mediastinum->GetLargestPossibleRegion().GetIndex(); - // MaskImagePointType p_Mediast; - // m_Mediastinum->TransformIndexToPhysicalPoint(z_Mediast, p_Mediast); - // z_Support = m_Working_Support->GetLargestPossibleRegion().GetIndex(); - // m_Working_Support->TransformIndexToPhysicalPoint(z_Support, p_Support); - // while (p_Mediast[2] < p_Support[2]) { - // z_Mediast[2] ++; - // m_Mediastinum->TransformIndexToPhysicalPoint(z_Mediast, p_Mediast); - // } - // sm = z_Mediast[2] - m_Mediastinum->GetLargestPossibleRegion().GetIndex()[2]; - // DD(sm); - - DD(EsophagusForSlice->GetLargestPossibleRegion().GetSize()[2]); - DD(m_Mediastinum->GetLargestPossibleRegion().GetSize()[2]); - DD(slices.size()); - DD(vert_slices.size()); - DD(eso_slices.size()); - DD(mediast_slices.size()); - // uint max = std::min(slices.size(), std::min(eso_slices.size()-s, mediast_slices.size()-sm)); - // DD(max); + // List of points + std::vector p_RightMostAnt; + std::vector p_RightMostPost; + std::vector p_LeftMostAnt; + std::vector p_LeftMostPost; + std::vector p_AllPoints; + std::vector p_LeftAorta; + std::vector p_LeftEso; - // DEBUG POINTS - std::ofstream osp; - openFileForWriting(osp, "point_s8.txt"); - osp << "LANDMARKS1" << std::endl; + /* + In the following, we search for the LeftRight limits. We search + for the most Right points in Esophagus and in VertebralBody and + consider a line between those to most right points. All points in + the support which are most right to this line are discarded. Same + for the left part. The underlying assumption is that the support + is concave between Eso/VertebralBody. Esophagus is a bit + dilatated. On VertebralBody we go right (or left) until we reach + the lung (but no more 20 mm). + */ // Loop slices MaskImagePointType p; - int pi=0; + MaskImagePointType pp; for(uint i=0; i(eso_slices[i], GetBackgroundValue(), - 1, false, p_maxRight_Eso); - // Debug point - clitk::PointsUtils::Convert2DTo3D(p_maxRight_Eso, EsophagusForSlice, i, p); - osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl; - ++pi; - - // Find right limit of Vertebra - typename MaskSliceType::PointType p_maxRight_Vertebra; - clitk::FindExtremaPointInAGivenDirection(vert_slices[i], GetBackgroundValue(), - 1, false, p_maxRight_Vertebra); - // Debug point - clitk::PointsUtils::Convert2DTo3D(p_maxRight_Vertebra, VertebralBody, i, p); - osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl; - ++pi; + // Find limit of Vertebral -> only at most Post part of current + // slice support. First found most ant point in VertebralBody + typedef MaskSliceType SliceType; + typename SliceType::PointType p_slice_ant; + bool found = clitk::FindExtremaPointInAGivenDirection(vert_slices[i], GetBackgroundValue(), 1, true, p_slice_ant); + if (!found) { + // 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; + int j=i++; + bool found = false; + while (!found) { + found = clitk::FindExtremaPointInAGivenDirection(vert_slices[j], GetBackgroundValue(), 1, true, p_slice_ant); + //clitkExceptionMacro("No foreground pixels in this VertebralBody slices ??"); + j++; + } + } + p_slice_ant[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset + + // The, find most Right and Left points on that AP position + typename SliceType::IndexType indexR; + typename SliceType::IndexType indexL; + vert_slices[i]->TransformPhysicalPointToIndex(p_slice_ant, indexR); + indexL = indexR; + // Check that is inside the mask + indexR[1] = std::min(indexR[1], (long)vert_slices[i]->GetLargestPossibleRegion().GetSize()[1]-1); + indexL[1] = indexR[1]; + while (vert_slices[i]->GetPixel(indexR) != GetBackgroundValue()) { + indexR[0] --; // Go to the right + } + while (vert_slices[i]->GetPixel(indexL) != GetBackgroundValue()) { + indexL[0] ++; // Go to the left + } + vert_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra); + clitk::PointsUtils::Convert2DTo3D(sp_maxRight_Vertebra, VertebralBody, i, p); + p_AllPoints.push_back(p); + vert_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra); + clitk::PointsUtils::Convert2DTo3D(sp_maxLeft_Vertebra, VertebralBody, i, p); + p_AllPoints.push_back(p); - // Find last point out of the mediastinum on this line : - typename MaskSliceType::IndexType index; - mediast_slices[i]->TransformPhysicalPointToIndex(p_maxRight_Vertebra, index); + // Find last point out of the mediastinum on this line, Right : + mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxRight_Vertebra, indexR); double distance = 0.0; - while (mediast_slices[i]->GetPixel(index) != GetBackgroundValue()) { - index[0] --; + while (mediast_slices[i]->GetPixel(indexR) != GetBackgroundValue()) { + indexR[0] --; distance += mediast_slices[i]->GetSpacing()[0]; } - DD(distance); - if (distance < 20) { // Ok in this case, we found limit with lung - mediast_slices[i]->TransformIndexToPhysicalPoint(index, p_maxRight_Vertebra); - clitk::PointsUtils::Convert2DTo3D(p_maxRight_Vertebra, m_Mediastinum, i, p); - osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl; - ++pi; + if (distance < 30) { // Ok in this case, we found limit with lung + mediast_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra); + clitk::PointsUtils::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p); } - - // Choose the most extrema one - if (p_maxRight_Vertebra[0] < p_maxRight_Eso[0]) { - p_maxRight = p_maxRight_Vertebra; + else { // in that case, we are probably below the diaphragm, so we + // add aribtrarly few mm + sp_maxRight_Vertebra[0] -= 2; // Leave 2 mm around the VertebralBody + clitk::PointsUtils::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p); + } + p_RightMostPost.push_back(p); + p_AllPoints.push_back(p); + + // Find last point out of the mediastinum on this line, Left : + mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxLeft_Vertebra, indexL); + distance = 0.0; + while (mediast_slices[i]->GetPixel(indexL) != GetBackgroundValue()) { + indexL[0] ++; + distance += mediast_slices[i]->GetSpacing()[0]; + } + if (distance < 30) { // Ok in this case, we found limit with lung + mediast_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra); + clitk::PointsUtils::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p); + } + else { // in that case, we are probably below the diaphragm, so we + // add aribtrarly few mm + sp_maxLeft_Vertebra[0] += 2; // Leave 2 mm around the VertebralBody + clitk::PointsUtils::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p); + } + p_LeftMostPost.push_back(p); + p_AllPoints.push_back(p); + + // Find Eso slice centroid and do not consider what is post to + // this centroid. + std::vector c; + clitk::ComputeCentroids(eso_slices[i], GetBackgroundValue(), c); + if (c.size() >1) { + eso_slices[i] = + clitk::CropImageRemoveGreaterThan(eso_slices[i], 1, c[1][1], false, GetBackgroundValue()); + eso_slices[i] = + clitk::ResizeImageLike(eso_slices[i], aorta_slices[i], GetBackgroundValue()); + // writeImage(eso_slices[i], "eso-slice-"+toString(i)+".mhd"); + } + + // Find right limit of Esophagus and Aorta + bool f = + clitk::FindExtremaPointInAGivenDirection(eso_slices[i], GetBackgroundValue(), + 0, true, sp_maxRight_Eso); + f = f && + clitk::FindExtremaPointInAGivenDirection(aorta_slices[i], GetBackgroundValue(), + 0, true, sp_maxRight_Aorta); + clitk::PointsUtils::Convert2DTo3D(sp_maxRight_Eso, EsophagusForSlice, i, p); + clitk::PointsUtils::Convert2DTo3D(sp_maxRight_Aorta, Aorta, i, pp); + pp[0] -= 2; // Add a margin of 2 mm to include the Aorta 'wall' + p_AllPoints.push_back(p); + if (f) { + p_AllPoints.push_back(pp); + MaskImagePointType A = p_RightMostPost.back(); + MaskImagePointType B = p; + MaskImagePointType C = pp; + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + if (s>0) p_RightMostAnt.push_back(p); + else p_RightMostAnt.push_back(pp); + } + else { // No more Esophagus in this slice : do nothing + // p_RightMostAnt.push_back(p); + p_RightMostPost.pop_back(); } - else p_maxRight = p_maxRight_Eso; // -------------------------------------------------------------------------- + // Find the limit on the Left: most left point between Eso and + // Eso. (Left is left on screen, coordinate increase) + + // Find left limit of Esophagus + clitk::FindExtremaPointInAGivenDirection(eso_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Eso); + f = clitk::FindExtremaPointInAGivenDirection(aorta_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Aorta); + clitk::PointsUtils::Convert2DTo3D(sp_maxLeft_Eso, EsophagusForSlice, i, p); + clitk::PointsUtils::Convert2DTo3D(sp_maxLeft_Aorta, Aorta, i, pp); + p_AllPoints.push_back(p); + pp[0] += 2; // Add a margin of 2 mm to include the 'wall' + if (f) { // not below Aorta + p_AllPoints.push_back(pp); + MaskImagePointType A = p_LeftMostPost.back(); + MaskImagePointType B = p; + MaskImagePointType C = pp; + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + if (s<0) { + p_LeftMostAnt.push_back(p); // Insert point most at Left, Eso + } + else { + // in this case -> two lines ! + p_LeftMostAnt.push_back(pp); // Insert point most at Left, Aorta (Vert to Aorta) + // but also consider Aorta to Eso + p_LeftAorta.push_back(pp); + p_LeftEso.push_back(p); + } + } + else { // No more Esophagus in this slice : do nothing + p_LeftMostPost.pop_back(); + //p_LeftMostAnt.push_back(p); + } + } // End of slice loop + + clitk::WriteListOfLandmarks(p_AllPoints, "S8-LR-Eso-Vert.txt"); + clitk::WriteListOfLandmarks(p_LeftEso, "S8-Left-Eso.txt"); + clitk::WriteListOfLandmarks(p_LeftAorta, "S8-Left-Aorta.txt"); + // Now uses these points to limit, slice by slice + // Line is mainly vertical, so mainDirection=0 + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + p_RightMostAnt, p_RightMostPost, + GetBackgroundValue(), 0, 10); + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + p_LeftMostAnt, p_LeftMostPost, + GetBackgroundValue(), 0, -10); + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + p_LeftEso,p_LeftAorta, + GetBackgroundValue(), 0, -10); + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; + return; +} +//-------------------------------------------------------------------- - /* - // Get most post of dilated Esophagus - typename MaskSliceType::PointType p_post; - bool f = clitk::FindExtremaPointInAGivenDirection(eso_slices[i], GetBackgroundValue(), 1, false, p_post); - clitk::PointsUtils::Convert2DTo3D(p_post, EsophagusForSlice, i, p); - osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl; - ++pi; - // DD(f); - // DD(p); - - // Get most left of the vertebral body - typename MaskSliceType::PointType s_left; - f = f && clitk::FindExtremaPointInAGivenDirection(vert_slices[i], GetBackgroundValue(), 0, false, s_left); - clitk::PointsUtils::Convert2DTo3D(s_left, VertebralBody, i, p); - osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl; - ++pi; - // DD(f); - // DD(p); - - // Find last point out of the mediastinum on this line : - typename MaskSliceType::IndexType index_left; - mediast_slices[i]->TransformPhysicalPointToIndex(s_left, index_left); - index_left[0] ++; // on more left to be inside the support - while (mediast_slices[i]->GetPixel(index_left) != GetBackgroundValue()) { - index_left[0] ++; - } - mediast_slices[i]->TransformIndexToPhysicalPoint(index_left, s_left); - clitk::PointsUtils::Convert2DTo3D(s_left, m_Mediastinum, i, p); - osp << pi << " " << p[0] << " " << p[1] << " " << p[2] << " 0 0 " << std::endl; - ++pi; - // DD(f); - // DD(p); - */ +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_LR_Limits() +{ - // Loop to suppress - // if (f) { - typename MaskSliceType::PointType p; - typedef itk::ImageRegionIteratorWithIndex IteratorType; - IteratorType iter(slices[i], slices[i]->GetLargestPossibleRegion()); - iter.GoToBegin(); - while (!iter.IsAtEnd()) { - if (iter.Get() != GetBackgroundValue()) { - slices[i]->TransformIndexToPhysicalPoint(iter.GetIndex(), p); + //-------------------------------------------------------------------- + StartNewStep("[Station8] Left and Right limits arround esophagus with lines from VertebralBody-Aorta-Esophagus"); - // Remove point too at RIGHT - if (p[0] < p_maxRight[0]) { - iter.Set(GetBackgroundValue()); - } + // Estract slices for current support for slice by slice processing + std::vector slices; + clitk::ExtractSlices(m_Working_Support, 2, slices); + + // Dilate the Esophagus to consider a margins around + MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt(); + m_Esophagus = clitk::Dilate(m_Esophagus, radiusInMM, + GetBackgroundValue(), GetForegroundValue(), true); + // m_Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(m_Esophagus); - /* - // Remove point from foreground if too right or to high - if ((p[1] > p_post[1]) && (p[0] > s_left[0])) { - iter.Set(GetBackgroundValue()); - } - */ + // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select + // 'better' extrema points (not too post). + MaskImagePointer Lungs = GetAFDB()->template GetImage("Lungs"); + clitk::AndNot(m_Esophagus, Lungs, GetBackgroundValue()); + GetAFDB()->template ReleaseImage("Lungs"); + + // Estract slices of Esophagus (resize like support before to have the same set of slices) + MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike(m_Esophagus, m_Working_Support, GetBackgroundValue()); + std::vector eso_slices; + clitk::ExtractSlices(EsophagusForSlice, 2, eso_slices); + + // Estract slices of Vertebral (resize like support before to have the same set of slices) + MaskImagePointer VertebralBody = GetAFDB()->template GetImage("VertebralBody"); + VertebralBody = clitk::ResizeImageLike(VertebralBody, m_Working_Support, GetBackgroundValue()); + std::vector vert_slices; + clitk::ExtractSlices(VertebralBody, 2, vert_slices); + + // Estract slices of Aorta (resize like support before to have the same set of slices) + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + Aorta = clitk::ResizeImageLike(Aorta, m_Working_Support, GetBackgroundValue()); + std::vector aorta_slices; + clitk::ExtractSlices(Aorta, 2, aorta_slices); + + // Extract slices of Mediastinum (resize like support before to have the same set of slices) + m_Mediastinum = GetAFDB()->template GetImage("Mediastinum"); + m_Mediastinum = clitk::ResizeImageLike(m_Mediastinum, m_Working_Support, GetBackgroundValue()); + std::vector mediast_slices; + clitk::ExtractSlices(m_Mediastinum, 2, mediast_slices); + + // List of points + std::vector p_MostLeftVertebralBody; + std::vector p_MostRightVertebralBody; + std::vector p_MostLeftAorta; + std::vector p_MostLeftEsophagus; + + /* + In the following, we search for the LeftRight limits. We search + for the most Right points in Esophagus and in VertebralBody and + consider a line between those to most right points. All points in + the support which are most right to this line are discarded. Same + for the left part. The underlying assumption is that the support + is concave between Eso/VertebralBody. Esophagus is a bit + dilatated. On VertebralBody we go right (or left) until we reach + the lung (but no more 20 mm). + */ + + // Temporary 3D point + MaskImagePointType p; + // Loop slices + for(uint i=0; i(vert_slices[j], GetBackgroundValue(), 1, true, sp_MostAntVertebralBody); + if (!found) { + // 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; + j++; } - ++iter; } - // } // end if f - // s++; - // sm++; - } - osp.close(); + sp_MostAntVertebralBody[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset + + // Crop the vertebralbody below this most post line + vert_slices[j] = + clitk::CropImageRemoveGreaterThan(vert_slices[j], 1, sp_MostAntVertebralBody[1], false, GetBackgroundValue()); + vert_slices[j] = + clitk::ResizeImageLike(vert_slices[j], aorta_slices[i], GetBackgroundValue()); + // writeImage(vert_slices[i], "vert-slice-"+toString(i)+".mhd"); + + // Find first point not in mediastinum + clitk::FindExtremaPointInAGivenDirection(vert_slices[j], GetBackgroundValue(), 0, false, sp_MostLeftVertebralBody); + sp_MostLeftVertebralBody = clitk::FindExtremaPointInAGivenLine(mediast_slices[i], 0, false, sp_MostLeftVertebralBody, GetBackgroundValue(), 30); + sp_MostLeftVertebralBody[0] += 2; // 2mm margin + clitk::FindExtremaPointInAGivenDirection(vert_slices[j], GetBackgroundValue(), 0, true, sp_MostRightVertebralBody); + sp_MostRightVertebralBody = clitk::FindExtremaPointInAGivenLine(mediast_slices[i], 0, true, sp_MostRightVertebralBody, GetBackgroundValue(),30); + sp_MostRightVertebralBody[0] -= 2; // 2 mm margin + + // Convert 2D points into 3D + clitk::PointsUtils::Convert2DTo3D(sp_MostRightVertebralBody, VertebralBody, i, p); + p_MostRightVertebralBody.push_back(p); + clitk::PointsUtils::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p); + p_MostLeftVertebralBody.push_back(p); + + /* ------------------------------------------------------------------------------------- + Find most Left point in Esophagus + */ + found = clitk::FindExtremaPointInAGivenDirection(eso_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftEsophagus); + if (!found) { // No more Esophagus, I remove the previous point + + // if (p_MostLeftEsophagus.size() < 1) { + p_MostLeftVertebralBody.pop_back(); + // } + // else { + // // Consider the previous point + // p = p_MostLeftEsophagus.back(); + // p_MostLeftEsophagus.push_back(p); + // sp_MostLeftEsophagus = sp_temp; // Retrieve previous 2D position + // found = true; + // } + } + else { + clitk::PointsUtils::Convert2DTo3D(sp_MostLeftEsophagus, EsophagusForSlice, i, p); + p_MostLeftEsophagus.push_back(p); + // sp_temp = sp_MostLeftEsophagus; // Store previous 2D position + } + + /* ------------------------------------------------------------------------------------- + Find most Left point in Aorta + */ + if (found) { + clitk::FindExtremaPointInAGivenDirection(aorta_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftAorta); + sp_MostLeftAorta = clitk::FindExtremaPointInAGivenLine(mediast_slices[i], 0, false, sp_MostLeftAorta, GetBackgroundValue(), 10); + typename MaskSliceType::PointType temp=sp_MostLeftEsophagus; + temp[0] -= 10; + if (clitk::IsOnTheSameLineSide(sp_MostLeftAorta, sp_MostLeftVertebralBody, sp_MostLeftEsophagus, temp)) { + // sp_MostLeftAorta is on the same side than temp (at Right) -> ignore Aorta + sp_MostLeftAorta = sp_MostLeftEsophagus; + } + clitk::PointsUtils::Convert2DTo3D(sp_MostLeftAorta, Aorta, i, p); + p_MostLeftAorta.push_back(p); + } + + } // End of slice loop + clitk::WriteListOfLandmarks(p_MostLeftVertebralBody, "S8-MostLeft-VB-points.txt"); + clitk::WriteListOfLandmarks(p_MostRightVertebralBody, "S8-MostRight-VB-points.txt"); + clitk::WriteListOfLandmarks(p_MostLeftAorta, "S8-MostLeft-Aorta-points.txt"); + clitk::WriteListOfLandmarks(p_MostLeftEsophagus, "S8-MostLeft-eso-points.txt"); - // Joint slices - DD("after loop"); - m_Working_Support = clitk::JoinSlices(slices, m_Working_Support, 2); - clitk::PrintMemory(GetVerboseMemoryFlag(), "after JoinSlices"); + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + p_MostLeftVertebralBody, p_MostLeftAorta, + GetBackgroundValue(), 0, -10); + + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + p_MostLeftAorta, p_MostLeftEsophagus, + GetBackgroundValue(), 0, -10); - // DEBUG + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; + return; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_Remove_Structures() +{ + + //-------------------------------------------------------------------- + StartNewStep("[Station8] remove some structures"); + + Remove_Structures("8", "Aorta"); + Remove_Structures("8", "Esophagus"); + + // END + StopCurrentStep(m_Working_Support); m_ListOfStations["8"] = m_Working_Support; return; } @@ -705,12 +1148,25 @@ typename clitk::ExtractLymphStationsFilter::MaskImagePointer clitk::ExtractLymphStationsFilter:: EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus) { - // Check if Esophagus is delineated at least until GOjunction. Now, + // Check if Esophagus is delineated at least until Diaphragm. Now, // because we use AutoCrop, Origin[2] gives this max inferior // position. - // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2); + + DD("BUGGY DONT USE"); + exit(0); if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) { + // crop first slice without mask + MaskImagePointType pt; + clitk::FindExtremaPointInAGivenDirection(Esophagus, GetBackgroundValue(), 2, true, pt); + DD(pt); + Esophagus = + clitk::CropImageRemoveLowerThan(Esophagus, 2, + pt[2], + false, // AutoCrop + GetBackgroundValue()); + writeImage(Esophagus, "crop-eso.mhd"); + std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it." << std::endl; double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit; @@ -734,7 +1190,7 @@ EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus) template void clitk::ExtractLymphStationsFilter:: -ExtractStation_8_LR_Limits() +ExtractStation_8_LR_Limits_old() { /* Station 8: paraeosphageal nodes @@ -774,12 +1230,12 @@ ExtractStation_8_LR_Limits() relPosFilter->WriteStepFlagOff(); relPosFilter->SetInput(m_Working_Support); relPosFilter->SetInputObject(Esophagus); - relPosFilter->AddOrientationTypeString("L"); - relPosFilter->InverseOrientationFlagOn(); // Not Left to + relPosFilter->AddOrientationTypeString("NotLeftTo"); + // relPosFilter->InverseOrientationFlagOn(); // Not Left to relPosFilter->SetDirection(2); // Z axis relPosFilter->UniqueConnectedComponentBySliceOff(); // important relPosFilter->SetIntermediateSpacing(4); - relPosFilter->ResampleBeforeRelativePositionFilterOn(); + relPosFilter->IntermediateSpacingFlagOn(); relPosFilter->SetFuzzyThreshold(0.9); // remove few part only relPosFilter->RemoveObjectFlagOff(); relPosFilter->Update(); @@ -809,7 +1265,7 @@ ExtractStation_8_LR_Limits() relPosFilter->SetDirection(2); // Z axis relPosFilter->UniqueConnectedComponentBySliceOff(); // important relPosFilter->SetIntermediateSpacing(4); - relPosFilter->ResampleBeforeRelativePositionFilterOn(); + relPosFilter->IntermediateSpacingFlagOn(); relPosFilter->CombineWithOrFlagOn(); relPosFilter->SetFuzzyThreshold(0.9); // remove few part only relPosFilter->RemoveObjectFlagOff(); @@ -821,3 +1277,73 @@ ExtractStation_8_LR_Limits() m_ListOfStations["8"] = m_Working_Support; } //-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +FindExtremaPointsInBronchus(MaskImagePointer input, + int direction, + double distance_max_from_center_point, + ListOfPointsType & LR, + ListOfPointsType & Ant, + ListOfPointsType & Post) +{ + + // Other solution ==> with auto bounding box ! (but pb to prevent to + // be too distant from the center point + + // Extract slices + std::vector slices; + clitk::ExtractSlices(input, 2, slices); + + // Loop on slices + bool found; + for(uint i=0; i(slices[i], 0, true, 10); + slices[i] = KeepLabels(slices[i], + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + */ + + // ------- Find rightmost or leftmost point ------- + MaskSliceType::PointType LRMost; + found = + clitk::FindExtremaPointInAGivenDirection(slices[i], + GetBackgroundValue(), + 0, // axis XY + (direction==0?false:true), // right or left according to direction + LRMost); + // ------- Find postmost point ------- + MaskSliceType::PointType postMost; + found = + clitk::FindExtremaPointInAGivenDirection(slices[i], + GetBackgroundValue(), + 1, false, LRMost, + distance_max_from_center_point, + postMost); + // ------- Find antmost point ------- + MaskSliceType::PointType antMost; + found = + clitk::FindExtremaPointInAGivenDirection(slices[i], + GetBackgroundValue(), + 1, true, LRMost, + distance_max_from_center_point, + antMost); + // Only add point if found + if (found) { + // ------- Convert 2D to 3D points -------- + MaskImageType::PointType p; + clitk::PointsUtils::Convert2DTo3D(LRMost, input, i, p); + LR.push_back(p); + clitk::PointsUtils::Convert2DTo3D(antMost, input, i, p); + Ant.push_back(p); + clitk::PointsUtils::Convert2DTo3D(postMost, input, i, p); + Post.push_back(p); + } + } +} +//--------------------------------------------------------------------