X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLymphStationsFilter.txx;h=5f9bef0c1865923b1905d2cadf641e5a467d4838;hb=2300716104b55267692efc2c0e506b3c5c45df38;hp=aec05d782df8553d4a0bdd861f7e5d40fdf8eedb;hpb=907b0bad00cbf772fbf362879c74d673253f97bb;p=clitk.git diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index aec05d7..5f9bef0 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,10 +57,12 @@ ExtractLymphStationsFilter(): this->SetNumberOfRequiredInputs(1); SetBackgroundValue(0); SetForegroundValue(1); + ComputeStationsSupportsFlagOn(); // Default values ExtractStation_8_SetDefaultValues(); ExtractStation_3P_SetDefaultValues(); + ExtractStation_2RL_SetDefaultValues(); ExtractStation_3A_SetDefaultValues(); ExtractStation_7_SetDefaultValues(); } @@ -71,22 +79,59 @@ 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 + 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"); - StartSubStep(); ExtractStation_8(); - StopSubStep(); // Extract Station3P - StartNewStep("Station 3P"); - StartSubStep(); ExtractStation_3P(); - StopSubStep(); // Extract Station3A - StartNewStep("Station 3A"); - StartSubStep(); ExtractStation_3A(); + + // HERE + + // Extract Station2RL + StartNewStep("Station 2RL"); + StartSubStep(); + ExtractStation_2RL(); StopSubStep(); // Extract Station7 @@ -102,20 +147,6 @@ GenerateOutputInformation() { //ExtractStation_4RL(); StopSubStep(); } - - - // - // typedef clitk::BooleanOperatorLabelImageFilter BFilter; - //BFilter::Pointer merge = BFilter::New(); - // writeImage(m_Output, "ouput.mhd"); - //writeImage(m_Working_Support, "ws.mhd"); - /*merge->SetInput1(m_Station7); - merge->SetInput2(m_Station4RL); // support - merge->SetOperationType(BFilter::AndNot); CHANGE OPERATOR - merge->SetForegroundValue(4); - merge->Update(); - m_Output = merge->GetOutput(); - */ } //-------------------------------------------------------------------- @@ -196,27 +227,114 @@ AddComputeStation(std::string station) //-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +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 clitk::ExtractLymphStationsFilter:: -ExtractStation_8() +ExtractStation_4RL() { + DD("TODO"); + exit(0); + /* + WARNING ONLY 4R FIRST !!! (not same inf limits) + */ + ExtractStation_4RL_SI_Limits(); + ExtractStation_4RL_LR_Limits(); + ExtractStation_4RL_AP_Limits(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Remove_Structures(std::string station, std::string s) { - if (CheckForStation("8")) { - ExtractStation_8_SI_Limits(); - ExtractStation_8_Post_Limits(); - ExtractStation_8_Ant_Sup_Limits(); - ExtractStation_8_Ant_Inf_Limits(); - ExtractStation_8_LR_1_Limits(); - ExtractStation_8_LR_2_Limits(); - ExtractStation_8_LR_Limits(); - ExtractStation_8_Ant_Injected_Limits(); - ExtractStation_8_Single_CCL_Limits(); - ExtractStation_8_Remove_Structures(); - // Store image filenames into AFDB - writeImage(m_ListOfStations["8"], "seg/Station8.mhd"); - GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd"); - WriteAFDB(); + try { + StartNewStep("[Station"+station+"] Remove "+s); + MaskImagePointer Structure = GetAFDB()->template GetImage(s); + clitk::AndNot(m_Working_Support, Structure, GetBackgroundValue()); + } + catch(clitk::ExceptionObject e) { + std::cout << s << " not found, skip." << std::endl; } } //-------------------------------------------------------------------- @@ -226,22 +344,9 @@ ExtractStation_8() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_3P() +SetFuzzyThreshold(std::string station, std::string tag, double value) { - 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(); - } + m_FuzzyThreshold[station][tag] = value; } //-------------------------------------------------------------------- @@ -250,71 +355,1064 @@ ExtractStation_3P() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_7() { - if (CheckForStation("7")) { - ExtractStation_7_SI_Limits(); - ExtractStation_7_RL_Limits(); - ExtractStation_7_Posterior_Limits(); - // ExtractStation_8_Single_CCL_Limits(); // Yes the same than 8 - ExtractStation_7_Remove_Structures(); - // Store image filenames into AFDB - writeImage(m_ListOfStations["7"], "seg/Station7.mhd"); - GetAFDB()->SetImageFilename("Station7", "seg/Station7.mhd"); - WriteAFDB(); +SetThreshold(std::string station, std::string tag, double value) +{ + m_Threshold[station][tag] = value; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +GetFuzzyThreshold(std::string station, std::string tag) +{ + if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) { + clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+")."); + return 0.0; } + + if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) { + clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+"."); + return 0.0; + } + + return m_FuzzyThreshold[station][tag]; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template -void +double clitk::ExtractLymphStationsFilter:: -ExtractStation_3A() +GetThreshold(std::string station, std::string tag) { - if (CheckForStation("3A")) { - ExtractStation_3A_SI_Limits(); - ExtractStation_3A_Ant_Limits(); - // Store image filenames into AFDB - writeImage(m_ListOfStations["3A"], "seg/Station3A.mhd"); - GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); - WriteAFDB(); + if (m_Threshold.find(station) == m_Threshold.end()) { + clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+")."); + return 0.0; } + + if (m_Threshold[station].find(tag) == m_Threshold[station].end()) { + clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+"."); + return 0.0; + } + + return m_Threshold[station][tag]; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- template -void +void clitk::ExtractLymphStationsFilter:: -ExtractStation_4RL() { - DD("TODO"); - exit(0); +FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B) +{ + // Create line from A to B with + // A = upper border of LLL at left + // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus + + try { + GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A); + GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B); + } + catch(clitk::ExceptionObject & o) { + + DD("FindLineForS7S8Separation"); + // Load LeftLowerLobeBronchus and get centroid point + MaskImagePointer LeftLowerLobeBronchus = + GetAFDB()->template GetImage ("LeftLowerLobeBronchus"); + std::vector c; + clitk::ComputeCentroids(LeftLowerLobeBronchus, GetBackgroundValue(), c); + A = c[1]; + + // Load RightMiddleLobeBronchus and get superior point (not centroid here) + MaskImagePointer RightMiddleLobeBronchus = + GetAFDB()->template GetImage ("RightMiddleLobeBronchus"); + bool b = FindExtremaPointInAGivenDirection(RightMiddleLobeBronchus, + GetBackgroundValue(), + 2, false, B); + if (!b) { + clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort"); + } + + // Insert into the DB + GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A); + GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B); + } +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindCarina() +{ + double z; + try { + z = GetAFDB()->GetDouble("CarinaZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindCarinaSlicePosition"); + // Get Carina + MaskImagePointer Carina = GetAFDB()->template GetImage("Carina"); + + // Get Centroid and Z value + std::vector centroids; + clitk::ComputeCentroids(Carina, GetBackgroundValue(), centroids); + + // We dont need Carina structure from now + GetAFDB()->template ReleaseImage("Carina"); + + // Put inside the AFDB + GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]); + GetAFDB()->SetDouble("CarinaZ", centroids[1][2]); + WriteAFDB(); + z = centroids[1][2]; + } + return z; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +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 +clitk::ExtractLymphStationsFilter:: +FindLeftAndRightBronchi() +{ + try { + m_RightBronchus = GetAFDB()->template GetImage ("RightBronchus"); + m_LeftBronchus = GetAFDB()->template GetImage ("LeftBronchus"); + } + catch(clitk::ExceptionObject & o) { + + DD("FindLeftAndRightBronchi"); + // The goal is to separate the trachea inferiorly to the carina into + // a Left and Right bronchus. + + // Get the trachea + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + + // Get the Carina position + double m_CarinaZ = FindCarina(); + + // Consider only inferiorly to the Carina + MaskImagePointer m_Working_Trachea = + clitk::CropImageRemoveGreaterThan(Trachea, 2, m_CarinaZ, true, // AutoCrop + GetBackgroundValue()); + + // Labelize the trachea + m_Working_Trachea = Labelize(m_Working_Trachea, 0, true, 1); + + // Carina position must at the first slice that separate the two + // main bronchus (not superiorly). We thus first check that the + // upper slice is composed of at least two labels + MaskImagePointer RightBronchus; + MaskImagePointer LeftBronchus; + typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion()); + iter.SetFirstDirection(0); + iter.SetSecondDirection(1); + iter.GoToReverseBegin(); // Start from the end (because image is IS not SI) + int maxLabel=0; + while (!iter.IsAtReverseEndOfSlice()) { + while (!iter.IsAtReverseEndOfLine()) { + if (iter.Get() > maxLabel) maxLabel = iter.Get(); + --iter; + } + iter.PreviousLine(); + } + if (maxLabel < 2) { + clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort"); + } + + // Compute 3D centroids of both parts to identify the left from the + // right bronchus + std::vector c; + clitk::ComputeCentroids(m_Working_Trachea, GetBackgroundValue(), c); + ImagePointType C1 = c[1]; + ImagePointType C2 = c[2]; + + ImagePixelType rightLabel; + ImagePixelType leftLabel; + if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; } + else { rightLabel = 2; leftLabel = 1; } + + // Select LeftLabel (set one label to Backgroundvalue) + RightBronchus = + clitk::Binarize(m_Working_Trachea, rightLabel, rightLabel, + GetBackgroundValue(), GetForegroundValue()); + /* + SetBackground(m_Working_Trachea, m_Working_Trachea, + leftLabel, GetBackgroundValue(), false); + */ + LeftBronchus = clitk::Binarize(m_Working_Trachea, leftLabel, leftLabel, + GetBackgroundValue(), GetForegroundValue()); + /* + SetBackground(m_Working_Trachea, m_Working_Trachea, + rightLabel, GetBackgroundValue(), false); + */ + + // Crop images + RightBronchus = clitk::AutoCrop(RightBronchus, GetBackgroundValue()); + LeftBronchus = clitk::AutoCrop(LeftBronchus, GetBackgroundValue()); + + // Insert int AFDB if need after + GetAFDB()->template SetImage ("RightBronchus", "seg/rightBronchus.mhd", + RightBronchus, true); + GetAFDB()->template SetImage ("LeftBronchus", "seg/leftBronchus.mhd", + LeftBronchus, true); + } +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +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 /* - WARNING ONLY 4R FIRST !!! (not same inf limits) - */ - ExtractStation_4RL_SI_Limits(); - ExtractStation_4RL_LR_Limits(); - ExtractStation_4RL_AP_Limits(); + 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 -void +typename clitk::ExtractLymphStationsFilter::MaskImagePointer clitk::ExtractLymphStationsFilter:: -Remove_Structures(std::string s) +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 { - StartNewStep("[Station7] Remove "+s); - MaskImagePointer Structure = GetAFDB()->template GetImage(s); - clitk::AndNot(m_Working_Support, Structure, GetBackgroundValue()); + DD("FindAntPostVessels try to get"); + binarizedContour = GetAFDB()->template GetImage ("AntPostVesselsSeparation"); } catch(clitk::ExceptionObject e) { - std::cout << s << " not found, skip." << std::endl; + 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; } //--------------------------------------------------------------------