From: David Sarrut Date: Mon, 25 Jul 2011 12:49:25 +0000 (+0200) Subject: Update supports & stations 3A,3P X-Git-Tag: v1.3.0~265 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=2061b9e1d00ffe188af506917c2d3fa4f65121ca;p=clitk.git Update supports & stations 3A,3P --- diff --git a/segmentation/clitkExtractLymphStation_2RL.txx b/segmentation/clitkExtractLymphStation_2RL.txx index 71298eb..ee9d2f3 100644 --- a/segmentation/clitkExtractLymphStation_2RL.txx +++ b/segmentation/clitkExtractLymphStation_2RL.txx @@ -11,83 +11,6 @@ // itk #include -//-------------------------------------------------------------------- -template -class comparePointsX { -public: - bool operator() (PointType i, PointType j) { return (i[0] -class comparePointsWithAngle { -public: - bool operator() (PairType i, PairType j) { return (i.second < j.second); } -}; -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void HypercubeCorners(std::vector > & out) { - std::vector > previous; - HypercubeCorners(previous); - out.resize(previous.size()*2); - for(unsigned int i=0; i p; - if (i -void HypercubeCorners<1>(std::vector > & out) { - out.resize(2); - out[0][0] = 0; - out[1][0] = 1; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image, - std::vector & bounds) -{ - // Get image max/min coordinates - const unsigned int dim=ImageType::ImageDimension; - typedef typename ImageType::PointType PointType; - typedef typename ImageType::IndexType IndexType; - PointType min_c, max_c; - IndexType min_i, max_i; - min_i = image->GetLargestPossibleRegion().GetIndex(); - for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; - image->TransformIndexToPhysicalPoint(min_i, min_c); - image->TransformIndexToPhysicalPoint(max_i, max_c); - - // Get corners coordinates - HypercubeCorners(bounds); - for(unsigned int i=0; i void @@ -104,6 +27,34 @@ ExtractStation_2RL_SetDefaultValues() //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL() +{ + if (CheckForStation("2RL")) { + ExtractStation_2RL_SI_Limits(); + ExtractStation_2RL_Post_Limits(); + + ExtractStation_2RL_Ant_Limits2(); + //ExtractStation_2RL_Ant_Limits(); + + ExtractStation_2RL_LR_Limits(); + ExtractStation_2RL_Remove_Structures(); + ExtractStation_2RL_SeparateRL(); + + // Store image filenames into AFDB + writeImage(m_ListOfStations["2R"], "seg/Station2R.mhd"); + writeImage(m_ListOfStations["2L"], "seg/Station2L.mhd"); + GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd"); + GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd"); + WriteAFDB(); + } +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -300,299 +251,23 @@ ExtractStation_2RL_Ant_Limits() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_Ant_Limits2() +ExtractStation_2RL_Ant_Limits2() { - // ----------------------------------------------------- - /* Rod says: "The anterior border, as with the Atlas – UM, is - posterior to the vessels (right subclavian vein, left - brachiocephalic vein, right brachiocephalic vein, left subclavian - artery, left common carotid artery and brachiocephalic trunk). - These vessels are not included in the nodal station. The anterior - border is drawn to the midpoint of the vessel and an imaginary - line joins the middle of these vessels. Between the vessels, - station 2 is in contact with station 3a." */ - // ----------------------------------------------------- StartNewStep("[Station 2RL] Ant limits with vessels centroids"); - - /* Here, we consider the vessels form a kind of anterior barrier. We - link all vessels centroids and remove what is post to it. - - select the list of structure - vessel1 = BrachioCephalicArtery - vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) - vessel3 = CommonCarotidArtery - vessel4 = SubclavianArtery - other = Thyroid - other = Aorta - - crop images as needed - - slice by slice, choose the CCL and extract centroids - - slice by slice, sort according to polar angle wrt Trachea centroid. - - slice by slice, link centroids and close contour - - remove outside this contour - - merge with support - */ - - // Read structures - MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage("BrachioCephalicArtery"); - MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); - MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); - MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); - MaskImagePointer Thyroid = GetAFDB()->template GetImage("Thyroid"); - MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); - - // Resize all structures like support - BrachioCephalicArtery = - clitk::ResizeImageLike(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue()); - CommonCarotidArtery = - clitk::ResizeImageLike(CommonCarotidArtery, m_Working_Support, GetBackgroundValue()); - SubclavianArtery = - clitk::ResizeImageLike(SubclavianArtery, m_Working_Support, GetBackgroundValue()); - Thyroid = - clitk::ResizeImageLike(Thyroid, m_Working_Support, GetBackgroundValue()); - Aorta = - clitk::ResizeImageLike(Aorta, m_Working_Support, GetBackgroundValue()); - BrachioCephalicVein = - clitk::ResizeImageLike(BrachioCephalicVein, m_Working_Support, GetBackgroundValue()); - Trachea = - clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); - - // Extract slices - std::vector slices_BrachioCephalicArtery; - clitk::ExtractSlices(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery); - std::vector slices_BrachioCephalicVein; - clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BrachioCephalicVein); - std::vector slices_CommonCarotidArtery; - clitk::ExtractSlices(CommonCarotidArtery, 2, slices_CommonCarotidArtery); - std::vector slices_SubclavianArtery; - clitk::ExtractSlices(SubclavianArtery, 2, slices_SubclavianArtery); - std::vector slices_Thyroid; - clitk::ExtractSlices(Thyroid, 2, slices_Thyroid); - std::vector slices_Aorta; - clitk::ExtractSlices(Aorta, 2, slices_Aorta); - std::vector slices_Trachea; - clitk::ExtractSlices(Trachea, 2, slices_Trachea); - unsigned int n= slices_BrachioCephalicArtery.size(); - // Get the boundaries of one slice - std::vector bounds; - ComputeImageBoundariesCoordinates(slices_BrachioCephalicArtery[0], bounds); - - // For all slices, for all structures, find the centroid and build the contour - // List of 3D points (for debug) - std::vector p3D; - - vtkSmartPointer append = vtkSmartPointer::New(); - for(unsigned int i=0; i(slices_CommonCarotidArtery[i], - GetBackgroundValue(), true, 1); - slices_SubclavianArtery[i] = Labelize(slices_SubclavianArtery[i], - GetBackgroundValue(), true, 1); - slices_BrachioCephalicArtery[i] = Labelize(slices_BrachioCephalicArtery[i], - GetBackgroundValue(), true, 1); - slices_BrachioCephalicVein[i] = Labelize(slices_BrachioCephalicVein[i], - GetBackgroundValue(), true, 1); - slices_Thyroid[i] = Labelize(slices_Thyroid[i], - GetBackgroundValue(), true, 1); - slices_Aorta[i] = Labelize(slices_Aorta[i], - GetBackgroundValue(), true, 1); - - // Search centroids - std::vector points2D; - std::vector centroids1; - std::vector centroids2; - std::vector centroids3; - std::vector centroids4; - std::vector centroids5; - std::vector centroids6; - ComputeCentroids(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1); - ComputeCentroids(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2); - ComputeCentroids(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3); - ComputeCentroids(slices_Thyroid[i], GetBackgroundValue(), centroids4); - ComputeCentroids(slices_Aorta[i], GetBackgroundValue(), centroids5); - ComputeCentroids(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6); - - // BrachioCephalicVein -> when it is separated into two CCL, we - // only consider the most at Right one - if (centroids6.size() > 2) { - if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2); - else centroids6.erase(centroids6.begin()+1); - } - - // BrachioCephalicVein -> when SubclavianArtery has 2 CCL - // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein - if ((centroids3.size() ==1) && (centroids2.size() > 2)) { - centroids6.clear(); - } - - for(unsigned int j=1; j centroids_trachea; - ComputeCentroids(slices_Trachea[i], GetBackgroundValue(), centroids_trachea); - typedef std::pair PointAngleType; - std::vector angles; - for(unsigned int j=0; j0) angle = atan(y/x); - if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; - if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; - if (x==0) { - if (y>0) angle = M_PI/2.0; - if (y<0) angle = -M_PI/2.0; - if (y==0) angle = 0; - } - angle = clitk::rad2deg(angle); - // Angle is [-180;180] wrt the X axis. We change the X axis to - // be the vertical line, because we want to sort from Right to - // Left from Post to Ant. - if (angle>0) - angle = (270-angle); - if (angle<0) { - angle = -angle-90; - if (angle<0) angle = 360-angle; - } - PointAngleType a; - a.first = points2D[j]; - a.second = angle; - angles.push_back(a); - } - - // Do nothing if less than 2 points --> n - if (points2D.size() < 3) { //continue; - continue; - } - - // Sort points2D according to polar angles - std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); - for(unsigned int j=0; j toadd; - unsigned int index = 0; - double dmax = 5; - while (indexdmax) { - - MaskSlicePointType b; - b[0] = a[0]+(c[0]-a[0])/2.0; - b[1] = a[1]+(c[1]-a[1])/2.0; - - // Compute distance to trachea centroid - MaskSlicePointType m = centroids_trachea[1]; - double da = m.EuclideanDistanceTo(a); - double db = m.EuclideanDistanceTo(b); - //double dc = m.EuclideanDistanceTo(c); - - // Mean distance, find point on the line from trachea centroid - // to b - double alpha = (da+db)/2.0; - MaskSlicePointType v; - double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); - v[0] = (b[0]-m[0])/n; - v[1] = (b[1]-m[1])/n; - MaskSlicePointType r; - r[0] = m[0]+alpha*v[0]; - r[1] = m[1]+alpha*v[1]; - points2D.insert(points2D.begin()+index+1, r); - } - else { - index++; - } - } - // DDV(points2D, points2D.size()); - - // Add some points to close the contour - // - H line towards Right - MaskSlicePointType p = points2D[0]; - p[0] = bounds[0][0]; - points2D.insert(points2D.begin(), p); - // - corner Right/Post - p = bounds[0]; - points2D.insert(points2D.begin(), p); - // - H line towards Left - p = points2D.back(); - p[0] = bounds[2][0]; - points2D.push_back(p); - // - corner Right/Post - p = bounds[2]; - points2D.push_back(p); - // Close contour with the first point - points2D.push_back(points2D[0]); - // DDV(points2D, points2D.size()); - - // Build 3D points from the 2D points - std::vector points3D; - clitk::PointsUtils::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D); - for(unsigned int x=0; x mesh = Build3DMeshFrom2DContour(points3D); - append->AddInput(mesh); - } - - // DEBUG: write points3D - clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); - - // Build the final 3D mesh form the list 2D mesh - append->Update(); - vtkSmartPointer mesh = append->GetOutput(); - - // Debug, write the mesh - /* - vtkSmartPointer w = vtkSmartPointer::New(); - w->SetInput(mesh); - w->SetFileName("bidon.vtk"); - w->Write(); - */ + // WARNING, as I used "And" after, empty slice in binarizedContour + // lead to remove part of the support, although we want to keep + // unchanged. So we decide to ResizeImageLike but pad with + // ForegroundValue instead of BG + + // Get or compute the binary mask that separate Ant/Post part + // according to vessels + MaskImagePointer binarizedContour = FindAntPostVessels2(); + binarizedContour = clitk::ResizeImageLike(binarizedContour, + m_Working_Support, + GetForegroundValue()); - // Compute a single binary 3D image from the list of contours - clitk::MeshToBinaryImageFilter::Pointer filter = - clitk::MeshToBinaryImageFilter::New(); - filter->SetMesh(mesh); - filter->SetLikeImage(m_Working_Support); - filter->Update(); - MaskImagePointer binarizedContour = filter->GetOutput(); - - // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse - ImagePointType p = p3D[2]; // This is the first centroid of the first slice - p[1] += 50; // 50 mm Post from this point must be kept - ImageIndexType index; - binarizedContour->TransformPhysicalPointToIndex(p, index); - bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); - // remove from support typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); @@ -601,10 +276,7 @@ ExtractStation_2RL_Ant_Limits2() boolFilter->SetInput2(binarizedContour); boolFilter->SetBackgroundValue1(GetBackgroundValue()); boolFilter->SetBackgroundValue2(GetBackgroundValue()); - if (isInside) - boolFilter->SetOperationType(BoolFilterType::And); - else - boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->SetOperationType(BoolFilterType::And); boolFilter->Update(); m_Working_Support = boolFilter->GetOutput(); diff --git a/segmentation/clitkExtractLymphStation_3A.txx b/segmentation/clitkExtractLymphStation_3A.txx index a55071c..e5529bf 100644 --- a/segmentation/clitkExtractLymphStation_3A.txx +++ b/segmentation/clitkExtractLymphStation_3A.txx @@ -12,33 +12,42 @@ ExtractStation_3A_SetDefaultValues() //-------------------------------------------------------------------- -template +template void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3A_SI_Limits() +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A() { - // Apex of the chest or Sternum & Carina. - StartNewStep("[Station 3A] Inf/Sup limits with Sternum and Carina"); - - // Get Carina position (has been determined in Station8) - m_CarinaZ = GetAFDB()->GetDouble("CarinaZ"); + if (!CheckForStation("3A")) return; + + // Get the current support + StartNewStep("[Station 3A] Get the current 3A suppport"); + m_Working_Support = m_ListOfSupports["S3A"]; + m_ListOfStations["3A"] = m_Working_Support; + StopCurrentStep(m_Working_Support); - // Get Sternum and search for the upper position - MaskImagePointer Sternum = GetAFDB()->template GetImage("Sternum"); + ExtractStation_3A_AntPost_S5(); + ExtractStation_3A_AntPost_S6(); + ExtractStation_3A_AntPost_Superiorly(); - // Search most sup point - MaskImagePointType ps = Sternum->GetOrigin(); // initialise to avoid warning - clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, ps); - double m_SternumZ = ps[2]+Sternum->GetSpacing()[2]; // One more slice, because it is below this point + ExtractStation_3A_Remove_Structures(); - //* Crop support : - m_Working_Support = - clitk::CropImageAlongOneAxis(m_Working_Support, 2, - m_CarinaZ, m_SternumZ, true, - GetBackgroundValue()); + Remove_Structures("3A", "Aorta"); + Remove_Structures("3A", "SubclavianArteryLeft"); + Remove_Structures("3A", "SubclavianArteryRight"); + Remove_Structures("3A", "Thyroid"); + Remove_Structures("3A", "CommonCarotidArteryLeft"); + Remove_Structures("3A", "CommonCarotidArteryRight"); + Remove_Structures("3A", "BrachioCephalicArtery"); + // Remove_Structures("3A", "BrachioCephalicVein"); ? - StopCurrentStep(m_Working_Support); - m_ListOfStations["3A"] = m_Working_Support; + + // ExtractStation_3A_Ant_Limits(); --> No, already in support; to remove + // ExtractStation_3A_Post_Limits(); --> No, more complex, see Vessels etc + + // Store image filenames into AFDB + writeImage(m_ListOfStations["3A"], "seg/Station3A.mhd"); + GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); + WriteAFDB(); } //-------------------------------------------------------------------- @@ -72,14 +81,270 @@ ExtractStation_3A_Post_Limits() StartNewStep("[Station 3A] Post limits with SubclavianArtery"); // Get Sternum, keep posterior part. - MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); + MaskImagePointer SubclavianArteryLeft = + GetAFDB()->template GetImage("SubclavianArteryLeft"); + MaskImagePointer SubclavianArteryRight = + GetAFDB()->template GetImage("SubclavianArteryRight"); + m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, SubclavianArtery, 2, + clitk::SliceBySliceRelativePosition(m_Working_Support, SubclavianArteryLeft, 2, GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo", false, 3, true, false); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, SubclavianArteryRight, 2, + GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo", + false, 3, true, false); + StopCurrentStep(m_Working_Support); + m_ListOfStations["3A"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A_AntPost_S5() +{ + StartNewStep("[Station 3A] Post limits around S5"); + + // First remove post to SVC + MaskImagePointer SVC = GetAFDB()->template GetImage ("SVC"); + + // Trial in 3D -> difficulties superiorly. Stay slice by slice. + /* + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); + relPosFilter->VerboseStepFlagOff(); + relPosFilter->WriteStepFlagOff(); + relPosFilter->SetBackgroundValue(GetBackgroundValue()); + relPosFilter->SetInput(m_Working_Support); + relPosFilter->SetInputObject(SVC); + relPosFilter->AddOrientationTypeString("PostTo"); + relPosFilter->SetInverseOrientationFlag(true); + relPosFilter->SetIntermediateSpacing(4); + relPosFilter->SetIntermediateSpacingFlag(false); + relPosFilter->SetFuzzyThreshold(0.5); + // relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop + relPosFilter->Update(); + m_Working_Support = relPosFilter->GetOutput(); + */ + + // Slice by slice not post to SVC. Use initial spacing + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, SVC, 2, + GetFuzzyThreshold("3A", "SVC"), + "NotPostTo", true, + SVC->GetSpacing()[0], false, false); + + // Consider Aorta, remove Left/Post part ; only around S5 + // Get S5 support and Aorta + MaskImagePointer S5 = m_ListOfSupports["S5"]; + MaskImagePointer Aorta = GetAFDB()->template GetImage ("Aorta"); + Aorta = clitk::ResizeImageLike(Aorta, S5, GetBackgroundValue()); + + // Inferiorly, Aorta has two CCL that merge into a single one when + // S6 appears. Loop on Aorta slices, select the most ant one, detect + // the most ant point. + std::vector slices; + clitk::ExtractSlices(Aorta, 2, slices); + std::vector points; + for(uint i=0; i(slices[i], GetBackgroundValue(), false, 1); + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), c); + assert(c.size() == 3); // only 2 CCL + typename MaskSliceType::PixelType l; + if (c[1][1] > c[2][1]) { // We will remove the label=1 + l = 1; + } + else { + l = 2;// We will remove the label=2 + } + slices[i] = clitk::SetBackground(slices[i], slices[i], l, + GetBackgroundValue(), true); + + // Detect the most ant point + MaskSlicePointType p; + MaskImagePointType pA; + clitk::FindExtremaPointInAGivenDirection(slices[i], GetBackgroundValue(), 1, true, p); + // Set the X coordinate to the X coordinate of the centroid + if (l==1) p[0] = c[2][0]; + else p[0] = c[1][0]; + + // Convert in 3D and store + clitk::PointsUtils::Convert2DTo3D(p, Aorta, i, pA); + points.push_back(pA); + } + + // DEBUG + // MaskImagePointer o = clitk::JoinSlices(slices, Aorta, 2); + // writeImage(o, "o.mhd"); + // clitk::WriteListOfLandmarks(points, "Ant-Aorta.txt"); + + // Remove Post/Left to this point + m_Working_Support = + clitk::SliceBySliceSetBackgroundFromPoints(m_Working_Support, + GetBackgroundValue(), 2, + points, + true, // Set BG if X greater than point[x], and + true); // if Y greater than point[y] + + StopCurrentStep(m_Working_Support); + m_ListOfStations["3A"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A_AntPost_S6() +{ + StartNewStep("[Station 3A] Post limits around S6"); + + // Consider Aorta + MaskImagePointer Aorta = GetAFDB()->template GetImage ("Aorta"); + + // Limits the support to S6 + MaskImagePointer S6 = m_ListOfSupports["S6"]; + Aorta = clitk::ResizeImageLike(Aorta, S6, GetBackgroundValue()); + + // Extend 1cm anteriorly + MaskImagePointType radius; // in mm + radius[0] = 10; + radius[1] = 10; + radius[2] = 0; // required + Aorta = clitk::Dilate(Aorta, radius, GetBackgroundValue(), GetForegroundValue(), false); + + // Not Post to + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, Aorta, 2, + GetFuzzyThreshold("3A", "Aorta"), + "NotPostTo", true, + Aorta->GetSpacing()[0], false, false); + StopCurrentStep(m_Working_Support); m_ListOfStations["3A"] = m_Working_Support; } //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A_AntPost_Superiorly() +{ + StartNewStep("[Station 3A] Post limits superiorly"); + + /* + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage ("BrachioCephalicVein"); + MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage ("BrachioCephalicArtery"); + MaskImagePointer CommonCarotidArteryLeft = GetAFDB()->template GetImage ("CommonCarotidArteryLeft"); + MaskImagePointer CommonCarotidArteryRight = GetAFDB()->template GetImage ("CommonCarotidArteryRight"); + MaskImagePointer SubclavianArteryLeft = GetAFDB()->template GetImage ("SubclavianArteryLeft"); + MaskImagePointer SubclavianArteryRight = GetAFDB()->template GetImage ("SubclavianArteryRight"); + + // Not Post to +#define RP(STRUCTURE) \ + m_Working_Support = \ + clitk::SliceBySliceRelativePosition(m_Working_Support, STRUCTURE, 2, \ + 0.5, \ + "NotPostTo", true, \ + STRUCTURE->GetSpacing()[0], false, false); + + // RP(BrachioCephalicVein); + RP(BrachioCephalicArtery); + RP(CommonCarotidArteryRight); + RP(CommonCarotidArteryLeft); + RP(SubclavianArteryRight); + RP(SubclavianArteryLeft); + */ + + // Get or compute the binary mask that separate Ant/Post part + // according to vessels + MaskImagePointer binarizedContour = FindAntPostVessels2(); + binarizedContour = clitk::ResizeImageLike(binarizedContour, + m_Working_Support, + GetBackgroundValue()); + + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["3A"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A_Remove_Structures() +{ + Remove_Structures("3A", "Aorta"); + Remove_Structures("3A", "SubclavianArteryLeft"); + Remove_Structures("3A", "SubclavianArteryRight"); + Remove_Structures("3A", "Thyroid"); + Remove_Structures("3A", "CommonCarotidArteryLeft"); + Remove_Structures("3A", "CommonCarotidArteryRight"); + Remove_Structures("3A", "BrachioCephalicArtery"); + // Remove_Structures("3A", "BrachioCephalicVein"); ? + + StartNewStep("[Station 3A] Remove part of BrachioCephalicVein"); + // resize like support, extract slices + // while single CCL -> remove + // when two remove only the most post + BrachioCephalicVein = clitk::ResizeImageLike(BrachioCephalicVein, + m_Working_Support, + GetBackgroundValue()); + std::vector slices; + std::vector slices_BCV; + clitk::ExtractSlices(m_Working_Support, 2, slices); + clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BCV); + for(uint i=0; i(slices_BCV[i], 0, true, 1); + // Compute centroids + std::vector centroids; + ComputeCentroids(slices_BCV[i], GetBackgroundValue(), centroids); + if (centroids.size() > 1) { + // Only keep the one most post + typename MaskSliceType::PixelType label; + if (centroids[1][1] > centroids[2][1]) { + label = 1; + } + else { + label = 2; + } + + HERE + + slices_BCV[i] = clitk::SetBackground(slices_BCV[i], slices_BCV[i], + label, + GetBackgroundValue(), true); + } + // Remove from the support + slices[i] = clitk::AndNot(slices[i], slices_BCV[i], GetBackgroundValue()); + } + + // Joint + m_Working_Support = clitk::JoinSlices(slices, m_Working_Support, 2); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["3A"] = m_Working_Support; +} +//-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_3P.txx b/segmentation/clitkExtractLymphStation_3P.txx index 1c9fdc0..58549df 100644 --- a/segmentation/clitkExtractLymphStation_3P.txx +++ b/segmentation/clitkExtractLymphStation_3P.txx @@ -10,40 +10,34 @@ ExtractStation_3P_SetDefaultValues() //-------------------------------------------------------------------- -template +template void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_SI_Limits() +clitk::ExtractLymphStationsFilter:: +ExtractStation_3P() { - /* - Apex of the chest & Carina. - */ - StartNewStep("[Station 3P] Inf/Sup limits with apex of the chest and carina"); + if (!CheckForStation("3P")) return; - // Get Carina position (has been determined in Station8) - m_CarinaZ = GetAFDB()->GetDouble("CarinaZ"); - - // Get Apex of the Chest. The "lungs" structure is autocroped so we - // consider the most superior point. - MaskImagePointer Lungs = GetAFDB()->template GetImage("Lungs"); - MaskImageIndexType index = Lungs->GetLargestPossibleRegion().GetIndex(); - index += Lungs->GetLargestPossibleRegion().GetSize(); - MaskImagePointType p; - Lungs->TransformIndexToPhysicalPoint(index, p); - p[2] -= 5; // 5 mm below because the autocrop is slightly above the apex - double m_ApexOfTheChest = p[2]; - - /* Crop support : - Superior limit = carina - Inferior limit = Apex of the chest */ - m_Working_Support = - clitk::CropImageAlongOneAxis(m_Working_Support, 2, - m_CarinaZ, - m_ApexOfTheChest, true, - GetBackgroundValue()); - - StopCurrentStep(m_Working_Support); + // Get the current support + StartNewStep("[Station 3P] Get the current 3P suppport"); + m_Working_Support = m_ListOfSupports["S3P"]; m_ListOfStations["3P"] = m_Working_Support; + StopCurrentStep(m_Working_Support); + + // LR limits inferiorly + ExtractStation_3P_LR_inf_Limits(); + + // LR limits superiorly => not here for the moment because not + // clear in the def + // ExtractStation_3P_LR_sup_Limits_2(); //TODO + // ExtractStation_3P_LR_sup_Limits(); // old version to change + + ExtractStation_8_Single_CCL_Limits(); // YES 8 ! + ExtractStation_3P_Remove_Structures(); // after CCL + + // Store image filenames into AFDB + writeImage(m_ListOfStations["3P"], "seg/Station3P.mhd"); + GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); + WriteAFDB(); } //-------------------------------------------------------------------- @@ -54,19 +48,14 @@ void clitk::ExtractLymphStationsFilter:: ExtractStation_3P_Remove_Structures() { - /* - 3 - (sup) remove Aorta, VB, subcl - not LR subcl ! -> a séparer LR ? - (inf) remove Eso, Aorta, Azygousvein - */ - StartNewStep("[Station 3P] Remove some structures."); Remove_Structures("3P", "Aorta"); Remove_Structures("3P", "VertebralBody"); - Remove_Structures("3P", "SubclavianArtery"); + Remove_Structures("3P", "SubclavianArteryLeft"); + Remove_Structures("3P", "SubclavianArteryRight"); Remove_Structures("3P", "Esophagus"); - Remove_Structures("3P", "Azygousvein"); + Remove_Structures("3P", "AzygousVein"); Remove_Structures("3P", "Thyroid"); Remove_Structures("3P", "VertebralArtery"); @@ -76,156 +65,6 @@ ExtractStation_3P_Remove_Structures() //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_Ant_Limits() -{ - /* - Ant Post limit : - - Anteriorly, it is in contact with the posterior aspect of Stations - 1–2 superiorly (Fig. 2A–C) and with Stations 4R and 4L inferiorly - (Fig. 2D–I and 3A–C). The anterior limit of Station 3P is kept - posterior to the trachea, which is defined by an imaginary - horizontal line running along the posterior wall of the trachea - (Fig. 2B,E, red line). Posteriorly, it is delineated along the - anterior and lateral borders of the vertebral body until an - imaginary horizontal line running 1 cm posteriorly to the - anterior border of the vertebral body (Fig. 2D). - - 1 - post to the trachea : define an imaginary line just post the - Trachea and remove what is anterior. - - */ - StartNewStep("[Station 3P] Ant limits with Trachea "); - - // Load Trachea - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); - - // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after) - Trachea = - clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); - - // Slice by slice, determine the most post point of the trachea (A) - // and consider a second point on the same horizontal line (B) - std::vector p_A; - std::vector p_B; - std::vector slices; - clitk::ExtractSlices(Trachea, 2, slices); - MaskImagePointType p; - typename MaskSliceType::PointType sp; - for(uint i=0; i(slices[i], 0, true, 100); - slices[i] = KeepLabels(slices[i], GetBackgroundValue(), - GetForegroundValue(), 1, 1, true); - // Find most posterior point - clitk::FindExtremaPointInAGivenDirection(slices[i], GetBackgroundValue(), - 1, false, sp); - clitk::PointsUtils::Convert2DTo3D(sp, Trachea, i, p); - p_A.push_back(p); - p[0] -= 20; - p_B.push_back(p); - } - clitk::WriteListOfLandmarks(p_A, "trachea-post.txt"); - - // Remove Ant part above those lines - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - p_A, p_B, - GetBackgroundValue(), - 1, 10); - // End, release memory - GetAFDB()->template ReleaseImage("Trachea"); - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_Post_Limits() -{ - /* - Ant Post limit : - - Anteriorly, it is in contact with the posterior aspect of Stations - 1–2 superiorly (Fig. 2A–C) and with Stations 4R and 4L inferiorly - (Fig. 2D–I and 3A–C). The anterior limit of Station 3P is kept - posterior to the trachea, which is defined by an imaginary - horizontal line running along the posterior wall of the trachea - (Fig. 2B,E, red line). Posteriorly, it is delineated along the - anterior and lateral borders of the vertebral body until an - imaginary horizontal line running 1 cm posteriorly to the - anterior border of the vertebral body (Fig. 2D). - - 2 - post to the trachea : define an imaginary line just post the - Trachea and remove what is anterior. - - */ - StartNewStep("[Station 3P] Post limits with VertebralBody "); - - // Load VertebralBody - MaskImagePointer VertebralBody = GetAFDB()->template GetImage("VertebralBody"); - - // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after) - VertebralBody = clitk::ResizeImageLike(VertebralBody, m_Working_Support, GetBackgroundValue()); - - writeImage(VertebralBody, "vb.mhd"); - - // Compute VertebralBody most Ant position (again because slices - // changes). Slice by slice, determine the most post point of the - // trachea (A) and consider a second point on the same horizontal - // line (B) - std::vector p_A; - std::vector p_B; - std::vector slices; - clitk::ExtractSlices(VertebralBody, 2, slices); - MaskImagePointType p; - typename MaskSliceType::PointType sp; - for(uint i=0; i(slices[i], GetBackgroundValue(), - 1, true, sp); - - // If the VertebralBody stop superiorly before the end of - // m_Working_Support, we consider the same previous point. - if (!found) { - p = p_A.back(); - p[2] += VertebralBody->GetSpacing()[2]; - p_A.push_back(p); - p = p_B.back(); - p[2] += VertebralBody->GetSpacing()[2]; - p_B.push_back(p); - } - else { - clitk::PointsUtils::Convert2DTo3D(sp, VertebralBody, i, p); - p[1] += 10; // Consider 10 mm more post - p_A.push_back(p); - p[0] -= 20; - p_B.push_back(p); - } - } - clitk::WriteListOfLandmarks(p_A, "vb-ant.txt"); - - - // Remove Ant part above those lines - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - p_A, p_B, - GetBackgroundValue(), - 1, -10); - writeImage(m_Working_Support, "a.mhd"); - - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; -} -//-------------------------------------------------------------------- - - //-------------------------------------------------------------------- template void @@ -424,19 +263,21 @@ clitk::ExtractLymphStationsFilter:: ExtractStation_3P_LR_sup_Limits_2() { /* - Use VertebralArtery to limit. - - - StartNewStep("[Station 3P] Left/Right limits with VertebralArtery"); - - // Load structures - MaskImagePointer VertebralArtery = GetAFDB()->template GetImage("VertebralArtery"); - - clitk::AndNot(m_Working_Support, VertebralArtery); + "On the right side, the limit is defined by the air–soft-tissue + interface. On the left side, it is defined by the air–tissue + interface superiorly (Fig. 2A–C) and the aorta inferiorly + (Figs. 2D–I and 3A–C)." - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; + sup : + Resizelike support : Trachea, SubclavianArtery + Trachea, slice by slice, get centroid trachea + SubclavianArtery, slice by slice, CCL + prendre la 1ère à L et R, not at Left + */ + // StartNewStep("[Station 3P] Left/Right limits (superior part) "); + + } //-------------------------------------------------------------------- @@ -460,6 +301,9 @@ ExtractStation_3P_LR_inf_Limits() MaskImagePointer AzygousVein = GetAFDB()->template GetImage("AzygousVein"); MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + DD("ici"); + writeImage(m_Working_Support, "ws.mhd"); + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); relPosFilter->VerboseStepFlagOff(); @@ -469,19 +313,20 @@ ExtractStation_3P_LR_inf_Limits() relPosFilter->SetInputObject(AzygousVein); relPosFilter->AddOrientationTypeString("R"); relPosFilter->SetInverseOrientationFlag(true); - // relPosFilter->SetIntermediateSpacing(3); + relPosFilter->SetIntermediateSpacing(3); relPosFilter->SetIntermediateSpacingFlag(false); relPosFilter->SetFuzzyThreshold(0.7); relPosFilter->AutoCropFlagOn(); relPosFilter->Update(); m_Working_Support = relPosFilter->GetOutput(); + DD("la"); writeImage(m_Working_Support, "before-L-aorta.mhd"); relPosFilter->SetInput(m_Working_Support); relPosFilter->SetInputObject(Aorta); relPosFilter->AddOrientationTypeString("L"); relPosFilter->SetInverseOrientationFlag(true); - // relPosFilter->SetIntermediateSpacing(3); + relPosFilter->SetIntermediateSpacing(3); relPosFilter->SetIntermediateSpacingFlag(false); relPosFilter->SetFuzzyThreshold(0.7); relPosFilter->AutoCropFlagOn(); diff --git a/segmentation/clitkExtractLymphStation_7.txx b/segmentation/clitkExtractLymphStation_7.txx index b12fa52..bcf2ccb 100644 --- a/segmentation/clitkExtractLymphStation_7.txx +++ b/segmentation/clitkExtractLymphStation_7.txx @@ -65,7 +65,7 @@ ExtractStation_7_SI_Limits() A, B, 2, 0, false); // Get the CarinaZ position - m_CarinaZ = FindCarinaSlicePosition(); + double m_CarinaZ = FindCarina(); // Crop support m_Working_Support = diff --git a/segmentation/clitkExtractLymphStation_Supports.txx b/segmentation/clitkExtractLymphStation_Supports.txx index 903948d..1008ff4 100644 --- a/segmentation/clitkExtractLymphStation_Supports.txx +++ b/segmentation/clitkExtractLymphStation_Supports.txx @@ -13,62 +13,536 @@ ExtractStationSupports() // Get initial Mediastinum m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage("Mediastinum", true); - // Superior limits: CricoidCartilag - // Inferior limits: lung - // (the Mediastinum support already stop at this limit) - - // Consider above Carina - m_CarinaZ = FindCarinaSlicePosition(); + // Consider sup/inf to Carina + double m_CarinaZ = FindCarina(); MaskImagePointer m_Support_Superior_to_Carina = clitk::CropImageRemoveLowerThan(m_Working_Support, 2, m_CarinaZ, true, GetBackgroundValue()); MaskImagePointer m_Support_Inferior_to_Carina = clitk::CropImageRemoveGreaterThan(m_Working_Support, 2, m_CarinaZ, true, GetBackgroundValue()); + m_ListOfSupports["Support_Superior_to_Carina"] = m_Support_Superior_to_Carina; + m_ListOfSupports["Support_Inferior_to_Carina"] = m_Support_Inferior_to_Carina; + writeImage(m_Support_Inferior_to_Carina, "seg/Support_Inf_Carina.mhd"); + GetAFDB()->SetImageFilename("Support_Inf_Carina", "seg/Support_Inf_Carina.mhd"); + writeImage(m_Support_Superior_to_Carina, "seg/Support_Sup_Carina.mhd"); + GetAFDB()->SetImageFilename("Support_Sup_Carina", "seg/Support_Sup_Carina.mhd"); + + // S1RL + Support_SupInf_S1RL(); + Support_LeftRight_S1R_S1L(); + + // S2RL + Support_SupInf_S2R_S2L(); + Support_LeftRight_S2R_S2L(); + + // S4RL + Support_SupInf_S4R_S4L(); + Support_LeftRight_S4R_S4L(); + + // Post limits of S1,S2,S4 + Support_Post_S1S2S4(); + + // S3AP + Support_S3P(); + Support_S3A(); + + // S5, S6 + Support_S5(); + Support_S6(); + + // Below Carina S7,8,9,10 + m_ListOfSupports["S7"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S8"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S9"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S10"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S11"] = clitk::Clone(m_Support_Inferior_to_Carina); + + // Store image filenames into AFDB + writeImage(m_ListOfSupports["S1R"], "seg/Support_S1R.mhd"); + GetAFDB()->SetImageFilename("Support_S1R", "seg/Support_S1R.mhd"); + writeImage(m_ListOfSupports["S1L"], "seg/Support_S1L.mhd"); + GetAFDB()->SetImageFilename("Support_S1L", "seg/Support_S1L.mhd"); + + writeImage(m_ListOfSupports["S2L"], "seg/Support_S2L.mhd"); + GetAFDB()->SetImageFilename("Support_S2L", "seg/Support_S2L.mhd"); + writeImage(m_ListOfSupports["S2R"], "seg/Support_S2R.mhd"); + GetAFDB()->SetImageFilename("Support_S2R", "seg/Support_S2R.mhd"); + + writeImage(m_ListOfSupports["S3P"], "seg/Support_S3P.mhd"); + GetAFDB()->SetImageFilename("Support_S3P", "seg/Support_S3P.mhd"); + writeImage(m_ListOfSupports["S3A"], "seg/Support_S3A.mhd"); + GetAFDB()->SetImageFilename("Support_S3A", "seg/Support_S3A.mhd"); + + writeImage(m_ListOfSupports["S4L"], "seg/Support_S4L.mhd"); + GetAFDB()->SetImageFilename("Support_S4L", "seg/Support_S4L.mhd"); + writeImage(m_ListOfSupports["S4R"], "seg/Support_S4R.mhd"); + GetAFDB()->SetImageFilename("Support_S4R", "seg/Support_S4R.mhd"); + + writeImage(m_ListOfSupports["S5"], "seg/Support_S5.mhd"); + GetAFDB()->SetImageFilename("Support_S5", "seg/Support_S5.mhd"); + writeImage(m_ListOfSupports["S6"], "seg/Support_S6.mhd"); + GetAFDB()->SetImageFilename("Support_S6", "seg/Support_S6.mhd"); + + writeImage(m_ListOfSupports["S7"], "seg/Support_S7.mhd"); + GetAFDB()->SetImageFilename("Support_S7", "seg/Support_S7.mhd"); + + writeImage(m_ListOfSupports["S8"], "seg/Support_S8.mhd"); + GetAFDB()->SetImageFilename("Support_S8", "seg/Support_S8.mhd"); + + writeImage(m_ListOfSupports["S9"], "seg/Support_S9.mhd"); + GetAFDB()->SetImageFilename("Support_S9", "seg/Support_S9.mhd"); + + writeImage(m_ListOfSupports["S10"], "seg/Support_S10.mhd"); + GetAFDB()->SetImageFilename("Support_S10", "seg/Support_S10.mhd"); + + writeImage(m_ListOfSupports["S11"], "seg/Support_S11.mhd"); + GetAFDB()->SetImageFilename("Support_S11", "seg/Support_S11.mhd"); + WriteAFDB(); +} +//-------------------------------------------------------------------- - // Consider only Superior to Carina - m_Working_Support = m_Support_Superior_to_Carina; +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_SupInf_S1RL() +{ // Step : S1RL - StartNewStep("[Support] sup-inf S1RL"); + StartNewStep("[Support] Sup-Inf S1RL"); /* - Lower border: clavicles bilaterally and, in the midline, the upper - border of the manubrium, 1R designates right-sided nodes, 1L, - left-sided nodes in this region - 2R: Upper border: apex of the right lung and pleural space, and in the midline, the upper border of the manubrium 2L: Upper border: apex of the left lung and pleural space, and in the midline, the upper border of the manubrium + + => apex / manubrium = up Sternum */ + m_Working_Support = m_ListOfSupports["Support_Superior_to_Carina"]; + MaskImagePointer Sternum = GetAFDB()->template GetImage ("Sternum"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); + MaskImagePointer S1RL = + clitk::CropImageRemoveLowerThan(m_Working_Support, 2, + p[2], true, GetBackgroundValue()); + m_Working_Support = + clitk::CropImageRemoveGreaterThan(m_Working_Support, 2, + p[2], true, GetBackgroundValue()); + m_ListOfSupports["S1RL"] = S1RL; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_LeftRight_S1R_S1L() +{ + // Step S1RL : Left-Right + StartNewStep("[Support] Left-Right S1R S1L"); + std::vector A; + std::vector B; + // Search for centroid positions of trachea + MaskImagePointer Trachea = GetAFDB()->template GetImage ("Trachea"); + MaskImagePointer S1RL = m_ListOfSupports["S1RL"]; + Trachea = clitk::ResizeImageLike(Trachea, S1RL, GetBackgroundValue()); + std::vector slices; + clitk::ExtractSlices(Trachea, 2, slices); + for(uint i=0; i(slices[i], 0, false, 10); + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), c); + ImagePointType a,b; + clitk::PointsUtils::Convert2DTo3D(c[1], Trachea, i, a); + A.push_back(a); + b = a; + b[1] += 50; + B.push_back(b); + } + clitk::WriteListOfLandmarks(A, "S1LR_A.txt"); + clitk::WriteListOfLandmarks(B, "S1LR_B.txt"); + // Clone support + MaskImagePointer S1R = clitk::Clone(S1RL); + MaskImagePointer S1L = clitk::Clone(S1RL); + // Right part + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1R, A, B, + GetBackgroundValue(), 0, 10); + S1R = clitk::AutoCrop(S1R, GetBackgroundValue()); + m_ListOfSupports["S1R"] = S1R; + // Left part + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1L, A, B, + GetBackgroundValue(), 0, -10); + S1L = clitk::AutoCrop(S1L, GetBackgroundValue()); + m_ListOfSupports["S1L"] = S1L; +} +//-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_SupInf_S2R_S2L() +{ + // Step : S2RL Sup-Inf limits + /* + 2R Lower border: intersection of caudal margin of innominate vein with + the trachea + 2L Lower border: superior border of the aortic arch + */ + StartNewStep("[Support] Sup-Inf S2RL"); + m_Working_Support = m_ListOfSupports["Support_Superior_to_Carina"]; + + // S2R Caudal Margin Of Left BrachiocephalicVein + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); + MaskImagePointType p; + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicVein, GetBackgroundValue(), 2, true, p); + // I add slightly more than a slice + double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2]+ 1.1*m_Working_Support->GetSpacing()[2]; + GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ); + MaskImagePointer S2R = + clitk::CropImageRemoveLowerThan(m_Working_Support, 2, + CaudalMarginOfLeftBrachiocephalicVeinZ, true, + GetBackgroundValue()); + m_ListOfSupports["S2R"] = S2R; + + // S2L : Top Of Aortic Arch + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + clitk::FindExtremaPointInAGivenDirection(Aorta, GetBackgroundValue(), 2, false, p); + // I add slightly more than a slice + double TopOfAorticArchZ=p[2]+ 1.1*m_Working_Support->GetSpacing()[2]; + GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ); + MaskImagePointer S2L = + clitk::CropImageRemoveLowerThan(m_Working_Support, 2, + TopOfAorticArchZ, true, + GetBackgroundValue()); + m_ListOfSupports["S2L"] = S2L; +} +//-------------------------------------------------------------------- - // // LeftRight cut along trachea - // MaskImagePointer Trachea = GetAFDB()->GetImage("Trachea"); - // // build a ant-post line for each slice - // MaskImagePointer m_Support_SupRight = - // clitk::CropImageRemoveLowerThan(m_Working_Support, 2, - // m_CarinaZ, true, GetBackgroundValue()); +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_LeftRight_S2R_S2L() +{ + // --------------------------------------------------------------------------- + /* Step : S2RL LeftRight + As for lymph node station 4R, 2R includes nodes extending to the + left lateral border of the trachea + Rod says: "For station 2 there is a shift, dividing 2R from 2L, from midline + to the left paratracheal border." + */ + StartNewStep("[Support] Separate 2R/2L according to Trachea"); + MaskImagePointer S2R = m_ListOfSupports["S2R"]; + MaskImagePointer S2L = m_ListOfSupports["S2L"]; + S2R = LimitsWithTrachea(S2R, 0, 1, -10); + S2L = LimitsWithTrachea(S2L, 0, 1, 10); + m_ListOfSupports["S2R"] = S2R; + m_ListOfSupports["S2L"] = S2L; + GetAFDB()->template ReleaseImage("Trachea"); +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_SupInf_S4R_S4L() +{ + // --------------------------------------------------------------------------- + /* Step : S4RL Sup-Inf + - start at the end of 2R and 2L + - stop ? + - 4R + Rod says : "The inferior border is at the lower border of the azygous vein." + Rod says : difficulties + (was : "ends at the upper lobe bronchus or where the right pulmonary artery + crosses the midline of the mediastinum ") + - 4L + Rod says : "The lower border is to upper margin of the left main pulmonary artery." + (was LLL bronchus) + */ + StartNewStep("[Support] Sup-Inf limits of 4R/4L"); + // Start from the support, remove 2R and 2L + MaskImagePointer S4RL = clitk::Clone(m_Working_Support); + MaskImagePointer S2R = m_ListOfSupports["S2R"]; + MaskImagePointer S2L = m_ListOfSupports["S2L"]; + clitk::AndNot(S4RL, S2R, GetBackgroundValue()); + clitk::AndNot(S4RL, S2L, GetBackgroundValue()); + S4RL = clitk::AutoCrop(S4RL, GetBackgroundValue()); + // Copy, stop 4R at AzygousVein and 4L at LeftPulmonaryArtery + MaskImagePointer S4R = clitk::Clone(S4RL); + MaskImagePointer S4L = clitk::Clone(S4RL); + + // Get AzygousVein and limit according to LowerBorderAzygousVein + MaskImagePointer LowerBorderAzygousVein + = GetAFDB()->template GetImage("LowerBorderAzygousVein"); + std::vector c; + clitk::ComputeCentroids(LowerBorderAzygousVein, GetBackgroundValue(), c); + S4R = clitk::CropImageRemoveLowerThan(S4RL, 2, + c[1][2], true, GetBackgroundValue()); + S4R = clitk::AutoCrop(S4R, GetBackgroundValue()); + m_ListOfSupports["S4R"] = S4R; - // Store image filenames into AFDB - m_ListOfSupports["S1R"] = m_Working_Support; - writeImage(m_ListOfSupports["S1R"], "seg/Support_S1R.mhd"); - GetAFDB()->SetImageFilename("Support_S1R", "seg/Support_S1R.mhd"); - WriteAFDB(); + // Limit according to LeftPulmonaryArtery + MaskImagePointer LeftPulmonaryArtery + = GetAFDB()->template GetImage("LeftPulmonaryArtery"); + MaskImagePointType p; + clitk::FindExtremaPointInAGivenDirection(LeftPulmonaryArtery, GetBackgroundValue(), + 2, false, p); + S4L = clitk::CropImageRemoveLowerThan(S4RL, 2, + p[2], true, GetBackgroundValue()); + S4L = clitk::AutoCrop(S4L, GetBackgroundValue()); + m_ListOfSupports["S4L"] = S4L; } //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_LeftRight_S4R_S4L() +{ + // --------------------------------------------------------------------------- + /* Step : S4RL LeftRight + + - 4R: includes right paratracheal nodes, and pretracheal nodes + extending to the left lateral border of trachea + + - 4L: includes nodes to the left of the left lateral border of + the trachea, medial to the ligamentum arteriosum + + => same than 2RL + */ + StartNewStep("[Support] Left Right separation of 4R/4L"); + + MaskImagePointer S4R = m_ListOfSupports["S4R"]; + MaskImagePointer S4L = m_ListOfSupports["S4L"]; + S4R = LimitsWithTrachea(S4R, 0, 1, -10); + S4L = LimitsWithTrachea(S4L, 0, 1, 10); + m_ListOfSupports["S4R"] = S4R; + m_ListOfSupports["S4L"] = S4L; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection, + double offset) +{ + MaskImagePointType min, max; + GetMinMaxBoundary(input, min, max); + return LimitsWithTrachea(input, extremaDirection, lineDirection, offset, max[2]); +} +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection, + double offset, double maxSupPosition) +{ + /* + Take the input mask, consider the trachea and limit according to + Left border of the trachea. Keep at Left or at Right according to + the offset + */ + // Read the trachea + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + + // Find extrema post positions + std::vector tracheaLeftPositionsA; + std::vector tracheaLeftPositionsB; + + // Crop Trachea only on the Sup-Inf axes, without autocrop + // Trachea = clitk::ResizeImageLike(Trachea, input, GetBackgroundValue()); + MaskImagePointType min, max; + GetMinMaxBoundary(input, min, max); + Trachea = clitk::CropImageAlongOneAxis(Trachea, 2, min[2], max[2], + false, GetBackgroundValue()); + + // Select the main CCL (because of bronchus) + Trachea = SliceBySliceKeepMainCCL(Trachea, GetBackgroundValue(), GetForegroundValue()); + + // Slice by slice, build the separation line + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(Trachea, + GetBackgroundValue(), 2, + extremaDirection, false, // Left + lineDirection, // Vertical line + 1, // margins + tracheaLeftPositionsA, + tracheaLeftPositionsB); + // Do not consider trachea above the limit + int indexMax=tracheaLeftPositionsA.size(); + for(uint i=0; i maxSupPosition) { + indexMax = i; + i = tracheaLeftPositionsA.size(); // stop loop + } + } + tracheaLeftPositionsA.erase(tracheaLeftPositionsA.begin()+indexMax, tracheaLeftPositionsA.end()); + tracheaLeftPositionsB.erase(tracheaLeftPositionsB.begin()+indexMax, tracheaLeftPositionsB.end()); + + // Cut post to this line + clitk::SliceBySliceSetBackgroundFromLineSeparation(input, + tracheaLeftPositionsA, + tracheaLeftPositionsB, + GetBackgroundValue(), + extremaDirection, offset); + MaskImagePointer output = clitk::AutoCrop(input, GetBackgroundValue()); + return output; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_Post_S1S2S4() +{ + StartNewStep("[Support] Post limits of S1RL, S2RL, S4RL"); + + double m_ApexOfTheChest = FindApexOfTheChest(); + + // Post limits with Trachea + MaskImagePointer S1R = m_ListOfSupports["S1R"]; + MaskImagePointer S1L = m_ListOfSupports["S1L"]; + MaskImagePointer S2R = m_ListOfSupports["S2R"]; + MaskImagePointer S2L = m_ListOfSupports["S2L"]; + MaskImagePointer S4R = m_ListOfSupports["S4R"]; + MaskImagePointer S4L = m_ListOfSupports["S4L"]; + S1L = LimitsWithTrachea(S1L, 1, 0, -10, m_ApexOfTheChest); + S1R = LimitsWithTrachea(S1R, 1, 0, -10, m_ApexOfTheChest); + S2R = LimitsWithTrachea(S2R, 1, 0, -10, m_ApexOfTheChest); + S2L = LimitsWithTrachea(S2L, 1, 0, -10, m_ApexOfTheChest); + S4R = LimitsWithTrachea(S4R, 1, 0, -10, m_ApexOfTheChest); + S4L = LimitsWithTrachea(S4L, 1, 0, -10, m_ApexOfTheChest); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S3P() +{ + StartNewStep("[Support] Ant limits of S3P and Post limits of S1RL, S2RL, S4RL"); + + // Initial S3P support + MaskImagePointer S3P = clitk::Clone(m_ListOfSupports["Support_Superior_to_Carina"]); + + // Stop at Lung Apex + double m_ApexOfTheChest = FindApexOfTheChest(); + S3P = + clitk::CropImageRemoveGreaterThan(S3P, 2, + m_ApexOfTheChest, true, + GetBackgroundValue()); + // Ant limits with Trachea + S3P = LimitsWithTrachea(S3P, 1, 0, 10); + m_ListOfSupports["S3P"] = S3P; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S3A() +{ + StartNewStep("[Support] Sup-Inf and Post limits for S3A"); + + // Initial S3A support + MaskImagePointer S3A = clitk::Clone(m_ListOfSupports["Support_Superior_to_Carina"]); + + // Stop at Lung Apex or like S2/S1 (upper border Sternum - manubrium) ? + + //double m_ApexOfTheChest = FindApexOfTheChest(); + + MaskImagePointer Sternum = GetAFDB()->template GetImage ("Sternum"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); + + S3A = + clitk::CropImageRemoveGreaterThan(S3A, 2, + //m_ApexOfTheChest + p[2], true, + GetBackgroundValue()); + // Ant limits with Trachea + S3A = LimitsWithTrachea(S3A, 1, 0, -10); + m_ListOfSupports["S3A"] = S3A; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S5() +{ + StartNewStep("[Support] Sup-Inf limits S5 with aorta"); + + // Initial S5 support + MaskImagePointer S5 = clitk::Clone(GetAFDB()->template GetImage("Mediastinum", true)); + + // Sup limits with Aorta + double sup = FindInferiorBorderOfAorticArch(); + + // Inf limits with "upper rim of the left main pulmonary artery" + // For the moment only, it will change. + MaskImagePointer PulmonaryTrunk = GetAFDB()->template GetImage("PulmonaryTrunk"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(PulmonaryTrunk, GetBackgroundValue(), 2, false, p); + + // Cut Sup/Inf + S5 = clitk::CropImageAlongOneAxis(S5, 2, p[2], sup, true, GetBackgroundValue()); + + m_ListOfSupports["S5"] = S5; +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S6() +{ + StartNewStep("[Support] Sup-Inf limits S6 with aorta"); + + // Initial S6 support like S3A + MaskImagePointer S6 = clitk::Clone(m_ListOfSupports["S3A"]); + + // Inf Sup limits with Aorta + double sup = FindSuperiorBorderOfAorticArch(); + double inf = FindInferiorBorderOfAorticArch(); + + // Cut Sup/Inf + S6 = clitk::CropImageAlongOneAxis(S6, 2, inf, sup, true, GetBackgroundValue()); + + m_ListOfSupports["S6"] = S6; +} +//-------------------------------------------------------------------- + diff --git a/segmentation/clitkExtractLymphStations.ggo b/segmentation/clitkExtractLymphStations.ggo index 96ac994..3cce000 100644 --- a/segmentation/clitkExtractLymphStations.ggo +++ b/segmentation/clitkExtractLymphStations.ggo @@ -18,6 +18,14 @@ option "afdb" a "Input Anatomical Feature DB" string no option "input" i "Input filename" string no option "output" o "Output lungs mask filename" string no option "station" - "Force to compute station even if already exist in the DB" string no multiple +option "nosupport" - "Do not recompute the station supports if already available" flag off + + +section "Options for Station 3A" +option "S3A_ft_SVC" - "Fuzzy Threshold for NotPost to SVC" double default="0.5" no +option "S3A_ft_Aorta" - "Fuzzy Threshold for NotPost to Aorta" double default="0.5" no +option "S3A_ft_SubclavianArtery" - "Fuzzy Threshold for Ant to SubclavianArtery" double default="0.2" no + section "Options for Station 8" option "S8_esophagusDilatationForAnt" - "Dilatation of esophagus, in mm, for 'anterior' limits (default=15,2,1)" double no multiple @@ -33,10 +41,6 @@ option "S7_ft_LeftPulmonaryArtery" - "Fuzzy Threshold for LeftPulmonaryAr option "S7_ft_SVC" - "Fuzzy Threshold for SVC" double default="0.2" no option "S7_UseMostInferiorPartOnly" - "Inferior separation with S8." flag off -section "Options for Station 3A" -option "S3A_ft_Sternum" - "Fuzzy Threshold for Post to Sternum" double default="0.5" no -option "S3A_ft_SubclavianArtery" - "Fuzzy Threshold for Ant to SubclavianArtery" double default="0.2" no - section "Options for Station 2RL" option "S2RL_ft_CommonCarotidArtery" - "Threshold for NotAntTo to CommonCarotidArtery" double default="0.7" no option "S2RL_ft_BrachioCephalicTrunk" - "Threshold for NotAntTo to BrachioCephalicTrunk" double default="0.7" no diff --git a/segmentation/clitkExtractLymphStationsFilter.h b/segmentation/clitkExtractLymphStationsFilter.h index dd99c5f..8bbc301 100644 --- a/segmentation/clitkExtractLymphStationsFilter.h +++ b/segmentation/clitkExtractLymphStationsFilter.h @@ -109,6 +109,9 @@ namespace clitk { void AddComputeStation(std::string station) ; void SetFuzzyThreshold(std::string station, std::string tag, double value); double GetFuzzyThreshold(std::string station, std::string tag); + itkGetConstMacro(ComputeStationsSupportsFlag, bool); + itkSetMacro(ComputeStationsSupportsFlag, bool); + itkBooleanMacro(ComputeStationsSupportsFlag); protected: ExtractLymphStationsFilter(); @@ -131,9 +134,14 @@ namespace clitk { void Remove_Structures(std::string station, std::string s); // Functions common to several stations - void FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B); - double FindCarinaSlicePosition(); + double FindCarina(); + double FindApexOfTheChest(); + double FindSuperiorBorderOfAorticArch(); + double FindInferiorBorderOfAorticArch(); void FindLeftAndRightBronchi(); + void FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B); + MaskImagePointer FindAntPostVessels(); + MaskImagePointer FindAntPostVessels2(); // Global parameters typedef std::map FuzzyThresholdByStructureType; @@ -141,11 +149,28 @@ namespace clitk { // Station's supports void ExtractStationSupports(); + void Support_SupInf_S1RL(); + void Support_LeftRight_S1R_S1L(); + void Support_SupInf_S2R_S2L(); + void Support_LeftRight_S2R_S2L(); + void Support_SupInf_S4R_S4L(); + void Support_LeftRight_S4R_S4L(); + void Support_Post_S1S2S4(); + void Support_S3P(); + void Support_S3A(); + void Support_S5(); + void Support_S6(); + + MaskImagePointer LimitsWithTrachea(MaskImageType * input, + int extremaDirection, int lineDirection, + double offset, double maxSupPosition); + MaskImagePointer LimitsWithTrachea(MaskImageType * input, + int extremaDirection, int lineDirection, + double offset); // Station 8 // double m_DistanceMaxToAnteriorPartOfTheSpine; double m_DiaphragmInferiorLimit; - double m_CarinaZ; double m_OriginOfRightMiddleLobeBronchusZ; double m_InjectedThresholdForS8; MaskImagePointer m_Esophagus; @@ -164,13 +189,20 @@ namespace clitk { // Station 3P void ExtractStation_3P(); void ExtractStation_3P_SetDefaultValues(); - void ExtractStation_3P_SI_Limits(); + void ExtractStation_3P_LR_inf_Limits(); + void ExtractStation_3P_LR_sup_Limits_2(); void ExtractStation_3P_Remove_Structures(); - void ExtractStation_3P_Ant_Limits(); - void ExtractStation_3P_Post_Limits(); void ExtractStation_3P_LR_sup_Limits(); - void ExtractStation_3P_LR_sup_Limits_2(); - void ExtractStation_3P_LR_inf_Limits(); + + // Station 3A + void ExtractStation_3A(); + void ExtractStation_3A_SetDefaultValues(); + void ExtractStation_3A_AntPost_S5(); + void ExtractStation_3A_AntPost_S6(); + void ExtractStation_3A_AntPost_Superiorly(); + void ExtractStation_3A_SI_Limits(); + void ExtractStation_3A_Ant_Limits(); + void ExtractStation_3A_Post_Limits(); // Station 2RL void ExtractStation_2RL(); @@ -184,13 +216,6 @@ namespace clitk { void ExtractStation_2RL_SeparateRL(); vtkSmartPointer Build3DMeshFrom2DContour(const std::vector & points); - // Station 3A - void ExtractStation_3A(); - void ExtractStation_3A_SetDefaultValues(); - void ExtractStation_3A_SI_Limits(); - void ExtractStation_3A_Ant_Limits(); - void ExtractStation_3A_Post_Limits(); - // Station 7 void ExtractStation_7(); void ExtractStation_7_SetDefaultValues(); @@ -202,6 +227,7 @@ namespace clitk { void ExtractStation_7_Posterior_Limits(); void ExtractStation_7_Remove_Structures(); bool m_S7_UseMostInferiorPartOnlyFlag; + bool m_ComputeStationsSupportsFlag; MaskImagePointer m_Working_Trachea; MaskImagePointer m_LeftBronchus; MaskImagePointer m_RightBronchus; diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index 58f2991..904bfc6 100644 --- a/segmentation/clitkExtractLymphStationsFilter.txx +++ b/segmentation/clitkExtractLymphStationsFilter.txx @@ -27,6 +27,7 @@ #include "clitkAutoCropFilter.h" #include "clitkSegmentationUtils.h" #include "clitkSliceBySliceRelativePositionFilter.h" +#include "clitkMeshToBinaryImageFilter.h" // itk #include @@ -40,6 +41,11 @@ // itk ENST #include "RelativePositionPropImageFilter.h" +// vtk +#include +#include +#include + //-------------------------------------------------------------------- template clitk::ExtractLymphStationsFilter:: @@ -51,6 +57,7 @@ ExtractLymphStationsFilter(): this->SetNumberOfRequiredInputs(1); SetBackgroundValue(0); SetForegroundValue(1); + ComputeStationsSupportsFlagOn(); // Default values ExtractStation_8_SetDefaultValues(); @@ -72,11 +79,43 @@ GenerateOutputInformation() { m_Input = dynamic_cast(itk::ProcessObject::GetInput(0)); m_Mediastinum = GetAFDB()->template GetImage ("Mediastinum"); + // Clean some computer landmarks to force the recomputation + GetAFDB()->RemoveTag("AntPostVesselsSeparation"); + // Global supports for stations - StartNewStep("Supports for stations"); - StartSubStep(); - ExtractStationSupports(); - StopSubStep(); + if (GetComputeStationsSupportsFlag()) { + StartNewStep("Supports for stations"); + StartSubStep(); + GetAFDB()->RemoveTag("CarinaZ"); + GetAFDB()->RemoveTag("ApexOfTheChestZ"); + GetAFDB()->RemoveTag("ApexOfTheChest"); + GetAFDB()->RemoveTag("RightBronchus"); + GetAFDB()->RemoveTag("LeftBronchus"); + GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ"); + GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch"); + GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ"); + GetAFDB()->RemoveTag("InferiorBorderOfAorticArch"); + ExtractStationSupports(); + StopSubStep(); + } + else { + m_ListOfSupports["S1R"] = GetAFDB()->template GetImage("Support_S1R"); + m_ListOfSupports["S1L"] = GetAFDB()->template GetImage("Support_S1L"); + m_ListOfSupports["S2R"] = GetAFDB()->template GetImage("Support_S2R"); + m_ListOfSupports["S2L"] = GetAFDB()->template GetImage("Support_S2L"); + m_ListOfSupports["S4R"] = GetAFDB()->template GetImage("Support_S4R"); + m_ListOfSupports["S4L"] = GetAFDB()->template GetImage("Support_S4L"); + + m_ListOfSupports["S3A"] = GetAFDB()->template GetImage("Support_S3A"); + m_ListOfSupports["S3P"] = GetAFDB()->template GetImage("Support_S3P"); + m_ListOfSupports["S5"] = GetAFDB()->template GetImage("Support_S5"); + m_ListOfSupports["S6"] = GetAFDB()->template GetImage("Support_S6"); + m_ListOfSupports["S7"] = GetAFDB()->template GetImage("Support_S7"); + m_ListOfSupports["S8"] = GetAFDB()->template GetImage("Support_S8"); + m_ListOfSupports["S9"] = GetAFDB()->template GetImage("Support_S9"); + m_ListOfSupports["S10"] = GetAFDB()->template GetImage("Support_S10"); + m_ListOfSupports["S11"] = GetAFDB()->template GetImage("Support_S11"); + } // Extract Station8 StartNewStep("Station 8"); @@ -90,20 +129,21 @@ GenerateOutputInformation() { ExtractStation_3P(); StopSubStep(); - // Extract Station2RL - /* - StartNewStep("Station 2RL"); - StartSubStep(); - ExtractStation_2RL(); - StopSubStep(); - */ - // Extract Station3A StartNewStep("Station 3A"); StartSubStep(); ExtractStation_3A(); StopSubStep(); + + // HERE + + // Extract Station2RL + StartNewStep("Station 2RL"); + StartSubStep(); + ExtractStation_2RL(); + StopSubStep(); + // Extract Station7 StartNewStep("Station 7"); StartSubStep(); @@ -197,70 +237,79 @@ AddComputeStation(std::string station) //-------------------------------------------------------------------- + //-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P() -{ - if (CheckForStation("3P")) { - ExtractStation_3P_SI_Limits(); - ExtractStation_3P_Ant_Limits(); - ExtractStation_3P_Post_Limits(); - ExtractStation_3P_LR_sup_Limits(); - // ExtractStation_3P_LR_sup_Limits_2(); - ExtractStation_3P_LR_inf_Limits(); - ExtractStation_8_Single_CCL_Limits(); // YES 8 ! - ExtractStation_3P_Remove_Structures(); // after CCL - // Store image filenames into AFDB - writeImage(m_ListOfStations["3P"], "seg/Station3P.mhd"); - GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); - WriteAFDB(); +template +class comparePointsX { +public: + bool operator() (PointType i, PointType j) { return (i[0] +class comparePointsWithAngle { +public: + bool operator() (PairType i, PairType j) { return (i.second < j.second); } +}; +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void HypercubeCorners(std::vector > & out) { + std::vector > previous; + HypercubeCorners(previous); + out.resize(previous.size()*2); + for(unsigned int i=0; i p; + if (i -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3A() -{ - if (CheckForStation("3A")) { - ExtractStation_3A_SI_Limits(); - ExtractStation_3A_Ant_Limits(); - ExtractStation_3A_Post_Limits(); - // Store image filenames into AFDB - writeImage(m_ListOfStations["3A"], "seg/Station3A.mhd"); - GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); - WriteAFDB(); - } +template<> +void HypercubeCorners<1>(std::vector > & out) { + out.resize(2); + out[0][0] = 0; + out[1][0] = 1; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL() +template +void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image, + std::vector & bounds) { - if (CheckForStation("2RL")) { - ExtractStation_2RL_SI_Limits(); - ExtractStation_2RL_Post_Limits(); - ExtractStation_2RL_Ant_Limits2(); - ExtractStation_2RL_Ant_Limits(); - ExtractStation_2RL_LR_Limits(); - ExtractStation_2RL_Remove_Structures(); - ExtractStation_2RL_SeparateRL(); - - // Store image filenames into AFDB - writeImage(m_ListOfStations["2R"], "seg/Station2R.mhd"); - writeImage(m_ListOfStations["2L"], "seg/Station2L.mhd"); - GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd"); - GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd"); - WriteAFDB(); + // Get image max/min coordinates + const unsigned int dim=ImageType::ImageDimension; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::IndexType IndexType; + PointType min_c, max_c; + IndexType min_i, max_i; + min_i = image->GetLargestPossibleRegion().GetIndex(); + for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; + image->TransformIndexToPhysicalPoint(min_i, min_c); + image->TransformIndexToPhysicalPoint(max_i, max_c); + + // Get corners coordinates + HypercubeCorners(bounds); + for(unsigned int i=0; i double clitk::ExtractLymphStationsFilter:: -FindCarinaSlicePosition() +FindCarina() { double z; try { @@ -395,7 +444,7 @@ FindCarinaSlicePosition() clitk::ComputeCentroids(Carina, GetBackgroundValue(), centroids); // We dont need Carina structure from now - Carina->Delete(); + GetAFDB()->template ReleaseImage("Carina"); // Put inside the AFDB GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]); @@ -408,6 +457,38 @@ FindCarinaSlicePosition() //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindApexOfTheChest() +{ + double z; + try { + z = GetAFDB()->GetDouble("ApexOfTheChestZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindApexOfTheChestPosition"); + MaskImagePointer Lungs = GetAFDB()->template GetImage("Lungs"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Lungs, GetBackgroundValue(), 2, false, p); + + // We dont need Lungs structure from now + GetAFDB()->template ReleaseImage("Lungs"); + + // Put inside the AFDB + GetAFDB()->SetPoint3D("ApexOfTheChest", p); + p[2] -= 5; // We consider 5 mm lower + GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]); + WriteAFDB(); + z = p[2]; + } + return z; +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -428,7 +509,7 @@ FindLeftAndRightBronchi() MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); // Get the Carina position - m_CarinaZ = FindCarinaSlicePosition(); + double m_CarinaZ = FindCarina(); // Consider only inferiorly to the Carina MaskImagePointer m_Working_Trachea = @@ -500,4 +581,819 @@ FindLeftAndRightBronchi() } //-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindSuperiorBorderOfAorticArch() +{ + double z; + try { + z = GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindSuperiorBorderOfAorticArch"); + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Aorta, GetBackgroundValue(), 2, false, p); + p[2] += Aorta->GetSpacing()[2]; // the slice above + + // We dont need Lungs structure from now + GetAFDB()->template ReleaseImage("Aorta"); + + // Put inside the AFDB + GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p); + GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]); + WriteAFDB(); + z = p[2]; + } + return z; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindInferiorBorderOfAorticArch() +{ + double z; + try { + z = GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindInferiorBorderOfAorticArch"); + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + std::vector slices; + clitk::ExtractSlices(Aorta, 2, slices); + bool found=false; + uint i = slices.size()-1; + while (!found) { + slices[i] = Labelize(slices[i], 0, false, 10); + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), c); + if (c.size()>2) { + found = true; + } + else { + i--; + } + } + MaskImageIndexType index; + index[0] = index[1] = 0.0; + index[2] = i+1; + MaskImagePointType lower; + Aorta->TransformIndexToPhysicalPoint(index, lower); + + // We dont need Lungs structure from now + GetAFDB()->template ReleaseImage("Aorta"); + + // Put inside the AFDB + GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower); + GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]); + WriteAFDB(); + z = lower[2]; + } + return z; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +FindAntPostVessels() +{ + // ----------------------------------------------------- + /* Rod says: "The anterior border, as with the Atlas – UM, is + posterior to the vessels (right subclavian vein, left + brachiocephalic vein, right brachiocephalic vein, left subclavian + artery, left common carotid artery and brachiocephalic trunk). + These vessels are not included in the nodal station. The anterior + border is drawn to the midpoint of the vessel and an imaginary + line joins the middle of these vessels. Between the vessels, + station 2 is in contact with station 3a." */ + + // Check if no already done + bool found = true; + MaskImagePointer binarizedContour; + try { + DD("FindAntPostVessels try to get"); + binarizedContour = GetAFDB()->template GetImage ("AntPostVesselsSeparation"); + } + catch(clitk::ExceptionObject e) { + DD("not found"); + found = false; + } + if (found) { + return binarizedContour; + } + + /* Here, we consider the vessels form a kind of anterior barrier. We + link all vessels centroids and remove what is post to it. + - select the list of structure + vessel1 = BrachioCephalicArtery + vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) + vessel3 = CommonCarotidArtery + vessel4 = SubclavianArtery + other = Thyroid + other = Aorta + - crop images as needed + - slice by slice, choose the CCL and extract centroids + - slice by slice, sort according to polar angle wrt Trachea centroid. + - slice by slice, link centroids and close contour + - remove outside this contour + - merge with support + */ + + // Read structures + MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage("BrachioCephalicArtery"); + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); + MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); + MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); + MaskImagePointer Thyroid = GetAFDB()->template GetImage("Thyroid"); + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + + // Create a temporay support + // From first slice of BrachioCephalicVein to end of 3A + MaskImagePointer support = GetAFDB()->template GetImage("Support_Sup_Carina"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicVein, GetBackgroundValue(), 2, true, p); + double inf = p [2]; + clitk::FindExtremaPointInAGivenDirection(GetAFDB()->template GetImage("Support_S3A"), + GetBackgroundValue(), 2, false, p); + double sup = p [2]; + support = clitk::CropImageAlongOneAxis(support, 2, inf, sup, + false, GetBackgroundValue()); + + // Resize all structures like support + BrachioCephalicArtery = + clitk::ResizeImageLike(BrachioCephalicArtery, support, GetBackgroundValue()); + CommonCarotidArtery = + clitk::ResizeImageLike(CommonCarotidArtery, support, GetBackgroundValue()); + SubclavianArtery = + clitk::ResizeImageLike(SubclavianArtery, support, GetBackgroundValue()); + Thyroid = + clitk::ResizeImageLike(Thyroid, support, GetBackgroundValue()); + Aorta = + clitk::ResizeImageLike(Aorta, support, GetBackgroundValue()); + BrachioCephalicVein = + clitk::ResizeImageLike(BrachioCephalicVein, support, GetBackgroundValue()); + Trachea = + clitk::ResizeImageLike(Trachea, support, GetBackgroundValue()); + + // Extract slices + std::vector slices_BrachioCephalicArtery; + clitk::ExtractSlices(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery); + std::vector slices_BrachioCephalicVein; + clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BrachioCephalicVein); + std::vector slices_CommonCarotidArtery; + clitk::ExtractSlices(CommonCarotidArtery, 2, slices_CommonCarotidArtery); + std::vector slices_SubclavianArtery; + clitk::ExtractSlices(SubclavianArtery, 2, slices_SubclavianArtery); + std::vector slices_Thyroid; + clitk::ExtractSlices(Thyroid, 2, slices_Thyroid); + std::vector slices_Aorta; + clitk::ExtractSlices(Aorta, 2, slices_Aorta); + std::vector slices_Trachea; + clitk::ExtractSlices(Trachea, 2, slices_Trachea); + unsigned int n= slices_BrachioCephalicArtery.size(); + + // Get the boundaries of one slice + std::vector bounds; + ComputeImageBoundariesCoordinates(slices_BrachioCephalicArtery[0], bounds); + + // For all slices, for all structures, find the centroid and build the contour + // List of 3D points (for debug) + std::vector p3D; + + vtkSmartPointer append = vtkSmartPointer::New(); + for(unsigned int i=0; i(slices_CommonCarotidArtery[i], + GetBackgroundValue(), true, 1); + slices_SubclavianArtery[i] = Labelize(slices_SubclavianArtery[i], + GetBackgroundValue(), true, 1); + slices_BrachioCephalicArtery[i] = Labelize(slices_BrachioCephalicArtery[i], + GetBackgroundValue(), true, 1); + slices_BrachioCephalicVein[i] = Labelize(slices_BrachioCephalicVein[i], + GetBackgroundValue(), true, 1); + slices_Thyroid[i] = Labelize(slices_Thyroid[i], + GetBackgroundValue(), true, 1); + slices_Aorta[i] = Labelize(slices_Aorta[i], + GetBackgroundValue(), true, 1); + + // Search centroids + std::vector points2D; + std::vector centroids1; + std::vector centroids2; + std::vector centroids3; + std::vector centroids4; + std::vector centroids5; + std::vector centroids6; + ComputeCentroids(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1); + ComputeCentroids(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2); + ComputeCentroids(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3); + ComputeCentroids(slices_Thyroid[i], GetBackgroundValue(), centroids4); + ComputeCentroids(slices_Aorta[i], GetBackgroundValue(), centroids5); + ComputeCentroids(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6); + + // BrachioCephalicVein -> when it is separated into two CCL, we + // only consider the most at Right one + if (centroids6.size() > 2) { + if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2); + else centroids6.erase(centroids6.begin()+1); + } + + // BrachioCephalicVein -> when SubclavianArtery has 2 CCL + // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein + if ((centroids3.size() ==1) && (centroids2.size() > 2)) { + centroids6.clear(); + } + + for(unsigned int j=1; j centroids_trachea; + ComputeCentroids(slices_Trachea[i], GetBackgroundValue(), centroids_trachea); + typedef std::pair PointAngleType; + std::vector angles; + for(unsigned int j=0; j0) angle = atan(y/x); + if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; + if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; + if (x==0) { + if (y>0) angle = M_PI/2.0; + if (y<0) angle = -M_PI/2.0; + if (y==0) angle = 0; + } + angle = clitk::rad2deg(angle); + // Angle is [-180;180] wrt the X axis. We change the X axis to + // be the vertical line, because we want to sort from Right to + // Left from Post to Ant. + if (angle>0) + angle = (270-angle); + if (angle<0) { + angle = -angle-90; + if (angle<0) angle = 360-angle; + } + PointAngleType a; + a.first = points2D[j]; + a.second = angle; + angles.push_back(a); + } + + // Do nothing if less than 2 points --> n + if (points2D.size() < 3) { //continue; + continue; + } + + // Sort points2D according to polar angles + std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); + for(unsigned int j=0; j toadd; + unsigned int index = 0; + double dmax = 5; + while (indexdmax) { + + MaskSlicePointType b; + b[0] = a[0]+(c[0]-a[0])/2.0; + b[1] = a[1]+(c[1]-a[1])/2.0; + + // Compute distance to trachea centroid + MaskSlicePointType m = centroids_trachea[1]; + double da = m.EuclideanDistanceTo(a); + double db = m.EuclideanDistanceTo(b); + //double dc = m.EuclideanDistanceTo(c); + + // Mean distance, find point on the line from trachea centroid + // to b + double alpha = (da+db)/2.0; + MaskSlicePointType v; + double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); + v[0] = (b[0]-m[0])/n; + v[1] = (b[1]-m[1])/n; + MaskSlicePointType r; + r[0] = m[0]+alpha*v[0]; + r[1] = m[1]+alpha*v[1]; + points2D.insert(points2D.begin()+index+1, r); + } + else { + index++; + } + } + // DDV(points2D, points2D.size()); + + // Add some points to close the contour + // - H line towards Right + MaskSlicePointType p = points2D[0]; + p[0] = bounds[0][0]; + points2D.insert(points2D.begin(), p); + // - corner Right/Post + p = bounds[0]; + points2D.insert(points2D.begin(), p); + // - H line towards Left + p = points2D.back(); + p[0] = bounds[2][0]; + points2D.push_back(p); + // - corner Right/Post + p = bounds[2]; + points2D.push_back(p); + // Close contour with the first point + points2D.push_back(points2D[0]); + // DDV(points2D, points2D.size()); + + // Build 3D points from the 2D points + std::vector points3D; + clitk::PointsUtils::Convert2DListTo3DList(points2D, i, support, points3D); + for(unsigned int x=0; x mesh = Build3DMeshFrom2DContour(points3D); + append->AddInput(mesh); + } + + // DEBUG: write points3D + clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); + + // Build the final 3D mesh form the list 2D mesh + append->Update(); + vtkSmartPointer mesh = append->GetOutput(); + + // Debug, write the mesh + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(mesh); + w->SetFileName("bidon.vtk"); + w->Write(); + */ + + // Compute a single binary 3D image from the list of contours + clitk::MeshToBinaryImageFilter::Pointer filter = + clitk::MeshToBinaryImageFilter::New(); + filter->SetMesh(mesh); + filter->SetLikeImage(support); + filter->Update(); + binarizedContour = filter->GetOutput(); + + // Crop + clitk::FindExtremaPointInAGivenDirection(binarizedContour, GetForegroundValue(), 2, true, p); + inf = p[2]; + DD(p); + clitk::FindExtremaPointInAGivenDirection(binarizedContour, GetForegroundValue(), 2, false, p); + sup = p[2]; + DD(p); + binarizedContour = clitk::CropImageAlongOneAxis(binarizedContour, 2, inf, sup, + false, GetBackgroundValue()); + // Update the AFDB + writeImage(binarizedContour, "seg/AntPostVesselsSeparation.mhd"); + GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd"); + WriteAFDB(); + return binarizedContour; + + /* + // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse + ImagePointType p = p3D[2]; // This is the first centroid of the first slice + p[1] += 50; // 50 mm Post from this point must be kept + ImageIndexType index; + binarizedContour->TransformPhysicalPointToIndex(p, index); + bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); + + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + if (isInside) + boolFilter->SetOperationType(BoolFilterType::And); + else + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + */ + + // End + //StopCurrentStep(m_Working_Support); + //m_ListOfStations["2R"] = m_Working_Support; + //m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + + + + + + + + + + + + + + + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +FindAntPostVessels2() +{ + // ----------------------------------------------------- + /* Rod says: "The anterior border, as with the Atlas – UM, is + posterior to the vessels (right subclavian vein, left + brachiocephalic vein, right brachiocephalic vein, left subclavian + artery, left common carotid artery and brachiocephalic trunk). + These vessels are not included in the nodal station. The anterior + border is drawn to the midpoint of the vessel and an imaginary + line joins the middle of these vessels. Between the vessels, + station 2 is in contact with station 3a." */ + + // Check if no already done + bool found = true; + MaskImagePointer binarizedContour; + try { + DD("FindAntPostVessels try to get"); + binarizedContour = GetAFDB()->template GetImage ("AntPostVesselsSeparation"); + } + catch(clitk::ExceptionObject e) { + DD("not found"); + found = false; + } + if (found) { + return binarizedContour; + } + + /* Here, we consider the vessels form a kind of anterior barrier. We + link all vessels centroids and remove what is post to it. + - select the list of structure + vessel1 = BrachioCephalicArtery + vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) + vessel3 = CommonCarotidArtery + vessel4 = SubclavianArtery + other = Thyroid + other = Aorta + - crop images as needed + - slice by slice, choose the CCL and extract centroids + - slice by slice, sort according to polar angle wrt Trachea centroid. + - slice by slice, link centroids and close contour + - remove outside this contour + - merge with support + */ + + // Read structures + std::map MapOfStructures; + typedef std::map::iterator MapIter; + MapOfStructures["BrachioCephalicArtery"] = GetAFDB()->template GetImage("BrachioCephalicArtery"); + MapOfStructures["BrachioCephalicVein"] = GetAFDB()->template GetImage("BrachioCephalicVein"); + MapOfStructures["CommonCarotidArteryLeft"] = GetAFDB()->template GetImage("CommonCarotidArteryLeft"); + MapOfStructures["CommonCarotidArteryRight"] = GetAFDB()->template GetImage("CommonCarotidArteryRight"); + MapOfStructures["SubclavianArteryLeft"] = GetAFDB()->template GetImage("SubclavianArteryLeft"); + MapOfStructures["SubclavianArteryRight"] = GetAFDB()->template GetImage("SubclavianArteryRight"); + MapOfStructures["Thyroid"] = GetAFDB()->template GetImage("Thyroid"); + MapOfStructures["Aorta"] = GetAFDB()->template GetImage("Aorta"); + MapOfStructures["Trachea"] = GetAFDB()->template GetImage("Trachea"); + + std::vector ListOfStructuresNames; + + // Create a temporay support + // From first slice of BrachioCephalicVein to end of 3A + MaskImagePointer support = GetAFDB()->template GetImage("Support_Sup_Carina"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(MapOfStructures["BrachioCephalicVein"], + GetBackgroundValue(), 2, true, p); + double inf = p[2]; + clitk::FindExtremaPointInAGivenDirection(GetAFDB()->template GetImage("Support_S3A"), + GetBackgroundValue(), 2, false, p); + double sup = p[2]+support->GetSpacing()[2];//one slice more ? + support = clitk::CropImageAlongOneAxis(support, 2, inf, sup, + false, GetBackgroundValue()); + + // Resize all structures like support + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + it->second = clitk::ResizeImageLike(it->second, support, GetBackgroundValue()); + } + + // Extract slices + typedef std::vector SlicesType; + std::map MapOfSlices; + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + SlicesType s; + clitk::ExtractSlices(it->second, 2, s); + MapOfSlices[it->first] = s; + } + + unsigned int n= MapOfSlices["Trachea"].size(); + + // Get the boundaries of one slice + std::vector bounds; + ComputeImageBoundariesCoordinates(MapOfSlices["Trachea"][0], bounds); + + // LOOP ON SLICES + // For all slices, for all structures, find the centroid and build the contour + // List of 3D points (for debug) + std::vector p3D; + vtkSmartPointer append = vtkSmartPointer::New(); + for(unsigned int i=0; ifirst][i]; + s = clitk::Labelize(s, GetBackgroundValue(), true, 1); + } + + // Search centroids + std::vector points2D; + typedef std::vector CentroidsType; + std::map MapOfCentroids; + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + std::string structure = it->first; + MaskSlicePointer & s = MapOfSlices[structure][i]; + CentroidsType c; + clitk::ComputeCentroids(s, GetBackgroundValue(), c); + MapOfCentroids[structure] = c; + } + + + // BrachioCephalicVein -> when it is separated into two CCL, we + // only consider the most at Right one + CentroidsType & c = MapOfCentroids["BrachioCephalicVein"]; + if (c.size() > 2) { + if (c[1][0] < c[2][0]) c.erase(c.begin()+2); + else c.erase(c.begin()+1); + } + + /* + // BrachioCephalicVein -> when SubclavianArtery has 2 CCL + // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein + if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1) + && (MapOfCentroids["SubclavianArtery"].size() > 2)) { + MapOfCentroids["BrachioCephalicVein"].clear(); + } + */ + + // Add all 2D points + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + std::string structure = it->first; + if (structure != "Trachea") { // do not add centroid of trachea + CentroidsType & c = MapOfCentroids[structure]; + for(unsigned int j=1; j centroids_trachea; + //ComputeCentroids(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea); + CentroidsType centroids_trachea = MapOfCentroids["Trachea"]; + typedef std::pair PointAngleType; + std::vector angles; + for(unsigned int j=0; j0) angle = atan(y/x); + if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; + if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; + if (x==0) { + if (y>0) angle = M_PI/2.0; + if (y<0) angle = -M_PI/2.0; + if (y==0) angle = 0; + } + angle = clitk::rad2deg(angle); + // Angle is [-180;180] wrt the X axis. We change the X axis to + // be the vertical line, because we want to sort from Right to + // Left from Post to Ant. + if (angle>0) + angle = (270-angle); + if (angle<0) { + angle = -angle-90; + if (angle<0) angle = 360-angle; + } + PointAngleType a; + a.first = points2D[j]; + a.second = angle; + angles.push_back(a); + } + + // Do nothing if less than 2 points --> n + if (points2D.size() < 3) { //continue; + continue; + } + + // Sort points2D according to polar angles + std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); + for(unsigned int j=0; j toadd; + unsigned int index = 0; + double dmax = 5; + while (indexdmax) { + + MaskSlicePointType b; + b[0] = a[0]+(c[0]-a[0])/2.0; + b[1] = a[1]+(c[1]-a[1])/2.0; + + // Compute distance to trachea centroid + MaskSlicePointType m = centroids_trachea[1]; + double da = m.EuclideanDistanceTo(a); + double db = m.EuclideanDistanceTo(b); + //double dc = m.EuclideanDistanceTo(c); + + // Mean distance, find point on the line from trachea centroid + // to b + double alpha = (da+db)/2.0; + MaskSlicePointType v; + double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); + v[0] = (b[0]-m[0])/n; + v[1] = (b[1]-m[1])/n; + MaskSlicePointType r; + r[0] = m[0]+alpha*v[0]; + r[1] = m[1]+alpha*v[1]; + points2D.insert(points2D.begin()+index+1, r); + } + else { + index++; + } + } + // DDV(points2D, points2D.size()); + + // Add some points to close the contour + // - H line towards Right + MaskSlicePointType p = points2D[0]; + p[0] = bounds[0][0]; + points2D.insert(points2D.begin(), p); + // - corner Right/Post + p = bounds[0]; + points2D.insert(points2D.begin(), p); + // - H line towards Left + p = points2D.back(); + p[0] = bounds[2][0]; + points2D.push_back(p); + // - corner Right/Post + p = bounds[2]; + points2D.push_back(p); + // Close contour with the first point + points2D.push_back(points2D[0]); + // DDV(points2D, points2D.size()); + + // Build 3D points from the 2D points + std::vector points3D; + clitk::PointsUtils::Convert2DListTo3DList(points2D, i, support, points3D); + for(unsigned int x=0; x mesh = Build3DMeshFrom2DContour(points3D); + append->AddInput(mesh); + // if (i ==n-1) { // last slice + // clitk::PointsUtils::Convert2DListTo3DList(points2D, i+1, support, points3D); + // vtkSmartPointer mesh = Build3DMeshFrom2DContour(points3D); + // append->AddInput(mesh); + // } + } + + // DEBUG: write points3D + clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); + + // Build the final 3D mesh form the list 2D mesh + append->Update(); + vtkSmartPointer mesh = append->GetOutput(); + + // Debug, write the mesh + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(mesh); + w->SetFileName("bidon.vtk"); + w->Write(); + */ + + // Compute a single binary 3D image from the list of contours + clitk::MeshToBinaryImageFilter::Pointer filter = + clitk::MeshToBinaryImageFilter::New(); + filter->SetMesh(mesh); + filter->SetLikeImage(support); + filter->Update(); + binarizedContour = filter->GetOutput(); + + // Crop + clitk::FindExtremaPointInAGivenDirection(binarizedContour, + GetForegroundValue(), 2, true, p); + inf = p[2]; + clitk::FindExtremaPointInAGivenDirection(binarizedContour, + GetForegroundValue(), 2, false, p); + sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice + binarizedContour = clitk::CropImageAlongOneAxis(binarizedContour, 2, inf, sup, + false, GetBackgroundValue()); + + // Update the AFDB + writeImage(binarizedContour, "seg/AntPostVesselsSeparation.mhd"); + GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd"); + WriteAFDB(); + return binarizedContour; + + /* + // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse + ImagePointType p = p3D[2]; // This is the first centroid of the first slice + p[1] += 50; // 50 mm Post from this point must be kept + ImageIndexType index; + binarizedContour->TransformPhysicalPointToIndex(p, index); + bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); + + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + if (isInside) + boolFilter->SetOperationType(BoolFilterType::And); + else + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + */ + + // End + //StopCurrentStep(m_Working_Support); + //m_ListOfStations["2R"] = m_Working_Support; + //m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + + #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX diff --git a/segmentation/clitkExtractLymphStationsGenericFilter.txx b/segmentation/clitkExtractLymphStationsGenericFilter.txx index 8b6167e..03fdb29 100644 --- a/segmentation/clitkExtractLymphStationsGenericFilter.txx +++ b/segmentation/clitkExtractLymphStationsGenericFilter.txx @@ -69,6 +69,8 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetVerboseMemoryFlag(mArgsInfo.verboseMemory_flag); f->SetAFDBFilename(mArgsInfo.afdb_arg); + f->SetComputeStationsSupportsFlag(!mArgsInfo.nosupport_flag); + // Station 8 //f->SetDistanceMaxToAnteriorPartOfTheSpine(mArgsInfo.S8_maxAntSpine_arg); f->SetFuzzyThreshold("8", "Esophagus", mArgsInfo.S8_ft_Esophagus_arg); @@ -111,6 +113,11 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) for(uint i=0; iAddComputeStation(mArgsInfo.station_arg[i]); + // Station 3A + f->SetFuzzyThreshold("3A", "SVC", mArgsInfo.S3A_ft_SVC_arg); + f->SetFuzzyThreshold("3A", "Aorta", mArgsInfo.S3A_ft_Aorta_arg); + f->SetFuzzyThreshold("3A", "SubclavianArtery", mArgsInfo.S3A_ft_SubclavianArtery_arg); + // Station 7 f->SetFuzzyThreshold("7", "Bronchi", mArgsInfo.S7_ft_Bronchi_arg); f->SetFuzzyThreshold("7", "LeftSuperiorPulmonaryVein", mArgsInfo.S7_ft_LeftSuperiorPulmonaryVein_arg); @@ -120,10 +127,6 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetFuzzyThreshold("7", "SVC", mArgsInfo.S7_ft_SVC_arg); f->SetS7_UseMostInferiorPartOnlyFlag(mArgsInfo.S7_UseMostInferiorPartOnly_flag); - // Station 3A - f->SetFuzzyThreshold("3A", "Sternum", mArgsInfo.S3A_ft_Sternum_arg); - f->SetFuzzyThreshold("3A", "SubclavianArtery", mArgsInfo.S3A_ft_SubclavianArtery_arg); - // Station 2RL f->SetFuzzyThreshold("2RL", "CommonCarotidArtery", mArgsInfo.S2RL_ft_CommonCarotidArtery_arg); f->SetFuzzyThreshold("2RL", "BrachioCephalicTrunk", mArgsInfo.S2RL_ft_BrachioCephalicTrunk_arg);