Fraxinus  22.04-rc5
An IGT application
AngleCorrection.cpp
Go to the documentation of this file.
1 #include "AngleCorrection.h"
2 
3 #include "spline3d.hpp"
4 #include "ErrorHandler.hpp"
5 #include <vtkDoubleArray.h>
6 #include <vtkSmartPointer.h>
7 
8 #include <vtkPolyDataReader.h>
9 #include <vtkPolyData.h>
10 #include <vtkPolyDataWriter.h>
11 #include <vtkPointData.h>
12 
13 using namespace std;
14 
15 
25  mOutput = NULL;
26  mValidInput= false;
27  mClSplinesPtr = new vectorSpline3dDouble();
28  mVelDataPtr = new vector<MetaImage<inData_t>>();
29  mClData=vtkSmartPointer<vtkPolyData>::New();
30  mVelImagePrefix="";
31  mIntersections = 0;
32  mBloodVessels = 0;
33  mNumOfStepsRan=0;
34  mVnyq=0;
35  mCutoff=0;
36  mnConvolutions=0;
37  mUncertainty_limit=0;
38  mMinArrowDist=0;
39  mBloodVesselsRemoved = 0;
40  mUpdate1=true;
41  mUpdate2=true;
42 }
43 
44 
46  mClSplinesPtr->clear();
47  delete mClSplinesPtr;
48 
49  mVelDataPtr->clear();
50  delete mVelDataPtr;
51 
52  mOutput = NULL;
53  mValidInput= false;
54  mClSplinesPtr = new vectorSpline3dDouble();
55  mVelDataPtr = new vector<MetaImage<inData_t>>();
56  mClData=vtkSmartPointer<vtkPolyData>::New();
57  mVelImagePrefix="";
58  mIntersections = 0;
59  mBloodVessels = 0;
60  mNumOfStepsRan=0;
61 }
62 
74 void AngleCorrection::setInput(vtkSmartPointer<vtkPolyData> vpd_centerline, vector<MetaImage<inData_t> >* velData, double Vnyq, double cutoff, int nConvolutions, double uncertainty_limit, double minArrowDist)
75 {
76  mValidInput= false;
77  if (uncertainty_limit < 0.0) reportError("ERROR: uncertainty_limit must be positive ");
78  if (minArrowDist < 0.0) reportError("ERROR: minArrowDist must be positive ");
79  if (Vnyq < 0.0) reportError("ERROR: vNyquist must be positive ");
80  if (nConvolutions < 0.0) reportError("ERROR: nConvolutions must be positive ");
81  if (vpd_centerline->GetNumberOfPoints()<=0) reportError("ERROR: No points found in the center line ");
82  if (vpd_centerline->GetNumberOfLines()<=0) reportError("ERROR: No lines found in the center line, the center line must be a linked list ");
83 
84  if((velData->size() > 0 && mVelDataPtr!=velData) || mVelImagePrefix.size()==0)
85  {
86  mVelDataPtr->clear();
87  mVelDataPtr=velData;
88  mUpdate1=true;
89  }
90 
91  if(mVnyq!=Vnyq ||
92  mCutoff!=cutoff ||
93  mnConvolutions!=nConvolutions ||
94  !EqualVtkPolyData(mClData,vpd_centerline))
95  {
96  mClData->DeepCopy(vpd_centerline);
97  mVnyq=Vnyq;
98  mCutoff=cutoff;
99  mnConvolutions=nConvolutions;
100  mUpdate1=true;
101  }
102 
103  if(mUncertainty_limit!=uncertainty_limit ||
104  mMinArrowDist!=minArrowDist)
105  {
106  mUncertainty_limit=uncertainty_limit;
107  mMinArrowDist=minArrowDist;
108  mUpdate2=true;
109  }
110  mValidInput= true;
111 }
112 
113 
114 void AngleCorrection::setInput(vtkSmartPointer<vtkPolyData> vpd_centerline, const char* velImagePrefix , double Vnyq, double cutoff, int nConvolutions, double uncertainty_limit, double minArrowDist)
115 {
116  if(mVelImagePrefix!=std::string(velImagePrefix))
117  {
118  mVelImagePrefix=std::string(velImagePrefix);
119  mVelDataPtr->clear();
120  mUpdate1=true;
121  }
122 
123  std::string filename=std::string(velImagePrefix);
124  filename.append("0.mhd");
125 
126  vtkSmartPointer<vtkMetaImageReader> reader= vtkSmartPointer<vtkMetaImageReader>::New();
127  vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
128  reader->AddObserver(vtkCommand::ErrorEvent,errorObserver);
129 
130  reader->SetFileName(filename.c_str());
131  reader->Update();
132 
133  if(!reader->CanReadFile(filename.c_str())){
134  reportError("ERROR: Could not read velocity data \n");
135  }
136 
137  setInput(vpd_centerline, new vector<MetaImage<inData_t>>(), Vnyq, cutoff, nConvolutions, uncertainty_limit, minArrowDist);
138 }
139 
140 
141 void AngleCorrection::setInput(const char* centerline,const char* image_prefix, double Vnyq, double cutoff, int nConvolutions, double uncertainty_limit, double minArrowDist)
142 {
143 
144  cerr << "Input params: " << centerline<< " " << image_prefix << " " << Vnyq << " " << nConvolutions << " " << uncertainty_limit<< " " << minArrowDist <<endl;
145 
146  vtkSmartPointer<vtkPolyDataReader> clReader = vtkSmartPointer<vtkPolyDataReader>::New();
147 
148  vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
149  clReader->AddObserver(vtkCommand::ErrorEvent,errorObserver);
150 
151  clReader->SetFileName(centerline);
152  clReader->Update();
153  if (errorObserver->GetError())
154  {
155  reportError("ERROR: Could not read center line data \n"+ errorObserver->GetErrorMessage());
156 
157  }
158  if (errorObserver->GetWarning()){
159  cerr << "Caught warning while reading center line data \n! " << errorObserver->GetWarningMessage();
160  }
161 
162  vtkSmartPointer<vtkPolyData> vpd_centerline = clReader->GetOutput();
163 
164  setInput(vpd_centerline, image_prefix, Vnyq, cutoff, nConvolutions, uncertainty_limit, minArrowDist);
165 }
166 
167 
168 
169 
171 {
172  cerr << "params: " << mnConvolutions<< " " << mCutoff << " " << mUncertainty_limit << " " << mMinArrowDist << " " << mVnyq <<endl;
173 
174  if(!mValidInput)
175  {
176  cerr << "Invalid input " << endl;
177  mOutput = NULL;
178  return false;
179  }
180  mValidInput=false;
181 
182  if(mVelDataPtr->size() == 0)
183  {
184  cerr << "Loading data " << endl;
185  mVelDataPtr->clear();
186  mVelDataPtr = MetaImage<inData_t>::readImages(mVelImagePrefix);
187  }
188 
189  mNumOfStepsRan=0;
190  if(mUpdate1)
191  {
192  mNumOfStepsRan++;
193  cerr << "started step 1 of 2 "<< endl;
194  angle_correction_impl(mClData, mVelDataPtr, mVnyq, mCutoff, mnConvolutions);
195  }
196 
197  if(mUpdate1 || mUpdate2)
198  {
199  mNumOfStepsRan++;
200  cerr << "started step 2 of 2" << endl;
201  mOutput= computeVtkPolyData(mClSplinesPtr, mUncertainty_limit, mMinArrowDist);
202  }
203 
204  mUpdate1=false;
205  mUpdate2=false;
206  return true;
207 }
208 
209 vtkSmartPointer<vtkPolyData> AngleCorrection::getOutput()
210 {
211  return mOutput;
212 }
213 
214 
216 {
217  return *mClSplinesPtr;
218 }
219 
220 
221 
223 {
224 
225  vtkSmartPointer<vtkPolyData> polydata = getOutput();
226  if(!polydata) polydata=vtkSmartPointer<vtkPolyData>::New();
227 
228  vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
229  vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
230  writer->AddObserver(vtkCommand::ErrorEvent,errorObserver);
231 
232  writer->SetFileName(filename);
233 #if VTK_MAJOR_VERSION <= 5
234  writer->SetInput(polydata);
235 #else
236  writer->SetInputData(polydata);
237 #endif
238  writer->Write();
239 
240  if (errorObserver->GetError())
241  {
242  reportError("ERROR: Could not write file to disk \n"+ errorObserver->GetErrorMessage());
243  }
244  if (errorObserver->GetWarning()){
245  cerr << "Caught warning while not writing file to disk! \n " << errorObserver->GetWarningMessage();
246  }
247 }
248 
249 void AngleCorrection::angle_correction_impl(vtkSmartPointer<vtkPolyData> vpd_centerline, vector<MetaImage<inData_t> >* images , double Vnyq, double cutoff, int nConvolutions)
250 {
251  bool verbose = false;
252 
253  mClSplinesPtr->clear();
254  mClSplinesPtr = Spline3D<double>::build(vpd_centerline);
255 
256  for(auto &spline: *mClSplinesPtr)
257  {
258  mBloodVessels++;
259  // Smooth the splines
260  for(int j = 0; j < nConvolutions; j++)
261  {
262  spline.applyConvolution({0.25, 0.50, 0.25 });
263  }
264 
265  // Compute control points for splines
266  spline.compute();
267 
268  // Find all the intersections
269  spline.findAllIntersections(*images);
270  spline.getIntersections().setVelocityEstimationCutoff(cutoff,1.0);
271  mIntersections += spline.getIntersections().size();
272 
273  // Now that we know the intersection points,
274  // we can go through all the image planes and do the region growing.
275  for_each(spline.getIntersections().begin(), spline.getIntersections().end(),[](Intersection<double> &it){it.regionGrow();});
276 
277  // We may now do the direction vector estimation
278  // Using the default parameters set in IntersectionSet constructor
279  spline.getIntersections().estimateDirection();
280 
281  // Do aliasing correction
282  if (Vnyq > 0){
283  spline.getIntersections().correctAliasing(Vnyq);
284  }
285 
286  // Least squares velocity estimates
287  spline.getIntersections().estimateVelocityLS();
288 
289  // Output direction and LS velocity
290  if (verbose)
291  {
292  cerr << "Spline " << mBloodVessels << " gave direction "
293  << spline.getIntersections().getEstimatedDirection()
294  << " LS velocity " << spline.getIntersections().getEstimatedVelocity() << endl;
295  }
296  }
297  if (verbose)
298  {
299  cerr << "Found " << mIntersections << " intersections in " << mBloodVessels << " splines.\n";
300  }
301 }
302 
303 
304 
305 vtkSmartPointer<vtkPolyData> AngleCorrection::computeVtkPolyData( vectorSpline3dDoublePtr splines, double uncertainty_limit, double minArrowDist)
306 {
307  vtkSmartPointer<vtkPoints> pointarray = vtkSmartPointer<vtkPoints>::New();
308  vtkSmartPointer<vtkDoubleArray> flowdirection = vtkSmartPointer<vtkDoubleArray>::New();
309  flowdirection->SetName("Flow direction");
310  flowdirection->SetNumberOfComponents(3);
311  vtkSmartPointer<vtkDoubleArray> dir_uncertainty = vtkSmartPointer<vtkDoubleArray>::New();
312  dir_uncertainty->SetName("FLow direction uncertainty");
313  vtkSmartPointer<vtkDoubleArray> velocitydata = vtkSmartPointer<vtkDoubleArray>::New();
314  velocitydata->SetName("Vessel velocity");
315 
316  double p[3];
317  double p_prev[3];
318  double p_temp[3];
319  double flow_vector[3];
320  double flow_vector_n[3];
321  double abs_dir;
322  double abs_vessel_vel;
323  mBloodVesselsRemoved = 0;
324  for(auto &spline: *splines)
325  {
326 
327  double flow_direction = spline.getIntersections().getEstimatedDirection();
328  abs_dir = flow_direction*sgn(flow_direction);
329  abs_vessel_vel = spline.getIntersections().getEstimatedVelocity();
330  abs_vessel_vel = abs_vessel_vel*sgn(abs_vessel_vel);
331 
332 
333  if( abs_dir< uncertainty_limit){
334  mBloodVesselsRemoved++;
335  continue;
336  }
337  if(std::isnan(abs_dir) || std::isnan(abs_vessel_vel) ){
338  continue;
339  }
340 
341  spline.evaluateSingle(0,p_prev);
342  for(double t = 0; t < (spline.length()-1); t+=0.3){
343  spline.evaluateSingle(t,p);
344  if(std::isnan(p[0]) || std::isnan(p[1]) || std::isnan(p[2]) ){
345  continue;
346  }
347 
348  p_temp[0]=p[0]-p_prev[0];
349  p_temp[1]=p[2]-p_prev[1];
350  p_temp[1]=p[2]-p_prev[2];
351  if(length3d(p_temp)<minArrowDist){
352  continue;
353  }
354  p_prev[0]=p[0];
355  p_prev[1]=p[1];
356  p_prev[2]=p[2];
357 
358  spline.derivativeSingle(t,flow_vector);
359 
360  flow_vector[0]=flow_direction*flow_vector[0];
361  flow_vector[1]=flow_direction*flow_vector[1];
362  flow_vector[2]=flow_direction*flow_vector[2];
363 
364  normalize3d(flow_vector_n, flow_vector);
365  if(std::isnan(flow_vector_n[0]) || std::isnan(flow_vector_n[1]) || std::isnan(flow_vector_n[2]) ){
366  continue;
367  }
368 
369  pointarray->InsertNextPoint(p);
370  flowdirection->InsertNextTupleValue(flow_vector_n);
371  dir_uncertainty->InsertNextTupleValue(&abs_dir);
372  velocitydata->InsertNextTupleValue(&abs_vessel_vel);
373  }
374  }
375  // Put it all inside a vtkPolyData
376  vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
377  vtkSmartPointer<vtkPointData> pointdata = polydata->GetPointData();
378  polydata->SetPoints(pointarray);
379 
380  pointdata->AddArray(flowdirection);
381  pointdata->AddArray(dir_uncertainty);
382  pointdata->AddArray(velocitydata);
383  if (mBloodVesselsRemoved){
384  cerr << "Removed " << mBloodVesselsRemoved << " segment(s) due to an uncertainty limit of " << uncertainty_limit << "\n";
385  }
386  return polydata;
387 }
388 
389 
390 
391 bool AngleCorrection::EqualVtkPolyData( vtkSmartPointer<vtkPolyData> leftHandSide, vtkSmartPointer<vtkPolyData> rightHandSide)
392 {
393  if( leftHandSide->GetNumberOfCells()!=rightHandSide->GetNumberOfCells()) return false;
394  if( leftHandSide->GetNumberOfVerts()!=rightHandSide->GetNumberOfVerts()) return false;
395  if( leftHandSide->GetNumberOfLines()!=rightHandSide->GetNumberOfLines()) return false;
396  if( leftHandSide->GetNumberOfPolys()!=rightHandSide->GetNumberOfPolys()) return false;
397  if(leftHandSide->GetNumberOfStrips()!=rightHandSide->GetNumberOfStrips())return false;
398  unsigned int numberOfPointsRight = rightHandSide->GetNumberOfPoints();
399  unsigned int numberOfPointsLeft = leftHandSide->GetNumberOfPoints();
400  if(numberOfPointsLeft!=numberOfPointsRight) return false;
401 
402  double pointOne[3];
403  double pointTwo[3];
404  for( unsigned int i( 0 ); i < numberOfPointsRight; i++ )
405  {
406  rightHandSide->GetPoint(i, pointOne);
407  leftHandSide->GetPoint(i, pointTwo);
408  double x = pointOne[0] - pointTwo[0];
409  double y = pointOne[1] - pointTwo[1];
410  double z = pointOne[2] - pointTwo[2];
411  double distance = x*x + y*y + z*z;
412  if( distance > 0.001 ) return false;
413  }
414  return true;
415 }
void normalize3d(T normalized[3], T const toNormalize[3])
Definition: helpers.hpp:34
void setInput(vtkSmartPointer< vtkPolyData > vpd_centerline, const char *image_prefix, double Vnyq, double cutoff, int nConvolutions, double uncertainty_limit=0.0, double minArrowDist=1.0)
int sgn(T val)
Definition: helpers.hpp:9
T length3d(T const a[3])
Definition: helpers.hpp:19
vector< Spline3D< double > > vectorSpline3dDouble
static vector< Spline3D< T > > * build(vtkSmartPointer< vtkPolyData > data)
Definition: spline3d.hpp:323
void writeDirectionToVtkFile(const char *filename)
static vector< MetaImage > * readImages(const string &prefix)
Definition: metaimage.hpp:254
vtkSmartPointer< vtkPolyData > getOutput()
vectorSpline3dDouble getClSpline()
void reportError(std::string errMsg)
vectorSpline3dDouble * vectorSpline3dDoublePtr