Fraxinus  16.5.0-fx-rc7
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxRouteToTarget.cpp
Go to the documentation of this file.
1 
2 
3 #include "cxRouteToTarget.h"
4 #include <vtkPolyData.h>
5 //#include "cxDoubleDataAdapterXml.h"
6 #include "cxBranchList.h"
7 #include "cxBranch.h"
8 #include <vtkCellArray.h>
9 #include "vtkCardinalSpline.h"
10 
11 typedef vtkSmartPointer<class vtkCardinalSpline> vtkCardinalSplinePtr;
12 
13 namespace cx
14 {
15 
17  mBranchListPtr(new BranchList),
18  mProjectedIndex(0)
19 {
20 }
21 
22 
24 {
25 }
26 
27 /*
28 void RouteToTarget::setCenterline(vtkPolyDataPtr centerline)
29 {
30  mCLpoints = this->getCenterlinePositions(centerline);
31 }
32 */
33 
35 {
36 
37  int N = centerline_r->GetNumberOfPoints();
38  Eigen::MatrixXd CLpoints(3,N);
39  for(vtkIdType i = 0; i < N; i++)
40  {
41  double p[3];
42  centerline_r->GetPoint(i,p);
43  Eigen::Vector3d position;
44  position(0) = p[0]; position(1) = p[1]; position(2) = p[2];
45  CLpoints.block(0 , i , 3 , 1) = position;
46  }
47  return CLpoints;
48 }
49 
51 {
52  if (mBranchListPtr)
53  mBranchListPtr->deleteAllBranches();
54 
55  Eigen::MatrixXd CLpoints_r = getCenterlinePositions(centerline_r);
56 
57  mBranchListPtr->findBranchesInCenterline(CLpoints_r);
58 
59  mBranchListPtr->calculateOrientations();
60  mBranchListPtr->smoothOrientations();
61 
62  std::cout << "Number of branches in CT centerline: " << mBranchListPtr->getBranches().size() << std::endl;
63 }
64 
65 
67 {
68 
69  double minDistance = 100000;
70  int minDistancePositionIndex;
71  BranchPtr minDistanceBranch;
72  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
73  for (int i = 0; i < branches.size(); i++)
74  {
75  Eigen::MatrixXd positions = branches[i]->getPositions();
76  for (int j = 0; j < positions.cols(); j++)
77  {
78  double D = findDistance(positions.col(j), targetCoordinate_r);
79  if (D < minDistance)
80  {
81  minDistance = D;
82  minDistanceBranch = branches[i];
83  minDistancePositionIndex = j;
84  }
85  }
86  }
87 
88  mProjectedBranchPtr = minDistanceBranch;
89  mProjectedIndex = minDistancePositionIndex;
90 }
91 
92 
94 {
95  mRoutePositions.clear();
96 
97  searchBranchUp(mProjectedBranchPtr, mProjectedIndex);
98 }
99 
100 void RouteToTarget::searchBranchUp(BranchPtr searchBranchPtr, int startIndex)
101 {
102  Eigen::MatrixXd positions = searchBranchPtr->getPositions();
103 
104  for (int i = startIndex; i>=0; i--)
105  mRoutePositions.push_back(positions.col(i));
106 
107  BranchPtr parentBranchPtr = searchBranchPtr->getParentBranch();
108  if (parentBranchPtr)
109  searchBranchUp(parentBranchPtr, parentBranchPtr->getPositions().cols()-1);
110 }
111 
112 
114 {
115 
116  findClosestPointInBranches(targetCoordinate_r);
118 
119  smoothPositions();
120 
121  vtkPolyDataPtr retval = addVTKPoints();
122 
123  return retval;
124 }
125 
127 {
128  vtkPolyDataPtr retval = vtkPolyDataPtr::New();
129  vtkPointsPtr points = vtkPointsPtr::New();
130  vtkCellArrayPtr lines = vtkCellArrayPtr::New();
131  int numberOfPositions = mRoutePositions.size();
132  for (int j = numberOfPositions - 1; j >= 0; j--)
133  {
134  points->InsertNextPoint(mRoutePositions[j](0),mRoutePositions[j](1),mRoutePositions[j](2));
135  }
136  for (int j = 0; j < numberOfPositions-1; j++)
137  {
138  vtkIdType connection[2] = {j, j+1};
139  lines->InsertNextCell(2, connection);
140  }
141  retval->SetPoints(points);
142  retval->SetLines(lines);
143  return retval;
144 }
145 
146 void RouteToTarget::smoothPositions()
147 {
148  int numberOfInputPoints = mRoutePositions.size();
149  int controlPointFactor = 10;
150  int numberOfControlPoints = numberOfInputPoints / controlPointFactor;
151 
152  vtkCardinalSplinePtr splineX = vtkSmartPointer<vtkCardinalSpline>::New();
153  vtkCardinalSplinePtr splineY = vtkSmartPointer<vtkCardinalSpline>::New();
154  vtkCardinalSplinePtr splineZ = vtkSmartPointer<vtkCardinalSpline>::New();
155 
156  if (numberOfControlPoints >= 2)
157  {
158  //add control points to spline
159  for(int j=0; j<numberOfControlPoints; j++)
160  {
161  int indexP = (j*numberOfInputPoints)/numberOfControlPoints;
162  splineX->AddPoint(indexP,mRoutePositions[indexP](0));
163  splineY->AddPoint(indexP,mRoutePositions[indexP](1));
164  splineZ->AddPoint(indexP,mRoutePositions[indexP](2));
165  }
166  //Always add the last point to complete spline
167  splineX->AddPoint(numberOfInputPoints-1,mRoutePositions[numberOfInputPoints-1](0));
168  splineY->AddPoint(numberOfInputPoints-1,mRoutePositions[numberOfInputPoints-1](1));
169  splineZ->AddPoint(numberOfInputPoints-1,mRoutePositions[numberOfInputPoints-1](2));
170 
171 
172  //evaluate spline - get smoothed positions
173  std::vector< Eigen::Vector3d > smoothingResult;
174  for(int j=0; j<numberOfInputPoints; j++)
175  {
176  double splineParameter = j;
177  Eigen::Vector3d tempPoint;
178  tempPoint(0) = splineX->Evaluate(splineParameter);
179  tempPoint(1) = splineY->Evaluate(splineParameter);
180  tempPoint(2) = splineZ->Evaluate(splineParameter);
181  smoothingResult.push_back(tempPoint);
182  }
183  mRoutePositions.clear();
184  mRoutePositions = smoothingResult;
185  }
186 }
187 
188 
189 double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
190 {
191  double d0 = p1(0) - p2(0);
192  double d1 = p1(1) - p2(1);
193  double d2 = p1(2) - p2(2);
194 
195  double D = sqrt( d0*d0 + d1*d1 + d2*d2 );
196 
197  return D;
198 }
199 
200 
201 } /* namespace cx */
vtkSmartPointer< class vtkCellArray > vtkCellArrayPtr
boost::shared_ptr< class Branch > BranchPtr
double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
vtkSmartPointer< class vtkPolyData > vtkPolyDataPtr
void findClosestPointInBranches(Vector3D targetCoordinate_r)
void searchBranchUp(BranchPtr searchBranchPtr, int startIndex)
void processCenterline(vtkPolyDataPtr centerline_r)
vtkSmartPointer< class vtkCardinalSpline > vtkCardinalSplinePtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
vtkPolyDataPtr findRouteToTarget(Vector3D targetCoordinate_r)
vtkPolyDataPtr addVTKPoints()
vtkSmartPointer< class vtkPoints > vtkPointsPtr
Eigen::MatrixXd getCenterlinePositions(vtkPolyDataPtr centerline_r)