From 5ff46b355b778ae675ca60085430c062ea6a3945 Mon Sep 17 00:00:00 2001 From: dsarrut Date: Wed, 15 Dec 2010 09:05:00 +0000 Subject: [PATCH] various update --- segmentation/clitkExtractLymphStation_4RL.txx | 100 ++++++ segmentation/clitkExtractLymphStation_7.txx | 313 ++++++++++++++++++ 2 files changed, 413 insertions(+) create mode 100644 segmentation/clitkExtractLymphStation_4RL.txx create mode 100644 segmentation/clitkExtractLymphStation_7.txx diff --git a/segmentation/clitkExtractLymphStation_4RL.txx b/segmentation/clitkExtractLymphStation_4RL.txx new file mode 100644 index 0000000..a4478a6 --- /dev/null +++ b/segmentation/clitkExtractLymphStation_4RL.txx @@ -0,0 +1,100 @@ +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_4RL_SI_Limits() +{ + /* SupInf limits : + - top of aortic arch + - ends at the upper lobe bronchus or where the right pulmonary artery crosses the midline of the mediastinum + */ + + // Local variables + double m_TopOfAorticArchInMM; + double m_UpperLobeBronchusZPositionInMM; + double m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM; + + // Get Inputs + m_TopOfAorticArchInMM = GetAFDB()->GetPoint3D("topOfAorticArch", 2); + DD(m_TopOfAorticArchInMM); + m_UpperLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("rightUpperLobeBronchus", 2); + DD(m_UpperLobeBronchusZPositionInMM); + m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM = GetAFDB()->GetPoint3D("rightPulmoArteryCrossesMidMediastinum", 2); + DD(m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM); + + /* Crop support */ + StartNewStep("Inf/Sup mediastinum limits with aortic arch/upperLBronchus"); + double inf = std::max(m_UpperLobeBronchusZPositionInMM, m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM); + m_Working_Support = + clitk::CropImageAlongOneAxis(m_Support, 2, + inf, + m_TopOfAorticArchInMM, true, + GetBackgroundValue()); + StopCurrentStep(m_Working_Support); + + m_Station4RL = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_4RL_LR_Limits() +{ + // 4R first + + // Left : midline of the trachea + // Right : "- upper part : contained within the pleural envelope + //- intermediate section : medial to the superior vena cava and the arch of the azygos vein + // - very caudal part : right upper lobe pulmonary vein" + // AAV ?? + + // Constraint at right from the SVC + MaskImagePointer svc = GetAFDB()->template GetImage("SVC"); + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); + relPosFilter->VerboseStepOff(); + relPosFilter->WriteStepOff(); + relPosFilter->SetInput(m_Working_Support); + relPosFilter->SetInputObject(svc); + relPosFilter->SetOrientationType(RelPosFilterType::RightTo); + relPosFilter->SetIntermediateSpacing(2); + relPosFilter->SetFuzzyThreshold(0.3); + relPosFilter->Update(); + m_Working_Support = relPosFilter->GetOutput(); + m_Station4RL = m_Working_Support; + + // Left -> midline of the trachea + // slice by slice, find X coord of 2D centroid (?) + // check with previous line in order to not move too fast + + // skeleton ? -> need path description ? follow from slice to slice + // OR CENTROID at each slice ? + + // Crop trachea + // Extract list of slice from trachea + // Loop slice -> Get centroid crop along line (BB limit) -> two supports + + // Crop the trachea like the current support + MaskImagePointer crop_trachea = + clitk::ResizeImageLike(m_Trachea, m_Working_Support, GetBackgroundValue()); + writeImage(crop_trachea, "croptrachea.mhd"); + + // Extract all the slices + typedef clitk::ExtractSliceFilter ExtractSliceFilterType; + typedef typename ExtractSliceFilterType::SliceType SliceType; + typename ExtractSliceFilterType::Pointer + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(crop_trachea); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + std::vector trachea_slices; + extractSliceFilter->GetOutputSlices(trachea_slices); + + + +} +//-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_7.txx b/segmentation/clitkExtractLymphStation_7.txx new file mode 100644 index 0000000..e1cbe1f --- /dev/null +++ b/segmentation/clitkExtractLymphStation_7.txx @@ -0,0 +1,313 @@ +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_7_SI_Limits() +{ + // Local variables + double m_CarinaZPositionInMM; + double m_MiddleLobeBronchusZPositionInMM; + + // Get Inputs + m_Trachea = GetAFDB()->template GetImage ("trachea"); + m_CarinaZPositionInMM = GetAFDB()->GetPoint3D("carina", 2); + DD(m_CarinaZPositionInMM); + m_MiddleLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("rightMiddleLobeBronchus", 2); + DD(m_MiddleLobeBronchusZPositionInMM); + + /* Crop support : + Superior limit = carina + Inferior limit = origin right middle lobe bronchus */ + StartNewStep("Inf/Sup mediastinum limits with carina/bronchus"); + m_Working_Support = + clitk::CropImageAlongOneAxis(m_Support, 2, + m_MiddleLobeBronchusZPositionInMM, + m_CarinaZPositionInMM, true, + GetBackgroundValue()); + StopCurrentStep(m_Working_Support); + + /* Crop trachea + Superior limit = carina + Inferior limit = origin right middle lobe bronchus*/ + StartNewStep("Inf/Sup trachea limits with carina/bronchus"); + m_working_trachea = + clitk::CropImageAlongOneAxis(m_Trachea, 2, + m_MiddleLobeBronchusZPositionInMM, + m_CarinaZPositionInMM, true, + GetBackgroundValue()); + StopCurrentStep(m_working_trachea); + + m_Station7 = m_Working_Support; +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_7_RL_Limits() +{ + // ---------------------------------------------------------------- + // Separate trachea in two CCL + StartNewStep("Separate trachea under carina"); + + // Labelize and consider two main labels + m_working_trachea = Labelize(m_working_trachea, 0, true, 1); + + // Carina position must at the first slice that separate the two main bronchus (not superiorly) + // Check that 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()) { + // DD(iter.GetIndex()); + 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 centroid of both parts to identify the left from the right bronchus + typedef long LabelType; + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(GetBackgroundValue()); + imageToLabelFilter->SetInput(m_working_trachea); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->Update(); + typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); + + ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid(); + ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid(); + + ImagePixelType leftLabel; + ImagePixelType rightLabel; + if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; } + else { leftLabel = 2; rightLabel = 1; } + DD(leftLabel); + DD(rightLabel); + + StopCurrentStep(m_working_trachea); + + //----------------------------------------------------- + StartNewStep("Left limits with bronchus (slice by slice)"); + // Select LeftLabel (set one label to Backgroundvalue) + m_LeftBronchus = + SetBackground(m_working_trachea, m_working_trachea, + rightLabel, GetBackgroundValue()); + m_RightBronchus = + SetBackground(m_working_trachea, m_working_trachea, + leftLabel, GetBackgroundValue()); + writeImage(m_LeftBronchus, "left.mhd"); + writeImage(m_RightBronchus, "right.mhd"); + + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, + m_LeftBronchus, 2, + GetFuzzyThreshold(), "RightTo", + true, 4); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, + m_RightBronchus, + 2, GetFuzzyThreshold(), "LeftTo", + true, 4); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, + m_LeftBronchus, + 2, GetFuzzyThreshold(), "AntTo", + true, 4, true); // NOT + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, + m_RightBronchus, + 2, GetFuzzyThreshold(), "AntTo", + true, 4, true); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, + m_LeftBronchus, + 2, GetFuzzyThreshold(), "PostTo", + true, 4, true); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, + m_RightBronchus, + 2, GetFuzzyThreshold(), "PostTo", + true, 4, true); + m_Station7 = m_Working_Support; + StopCurrentStep(m_Station7); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_7_Posterior_Limits() +{ + StartNewStep("Posterior limits"); + + // Left Bronchus slices + typedef clitk::ExtractSliceFilter ExtractSliceFilterType; + typedef typename ExtractSliceFilterType::SliceType SliceType; + typename ExtractSliceFilterType::Pointer + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(m_LeftBronchus); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + std::vector leftBronchusSlices; + extractSliceFilter->GetOutputSlices(leftBronchusSlices); + + // Right Bronchus slices + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(m_RightBronchus); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + std::vector rightBronchusSlices ; + extractSliceFilter->GetOutputSlices(rightBronchusSlices); + + assert(leftBronchusSlices.size() == rightBronchusSlices.size()); + + std::vector leftPoints; + std::vector rightPoints; + for(uint i=0; i(leftBronchusSlices[i], 0, true, 10); + leftBronchusSlices[i] = KeepLabels(leftBronchusSlices[i], + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + rightBronchusSlices[i] = Labelize(rightBronchusSlices[i], 0, true, 10); + rightBronchusSlices[i] = KeepLabels(rightBronchusSlices[i], + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + double distance_max_from_center_point = 15; + + // ------- Find point in left Bronchus ------- + // find right most point in left = rightMost + SliceType::PointType a; + SliceType::PointType rightMost = + clitk::FindExtremaPointInAGivenDirection(leftBronchusSlices[i], + GetBackgroundValue(), + 0, false, a, 0); + // find post most point in left, not far away from rightMost + SliceType::PointType p = + clitk::FindExtremaPointInAGivenDirection(leftBronchusSlices[i], + GetBackgroundValue(), + 1, false, rightMost, + distance_max_from_center_point); + MaskImageType::PointType pp; + pp[0] = p[0]; pp[1] = p[1]; + pp[2] = i*m_LeftBronchus->GetSpacing()[2] + m_LeftBronchus->GetOrigin()[2]; + leftPoints.push_back(pp); + + // ------- Find point in right Bronchus ------- + // find left most point in right = leftMost + SliceType::PointType leftMost = + clitk::FindExtremaPointInAGivenDirection(rightBronchusSlices[i], + GetBackgroundValue(), + 0, true, a, 0); + // find post most point in left, not far away from leftMost + p = clitk::FindExtremaPointInAGivenDirection(rightBronchusSlices[i], + GetBackgroundValue(), + 1, false, leftMost, + distance_max_from_center_point); + pp[0] = p[0]; pp[1] = p[1]; + pp[2] = i*m_RightBronchus->GetSpacing()[2] + m_RightBronchus->GetOrigin()[2]; + rightPoints.push_back(pp); + } + + // DEBUG + std::ofstream osl; + openFileForWriting(osl, "left.txt"); + osl << "LANDMARKS1" << std::endl; + std::ofstream osr; + openFileForWriting(osr, "right.txt"); + osr << "LANDMARKS1" << std::endl; + for(uint i=0; i to accelerate, start with formula, when change sign -> stop and fill + */ + typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + SliceIteratorType iter = SliceIteratorType(m_Working_Support, + m_Working_Support->GetLargestPossibleRegion()); + iter.SetFirstDirection(0); + iter.SetSecondDirection(1); + iter.GoToBegin(); + int i=0; + MaskImageType::PointType A; + MaskImageType::PointType B; + MaskImageType::PointType C; + while (!iter.IsAtEnd()) { + A = leftPoints[i]; + B = rightPoints[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 (!iter.IsAtEndOfSlice()) { + while (!iter.IsAtEndOfLine()) { + // Very slow, I know ... but image should be very small + m_Working_Support->TransformIndexToPhysicalPoint(iter.GetIndex(), C); + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + if (s == 0) iter.Set(2); + if (isPositive) { + if (s > 0) iter.Set(GetBackgroundValue()); + } + else { + if (s < 0) iter.Set(GetBackgroundValue()); + } + ++iter; + } + iter.NextLine(); + } + iter.NextSlice(); + ++i; + } + + //----------------------------------------------------- + // StartNewStep("Anterior limits"); + + + // MISSING FROM NOW + + // Station 4R, Station 4L, the right pulmonary artery, and/or the + // left superior pulmonary vein + + + //----------------------------------------------------- + //----------------------------------------------------- + // ALSO SUBSTRACT ARTERY/VEIN (initially in the support) + + + // Set output + m_Station7 = m_Working_Support; +} +//-------------------------------------------------------------------- -- 2.47.1