5 #include <vtkDoubleArray.h> 6 #include <vtkSmartPointer.h> 8 #include <vtkPolyDataReader.h> 9 #include <vtkPolyData.h> 10 #include <vtkPolyDataWriter.h> 11 #include <vtkPointData.h> 28 mVelDataPtr =
new vector<MetaImage<inData_t>>();
29 mClData=vtkSmartPointer<vtkPolyData>::New();
39 mBloodVesselsRemoved = 0;
46 mClSplinesPtr->clear();
55 mVelDataPtr =
new vector<MetaImage<inData_t>>();
56 mClData=vtkSmartPointer<vtkPolyData>::New();
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 ");
84 if((velData->size() > 0 && mVelDataPtr!=velData) || mVelImagePrefix.size()==0)
93 mnConvolutions!=nConvolutions ||
94 !EqualVtkPolyData(mClData,vpd_centerline))
96 mClData->DeepCopy(vpd_centerline);
99 mnConvolutions=nConvolutions;
103 if(mUncertainty_limit!=uncertainty_limit ||
104 mMinArrowDist!=minArrowDist)
106 mUncertainty_limit=uncertainty_limit;
107 mMinArrowDist=minArrowDist;
114 void AngleCorrection::setInput(vtkSmartPointer<vtkPolyData> vpd_centerline,
const char* velImagePrefix ,
double Vnyq,
double cutoff,
int nConvolutions,
double uncertainty_limit,
double minArrowDist)
116 if(mVelImagePrefix!=std::string(velImagePrefix))
118 mVelImagePrefix=std::string(velImagePrefix);
119 mVelDataPtr->clear();
123 std::string filename=std::string(velImagePrefix);
124 filename.append(
"0.mhd");
126 vtkSmartPointer<vtkMetaImageReader> reader= vtkSmartPointer<vtkMetaImageReader>::New();
127 vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
128 reader->AddObserver(vtkCommand::ErrorEvent,errorObserver);
130 reader->SetFileName(filename.c_str());
133 if(!reader->CanReadFile(filename.c_str())){
134 reportError(
"ERROR: Could not read velocity data \n");
137 setInput(vpd_centerline,
new vector<
MetaImage<inData_t>>(), Vnyq, cutoff, nConvolutions, uncertainty_limit, minArrowDist);
141 void AngleCorrection::setInput(
const char* centerline,
const char* image_prefix,
double Vnyq,
double cutoff,
int nConvolutions,
double uncertainty_limit,
double minArrowDist)
144 cerr <<
"Input params: " << centerline<<
" " << image_prefix <<
" " << Vnyq <<
" " << nConvolutions <<
" " << uncertainty_limit<<
" " << minArrowDist <<endl;
146 vtkSmartPointer<vtkPolyDataReader> clReader = vtkSmartPointer<vtkPolyDataReader>::New();
148 vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
149 clReader->AddObserver(vtkCommand::ErrorEvent,errorObserver);
151 clReader->SetFileName(centerline);
153 if (errorObserver->GetError())
155 reportError(
"ERROR: Could not read center line data \n"+ errorObserver->GetErrorMessage());
158 if (errorObserver->GetWarning()){
159 cerr <<
"Caught warning while reading center line data \n! " << errorObserver->GetWarningMessage();
162 vtkSmartPointer<vtkPolyData> vpd_centerline = clReader->GetOutput();
164 setInput(vpd_centerline, image_prefix, Vnyq, cutoff, nConvolutions, uncertainty_limit, minArrowDist);
172 cerr <<
"params: " << mnConvolutions<<
" " << mCutoff <<
" " << mUncertainty_limit <<
" " << mMinArrowDist <<
" " << mVnyq <<endl;
176 cerr <<
"Invalid input " << endl;
182 if(mVelDataPtr->size() == 0)
184 cerr <<
"Loading data " << endl;
185 mVelDataPtr->clear();
193 cerr <<
"started step 1 of 2 "<< endl;
194 angle_correction_impl(mClData, mVelDataPtr, mVnyq, mCutoff, mnConvolutions);
197 if(mUpdate1 || mUpdate2)
200 cerr <<
"started step 2 of 2" << endl;
201 mOutput= computeVtkPolyData(mClSplinesPtr, mUncertainty_limit, mMinArrowDist);
217 return *mClSplinesPtr;
225 vtkSmartPointer<vtkPolyData> polydata = getOutput();
226 if(!polydata) polydata=vtkSmartPointer<vtkPolyData>::New();
228 vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
229 vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
230 writer->AddObserver(vtkCommand::ErrorEvent,errorObserver);
232 writer->SetFileName(filename);
233 #if VTK_MAJOR_VERSION <= 5 234 writer->SetInput(polydata);
236 writer->SetInputData(polydata);
240 if (errorObserver->GetError())
242 reportError(
"ERROR: Could not write file to disk \n"+ errorObserver->GetErrorMessage());
244 if (errorObserver->GetWarning()){
245 cerr <<
"Caught warning while not writing file to disk! \n " << errorObserver->GetWarningMessage();
249 void AngleCorrection::angle_correction_impl(vtkSmartPointer<vtkPolyData> vpd_centerline, vector<
MetaImage<inData_t> >* images ,
double Vnyq,
double cutoff,
int nConvolutions)
251 bool verbose =
false;
253 mClSplinesPtr->clear();
256 for(
auto &spline: *mClSplinesPtr)
260 for(
int j = 0; j < nConvolutions; j++)
262 spline.applyConvolution({0.25, 0.50, 0.25 });
269 spline.findAllIntersections(*images);
270 spline.getIntersections().setVelocityEstimationCutoff(cutoff,1.0);
271 mIntersections += spline.getIntersections().size();
275 for_each(spline.getIntersections().begin(), spline.getIntersections().end(),[](
Intersection<double> &it){it.regionGrow();});
279 spline.getIntersections().estimateDirection();
283 spline.getIntersections().correctAliasing(Vnyq);
287 spline.getIntersections().estimateVelocityLS();
292 cerr <<
"Spline " << mBloodVessels <<
" gave direction " 293 << spline.getIntersections().getEstimatedDirection()
294 <<
" LS velocity " << spline.getIntersections().getEstimatedVelocity() << endl;
299 cerr <<
"Found " << mIntersections <<
" intersections in " << mBloodVessels <<
" splines.\n";
305 vtkSmartPointer<vtkPolyData> AngleCorrection::computeVtkPolyData(
vectorSpline3dDoublePtr splines,
double uncertainty_limit,
double minArrowDist)
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");
319 double flow_vector[3];
320 double flow_vector_n[3];
322 double abs_vessel_vel;
323 mBloodVesselsRemoved = 0;
324 for(
auto &spline: *splines)
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);
333 if( abs_dir< uncertainty_limit){
334 mBloodVesselsRemoved++;
337 if(std::isnan(abs_dir) || std::isnan(abs_vessel_vel) ){
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]) ){
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];
358 spline.derivativeSingle(t,flow_vector);
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];
365 if(std::isnan(flow_vector_n[0]) || std::isnan(flow_vector_n[1]) || std::isnan(flow_vector_n[2]) ){
369 pointarray->InsertNextPoint(p);
370 flowdirection->InsertNextTupleValue(flow_vector_n);
371 dir_uncertainty->InsertNextTupleValue(&abs_dir);
372 velocitydata->InsertNextTupleValue(&abs_vessel_vel);
376 vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
377 vtkSmartPointer<vtkPointData> pointdata = polydata->GetPointData();
378 polydata->SetPoints(pointarray);
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";
391 bool AngleCorrection::EqualVtkPolyData( vtkSmartPointer<vtkPolyData> leftHandSide, vtkSmartPointer<vtkPolyData> rightHandSide)
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;
404 for(
unsigned int i( 0 ); i < numberOfPointsRight; i++ )
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;
void normalize3d(T normalized[3], T const toNormalize[3])
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)
vector< Spline3D< double > > vectorSpline3dDouble
static vector< Spline3D< T > > * build(vtkSmartPointer< vtkPolyData > data)
void writeDirectionToVtkFile(const char *filename)
vtkSmartPointer< vtkPolyData > getOutput()
vectorSpline3dDouble getClSpline()
void reportError(std::string errMsg)
vectorSpline3dDouble * vectorSpline3dDoublePtr