Fraxinus  17.12-rc3
An IGT application
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  //mBranchListPtr->smoothBranchPositions(40);
62 
63  std::cout << "Number of branches in CT centerline: " << mBranchListPtr->getBranches().size() << std::endl;
64 }
65 
66 
68 {
69 
70  double minDistance = 100000;
71  int minDistancePositionIndex;
72  BranchPtr minDistanceBranch;
73  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
74  for (int i = 0; i < branches.size(); i++)
75  {
76  Eigen::MatrixXd positions = branches[i]->getPositions();
77  for (int j = 0; j < positions.cols(); j++)
78  {
79  double D = findDistance(positions.col(j), targetCoordinate_r);
80  if (D < minDistance)
81  {
82  minDistance = D;
83  minDistanceBranch = branches[i];
84  minDistancePositionIndex = j;
85  }
86  }
87  }
88 
89  mProjectedBranchPtr = minDistanceBranch;
90  mProjectedIndex = minDistancePositionIndex;
91 }
92 
93 
95 {
96  mRoutePositions.clear();
97 
98  searchBranchUp(mProjectedBranchPtr, mProjectedIndex);
99 }
100 
101 /*
102  RouteToTarget::searchBranchUp is finding all positions from a given index on a branch and up
103  the airway tree to the top of trachea. All positions are added to mRoutePositions, which stores
104  all positions along the route-to-target.
105  Before the positions are added they are smoothed by RouteToTarget::smoothBranch.
106 */
107 void RouteToTarget::searchBranchUp(BranchPtr searchBranchPtr, int startIndex)
108 {
109  if (!searchBranchPtr)
110  return;
111 
112  std::vector< Eigen::Vector3d > positions = smoothBranch(searchBranchPtr, startIndex, searchBranchPtr->getPositions().col(startIndex));
113 
114  for (int i = 0; i<=startIndex && i<positions.size(); i++)
115  mRoutePositions.push_back(positions[i]);
116 
117  BranchPtr parentBranchPtr = searchBranchPtr->getParentBranch();
118  if (parentBranchPtr)
119  searchBranchUp(parentBranchPtr, parentBranchPtr->getPositions().cols()-1);
120 }
121 
122 
124 {
125 
126  findClosestPointInBranches(targetCoordinate_r);
128 
129  //smoothPositions();
130 
131  vtkPolyDataPtr retval = addVTKPoints(mRoutePositions);
132 
133  return retval;
134 }
135 
137 {
138  float extentionPointIncrement = 0.5; //mm
139  mExtendedRoutePositions.clear();
140  mExtendedRoutePositions = mRoutePositions;
141  if(mRoutePositions.size() > 0)
142  {
143  double extentionDistance = findDistance(mRoutePositions.front(),targetCoordinate_r);
144  Eigen::Vector3d extentionVector = ( targetCoordinate_r - mRoutePositions.front() ) / extentionDistance;
145  int numberOfextentionPoints = (int) extentionDistance * extentionPointIncrement;
146  Eigen::Vector3d extentionPointIncrementVector = extentionVector / extentionPointIncrement;
147 
148  for (int i = 1; i<= numberOfextentionPoints; i++)
149  {
150  mExtendedRoutePositions.insert(mExtendedRoutePositions.begin(), mRoutePositions.front() + extentionPointIncrementVector*i);
151  //std::cout << mRoutePositions.front() + extentionPointIncrementVector*i << std::endl;
152  }
153  }
154 
155  vtkPolyDataPtr retval = addVTKPoints(mExtendedRoutePositions);
156  return retval;
157 }
158 
159 vtkPolyDataPtr RouteToTarget::addVTKPoints(std::vector<Eigen::Vector3d> positions)
160 {
161  vtkPolyDataPtr retval = vtkPolyDataPtr::New();
162  vtkPointsPtr points = vtkPointsPtr::New();
163  vtkCellArrayPtr lines = vtkCellArrayPtr::New();
164  int numberOfPositions = positions.size();
165  for (int j = numberOfPositions - 1; j >= 0; j--)
166  {
167  points->InsertNextPoint(positions[j](0),positions[j](1),positions[j](2));
168  }
169  for (int j = 0; j < numberOfPositions-1; j++)
170  {
171  vtkIdType connection[2] = {j, j+1};
172  lines->InsertNextCell(2, connection);
173  }
174  retval->SetPoints(points);
175  retval->SetLines(lines);
176  return retval;
177 }
178 
179 /*
180  RouteToTarget::smoothBranch is smoothing the positions of a centerline branch by using vtkCardinalSpline.
181  The degree of smoothing is dependent on the branch radius and the shape of the branch.
182  First, the method tests if a straight line from start to end of the branch is sufficient by the condition of
183  all positions along the line being within the lumen of the airway (max distance from original centerline
184  is set to branch radius).
185  If this fails, one more control point is added to the spline at the time, until the condition is fulfilled.
186  The control point added for each iteration is the position with the larges deviation from the original/unfiltered
187  centerline.
188 */
189 std::vector< Eigen::Vector3d > RouteToTarget::smoothBranch(BranchPtr branchPtr, int startIndex, Eigen::MatrixXd startPosition)
190 {
191  vtkCardinalSplinePtr splineX = vtkSmartPointer<vtkCardinalSpline>::New();
192  vtkCardinalSplinePtr splineY = vtkSmartPointer<vtkCardinalSpline>::New();
193  vtkCardinalSplinePtr splineZ = vtkSmartPointer<vtkCardinalSpline>::New();
194 
195  double branchRadius = branchPtr->findBranchRadius();
196 
197  //add control points to spline
198 
199  //add first position
200  Eigen::MatrixXd positions = branchPtr->getPositions();
201  splineX->AddPoint(0,startPosition(0));
202  splineY->AddPoint(0,startPosition(1));
203  splineZ->AddPoint(0,startPosition(2));
204 
205 
206  // Add last position if no parent branch, else add parents position closest to current branch.
207  // Branch positions are stored in order from head to feet (e.g. first position is top of trachea),
208  // while route-to-target is generated from target to top of trachea.
209  if(!branchPtr->getParentBranch())
210  {
211  splineX->AddPoint(startIndex,positions(0,0));
212  splineY->AddPoint(startIndex,positions(1,0));
213  splineZ->AddPoint(startIndex,positions(2,0));
214  }
215  else
216  {
217  Eigen::MatrixXd parentPositions = branchPtr->getParentBranch()->getPositions();
218  splineX->AddPoint(startIndex,parentPositions(0,parentPositions.cols()-1));
219  splineY->AddPoint(startIndex,parentPositions(1,parentPositions.cols()-1));
220  splineZ->AddPoint(startIndex,parentPositions(2,parentPositions.cols()-1));
221 
222  }
223 
224  //Add points until all filtered/smoothed postions are within branch radius distance of the unfiltered branch centerline posiition.
225  //This is to make sure the smoothed centerline is within the lumen of the airways.
226  double maxDistanceToOriginalPosition = branchRadius;
227  int maxDistanceIndex = -1;
228  std::vector< Eigen::Vector3d > smoothingResult;
229 
230  //add positions to spline
231  while (maxDistanceToOriginalPosition >= branchRadius && splineX->GetNumberOfPoints() < startIndex)
232  {
233  if(maxDistanceIndex > 0)
234  {
235  //add to spline the position with largest distance to original position
236  splineX->AddPoint(maxDistanceIndex,positions(0,startIndex - maxDistanceIndex));
237  splineY->AddPoint(maxDistanceIndex,positions(1,startIndex - maxDistanceIndex));
238  splineZ->AddPoint(maxDistanceIndex,positions(2,startIndex - maxDistanceIndex));
239  }
240 
241  //evaluate spline - get smoothed positions
242  maxDistanceToOriginalPosition = 0.0;
243  smoothingResult.clear();
244  for(int j=0; j<=startIndex; j++)
245  {
246  double splineParameter = j;
247  Eigen::Vector3d tempPoint;
248  tempPoint(0) = splineX->Evaluate(splineParameter);
249  tempPoint(1) = splineY->Evaluate(splineParameter);
250  tempPoint(2) = splineZ->Evaluate(splineParameter);
251  smoothingResult.push_back(tempPoint);
252 
253  //calculate distance to original (non-filtered) position
254  double distance = findDistanceToLine(tempPoint, positions);
255  //finding the index with largest distance
256  if (distance > maxDistanceToOriginalPosition)
257  {
258  maxDistanceToOriginalPosition = distance;
259  maxDistanceIndex = j;
260  }
261  }
262  }
263 
264  return smoothingResult;
265 }
266 
267 double findDistanceToLine(Eigen::MatrixXd point, Eigen::MatrixXd line)
268 {
269  double minDistance = findDistance(point, line.col(0));
270  for (int i=1; i<line.cols(); i++)
271  if (minDistance > findDistance(point, line.col(i)))
272  minDistance = findDistance(point, line.col(i));
273 
274  return minDistance;
275 }
276 
277 double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
278 {
279  double d0 = p1(0) - p2(0);
280  double d1 = p1(1) - p2(1);
281  double d2 = p1(2) - p2(2);
282 
283  double D = sqrt( d0*d0 + d1*d1 + d2*d2 );
284 
285  return D;
286 }
287 
288 
289 } /* namespace cx */
vtkPolyDataPtr findExtendedRoute(Vector3D targetCoordinate_r)
double findDistanceToLine(Eigen::MatrixXd point, Eigen::MatrixXd line)
vtkSmartPointer< class vtkCellArray > vtkCellArrayPtr
vtkSmartPointer< vtkPoints > vtkPointsPtr
vtkPolyDataPtr addVTKPoints(std::vector< Eigen::Vector3d > positions)
boost::shared_ptr< class Branch > BranchPtr
double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
vtkSmartPointer< class vtkCardinalSpline > vtkCardinalSplinePtr
Definition: cxAccusurf.cpp:17
void findClosestPointInBranches(Vector3D targetCoordinate_r)
void searchBranchUp(BranchPtr searchBranchPtr, int startIndex)
void processCenterline(vtkPolyDataPtr centerline_r)
vtkSmartPointer< class vtkCardinalSpline > vtkCardinalSplinePtr
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
vtkPolyDataPtr findRouteToTarget(Vector3D targetCoordinate_r)
Eigen::MatrixXd getCenterlinePositions(vtkPolyDataPtr centerline_r)
Namespace for all CustusX production code.