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