- // Resize all structures like support
- BrachioCephalicArtery =
- clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue());
- CommonCarotidArtery =
- clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, m_Working_Support, GetBackgroundValue());
- SubclavianArtery =
- clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, m_Working_Support, GetBackgroundValue());
- Thyroid =
- clitk::ResizeImageLike<MaskImageType>(Thyroid, m_Working_Support, GetBackgroundValue());
- Aorta =
- clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
- BrachioCephalicVein =
- clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, m_Working_Support, GetBackgroundValue());
- Trachea =
- clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
-
- // Extract slices
- std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
- clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
- std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
- clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
- std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
- clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
- std::vector<MaskSlicePointer> slices_SubclavianArtery;
- clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
- std::vector<MaskSlicePointer> slices_Thyroid;
- clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
- std::vector<MaskSlicePointer> slices_Aorta;
- clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
- std::vector<MaskSlicePointer> slices_Trachea;
- clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
- uint n= slices_BrachioCephalicArtery.size();
-
- // 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(uint 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(uint j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
- for(uint j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
- for(uint j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
- for(uint j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
- for(uint j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
- for(uint 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(uint 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(uint 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;
- uint 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, m_Working_Support, points3D);
- for(uint 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(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());
-