+ // Get the boundaries of one slice
+ std::vector<MaskSlicePointType> bounds;
+ ComputeImageBoundariesCoordinates<MaskSliceType>(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<MaskImagePointType> p3D;
+
+ vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
+ for(unsigned int i=0; i<n; i++) {
+ // Labelize the slices
+ slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
+ GetBackgroundValue(), true, 1);
+ slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
+ GetBackgroundValue(), true, 1);
+ slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
+ GetBackgroundValue(), true, 1);
+ slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
+ GetBackgroundValue(), true, 1);
+ slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
+ GetBackgroundValue(), true, 1);
+ slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
+ GetBackgroundValue(), true, 1);
+
+ // Search centroids
+ std::vector<MaskSlicePointType> points2D;
+ std::vector<MaskSlicePointType> centroids1;
+ std::vector<MaskSlicePointType> centroids2;
+ std::vector<MaskSlicePointType> centroids3;
+ std::vector<MaskSlicePointType> centroids4;
+ std::vector<MaskSlicePointType> centroids5;
+ std::vector<MaskSlicePointType> centroids6;
+ ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
+ ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
+ ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
+ ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
+ ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
+ ComputeCentroids<MaskSliceType>(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<centroids1.size(); j++) points2D.push_back(centroids1[j]);
+ for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
+ for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
+ for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
+ for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
+ for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
+
+ // Sort by angle according to trachea centroid and vertical line,
+ // in polar coordinates :
+ // http://en.wikipedia.org/wiki/Polar_coordinate_system
+ std::vector<MaskSlicePointType> centroids_trachea;
+ ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
+ typedef std::pair<MaskSlicePointType, double> PointAngleType;
+ std::vector<PointAngleType> angles;
+ for(unsigned int j=0; j<points2D.size(); j++) {
+ //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
+ double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
+ double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
+ double angle = 0;
+ if (x>0) 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<PointAngleType>());
+ for(unsigned int j=0; j<angles.size(); j++) {
+ points2D[j] = angles[j].first;
+ }
+ // DDV(points2D, points2D.size());
+
+ /* When vessels are far away, we try to replace the line segment
+ with a curved line that join the two vessels but stay
+ approximately at the same distance from the trachea centroids
+ than the vessels.
+
+ For that:
+ - let a and c be two successive vessels centroids
+ - id distance(a,c) < threshold, next point
+
+ TODO HERE
+
+ - compute mid position between two successive points -
+ compute dist to trachea centroid for the 3 pts - if middle too
+ low, add one point
+ */
+ std::vector<MaskSlicePointType> toadd;
+ unsigned int index = 0;
+ double dmax = 5;
+ while (index<points2D.size()-1) {
+ MaskSlicePointType a = points2D[index];
+ MaskSlicePointType c = points2D[index+1];
+
+ double dac = a.EuclideanDistanceTo(c);
+ if (dac>dmax) {
+
+ 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<ImagePointType> points3D;
+ clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
+ for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
+
+ // Build the mesh from the contour's points
+ vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
+ append->AddInput(mesh);
+ }
+
+ // DEBUG: write points3D
+ clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
+
+ // Build the final 3D mesh form the list 2D mesh
+ append->Update();
+ vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
+
+ // Debug, write the mesh
+ /*
+ vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
+ w->SetInput(mesh);
+ w->SetFileName("bidon.vtk");
+ w->Write();
+ */
+
+ // Compute a single binary 3D image from the list of contours
+ clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
+ clitk::MeshToBinaryImageFilter<MaskImageType>::New();
+ filter->SetMesh(mesh);
+ filter->SetLikeImage(support);
+ filter->Update();
+ binarizedContour = filter->GetOutput();
+
+ // Crop
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
+ inf = p[2];
+ DD(p);
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
+ sup = p[2];
+ DD(p);
+ binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
+ false, GetBackgroundValue());
+ // Update the AFDB
+ writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
+ this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
+ this->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<MaskImageType> 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<MaskImageType>(m_Working_Support);
+ //m_ListOfStations["2R"] = m_Working_Support;
+ //m_ListOfStations["2L"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 {
+ binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
+ }
+ catch(clitk::ExceptionObject e) {
+ 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<std::string, MaskImagePointer> MapOfStructures;
+ typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
+ MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
+ MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
+ MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
+ MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
+ MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
+ MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
+ MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
+ MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
+
+ std::vector<std::string> ListOfStructuresNames;
+
+ // Create a temporay support
+ // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
+ MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
+ GetBackgroundValue(), 2, true, p);
+ double inf = p[2];
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
+ GetBackgroundValue(), 2, false, p);
+ MaskImagePointType p2;
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
+ GetBackgroundValue(), 2, false, p2);
+ if (p2[2] > p[2]) p = p2;
+
+ double sup = p[2]+support->GetSpacing()[2];//one slice more ?
+ support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
+ false, GetBackgroundValue());
+
+ // Resize all structures like support
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
+ }
+
+ // Extract slices
+ typedef std::vector<MaskSlicePointer> SlicesType;
+ std::map<std::string, SlicesType> MapOfSlices;
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ SlicesType s;
+ clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
+ MapOfSlices[it->first] = s;
+ }
+
+ unsigned int n= MapOfSlices["Trachea"].size();
+
+ // Get the boundaries of one slice
+ std::vector<MaskSlicePointType> bounds;
+ ComputeImageBoundariesCoordinates<MaskSliceType>(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<MaskImagePointType> p3D;
+ vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
+ for(unsigned int i=0; i<n; i++) {
+
+ // Labelize the slices
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ MaskSlicePointer & s = MapOfSlices[it->first][i];
+ s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
+ }
+
+ // Search centroids
+ std::vector<MaskSlicePointType> points2D;
+ typedef std::vector<MaskSlicePointType> CentroidsType;
+ std::map<std::string, CentroidsType> MapOfCentroids;
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ std::string structure = it->first;
+ MaskSlicePointer & s = MapOfSlices[structure][i];
+ CentroidsType c;
+ clitk::ComputeCentroids<MaskSliceType>(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);
+ }
+