Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
cxAccusurf.cpp
Go to the documentation of this file.
1 
2 
3 #include "cxAccusurf.h"
4 #include <vtkPolyData.h>
5 #include <vtkSmartPointer.h>
6 #include <vtkImageData.h>
7 //#include "cxDoubleDataAdapterXml.h"
8 #include <vtkCellArray.h>
9 #include "vtkCardinalSpline.h"
10 #include "cxVolumeHelpers.h"
11 #include <boost/math/special_functions/round.hpp>
12 #include "cxReporter.h"
13 #include "cxBoundingBox3D.h"
14 #include "cxImageAlgorithms.h"
15 #include "cxImage.h"
16 
17 typedef vtkSmartPointer<class vtkCardinalSpline> vtkCardinalSplinePtr;
18 
19 namespace cx
20 {
21 
23 {
24 }
25 
26 
28 {
29 }
30 
32 {
33  mRoutePositions.clear();
34  int N = route->GetNumberOfPoints();
35  for(vtkIdType i = 0; i < N; i++)
36  {
37  double p[3];
38  route->GetPoint(i,p);
39  Eigen::Vector3d position;
40  position(0) = p[0]; position(1) = p[1]; position(2) = p[2];
41  mRoutePositions.push_back(position);
42  }
43  smoothPositions();
44 }
45 
47 {
48  mInputImage = inputImage->getBaseVtkImageData();
49  mVtkScalarType = mInputImage->GetScalarType();
50  mMinVoxelValue = inputImage->getMin();
51 }
52 
53 void Accusurf::setThickness(int thicknessUp, int thicknessDown)
54 {
55  mThicknessUp = thicknessUp;
56  mThicknessDown = thicknessDown;
57 }
58 
59 vtkImageDataPtr Accusurf::createNewEmptyImage()
60 {
61  vtkImageDataPtr newImage = vtkImageDataPtr::New();
62  newImage->DeepCopy(mInputImage);
63 
64  switch (mVtkScalarType)
65  {
66  case VTK_CHAR:
67  newImage->AllocateScalars(VTK_CHAR, 1);
68  this->insertValuesAtInitialization(static_cast<char*> (newImage->GetScalarPointer()), newImage);
69  break;
70  case VTK_UNSIGNED_CHAR:
71  newImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
72  this->insertValuesAtInitialization(static_cast<unsigned char*> (newImage->GetScalarPointer()), newImage);
73  break;
74  case VTK_SIGNED_CHAR:
75  newImage->AllocateScalars(VTK_SIGNED_CHAR, 1);
76  this->insertValuesAtInitialization(static_cast<signed char*> (newImage->GetScalarPointer()), newImage);
77  break;
78  case VTK_UNSIGNED_SHORT:
79  newImage->AllocateScalars(VTK_UNSIGNED_SHORT, 1);
80  this->insertValuesAtInitialization(static_cast<unsigned short*> (newImage->GetScalarPointer()), newImage);
81  break;
82  case VTK_SHORT:
83  newImage->AllocateScalars(VTK_SHORT, 1);
84  this->insertValuesAtInitialization(static_cast<short*> (newImage->GetScalarPointer()), newImage);
85  break;
86  case VTK_UNSIGNED_INT:
87  newImage->AllocateScalars(VTK_UNSIGNED_INT, 1);
88  this->insertValuesAtInitialization(static_cast<unsigned int*> (newImage->GetScalarPointer()), newImage);
89  break;
90  case VTK_INT:
91  newImage->AllocateScalars(VTK_INT, 1);
92  this->insertValuesAtInitialization(static_cast<int*> (newImage->GetScalarPointer()), newImage);
93  break;
94  default:
95  reportError(QString("Unknown VTK ScalarType: %1").arg(mVtkScalarType));
96  return newImage;
97  break;
98  }
99 
100  setDeepModified(newImage);
101  return newImage;
102 }
103 
105 {
106  vtkImageDataPtr accusurfImage = createNewEmptyImage();
107  int* dim = accusurfImage->GetDimensions();
108  double* spacing = accusurfImage->GetSpacing();
109  std::vector<int> yIndexes (dim[2],0);
110 
111  for (int i = 0; i<mRoutePositions.size(); i++)
112  {
113  int y = (int) boost::math::round(mRoutePositions[i](1)/spacing[1]);
114  int z = (int) boost::math::round(mRoutePositions[i](2)/spacing[2]);
115  //std::cout << "y: " << y << " - z: " << z << " - size YIndexes: " << yIndexes.size() << std::endl;
116  if (z >= 0 && z < dim[2] && y >= 0 && y < dim[1])
117  yIndexes[z] = y;
118  }
119 
120  for (int i = 1; i<dim[2]; i++)
121  {
122  if (yIndexes[i] == 0)
123  yIndexes[i] = yIndexes[i-1];
124  }
125  for (int i = dim[2]-2; i>=0; i--)
126  {
127  if (yIndexes[i] == 0)
128  yIndexes[i] = yIndexes[i+1];
129  }
130 
131 
132  switch (mVtkScalarType)
133  {
134  case VTK_CHAR:
135  this->insertValuesFromOriginalImage<char>(accusurfImage, dim, yIndexes);
136  break;
137  case VTK_UNSIGNED_CHAR:
138  this->insertValuesFromOriginalImage<unsigned char>(accusurfImage, dim, yIndexes);
139  break;
140  case VTK_SIGNED_CHAR:
141  this->insertValuesFromOriginalImage<signed char>(accusurfImage, dim, yIndexes);
142  break;
143  case VTK_UNSIGNED_SHORT:
144  insertValuesFromOriginalImage<unsigned short>(accusurfImage, dim, yIndexes);
145  break;
146  case VTK_SHORT:
147  this->insertValuesFromOriginalImage<short>(accusurfImage, dim, yIndexes);
148  break;
149  case VTK_UNSIGNED_INT:
150  this->insertValuesFromOriginalImage<unsigned int>(accusurfImage, dim, yIndexes);
151  break;
152  case VTK_INT:
153  this->insertValuesFromOriginalImage<int>(accusurfImage, dim, yIndexes);
154  break;
155  default:
156  reportError(QString("Unknown VTK ScalarType: %1").arg(mVtkScalarType));
157  break;
158  }
159 
160  int yMin = *std::min_element(yIndexes.begin(), yIndexes.end());
161  int yMax = *std::max_element(yIndexes.begin(), yIndexes.end());
162  accusurfImage = crop(accusurfImage, yMin, yMax);
163  setDeepModified(accusurfImage);
164  return accusurfImage;
165 }
166 
167 void Accusurf::smoothPositions()
168 {
169  int numberOfInputPoints = mRoutePositions.size();
170  int controlPointFactor = 10;
171  int numberOfControlPoints = numberOfInputPoints / controlPointFactor;
172 
173  vtkCardinalSplinePtr splineX = vtkSmartPointer<vtkCardinalSpline>::New();
174  vtkCardinalSplinePtr splineY = vtkSmartPointer<vtkCardinalSpline>::New();
175  vtkCardinalSplinePtr splineZ = vtkSmartPointer<vtkCardinalSpline>::New();
176 
177  if (numberOfControlPoints >= 2)
178  {
179  //add control points to spline
180  for(int j=0; j<numberOfControlPoints; j++)
181  {
182  int indexP = (j*numberOfInputPoints)/numberOfControlPoints;
183  splineX->AddPoint(indexP,mRoutePositions[indexP](0));
184  splineY->AddPoint(indexP,mRoutePositions[indexP](1));
185  splineZ->AddPoint(indexP,mRoutePositions[indexP](2));
186  }
187  //Always add the last point to complete spline
188  splineX->AddPoint(numberOfInputPoints-1,mRoutePositions[numberOfInputPoints-1](0));
189  splineY->AddPoint(numberOfInputPoints-1,mRoutePositions[numberOfInputPoints-1](1));
190  splineZ->AddPoint(numberOfInputPoints-1,mRoutePositions[numberOfInputPoints-1](2));
191 
192  //evaluate spline - get smoothed positions
193  std::vector< Eigen::Vector3d > smoothingResult;
194  for(int j=0; j<numberOfInputPoints; j++)
195  {
196  double splineParameter = j;
197  Eigen::Vector3d tempPoint;
198  tempPoint(0) = splineX->Evaluate(splineParameter);
199  tempPoint(1) = splineY->Evaluate(splineParameter);
200  tempPoint(2) = splineZ->Evaluate(splineParameter);
201  smoothingResult.push_back(tempPoint);
202  }
203  mRoutePositions.clear();
204  mRoutePositions = smoothingResult;
205  }
206 }
207 
208 vtkImageDataPtr Accusurf::crop(vtkImageDataPtr image, int ymin, int ymax){
209 
210  int* dim = image->GetDimensions();
211  IntBoundingBox3D cropbox(0, dim[0], std::max(0, ymin-mThicknessUp-2) , std::min(dim[1]-1, ymax+mThicknessDown+2) , 0, dim[2]);
212  vtkImageDataPtr cropedImage = cropImage(image, cropbox);
213  return cropedImage;
214 }
215 
216 template <class TYPE>
217 void Accusurf::insertValuesAtInitialization(TYPE* dataPtr, vtkImageDataPtr image)
218 {
219  int* dim = image->GetDimensions();
220  int numVoxels = dim[0]*dim[1]*dim[2];
221  for (int i = 0; i < numVoxels; ++i)
222  {
223  dataPtr[i] = mMinVoxelValue;
224  }
225 }
226 
227 template <class TYPE>
228 void Accusurf::insertValuesFromOriginalImage(vtkImageDataPtr image, int* dim, std::vector<int> yIndexes)
229 {
230  for (int z = 0; z<dim[2]; z++)
231  {
232  for (int x = 0; x<dim[0]; x++)
233  {
234  int ymin = std::max(yIndexes[z]-mThicknessUp,0);
235  int ymax = std::min(yIndexes[z]+mThicknessDown,dim[1]-1);
236  for (int y = ymin; y<=ymax; y++)
237  {
238 
239  TYPE* dataPtrInputImage = static_cast<TYPE*> (mInputImage->GetScalarPointer(x,y,z));
240  TYPE* dataPtrAccusurfImage = static_cast<TYPE*> (image->GetScalarPointer(x,y,z));
241  dataPtrAccusurfImage[0] = dataPtrInputImage[0];
242  }
243  }
244  }
245 }
246 
247 } /* namespace cx */
void reportError(QString msg)
Definition: cxLogger.cpp:71
vtkImageDataPtr cropImage(vtkImageDataPtr input, IntBoundingBox3D cropbox)
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
void setInputImage(ImagePtr inputImage)
Definition: cxAccusurf.cpp:46
vtkSmartPointer< class vtkCardinalSpline > vtkCardinalSplinePtr
Definition: cxAccusurf.cpp:17
void setRoutePositions(vtkPolyDataPtr route)
Definition: cxAccusurf.cpp:31
virtual ~Accusurf()
Definition: cxAccusurf.cpp:27
Representation of an integer bounding box in 3D. The data are stored as {xmin,xmax,ymin,ymax,zmin,zmax}, in order to simplify communication with vtk.
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
void setThickness(int thicknessUp, int thicknessDown)
Definition: cxAccusurf.cpp:53
void setDeepModified(vtkImageDataPtr image)
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Vector3D round(const Vector3D &a)
Definition: cxVector3D.cpp:75
vtkImageDataPtr createAccusurfImage()
Definition: cxAccusurf.cpp:104
Namespace for all CustusX production code.