CustusX  22.04-rc5
An IGT application
cxAirwaysFromCenterline.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
11 
13 #include <boost/math/special_functions/round.hpp>
14 #include <vtkPolyData.h>
15 #include "cxBranchList.h"
16 #include "cxBranch.h"
17 #include "vtkCardinalSpline.h"
18 #include <cxImage.h>
19 #include "cxContourFilter.h"
20 #include <vtkImageData.h>
21 #include <vtkPointData.h>
22 #include "cxVolumeHelpers.h"
23 #include "vtkCardinalSpline.h"
24 #include "cxLogger.h"
25 #include <vtkImageResample.h>
26 
27 typedef vtkSmartPointer<class vtkCardinalSpline> vtkCardinalSplinePtr;
28 
29 
30 namespace cx
31 {
32 
34  mBranchListPtr(new BranchList),
35  mAirwaysVolumeBoundaryExtention(10),
36  mAirwaysVolumeBoundaryExtentionTracheaStart(2),
37  mAirwaysVolumeSpacing(0.5)
38 {
39 }
40 
42 {
43 }
44 
46 {
47  mBloodVessel = bloodVessel;
48 }
49 
51 {
52 
53  int N = centerline_r->GetNumberOfPoints();
54  Eigen::MatrixXd CLpoints(3,N);
55  for(vtkIdType i = 0; i < N; i++)
56  {
57  double p[3];
58  centerline_r->GetPoint(i,p);
59  Eigen::Vector3d position;
60  position(0) = p[0]; position(1) = p[1]; position(2) = p[2];
61  CLpoints.block(0 , i , 3 , 1) = position;
62  }
63  return CLpoints;
64 }
65 
67 {
68  mBranchListPtr = branches;
69 }
70 
72 {
73  mOriginalSegmentedVolume = segmentedVolume;
74 }
75 
77 {
78  if (mBranchListPtr)
79  mBranchListPtr->deleteAllBranches();
80 
81  Eigen::MatrixXd CLpoints_r = getCenterlinePositions(centerline_r);
82 
83  mBranchListPtr->findBranchesInCenterline(CLpoints_r);
84 
85  mBranchListPtr->smoothBranchPositions(40);
86  mBranchListPtr->interpolateBranchPositions(0.1);
87  this->smoothAllBranchesForVB();
88 
89  mBranchListPtr->smoothOrientations();
90 }
91 
93 {
94  return mBranchListPtr;
95 }
96 
97 /*
98  AirwaysFromCenterline::generateTubes makes artificial airway tubes around the input centerline. The radius
99  of the tubes is decided by the generation number, based on Weibel's model of airways. In contradiction to the model,
100  it is set a lower boundary for the tube radius (2 mm) making the peripheral airways larger than in reality,
101  which makes it possible to virtually navigate inside the tubes. The airways are generated by adding a sphere to
102  a volume (image) at each point along every branch. The output is a surface model generated from the volume.
103 */
104 vtkPolyDataPtr AirwaysFromCenterline::generateTubes(double staticRadius, bool mergeWithOriginalAirways) // if staticRadius == 0, radius is retrieved from branch generation number
105 {
106  mMergeWithOriginalAirways = mergeWithOriginalAirways;
107  vtkImageDataPtr airwaysVolumePtr;
108 
109  if (mergeWithOriginalAirways)
110  {
111  if (mOriginalSegmentedVolume)
112  {
113  airwaysVolumePtr = this->initializeAirwaysVolumeFromOriginalSegmentation();
114  }
115  else
116  {
117  CX_LOG_WARNING() << "AirwaysFromCenterline::generateTubes: Segmented airways volume not set. Creating pure artificaial tubes around centerlines.";
118  airwaysVolumePtr = this->initializeEmptyAirwaysVolume();
119  }
120  }
121  else
122  airwaysVolumePtr = this->initializeEmptyAirwaysVolume();
123 
124  airwaysVolumePtr = addSpheresAlongCenterlines(airwaysVolumePtr, staticRadius);
125 
126  //create contour from image
128  airwaysVolumePtr,
129  1, //treshold
130  false, // reduce resolution
131  true, // smoothing
132  true, // keep topology
133  0, // target decimation
134  30, // number of iterations smoothing
135  0.10 // band pass smoothing
136  );
137 
138  return rawContour;
139 }
140 
142 {
143  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
144  vtkPointsPtr pointsPtr = vtkPointsPtr::New();
145 
146  int numberOfPoints = 0;
147  for (int i = 0; i < branches.size(); i++)
148  numberOfPoints += branches[i]->getPositions().cols();
149 
150  pointsPtr->SetNumberOfPoints(numberOfPoints);
151 
152  int pointIndex = 0;
153  for (int i = 0; i < branches.size(); i++)
154  {
155  Eigen::MatrixXd positions = branches[i]->getPositions();
156  for (int j = 0; j < positions.cols(); j++)
157  {
158  pointsPtr->SetPoint(pointIndex, positions(0,j), positions(1,j), positions(2,j));
159  pointIndex += 1;
160  }
161  }
162 
163  pointsPtr->GetBounds(mBounds);
164 
165  //Extend bounds to make room for surface model extended from centerline
166  mBounds[0] -= mAirwaysVolumeBoundaryExtention;
167  mBounds[1] += mAirwaysVolumeBoundaryExtention;
168  mBounds[2] -= mAirwaysVolumeBoundaryExtention;
169  mBounds[3] += mAirwaysVolumeBoundaryExtention;
170  mBounds[4] -= mAirwaysVolumeBoundaryExtention;
171  mBounds[5] -= mAirwaysVolumeBoundaryExtentionTracheaStart; // to make top of trachea open
172  if (mBloodVessel)
173  mBounds[5] += mAirwaysVolumeBoundaryExtention;
174 
175  mSpacing[0] = mAirwaysVolumeSpacing; //Smaller spacing improves resolution but increases run-time
176  mSpacing[1] = mAirwaysVolumeSpacing;
177  mSpacing[2] = mAirwaysVolumeSpacing;
178 
179  // compute dimensions
180  for (int i = 0; i < 3; i++)
181  mDim[i] = static_cast<int>(std::ceil((mBounds[i * 2 + 1] - mBounds[i * 2]) / mSpacing[i]));
182 
183  mOrigin[0] = mBounds[0] + mSpacing[0] / 2;
184  mOrigin[1] = mBounds[2] + mSpacing[1] / 2;
185  mOrigin[2] = mBounds[4] + mSpacing[2] / 2;
186 
187  vtkImageDataPtr airwaysVolumePtr = generateVtkImageData(mDim, mSpacing, 0);
188  airwaysVolumePtr->SetOrigin(mOrigin);
189 
190  return airwaysVolumePtr;
191 }
192 
194 {
195  vtkImageDataPtr airwaysVolumePtr;
196  if (!mOriginalSegmentedVolume)
197  return airwaysVolumePtr;
198 
199  double magnificationFactor = mOriginalSegmentedVolume->GetSpacing()[0] / mAirwaysVolumeSpacing;
200  vtkImageResamplePtr resampler = vtkImageResamplePtr::New();
201  resampler->SetInterpolationModeToLinear();
202  resampler->SetAxisMagnificationFactor(0, magnificationFactor);
203  resampler->SetAxisMagnificationFactor(1, magnificationFactor);
204  resampler->SetAxisMagnificationFactor(2, magnificationFactor);
205  resampler->SetInputData(mOriginalSegmentedVolume);
206  resampler->Update();
207  airwaysVolumePtr = resampler->GetOutput();
208 
209  Vector3D origin(airwaysVolumePtr->GetOrigin());
210  mOrigin[0] = origin[0];
211  mOrigin[1] = origin[1];
212  mOrigin[2] = origin[2];
213 
214  Vector3D spacing(airwaysVolumePtr->GetSpacing());
215  mSpacing = spacing;
216 
217  airwaysVolumePtr->GetBounds(mBounds);
218 
219  // compute dimensions
220  for (int i = 0; i < 3; i++)
221  mDim[i] = static_cast<int>(std::ceil((mBounds[i * 2 + 1] - mBounds[i * 2]) / mSpacing[i]));
222 
223  return airwaysVolumePtr;
224 }
225 
226 
228 {
229  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
230 
231  for (int i = 0; i < branches.size(); i++)
232  {
233  Eigen::MatrixXd positions = branches[i]->getPositions();
234  vtkPointsPtr pointsPtr = vtkPointsPtr::New();
235  int numberOfPositionsInBranch = positions.cols();
236  pointsPtr->SetNumberOfPoints(numberOfPositionsInBranch);
237 
238  double radius = staticRadius;
239  if (similar(staticRadius, 0))
240  {
241  radius = branches[i]->findBranchRadius();
242  if (mMergeWithOriginalAirways)
243  radius = radius/2;
244  }
245 
246  for (int j = 0; j < numberOfPositionsInBranch; j++)
247  {
248  double spherePos[3];
249  spherePos[0] = positions(0,j);
250  spherePos[1] = positions(1,j);
251  spherePos[2] = positions(2,j);
252  airwaysVolumePtr = addSphereToImage(airwaysVolumePtr, spherePos, radius);
253  }
254  }
255  return airwaysVolumePtr;
256 }
257 
258 vtkImageDataPtr AirwaysFromCenterline::addSphereToImage(vtkImageDataPtr airwaysVolumePtr, double position[3], double radius)
259 {
260  int value = 1;
261  int centerIndex[3];
262  int sphereBoundingBoxIndex[6];
263 
264  for (int i=0; i<3; i++)
265  {
266  centerIndex[i] = static_cast<int>(boost::math::round( (position[i]-mOrigin[i]) / mSpacing[i] ));
267  sphereBoundingBoxIndex[2*i] = std::max(
268  static_cast<int>(boost::math::round( (position[i]-mOrigin[i] - radius) / mSpacing[i] )),
269  0);
270  sphereBoundingBoxIndex[2*i+1] = std::min(
271  static_cast<int>(boost::math::round( (position[i]-mOrigin[i] + radius) / mSpacing[i] )),
272  mDim[i]-1);
273  }
274 
275  for (int x = sphereBoundingBoxIndex[0]; x<=sphereBoundingBoxIndex[1]; x++)
276  for (int y = sphereBoundingBoxIndex[2]; y<=sphereBoundingBoxIndex[3]; y++)
277  for (int z = sphereBoundingBoxIndex[4]; z<=sphereBoundingBoxIndex[5]; z++)
278  {
279  double distanceFromCenter = sqrt((x-centerIndex[0])*mSpacing[0]*(x-centerIndex[0])*mSpacing[0] +
280  (y-centerIndex[1])*mSpacing[1]*(y-centerIndex[1])*mSpacing[1] +
281  (z-centerIndex[2])*mSpacing[2]*(z-centerIndex[2])*mSpacing[2]);
282 
283  if (distanceFromCenter < radius)
284  {
285  unsigned char* dataPtrImage = static_cast<unsigned char*>(airwaysVolumePtr->GetScalarPointer(x,y,z));
286  dataPtrImage[0] = value;
287  }
288  }
289 
290  return airwaysVolumePtr;
291 }
292 
294 {
295 
296  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
297  for (int i=0; i<branches.size(); i++)
298  {
299  Eigen::MatrixXd positions = branches[i]->getPositions();
300  std::vector< Eigen::Vector3d > smoothedPositions = smoothBranch(branches[i], positions.cols()-1, positions.col(positions.cols()-1));
301  for (int j=0; j<smoothedPositions.size(); j++)
302  {
303  positions(0,j) = smoothedPositions[smoothedPositions.size() - j - 1](0);
304  positions(1,j) = smoothedPositions[smoothedPositions.size() - j - 1](1);
305  positions(2,j) = smoothedPositions[smoothedPositions.size() - j - 1](2);
306  }
307  branches[i]->setPositions(positions);
308  }
309 }
310 
312 {
313  vtkPolyDataPtr retval = vtkPolyDataPtr::New();
314  vtkPointsPtr points = vtkPointsPtr::New();
315  vtkCellArrayPtr lines = vtkCellArrayPtr::New();
316 
317  if (!mBranchListPtr)
318  return retval;
319 
320  double minPointDistance = 0.5; //mm
321  mBranchListPtr->excludeClosePositionsInCTCenterline(minPointDistance); // to reduce number of positions in smoothed centerline
322 
323  std::vector<BranchPtr> branches = mBranchListPtr->getBranches();
324  int pointIndex = 0;
325 
326  for (int i = 0; i < branches.size(); i++)
327  {
328  Eigen::MatrixXd positions = branches[i]->getPositions();
329  int numberOfPositions = positions.cols();
330 
331  if (branches[i]->getParentBranch()) // Add parents last position to get connected centerlines
332  {
333  Eigen::MatrixXd parentPositions = branches[i]->getParentBranch()->getPositions();
334  points->InsertNextPoint(parentPositions(0,parentPositions.cols()-1),parentPositions(1,parentPositions.cols()-1),parentPositions(2,parentPositions.cols()-1));
335  pointIndex += 1;
336  }
337 
338  for (int j = 0; j < numberOfPositions; j++)
339  {
340  points->InsertNextPoint(positions(0,j),positions(1,j),positions(2,j));
341 
342  if (j>1 || branches[i]->getParentBranch())
343  {
344  vtkIdType connection[2] = {pointIndex-1, pointIndex};
345  lines->InsertNextCell(2, connection);
346  }
347  pointIndex += 1;
348  }
349  }
350 
351  retval->SetPoints(points);
352  retval->SetLines(lines);
353  return retval;
354 }
355 
356 std::pair<int, double> findDistanceToLine(Eigen::Vector3d point, Eigen::MatrixXd line)
357 {
358  int index = 0;
359  double minDistance = findDistance(point, line.col(0));
360  for (int i=1; i<line.cols(); i++)
361  if (minDistance > findDistance(point, line.col(i)))
362  {
363  minDistance = findDistance(point, line.col(i));
364  index = i;
365  }
366 
367  return std::make_pair(index , minDistance);
368 }
369 
370 double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
371 {
372  double d0 = p1(0) - p2(0);
373  double d1 = p1(1) - p2(1);
374  double d2 = p1(2) - p2(2);
375 
376  double D = sqrt( d0*d0 + d1*d1 + d2*d2 );
377 
378  return D;
379 }
380 
381 } /* namespace cx */
std::vector< Eigen::Vector3d > smoothBranch(BranchPtr branchPtr, int startIndex, Eigen::MatrixXd startPosition)
vtkImageDataPtr initializeAirwaysVolumeFromOriginalSegmentation()
Vector3D ceil(const Vector3D &a)
Definition: cxVector3D.cpp:84
vtkSmartPointer< class vtkCellArray > vtkCellArrayPtr
boost::shared_ptr< class BranchList > BranchListPtr
vtkImageDataPtr addSphereToImage(vtkImageDataPtr airwaysVolumePtr, double position[3], double radius)
void setSegmentedVolume(vtkImageDataPtr segmentedVolume)
vtkSmartPointer< vtkPoints > vtkPointsPtr
Eigen::MatrixXd getCenterlinePositions(vtkPolyDataPtr centerline_r)
void processCenterline(vtkPolyDataPtr centerline_r)
void setTypeToBloodVessel(bool bloodVessel)
double findDistance(Eigen::MatrixXd p1, Eigen::MatrixXd p2)
vtkImageDataPtr initializeEmptyAirwaysVolume()
vtkPolyDataPtr generateTubes(double staticRadius=0, bool mergeWithOriginalAirways=false)
virtual bool execute()
std::pair< int, double > findDistanceToLine(Eigen::Vector3d point, Eigen::MatrixXd line)
vtkImageDataPtr generateVtkImageData(Eigen::Array3i dim, Vector3D spacing, const unsigned char initValue, int components)
vtkSmartPointer< class vtkImageResample > vtkImageResamplePtr
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
void setBranches(BranchListPtr branches)
bool similar(const CameraInfo &lhs, const CameraInfo &rhs, double tol)
#define CX_LOG_WARNING
Definition: cxLogger.h:98
vtkSmartPointer< class vtkCardinalSpline > vtkCardinalSplinePtr
vtkImageDataPtr addSpheresAlongCenterlines(vtkImageDataPtr airwaysVolumePtr, double staticRadius=0)
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Vector3D round(const Vector3D &a)
Definition: cxVector3D.cpp:75
Namespace for all CustusX production code.