NorMIT-nav  22.09
An IGT application
cxRouteToTarget.cpp
Go to the documentation of this file.
1 
2 
3 #include "cxRouteToTarget.h"
4 #include <vtkPolyData.h>
5 #include "cxBranchList.h"
6 #include "cxBranch.h"
7 #include "cxAirwaysFromCenterline.h"
8 #include "cxPointMetric.h"
9 #include <vtkCellArray.h>
10 #include "vtkCardinalSpline.h"
11 #include "vtkImageData.h"
12 #include <boost/math/special_functions/round.hpp>
13 #include "cxLogger.h"
14 #include <QDir>
15 #include "cxTime.h"
16 #include "cxImage.h"
17 #include "cxVisServices.h"
18 #include <QTextStream>
19 #include <QJsonObject>
20 #include <QJsonArray>
21 #include <QList>
22 
23 #define PI 3.1415926535897
24 
25 typedef vtkSmartPointer<class vtkCardinalSpline> vtkCardinalSplinePtr;
26 
27 namespace cx
28 {
29 
30 RouteToTarget::RouteToTarget():
31  mBranchListPtr(new BranchList),
32  mProjectedIndex(0),
33  mBloodVesselBranchListPtr(new BranchList),
34  mProjectedBloodVesselIndex(0)
35 {
36 }
37 
38 
39 RouteToTarget::~RouteToTarget()
40 {
41 }
42 
43 void RouteToTarget::setBloodVesselVolume(ImagePtr bloodVesselVolume)
44 {
45  mBloodVesselVolume = bloodVesselVolume;
46 }
47 
48 Eigen::MatrixXd RouteToTarget::getCenterlinePositions(vtkPolyDataPtr centerline_r)
49 {
50 
51  int N = centerline_r->GetNumberOfPoints();
52  Eigen::MatrixXd CLpoints(3,N);
53  for(vtkIdType i = 0; i < N; i++)
54  {
55  double p[3];
56  centerline_r->GetPoint(i,p);
57  Eigen::Vector3d position;
58  position(0) = p[0]; position(1) = p[1]; position(2) = p[2];
59  CLpoints.block(0 , i , 3 , 1) = position;
60  }
61  return CLpoints;
62 }
63 
64 void RouteToTarget::setSmoothing(bool smoothing)
65 {
66  mSmoothing = smoothing; // default true
67 }
68 
69 void RouteToTarget::processCenterline(MeshPtr mesh)
70 {
71  if (mBranchListPtr)
72  mBranchListPtr->deleteAllBranches();
73 
74  vtkPolyDataPtr centerline_r = mesh->getTransformedPolyDataCopy(mesh->get_rMd());
75 
76  mCLpoints = getCenterlinePositions(centerline_r);
77 
78  mBranchListPtr->findBranchesInCenterline(mCLpoints);
79 
80  mBranchListPtr->smoothOrientations();
81  //mBranchListPtr->smoothBranchPositions(40);
82  mBranchListPtr->findBronchoscopeRotation();
83 
84  std::cout << "Number of branches in CT centerline: " << mBranchListPtr->getBranches().size() << std::endl;
85 }
86 
87 //Can be used instead of processCenterline if you wan to use an existing branchList.
88 void RouteToTarget::setBranchList(BranchListPtr branchList)
89 {
90  mBranchListPtr = branchList;
91 }
92 
93 void RouteToTarget::processBloodVesselCenterline(Eigen::MatrixXd positions)
94 {
95  if (mBloodVesselBranchListPtr)
96  mBloodVesselBranchListPtr->deleteAllBranches();
97 
98  mBloodVesselBranchListPtr->findBranchesInCenterline(positions, false);
99 
100  mBloodVesselBranchListPtr->smoothOrientations();
101  mBloodVesselBranchListPtr->smoothBranchPositions(40);
102  setBloodVesselRadius();
103  mBloodVesselBranchListPtr->smoothRadius();
104 
105  BranchPtr branchWithLargestRadius = mBloodVesselBranchListPtr->findBranchWithLargestRadius();
106 
107  if (branchWithLargestRadius->getParentBranch()) //If the largest branch has a parent, something is wrong. -> Reprocess with largest branch as root.
108  {
109  Eigen::MatrixXd positions = branchWithLargestRadius->getPositions();
110  Eigen::VectorXd radius = branchWithLargestRadius->getRadius();
111  int maxRadiusIndex;
112  radius.maxCoeff(&maxRadiusIndex);
113  positions.col(positions.cols()-1).swap(positions.col(maxRadiusIndex)); //setting largest radius position last (which is processed first)
114 
115  std::vector<BranchPtr> allBranches = mBloodVesselBranchListPtr->getBranches();
116  for (int i = 1; i < allBranches.size(); i++)
117  if (allBranches[i] != branchWithLargestRadius) //add positions from all other branches
118  {
119  Eigen::MatrixXd positionsToAppend = allBranches[i]->getPositions();
120  Eigen::MatrixXd resizedPositions(positions.rows(), positions.cols() + positionsToAppend.cols());
121  resizedPositions.rightCols(positions.cols()) = positions;
122  resizedPositions.leftCols(positionsToAppend.cols()) = positionsToAppend;
123  positions = resizedPositions;
124  }
125 
126  mBloodVesselBranchListPtr->deleteAllBranches();
127 
128  mBloodVesselBranchListPtr->findBranchesInCenterline(positions, false);
129  mBloodVesselBranchListPtr->smoothOrientations();
130  mBloodVesselBranchListPtr->smoothBranchPositions(40);
131  setBloodVesselRadius();
132  mBloodVesselBranchListPtr->smoothRadius();
133  }
134 
135 
136  CX_LOG_INFO() << "Number of branches in CT blood vessel centerline: " << mBloodVesselBranchListPtr->getBranches().size();
137 }
138 
139 void RouteToTarget::findClosestPointInBranches(Vector3D targetCoordinate_r)
140 {
141  double minDistance = 100000;
142  int minDistancePositionIndex;
143  BranchPtr minDistanceBranch;
144  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
145  for (int i = 0; i < branches.size(); i++)
146  {
147  Eigen::MatrixXd positions = branches[i]->getPositions();
148  for (int j = 0; j < positions.cols(); j++)
149  {
150  double D = findDistance(positions.col(j), targetCoordinate_r);
151  if (D < minDistance)
152  {
153  minDistance = D;
154  minDistanceBranch = branches[i];
155  minDistancePositionIndex = j;
156  }
157  }
158  }
159 
160  mProjectedBranchPtr = minDistanceBranch;
161  mProjectedIndex = minDistancePositionIndex;
162 }
163 
164 void RouteToTarget::findClosestPointInBloodVesselBranches(Vector3D targetCoordinate_r)
165 {
166  double minDistance = 100000;
167  int minDistancePositionIndex;
168  BranchPtr minDistanceBranch;
169  std::vector<BranchPtr> branches = mBloodVesselBranchListPtr->getBranches();
170  for (int i = 0; i < branches.size(); i++)
171  {
172  Eigen::MatrixXd positions = branches[i]->getPositions();
173  for (int j = 0; j < positions.cols(); j++)
174  {
175  double D = findDistance(positions.col(j), targetCoordinate_r);
176  if (D < minDistance)
177  {
178  minDistance = D;
179  minDistanceBranch = branches[i];
180  minDistancePositionIndex = j;
181  }
182  }
183  }
184 
185  mProjectedBloodVesselBranchPtr = minDistanceBranch;
186  mProjectedBloodVesselIndex = minDistancePositionIndex;
187 }
188 
189 
190 void RouteToTarget::findRoutePositions()
191 {
192  mRoutePositions.clear();
193 
194  searchBranchUp(mProjectedBranchPtr, mProjectedIndex);
195 }
196 
197 void RouteToTarget::findRoutePositionsInBloodVessels()
198 {
199  mBloodVesselRoutePositions.clear();
200 
201  searchBloodVesselBranchUp(mProjectedBloodVesselBranchPtr, mProjectedBloodVesselIndex);
202 }
203 
204 /*
205  RouteToTarget::searchBranchUp is finding all positions from a given index on a branch and up
206  the airway tree to the top of trachea. All positions are added to mRoutePositions, which stores
207  all positions along the route-to-target.
208  Before the positions are added they are smoothed by RouteToTarget::smoothBranch.
209 */
210 void RouteToTarget::searchBranchUp(BranchPtr searchBranchPtr, int startIndex)
211 {
212  if (!searchBranchPtr)
213  return;
214 
215  std::vector< Eigen::Vector3d > positions;
216  if (mSmoothing)
217  positions = smoothBranch(searchBranchPtr, startIndex, searchBranchPtr->getPositions().col(startIndex));
218  else
219  positions = getBranchPositions(searchBranchPtr, startIndex);
220 
221  double cameraRotation = searchBranchPtr->getBronchoscopeRotation();
222 
223  for (int i = 0; i<=startIndex && i<positions.size(); i++)
224  {
225  mRoutePositions.push_back(positions[i]);
226  mCameraRotation.push_back(cameraRotation);
227  }
228 
229  mBranchingIndex.push_back(mRoutePositions.size()-1);
230 
231  BranchPtr parentBranchPtr = searchBranchPtr->getParentBranch();
232  if (parentBranchPtr)
233  searchBranchUp(parentBranchPtr, parentBranchPtr->getPositions().cols()-1);
234 }
235 
236 void RouteToTarget::searchBloodVesselBranchUp(BranchPtr searchBranchPtr, int startIndex)
237 {
238  if (!searchBranchPtr)
239  return;
240 
241  //std::vector< Eigen::Vector3d > positions = smoothBranch(searchBranchPtr, startIndex, searchBranchPtr->getPositions().col(startIndex));
242  std::vector< Eigen::Vector3d > positions = getBranchPositions(searchBranchPtr, startIndex);
243 
244  for (int i = 0; i<=startIndex && i<positions.size(); i++)
245  {
246  mBloodVesselRoutePositions.push_back(positions[i]);
247  }
248 
249  BranchPtr parentBranchPtr = searchBranchPtr->getParentBranch();
250  if (parentBranchPtr)
251  {
252  if(parentBranchPtr->getAverageRadius() >= searchBranchPtr->getAverageRadius() && variance(parentBranchPtr->getRadius()) < 1.0 )
253  {
254  searchBloodVesselBranchUp(parentBranchPtr, parentBranchPtr->getPositions().cols()-1);
255  }
256  }
257 }
258 
259 
260 vtkPolyDataPtr RouteToTarget::findRouteToTarget(PointMetricPtr targetPoint)
261 {
262  mTargetPosition = targetPoint->getCoordinate();
263 
264  findClosestPointInBranches(mTargetPosition);
265  findRoutePositions();
266 
267  vtkPolyDataPtr retval = addVTKPoints(mRoutePositions);
268 
269  return retval;
270 }
271 
272 
273 vtkPolyDataPtr RouteToTarget::findExtendedRoute(PointMetricPtr targetPoint)
274 {
275  mTargetPosition = targetPoint->getCoordinate();
276  double extensionPointIncrement = 0.25; //mm
277  mExtendedRoutePositions.clear();
278  mExtendedRoutePositions = mRoutePositions;
279  mExtendedCameraRotation.clear();
280  mExtendedCameraRotation = mCameraRotation;
281  if(mRoutePositions.size() > 0)
282  {
283  double extensionDistance = findDistance(mRoutePositions.front(),mTargetPosition);
284  Eigen::Vector3d extensionVectorNormalized = ( mTargetPosition - mRoutePositions.front() ) / extensionDistance;
285  int numberOfextensionPoints = int(extensionDistance / extensionPointIncrement);
286  Eigen::Vector3d extensionPointIncrementVector = extensionVectorNormalized * extensionDistance / numberOfextensionPoints;
287 
288  for (int i = 1; i<= numberOfextensionPoints; i++)
289  {
290  mExtendedRoutePositions.insert(mExtendedRoutePositions.begin(), mRoutePositions.front() + extensionPointIncrementVector*i);
291  mExtendedCameraRotation.insert(mExtendedCameraRotation.begin(), mExtendedCameraRotation.front());
292  }
293  }
294 
295  vtkPolyDataPtr retval = addVTKPoints(mExtendedRoutePositions);
296  return retval;
297 }
298 
299 
300 vtkPolyDataPtr RouteToTarget::findRouteToTargetAlongBloodVesselCenterlines(MeshPtr bloodVesselCenterlineMesh, PointMetricPtr targetPoint)
301 {
302  vtkPolyDataPtr retval;
303 
304  vtkPolyDataPtr BVcenterline_r = bloodVesselCenterlineMesh->getTransformedPolyDataCopy(bloodVesselCenterlineMesh->get_rMd());
305  Eigen::MatrixXd BVCLpoints_r = getCenterlinePositions(BVcenterline_r);
306  mConnectedPointsInBVCL = findClosestBloodVesselSegments(BVCLpoints_r , mCLpoints, targetPoint->getCoordinate());
307 
308  if (mRoutePositions.empty())
309  return retval;
310 
311  Eigen::MatrixXd::Index closestPositionToEndOfAirwayRTTIndex = dsearch(mRoutePositions[0], mConnectedPointsInBVCL).first;
312 
313  //setting position closest to RTT from airways in first index, where RTT should continue
314  mConnectedPointsInBVCL.col(mConnectedPointsInBVCL.cols()-1).swap(mConnectedPointsInBVCL.col(closestPositionToEndOfAirwayRTTIndex));
315 
316  processBloodVesselCenterline(mConnectedPointsInBVCL);
317  findClosestPointInBloodVesselBranches(mTargetPosition);
318  findRoutePositionsInBloodVessels();
319  mPathToBloodVesselsFound = makeConnectedAirwayAndBloodVesselRoute();
320 
321  retval = addVTKPoints(mBloodVesselRoutePositions);
322 
323  return retval;
324 }
325 
326 vtkPolyDataPtr RouteToTarget::generateAirwaysFromBloodVesselCenterlines()
327 {
328  vtkPolyDataPtr airwayMesh;
329  if (mConnectedPointsInBVCL.cols() == 0 || !mPathToBloodVesselsFound)
330  return airwayMesh;
331 
332  AirwaysFromCenterlinePtr airwaysFromBVCenterlinePtr = AirwaysFromCenterlinePtr(new AirwaysFromCenterline());
333  airwaysFromBVCenterlinePtr->setTypeToBloodVessel(true);
334  mBloodVesselBranchListPtr->interpolateBranchPositions(0.1);
335  airwaysFromBVCenterlinePtr->setBranches(mBloodVesselBranchListPtr);
336 
337  airwayMesh = airwaysFromBVCenterlinePtr->generateTubes(2);
338 
339  return airwayMesh;
340 }
341 
342 bool RouteToTarget::makeConnectedAirwayAndBloodVesselRoute()
343 {
344  //vtkPolyDataPtr mergedRoute;
345 
346  if ( mRoutePositions.empty() )
347  return false;
348 
349  if ( mBloodVesselRoutePositions.empty() )
350  {
351  mMergedAirwayAndBloodVesselRoutePositions = mRoutePositions;
352  return false;
353  }
354 
355  if ( findDistance(mRoutePositions.front(),mTargetPosition) < findDistance(mBloodVesselRoutePositions.front(),mTargetPosition) )
356  {
357  CX_LOG_INFO() << "No improved route to target found along blood vessels";
358  mMergedAirwayAndBloodVesselRoutePositions = mRoutePositions; // do not add blood vessel route if that is not leading closer to target than airway route alone
359  return false;
360  }
361 
362  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
363  double minDistance = ( mRoutePositions[0] - mBloodVesselRoutePositions[mBloodVesselRoutePositions.size()-1] ).norm();
364  int connectionIndexBloodVesselRoute = mBloodVesselRoutePositions.size()-1;
365  bool closerAirwayFound = false;
366  Vector3D closestPosition;
367  for (int i = 0; i<branches.size(); i++) //check for closer airway branch to blood vessel route
368  {
369  if ( branches[i]->findGenerationNumber() > 2 )//Needs to be deeper branch
370  {
371  Eigen::MatrixXd branchPositinos = branches[i]->getPositions();
372  for (int j = 0; j<branchPositinos.cols(); j++)
373  {
374  double distance = findDistanceFromPointToLine(branchPositinos.col(j), mBloodVesselRoutePositions).second;
375  if (minDistance > distance)
376  {
377  minDistance = distance;
378  closestPosition = branchPositinos.col(j);
379  closerAirwayFound = true;
380  mProjectedBranchPtr = branches[i];
381  mProjectedIndex = j;
382  connectionIndexBloodVesselRoute = findDistanceFromPointToLine(branchPositinos.col(j), mBloodVesselRoutePositions).first;
383  }
384  }
385  }
386  }
387 
388  if (closerAirwayFound) //calculating new route
389  {
390  findClosestPointInBranches(closestPosition);
391  findRoutePositions();
392  }
393 
394  std::vector< Eigen::Vector3d > mergedPositions;
395  mergedPositions.insert( mergedPositions.end(), mBloodVesselRoutePositions.begin(), mBloodVesselRoutePositions.begin() + connectionIndexBloodVesselRoute );
396  mergedPositions.insert( mergedPositions.end(), mRoutePositions.begin(), mRoutePositions.end() );
397 
398  //make extension from blood vessel to target
399  double extensionPointIncrement = 0.25; //mm
400  double extensionDistance = findDistance(mergedPositions.front(),mTargetPosition);
401  Eigen::Vector3d extensionVectorNormalized = ( mTargetPosition - mergedPositions.front() ) / extensionDistance;
402  int numberOfextensionPoints = int(extensionDistance / extensionPointIncrement);
403  Eigen::Vector3d extensionPointIncrementVector = extensionVectorNormalized * extensionDistance / numberOfextensionPoints;
404 
405  for (int i = 1; i<= numberOfextensionPoints; i++)
406  {
407  mergedPositions.insert(mergedPositions.begin(), mBloodVesselRoutePositions.front() + extensionPointIncrementVector*i);
408  }
409 
410  mMergedAirwayAndBloodVesselRoutePositions = mergedPositions;
411 
412  return true;
413 }
414 
415 vtkPolyDataPtr RouteToTarget::getConnectedAirwayAndBloodVesselRoute()
416 {
417  return addVTKPoints(mMergedAirwayAndBloodVesselRoutePositions);
418 }
419 
420 
421 bool RouteToTarget::checkIfRouteToTargetEndsAtEndOfLastBranch() // remove if not in use?
422 {
423  if (!mProjectedBranchPtr)
424  return false;
425 
426  if (!mProjectedIndex)
427  return false;
428 
429  if (mProjectedBranchPtr->getChildBranches().empty())
430  if (mProjectedBranchPtr->getPositions().cols()-1 == mProjectedIndex)
431  return true;
432 
433  return false;
434 }
435 
436 
437 vtkPolyDataPtr RouteToTarget::addVTKPoints(std::vector<Eigen::Vector3d> positions)
438 {
439  vtkPolyDataPtr retval = vtkPolyDataPtr::New();
440  vtkPointsPtr points = vtkPointsPtr::New();
441  vtkCellArrayPtr lines = vtkCellArrayPtr::New();
442  int numberOfPositions = positions.size();
443  for (int j = numberOfPositions - 1; j >= 0; j--)
444  {
445  points->InsertNextPoint(positions[j](0),positions[j](1),positions[j](2));
446  }
447  for (int j = 0; j < numberOfPositions-1; j++)
448  {
449  vtkIdType connection[2] = {j, j+1};
450  lines->InsertNextCell(2, connection);
451  }
452  retval->SetPoints(points);
453  retval->SetLines(lines);
454  return retval;
455 }
456 
457 
458 
459 
460 void RouteToTarget::addRouteInformationToFile(VisServicesPtr services)
461 {
462  QString RTTpath = services->patient()->getActivePatientFolder() + "/RouteToTargetInformation/";
463  QDir RTTDirectory(RTTpath);
464  if (!RTTDirectory.exists()) // Creating RouteToTargetInformation folder if it does not exist
465  RTTDirectory.mkpath(RTTpath);
466 
467  QString format = timestampSecondsFormat();
468  QString filePath = RTTpath + QDateTime::currentDateTime().toString(format) + "_RouteToTargetInformation.txt";
469 
470  QFile outfile(filePath);
471  if (outfile.open(QIODevice::ReadWrite))
472  {
473  QTextStream stream(&outfile);
474 
475  stream << "#Target position:" << endl;
476  stream << mTargetPosition(0) << " " << mTargetPosition(1) << " " << mTargetPosition(2) << " " << endl;
477  if (mProjectedBranchPtr)
478  {
479  stream << "#Route to target generations:" << endl;
480  stream << mProjectedBranchPtr->findGenerationNumber() << endl;
481  }
482 
483  stream << "#Trachea length (mm):" << endl;
484  double tracheaLength = this->getTracheaLength();
485  stream << tracheaLength << endl;
486 
487  stream << "#Route to target length - from Carina (mm):" << endl;
488  stream << calculateRouteLength(mRoutePositions) - tracheaLength << endl;
489  stream << "#Extended route to target length - from Carina (mm):" << endl;
490  stream << calculateRouteLength(mExtendedRoutePositions) - tracheaLength << endl;
491  }
492 }
493 
494 double RouteToTarget::getTracheaLength()
495 {
496  BranchPtr trachea = mBranchListPtr->getBranches()[0];
497  int numberOfPositionsInTrachea = trachea->getPositions().cols();
498  double tracheaLength = calculateRouteLength(smoothBranch(trachea, numberOfPositionsInTrachea-1, trachea->getPositions().col(numberOfPositionsInTrachea-1)));
499  return tracheaLength;
500 }
501 
502 std::vector< Eigen::Vector3d > RouteToTarget::getRoutePositions(MeshPtr route)
503 {
504  vtkPolyDataPtr centerline_r = route->getTransformedPolyDataCopy(route->get_rMd());
505 
506  std::vector< Eigen::Vector3d > routePositions;
507 
508  //Used example from getCenterlinePositions
509  int N = centerline_r->GetNumberOfPoints();
510  for(vtkIdType i = 0; i < N; i++)
511  {
512  double p[3];
513  centerline_r->GetPoint(i,p);
514  Eigen::Vector3d position;
515  position(0) = p[0]; position(1) = p[1]; position(2) = p[2];
516  routePositions.push_back(position);
517  }
518  return routePositions;
519 }
520 
521 double RouteToTarget::calculateRouteLength(std::vector< Eigen::Vector3d > route)
522 {
523  double routeLenght = 0;
524  for (int i=0; i<route.size()-1; i++)
525  {
526  double d0 = route[i+1](0) - route[i](0);
527  double d1 = route[i+1](1) - route[i](1);
528  double d2 = route[i+1](2) - route[i](2);
529 
530  routeLenght += sqrt( d0*d0 +d1*d1 + d2*d2 );
531  }
532 
533  return routeLenght;
534 }
535 
536 void RouteToTarget::setBloodVesselRadius()
537 {
538  std::vector<BranchPtr> branches = mBloodVesselBranchListPtr->getBranches();
539 
540  for (int i = 0; i < branches.size(); i++)
541  {
542  Eigen::MatrixXd positions = branches[i]->getPositions();
543  Eigen::MatrixXd orientations = branches[i]->getOrientations();
544  Eigen::VectorXd radius(positions.cols());
545 
546  for (int j = 0; j < positions.cols(); j++)
547  {
548  radius(j) = calculateBloodVesselRadius(positions.col(j), orientations.col(j));
549  }
550 
551  branches[i]->setRadius(radius);
552  }
553 
554 }
555 
556 double RouteToTarget::calculateBloodVesselRadius(Eigen::Vector3d position, Eigen::Vector3d orientation)
557 {
558  double radius = 0;
559  if (!mBloodVesselVolume)
560  return radius;
561 
562  vtkImageDataPtr bloodVesselImage = mBloodVesselVolume->getBaseVtkImageData();
563  int* dim = bloodVesselImage->GetDimensions();
564  double* spacing = bloodVesselImage->GetSpacing();
565  Transform3D dMr = mBloodVesselVolume->get_rMd().inverse();
566  Eigen::Vector3d position_r = dMr.coord(position);
567  int x = (int) boost::math::round( position_r[0]/spacing[0] );
568  int y = (int) boost::math::round( position_r[1]/spacing[1] );
569  int z = (int) boost::math::round( position_r[2]/spacing[2] );
570  Eigen::Vector3i indexVector;
571  indexVector(0) = x;
572  indexVector(1) = y;
573  indexVector(2) = z;
574 
575  Eigen::MatrixXd maxRadius(3,2);
576  Eigen::Vector3d perpendicularX = orientation.cross(Eigen::Vector3d::UnitX());
577  maxRadius(0,0) = findDistanceToSegmentationEdge(bloodVesselImage, indexVector, perpendicularX, dim, spacing, 1);
578  maxRadius(0,1) = findDistanceToSegmentationEdge(bloodVesselImage, indexVector, perpendicularX, dim, spacing, -1);
579  Eigen::Vector3d perpendicularY = orientation.cross(Eigen::Vector3d::UnitY());
580  maxRadius(1,0) = findDistanceToSegmentationEdge(bloodVesselImage, indexVector, perpendicularY, dim, spacing, 1);
581  maxRadius(1,1) = findDistanceToSegmentationEdge(bloodVesselImage, indexVector, perpendicularY, dim, spacing, -1);
582  Eigen::Vector3d perpendicularZ = orientation.cross(Eigen::Vector3d::UnitZ());
583  maxRadius(2,0) = findDistanceToSegmentationEdge(bloodVesselImage, indexVector, perpendicularZ, dim, spacing, 1);
584  maxRadius(2,1) = findDistanceToSegmentationEdge(bloodVesselImage, indexVector, perpendicularZ, dim, spacing, -1);
585 
586  radius = maxRadius.rowwise().mean().minCoeff();
587 
588  if (std::isnan(radius))
589  radius = 0;
590 
591  return radius;
592 }
593 
594 double RouteToTarget::findDistanceToSegmentationEdge(vtkImageDataPtr bloodVesselImage, Eigen::Vector3i indexVector, Eigen::Vector3d perpendicularVector, int* dim, double* spacing, int direction)
595 {
596  double retval;
597  double maxValue = bloodVesselImage->GetScalarRange()[1];
598  for (int radiusVoxels=1; radiusVoxels<30; radiusVoxels++)
599  {
600  if (perpendicularVector.sum() != 0)
601  {
602  Eigen::Vector3d searchDirection = perpendicularVector.normalized() * radiusVoxels;
603  int xIndex = std::max(std::min(indexVector(0) + direction * (int) std::round(searchDirection(0)), dim[0]-1), 0);
604  int yIndex = std::max(std::min(indexVector(1) + direction * (int) std::round(searchDirection(1)), dim[1]-1), 0);
605  int zIndex = std::max(std::min(indexVector(2) + direction * (int) std::round(searchDirection(2)), dim[2]-1), 0);
606  if (bloodVesselImage->GetScalarComponentAsDouble(xIndex, yIndex, zIndex, 0) < maxValue)
607  {
608  searchDirection = perpendicularVector.normalized() * (radiusVoxels-1);
609  retval = std::sqrt( std::pow(searchDirection(0)*spacing[0],2) + std::pow(searchDirection(1)*spacing[1],2) + std::pow(searchDirection(2)*spacing[2],2) );
610  break;
611  }
612  }
613  }
614  return retval;
615 }
616 
617 std::vector< Eigen::Vector3d > RouteToTarget::getRoutePositions()
618 {
619  std::vector< Eigen::Vector3d > positions = mExtendedRoutePositions;
620  std::reverse(positions.begin(), positions.end());
621  return positions;
622 }
623 
624 std::vector< double > RouteToTarget::getCameraRotation()
625 {
626  std::vector< double > rotations = mExtendedCameraRotation;
627  std::reverse(rotations.begin(), rotations.end());
628  return rotations;
629 }
630 
631 void RouteToTarget::makeMarianaCenterlineFile(QString filename)
632 {
633  if (mExtendedRoutePositions.empty())
634  {
635  std::cout << "mExtendedRoutePositions is empty." << std::endl;
636  return;
637  }
638 
639  int numberOfExtendedPositions = mExtendedRoutePositions.size() - mRoutePositions.size();
640 
641  ofstream out(filename.toStdString().c_str());
642  out << "# [xPos yPos zPos BranchingPoint (0=Normal, 1=Branching position, 2=Extended from airway to target)] ";
643  out << "Number of positions: ";
644  out << mExtendedRoutePositions.size(); // write number of positions
645  out << "\n";
646 
647  for (int i = 1; i < mExtendedRoutePositions.size(); i++)
648  {
649  out << mExtendedRoutePositions[i](0) << " "; // write x coordinate
650  out << mExtendedRoutePositions[i](1) << " "; // write y coordinate
651  out << mExtendedRoutePositions[i](2) << " "; // write z coordinate
652 
653  if ( std::find(mBranchingIndex.begin(), mBranchingIndex.end(), i - numberOfExtendedPositions) != mBranchingIndex.end() )
654  out << "1 ";
655  else if (i <= numberOfExtendedPositions)
656  out << "2 ";
657  else
658  out << "0 ";
659 
660  out << "\n";
661  }
662 
663  out.close();
664 }
665 
666 QJsonArray RouteToTarget::makeMarianaCenterlineJSON()
667 {
668  QJsonArray array;
669  if ( mRoutePositions.empty() || mExtendedRoutePositions.empty() )
670  {
671  std::cout << "mRoutePositions is empty." << std::endl;
672  return array;
673  }
674 
675  int numberOfExtendedPositions = mExtendedRoutePositions.size() - mRoutePositions.size();
676 
677  for (int i = 1; i < mExtendedRoutePositions.size(); i++)
678  {
679  QJsonObject position;
680  position.insert( "x", mExtendedRoutePositions[i](0) );
681  position.insert( "y", mExtendedRoutePositions[i](1) );
682  position.insert( "z", mExtendedRoutePositions[i](2) );
683 
684  if ( std::find(mBranchingIndex.begin(), mBranchingIndex.end(), i - numberOfExtendedPositions) != mBranchingIndex.end() )
685  position.insert("Flag", 1);
686  else if (i <= numberOfExtendedPositions)
687  position.insert("Flag", 2);
688  else
689  position.insert("Flag", 0);
690 
691  array.append(position);
692  }
693 
694  return array;
695 }
696 
697 QJsonArray makeMarianaCenterlineOfFullBranchTreeJSON(BranchListPtr branchList)
698 {
699  QJsonArray array;
700 
701  std::vector<BranchPtr> branches = branchList->getBranches();
702  for (int i = 0; i < branches.size(); i++)
703  {
704  Eigen::MatrixXd positions = branches[i]->getPositions();
705  for (int j = 0; j < positions.cols(); j++)
706  {
707  QJsonObject JsonPosition;
708  JsonPosition.insert( "x", positions(0,j) );
709  JsonPosition.insert( "y", positions(1,j) );
710  JsonPosition.insert( "z", positions(2,j) );
711  array.append(JsonPosition);
712  }
713  }
714 
715  return array;
716 }
717 
718 /*
719  RouteToTarget::getBranchPositions is used to get positions of a branch without smoothing.
720  Equivalent to smoothBranch without smoothing.
721 */
722 std::vector< Eigen::Vector3d > getBranchPositions(BranchPtr branchPtr, int startIndex)
723 {
724  Eigen::MatrixXd branchPositions = branchPtr->getPositions();
725  std::vector< Eigen::Vector3d > positions;
726 
727  for (int i = startIndex; i >=0; i--)
728  positions.push_back(branchPositions.col(i));
729 
730  return positions;
731 }
732 
733 
734 
735 Eigen::MatrixXd findClosestBloodVesselSegments(Eigen::MatrixXd bloodVesselPositions , Eigen::MatrixXd airwayPositions, Vector3D targetPosition)
736 {
737  double maxDistanceToAirway = 10; //mm
738  int minNumberOfPositionsInSegment = 100; //to avoid small segments which are probably not true blood vessels
739 
740  Eigen::MatrixXd bloodVesselSegment;
741 
742  while (bloodVesselPositions.cols() > minNumberOfPositionsInSegment)
743  {
744  Eigen::MatrixXd::Index closestBloodVesselPositionToTarget = dsearch(targetPosition, bloodVesselPositions).first;
745  std::pair< Eigen::MatrixXd, Eigen::MatrixXd > localPositions = findLocalPointsInCT(closestBloodVesselPositionToTarget , bloodVesselPositions);
746  bloodVesselPositions = localPositions.second;
747 
748  if ( localPositions.first.cols() >= minNumberOfPositionsInSegment &&
749  dsearchn(airwayPositions, localPositions.first).second.minCoeff() <= maxDistanceToAirway )
750  {
751  bloodVesselSegment = localPositions.first;
752  break;
753  }
754  }
755 
756  return bloodVesselSegment;
757 }
758 
759 std::pair< Eigen::MatrixXd, Eigen::MatrixXd > findLocalPointsInCT(int closestCLIndex , Eigen::MatrixXd CLpositions)
760 {
761  Eigen::MatrixXd includedPositions;
762  Eigen::MatrixXd positionsNotIncluded = CLpositions;
763  int startIndex = closestCLIndex;
764 
765  bool closePositionFound = true;
766  while (closePositionFound)
767  {
768  closePositionFound = false;
769  std::pair<Eigen::MatrixXd,Eigen::MatrixXd> connectedPoints = findConnectedPointsInCT(startIndex , positionsNotIncluded);
770  positionsNotIncluded = connectedPoints.second;
771 
772  if (includedPositions.cols() > 0)
773  {
774  includedPositions.conservativeResize(Eigen::NoChange, includedPositions.cols() + connectedPoints.first.cols());
775  includedPositions.rightCols(connectedPoints.first.cols()) = connectedPoints.first;
776  }
777  else
778  includedPositions = connectedPoints.first;
779 
780  for (int i = 0; i < includedPositions.cols(); i++)
781  {
782  std::pair<Eigen::MatrixXd::Index, double> closePositionSearch = dsearch(includedPositions.col(i), positionsNotIncluded);
783  if (closePositionSearch.second < 3) //Invlude positions closer than 3 mm
784  {
785  closePositionFound = true;
786  startIndex = closePositionSearch.first;
787  break;
788  }
789  }
790  }
791 
792  return std::make_pair(includedPositions, positionsNotIncluded);
793 }
794 
795 std::pair<int, double> findDistanceFromPointToLine(Eigen::MatrixXd point, std::vector< Eigen::Vector3d > line)
796 {
797  int index = 0;
798  double minDistance = findDistance(point, line[0]);
799  for (int i=1; i<line.size(); i++)
800  if (minDistance > findDistance(point, line[i]))
801  {
802  minDistance = findDistance(point, line[i]);
803  index = i;
804  }
805 
806  return std::make_pair(index , minDistance);
807 }
808 
809 double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
810 {
811  double d0 = p1(0) - p2(0);
812  double d1 = p1(1) - p2(1);
813  double d2 = p1(2) - p2(2);
814 
815  double D = sqrt( d0*d0 + d1*d1 + d2*d2 );
816 
817  return D;
818 }
819 
820 Eigen::MatrixXd convertToEigenMatrix(std::vector< Eigen::Vector3d > positionsVector)
821 {
822  Eigen::MatrixXd positionsMatrix(3, positionsVector.size());
823  for (int i = 0; i < positionsVector.size(); i++)
824  {
825  positionsMatrix(0, i) = positionsVector[i](0);
826  positionsMatrix(1, i) = positionsVector[i](1);
827  positionsMatrix(2, i) = positionsVector[i](2);
828  }
829  return positionsMatrix;
830 }
831 
832 double variance(Eigen::VectorXd X)
833 {
834  double mean_X = X.mean();
835  double var = 0;
836  for (int i = 0; i < X.size(); i++)
837  var += ( X[i]-mean_X ) * ( X[i]-mean_X );
838 
839  var = var/X.size();
840  return var;
841 }
842 
843 } /* namespace cx */
cx::RouteToTarget::setBloodVesselRadius
void setBloodVesselRadius()
Definition: cxRouteToTarget.cpp:536
vtkCellArrayPtr
vtkSmartPointer< class vtkCellArray > vtkCellArrayPtr
Definition: vtkForwardDeclarations.h:41
cxLogger.h
cx::RouteToTarget::generateAirwaysFromBloodVesselCenterlines
vtkPolyDataPtr generateAirwaysFromBloodVesselCenterlines()
Definition: cxRouteToTarget.cpp:326
cx::RouteToTarget::calculateBloodVesselRadius
double calculateBloodVesselRadius(Eigen::Vector3d position, Eigen::Vector3d orientation)
Definition: cxRouteToTarget.cpp:556
vtkCardinalSplinePtr
vtkSmartPointer< class vtkCardinalSpline > vtkCardinalSplinePtr
Definition: cxRouteToTarget.cpp:25
cx::RouteToTarget::setBranchList
void setBranchList(BranchListPtr branchList)
Definition: cxRouteToTarget.cpp:88
cx::RouteToTarget::findDistanceToSegmentationEdge
double findDistanceToSegmentationEdge(vtkImageDataPtr bloodVesselImage, Eigen::Vector3i indexVector, Eigen::Vector3d perpendicularVector, int *dim, double *spacing, int direction)
Definition: cxRouteToTarget.cpp:594
cx::RouteToTarget::getCenterlinePositions
Eigen::MatrixXd getCenterlinePositions(vtkPolyDataPtr centerline_r)
Definition: cxRouteToTarget.cpp:48
cx
Namespace for all CustusX production code.
Definition: cx_dev_group_definitions.h:13
cx::RouteToTarget::findRoutePositions
void findRoutePositions()
Definition: cxRouteToTarget.cpp:190
cxImage.h
cx::dsearch
std::pair< Eigen::MatrixXd::Index, double > dsearch(Eigen::Vector3d p, Eigen::MatrixXd positions)
Definition: cxBranchList.cpp:636
cx::findDistance
double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
Definition: cxBronchoscopePositionProjection.cpp:372
cx::VisServicesPtr
boost::shared_ptr< class VisServices > VisServicesPtr
Definition: cxMainWindow.h:40
CX_LOG_INFO
#define CX_LOG_INFO
Definition: cxLogger.h:96
cx::RouteToTarget::findRoutePositionsInBloodVessels
void findRoutePositionsInBloodVessels()
Definition: cxRouteToTarget.cpp:197
vtkImageDataPtr
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Definition: cxVideoConnectionWidget.h:30
cxBranchList.h
cx::RouteToTarget::findRouteToTarget
vtkPolyDataPtr findRouteToTarget(PointMetricPtr targetPoint)
Definition: cxRouteToTarget.cpp:260
cx::RouteToTarget::searchBranchUp
void searchBranchUp(BranchPtr searchBranchPtr, int startIndex)
Definition: cxRouteToTarget.cpp:210
cx::variance
double variance(Eigen::VectorXd X)
Definition: cxRouteToTarget.cpp:832
cx::RouteToTarget::addVTKPoints
vtkPolyDataPtr addVTKPoints(std::vector< Eigen::Vector3d > positions)
Definition: cxRouteToTarget.cpp:437
cx::MeshPtr
boost::shared_ptr< class Mesh > MeshPtr
Definition: cxForwardDeclarations.h:48
cx::AirwaysFromCenterline
Definition: cxAirwaysFromCenterline.h:28
cx::Transform3D
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
Definition: cxLandmarkPatientRegistrationWidget.h:33
cx::RouteToTarget::~RouteToTarget
virtual ~RouteToTarget()
Definition: cxRouteToTarget.cpp:39
cx::findLocalPointsInCT
std::pair< Eigen::MatrixXd, Eigen::MatrixXd > findLocalPointsInCT(int closestCLIndex, Eigen::MatrixXd CLpositions)
Definition: cxRouteToTarget.cpp:759
cx::vtkPointsPtr
vtkSmartPointer< vtkPoints > vtkPointsPtr
Definition: cxCenterlineRegistration.h:41
cx::smoothBranch
std::vector< Eigen::Vector3d > smoothBranch(BranchPtr branchPtr, int startIndex, Eigen::MatrixXd startPosition)
Definition: cxBranchList.cpp:705
cx::timestampSecondsFormat
QString timestampSecondsFormat()
Definition: cxTime.cpp:18
cx::convertToEigenMatrix
Eigen::MatrixXd convertToEigenMatrix(std::vector< Eigen::Vector3d > positionsVector)
Definition: cxRouteToTarget.cpp:820
cx::ImagePtr
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
cx::dsearchn
std::pair< std::vector< Eigen::MatrixXd::Index >, Eigen::VectorXd > dsearchn(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
Definition: cxBranchList.cpp:646
cx::RouteToTarget::addRouteInformationToFile
void addRouteInformationToFile(VisServicesPtr services)
Definition: cxRouteToTarget.cpp:460
cxBranch.h
cx::RouteToTarget::getCameraRotation
std::vector< double > getCameraRotation()
Definition: cxRouteToTarget.cpp:624
cx::RouteToTarget::findExtendedRoute
vtkPolyDataPtr findExtendedRoute(PointMetricPtr targetPoint)
Definition: cxRouteToTarget.cpp:273
cx::findConnectedPointsInCT
std::pair< Eigen::MatrixXd, Eigen::MatrixXd > findConnectedPointsInCT(int startIndex, Eigen::MatrixXd positionsNotUsed)
Definition: cxBranchList.cpp:661
cx::RouteToTarget::processCenterline
void processCenterline(MeshPtr mesh)
Definition: cxRouteToTarget.cpp:69
cx::findDistanceFromPointToLine
std::pair< int, double > findDistanceFromPointToLine(Eigen::MatrixXd point, std::vector< Eigen::Vector3d > line)
Definition: cxRouteToTarget.cpp:795
cx::RouteToTarget::findClosestPointInBranches
void findClosestPointInBranches(Vector3D targetCoordinate_r)
Definition: cxRouteToTarget.cpp:139
cx::vtkPolyDataPtr
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
Definition: cxCenterlineRegistration.h:42
cxPointMetric.h
cx::AirwaysFromCenterlinePtr
boost::shared_ptr< AirwaysFromCenterline > AirwaysFromCenterlinePtr
Definition: cxAirwaysFromCenterline.h:67
cx::RouteToTarget::getConnectedAirwayAndBloodVesselRoute
vtkPolyDataPtr getConnectedAirwayAndBloodVesselRoute()
Definition: cxRouteToTarget.cpp:415
cx::round
Vector3D round(const Vector3D &a)
Definition: cxVector3D.cpp:75
cx::RouteToTarget::RouteToTarget
RouteToTarget()
Definition: cxRouteToTarget.cpp:30
cx::RouteToTarget::setSmoothing
void setSmoothing(bool smoothing)
Definition: cxRouteToTarget.cpp:64
cx::BranchList
Definition: cxBranchList.h:28
cxTime.h
cx::RouteToTarget::searchBloodVesselBranchUp
void searchBloodVesselBranchUp(BranchPtr searchBranchPtr, int startIndex)
Definition: cxRouteToTarget.cpp:236
cx::BranchListPtr
boost::shared_ptr< class BranchList > BranchListPtr
Definition: cxBronchoscopePositionProjection.h:14
cx::RouteToTarget::getTracheaLength
double getTracheaLength()
Definition: cxRouteToTarget.cpp:494
cx::findClosestBloodVesselSegments
Eigen::MatrixXd findClosestBloodVesselSegments(Eigen::MatrixXd bloodVesselPositions, Eigen::MatrixXd airwayPositions, Vector3D targetPosition)
Definition: cxRouteToTarget.cpp:735
cxAirwaysFromCenterline.h
cx::makeMarianaCenterlineOfFullBranchTreeJSON
QJsonArray makeMarianaCenterlineOfFullBranchTreeJSON(BranchListPtr branchList)
Definition: cxRouteToTarget.cpp:697
cx::RouteToTarget::findRouteToTargetAlongBloodVesselCenterlines
vtkPolyDataPtr findRouteToTargetAlongBloodVesselCenterlines(MeshPtr bloodVesselCenterlineMesh, PointMetricPtr targetPoint)
Definition: cxRouteToTarget.cpp:300
cx::RouteToTarget::calculateRouteLength
static double calculateRouteLength(std::vector< Eigen::Vector3d > route)
Definition: cxRouteToTarget.cpp:521
cxRouteToTarget.h
cx::BranchPtr
boost::shared_ptr< class Branch > BranchPtr
Definition: cxBronchoscopePositionProjection.h:15
cx::RouteToTarget::getRoutePositions
std::vector< Eigen::Vector3d > getRoutePositions()
Definition: cxRouteToTarget.cpp:617
cx::PointMetricPtr
boost::shared_ptr< class PointMetric > PointMetricPtr
Definition: cxForwardDeclarations.h:84
cx::RouteToTarget::makeMarianaCenterlineJSON
QJsonArray makeMarianaCenterlineJSON()
Definition: cxRouteToTarget.cpp:666
cx::RouteToTarget::findClosestPointInBloodVesselBranches
void findClosestPointInBloodVesselBranches(Vector3D targetCoordinate_r)
Definition: cxRouteToTarget.cpp:164
cx::Vector3D
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
cx::RouteToTarget::processBloodVesselCenterline
void processBloodVesselCenterline(Eigen::MatrixXd positions)
Definition: cxRouteToTarget.cpp:93
cx::RouteToTarget::setBloodVesselVolume
void setBloodVesselVolume(ImagePtr bloodVesselVolume)
Definition: cxRouteToTarget.cpp:43
cx::RouteToTarget::makeConnectedAirwayAndBloodVesselRoute
bool makeConnectedAirwayAndBloodVesselRoute()
Definition: cxRouteToTarget.cpp:342
cx::getBranchPositions
std::vector< Eigen::Vector3d > getBranchPositions(BranchPtr branchPtr, int startIndex)
Definition: cxRouteToTarget.cpp:722
cx::RouteToTarget::makeMarianaCenterlineFile
void makeMarianaCenterlineFile(QString filename)
Definition: cxRouteToTarget.cpp:631
cxVisServices.h