From: dsarrut Date: Tue, 15 Feb 2011 11:31:14 +0000 (+0000) Subject: Station8 X-Git-Tag: v1.2.0~244 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=9de203e39ab0d5c9779bc0827abdb8b5d8a58975;p=clitk.git Station8 --- diff --git a/segmentation/clitkExtractLymphStation_8.txx b/segmentation/clitkExtractLymphStation_8.txx new file mode 100644 index 0000000..469ff2f --- /dev/null +++ b/segmentation/clitkExtractLymphStation_8.txx @@ -0,0 +1,823 @@ + +#include +#include + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_SI_Limits() +{ + /* + Station 8: paraeosphageal nodes + + "Superiorly, Station 8 begins at the level of the carina and, + therefore, abuts Station 3P (Fig. 3C). Inferiorly, the lower limit + of Station 8 was not specified in the Mountain and Dresler + classification (1). We delineate this volume until it reaches the + gastroesphogeal junction. " + + => new classification IASCL 2009: + + "Lower border: the diaphragm" + + "Upper border: the upper border of the lower lobe bronchus on the + left; the lower border of the bronchus intermedius on the right" + + */ + StartNewStep("[Station8] Inf/Sup mediastinum limits with carina/diaphragm junction"); + + // 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); + // 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"); + + // 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"); + + /* Crop support : + Superior limit = carina + Inferior limit = DiaphragmInferiorLimit (old=Gastroesphogeal junction) */ + m_Working_Support = + clitk::CropImageAlongOneAxis(m_Working_Support, 2, + m_DiaphragmInferiorLimit, //GOjunctionZ, + m_CarinaZ, true, + GetBackgroundValue()); + + // Remove some structures that we know are excluded + MaskImagePointer VertebralBody = + GetAFDB()->template GetImage ("VertebralBody"); + MaskImagePointer Aorta = + GetAFDB()->template GetImage ("Aorta"); + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + 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(); + m_Working_Support = boolFilter->GetOutput(); + + // Done. + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_AP_Limits() +{ + /* + Station 8: paraeosphageal nodes + + Anteriorly, it is in contact with Station 7 and the + left main stem bronchus in its superior aspect (Fig. 3C–H) and + with the heart more inferiorly. 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). + + New classification IASCL 2009 : + + "Nodes lying adjacent to the wall of the esophagus and to the + right or left of the midline, excluding subcarinal nodes (S7) + Upper" + + */ + + StartNewStep("[Station8] Post limits with vertebral body"); + MaskImagePointer VertebralBody = + GetAFDB()->template GetImage ("VertebralBody"); + + // Consider vertebral body slice by slice + typedef clitk::ExtractSliceFilter ExtractSliceFilterType; + typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); + typedef typename ExtractSliceFilterType::SliceType SliceType; + std::vector vertebralSlices; + extractSliceFilter->SetInput(VertebralBody); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + extractSliceFilter->GetOutputSlices(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 { + // no Foreground in this slice (it appends generally at the + // beginning and the end of the volume). Do nothing in this + // case. + } + } + + // Convert 2D points in slice into 3D points + std::vector vertebralAntPositions; + clitk::PointsUtils::Convert2DTo3DList(vertebralAntPositionBySlice, + VertebralBody, + vertebralAntPositions); + + // DEBUG : write list of points + clitk::WriteListOfLandmarks(vertebralAntPositions, + "vertebralMostAntPositions.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] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[0]; + size[0] = m_Working_Support->GetLargestPossibleRegion().GetSize()[0]; + for(uint i=0; i IteratorType; + IteratorType iter = + IteratorType(m_Working_Support, m_Working_Support->GetLargestPossibleRegion()); + MaskImageIndexType index; + // Consider some cm posterior to most anterior positions (usually + // 1 cm). + vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheSpine(); + // Get index of this point + m_Working_Support->TransformPhysicalPointToIndex(vertebralAntPositions[i], index); + // Compute region (a single slice) + start[2] = index[2]; + start[1] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[1]+index[1]; + size[1] = m_Working_Support->GetLargestPossibleRegion().GetSize()[1]-start[1]; + region.SetSize(size); + region.SetIndex(start); + // Fill region + if (m_Working_Support->GetLargestPossibleRegion().IsInside(start)) { + itk::ImageRegionIterator it(m_Working_Support, region); + it.GoToBegin(); + while (!it.IsAtEnd()) { + it.Set(GetBackgroundValue()); + ++it; + } + } + } + + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; + + //-------------------------------------------------------------------- + StartNewStep("[Station8] Ant limits with S7 above Carina"); + /* + Anteriorly, it is in contact with Station 7 and the + left main stem bronchus in its superior aspect (Fig. 3C–H) and + with the heart more inferiorly. + + 1) line post wall bronchi (S7), till originRMLB + - LUL bronchus : to detect in trachea. But not needed here ?? + 2) heart ! -> not delineated. + ==> below S7, crop CT not to far away (ant), then try with threshold + + 1) ==> how to share with S7 ? Prepare both support at the same time ! + NEED vector of initial support for each station ? No -> map if it exist before, take it !! + + */ + + // 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 + GetBackgroundValue()); + clitk::PrintMemory(GetVerboseMemoryFlag(), "after crop"); + + // Seprate into two main bronchi + MaskImagePointer LeftBronchus; + MaskImagePointer RightBronchus; + + // 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 + // upper slice is composed of at least two labels + typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion()); + iter.SetFirstDirection(0); + iter.SetSecondDirection(1); + iter.GoToReverseBegin(); // Start from the end (because image is IS not SI) + int maxLabel=0; + while (!iter.IsAtReverseEndOfSlice()) { + while (!iter.IsAtReverseEndOfLine()) { + if (iter.Get() > maxLabel) maxLabel = iter.Get(); + --iter; + } + iter.PreviousLine(); + } + if (maxLabel < 2) { + clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort"); + } + + // Compute 3D centroids of both parts to identify the left from the + // right bronchus + std::vector c; + clitk::ComputeCentroids(m_Working_Trachea, GetBackgroundValue(), c); + 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; } + + // Select LeftLabel (set one label to Backgroundvalue) + LeftBronchus = + SetBackground(m_Working_Trachea, m_Working_Trachea, + rightLabel, GetBackgroundValue(), false); + RightBronchus = + SetBackground(m_Working_Trachea, m_Working_Trachea, + leftLabel, GetBackgroundValue(), false); + // Crop images + LeftBronchus = clitk::AutoCrop(LeftBronchus, GetBackgroundValue()); + RightBronchus = clitk::AutoCrop(RightBronchus, 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); + + + // 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()); + // We dont need Carina structure from now + OriginOfRightMiddleLobeBronchus->Delete(); + clitk::PrintMemory(GetVerboseMemoryFlag(), "after delete OriginOfRightMiddleLobeBronchus"); + + + LeftBronchus = + clitk::CropImageBelow(LeftBronchus, 2, + m_OriginOfRightMiddleLobeBronchusZ, + true, // AutoCrop + GetBackgroundValue()); + RightBronchus = + clitk::CropImageBelow(RightBronchus, 2, + m_OriginOfRightMiddleLobeBronchusZ, + true, // AutoCrop + GetBackgroundValue()); + + // 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, + m_AntMostInRightBronchus, m_PostMostInRightBronchus); + + // DEBUG : write the list of points + ListOfPointsType v; + 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"); + + 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.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() ); + clitk::WriteListOfLandmarks(v, "RightBronchusPoints.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(); + } + + // Keep main 3D CCL : + m_Working_Support = Labelize(m_Working_Support, 0, false, 10); + m_Working_Support = KeepLabels(m_Working_Support, + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + + // Autocrop + m_Working_Support = clitk::AutoCrop(m_Working_Support, GetBackgroundValue()); + + // End of step + StopCurrentStep(m_Working_Support); + // m_ListOfStations["8"] = m_Working_Support; + + //-------------------------------------------------------------------- + StartNewStep("[Station8] Ant limits arround esophagus below Carina"); + /* + Consider Esophagus, dilate it and remove ant part. It remains part + on L & R, than can be partly removed by cutting what remains at + right of vertebral body. + */ + + // Get Esophagus + DD("Esophagus"); + MaskImagePointer Esophagus = GetAFDB()->template GetImage("Esophagus"); + clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Esophagus"); + + // Crop Esophagus : keep only below the OriginOfRightMiddleLobeBronchusZ + Esophagus = + clitk::CropImageAbove(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 + typedef clitk::SliceBySliceRelativePositionFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->VerboseStepFlagOff(); + relPosFilter->SetInput(m_Working_Support); + relPosFilter->SetInputObject(Esophagus); + relPosFilter->AddOrientationTypeString("P"); + relPosFilter->InverseOrientationFlagOff(); + relPosFilter->SetDirection(2); // Z axis + relPosFilter->UniqueConnectedComponentBySliceOff(); + relPosFilter->SetIntermediateSpacing(3); + relPosFilter->ResampleBeforeRelativePositionFilterOn(); + relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8()); + relPosFilter->RemoveObjectFlagOff(); + 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"); + + clitk::PrintMemory(GetVerboseMemoryFlag(), "after autocrop"); + + // 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"); + + // 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"); + 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"); + 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"); + + // 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); + + // DEBUG POINTS + std::ofstream osp; + openFileForWriting(osp, "point_s8.txt"); + osp << "LANDMARKS1" << std::endl; + + // Loop slices + MaskImagePointType p; + int pi=0; + 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 last point out of the mediastinum on this line : + typename MaskSliceType::IndexType index; + mediast_slices[i]->TransformPhysicalPointToIndex(p_maxRight_Vertebra, index); + double distance = 0.0; + while (mediast_slices[i]->GetPixel(index) != GetBackgroundValue()) { + index[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; + } + + // Choose the most extrema one + if (p_maxRight_Vertebra[0] < p_maxRight_Eso[0]) { + p_maxRight = p_maxRight_Vertebra; + } + else p_maxRight = p_maxRight_Eso; + + // -------------------------------------------------------------------------- + + + /* + // 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); + */ + + + // 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); + + // Remove point too at RIGHT + if (p[0] < p_maxRight[0]) { + iter.Set(GetBackgroundValue()); + } + + /* + // Remove point from foreground if too right or to high + if ((p[1] > p_post[1]) && (p[0] > s_left[0])) { + iter.Set(GetBackgroundValue()); + } + */ + + } + ++iter; + } + // } // end if f + // s++; + // sm++; + } + osp.close(); + + + // Joint slices + DD("after loop"); + m_Working_Support = clitk::JoinSlices(slices, m_Working_Support, 2); + clitk::PrintMemory(GetVerboseMemoryFlag(), "after JoinSlices"); + + // DEBUG + m_ListOfStations["8"] = m_Working_Support; + return; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus) +{ + // Check if Esophagus is delineated at least until GOjunction. Now, + // because we use AutoCrop, Origin[2] gives this max inferior + // position. + // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2); + + if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) { + std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it." + << std::endl; + double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit; + + // Pad with few more slices + typedef itk::MirrorPadImageFilter PadFilterType; + typename PadFilterType::Pointer padFilter = PadFilterType::New(); + padFilter->SetInput(Esophagus); + MaskImageSizeType b; + b[0] = 0; b[1] = 0; b[2] = (uint)ceil(extraSize/Esophagus->GetSpacing()[2])+1; + padFilter->SetPadLowerBound(b); + padFilter->Update(); + Esophagus = padFilter->GetOutput(); + } + return Esophagus; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_8_LR_Limits() +{ + /* + Station 8: paraeosphageal nodes + + Laterally, it is within the pleural envelope and again abuts the + descending aorta on the left. Reasonably, the delineation of + Station 8 is limited to the soft tissue surrounding the esophagus + (Fig. 3C–H). + */ + + StartNewStep("[Station8] Right limits (around esophagus)"); + // Get Esophagus + MaskImagePointer Esophagus = GetAFDB()->template GetImage("Esophagus"); + + // Autocrop to get first slice with starting Esophagus + Esophagus = clitk::AutoCrop(Esophagus, GetBackgroundValue()); + + // Dilate + // LR dilatation -> large to keep point inside + // AP dilatation -> few mm + // SI dilatation -> enough to cover Diaphragm (old=GOjunctionZ) + MaskImagePointType radiusInMM = GetEsophagusDiltationForRight(); + Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus); + Esophagus = clitk::Dilate(Esophagus, + radiusInMM, + GetBackgroundValue(), + GetForegroundValue(), true); + writeImage(Esophagus, "dilateEso2.mhd"); + + writeImage(m_Working_Support, "before-relpos2.mhd"); + + // Remove Right (Left on the screen) part according to this + // dilatated esophagus + typedef clitk::SliceBySliceRelativePositionFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->VerboseStepFlagOff(); + relPosFilter->WriteStepFlagOff(); + relPosFilter->SetInput(m_Working_Support); + relPosFilter->SetInputObject(Esophagus); + relPosFilter->AddOrientationTypeString("L"); + relPosFilter->InverseOrientationFlagOn(); // Not Left to + relPosFilter->SetDirection(2); // Z axis + relPosFilter->UniqueConnectedComponentBySliceOff(); // important + relPosFilter->SetIntermediateSpacing(4); + relPosFilter->ResampleBeforeRelativePositionFilterOn(); + relPosFilter->SetFuzzyThreshold(0.9); // remove few part only + relPosFilter->RemoveObjectFlagOff(); + relPosFilter->Update(); + m_Working_Support = relPosFilter->GetOutput(); + + // Get a single 3D CCL + m_Working_Support = Labelize(m_Working_Support, 0, false, 10); + m_Working_Support = KeepLabels(m_Working_Support, + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + + + /* + // Re-Add post to Esophagus -> sometimes previous relpos remove some + // correct part below esophagus. + MaskImagePointer Esophagus = GetAFDB()->template GetImage("Esophagus"); + EnlargeEsophagusDilatationRadiusInferiorly(Esophagus); + writeImage(Esophagus, "e-again.mhd"); + typedef clitk::SliceBySliceRelativePositionFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->VerboseStepOff(); + relPosFilter->WriteStepOff(); + relPosFilter->SetInput(m_Working_Support); + relPosFilter->SetInputObject(Esophagus); + relPosFilter->SetOrientationTypeString("P"); + relPosFilter->InverseOrientationFlagOff(); + relPosFilter->SetDirection(2); // Z axis + relPosFilter->UniqueConnectedComponentBySliceOff(); // important + relPosFilter->SetIntermediateSpacing(4); + relPosFilter->ResampleBeforeRelativePositionFilterOn(); + relPosFilter->CombineWithOrFlagOn(); + relPosFilter->SetFuzzyThreshold(0.9); // remove few part only + relPosFilter->RemoveObjectFlagOff(); + relPosFilter->Update(); + m_Working_Support = relPosFilter->GetOutput(); + */ + + StopCurrentStep(m_Working_Support); + m_ListOfStations["8"] = m_Working_Support; +} +//--------------------------------------------------------------------