Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
SeansVesselReg.cxx
Go to the documentation of this file.
1 #include "SeansVesselReg.hxx"
2 #include "HackTPSTransform.hxx"
3 
4 #include <iostream>
5 #include <time.h>
6 #include <fstream>
7 
8 #include <QFileInfo>
9 
10 #include "cxImage.h"
11 #include "cxTypeConversions.h"
13 #include "cxReporter.h"
14 #include <boost/math/special_functions/fpclassify.hpp>
15 
16 #include "vtkClipPolyData.h"
17 #include "vtkPlanes.h"
18 #include "cxLogger.h"
19 #include "vtkPoints.h"
20 #include "vtkPolyData.h"
21 #include "vtkCellArray.h"
22 #include "vtkCellLocator.h"
23 #include "vtkMINCImageReader.h"
24 #include "vtkTransform.h"
25 #include "vtkImageData.h"
26 #include "vtkGeneralTransform.h"
27 #include "vtkMath.h"
28 #include "vtkSortDataArray.h"
29 #include "vtkMaskPoints.h"
30 #include "vtkPointData.h"
31 #include "vtkLandmarkTransform.h"
32 #include "vtkFloatArray.h"
33 #include "cxMesh.h"
34 #include "cxLogger.h"
35 
36 namespace cx
37 {
38 
39 SeansVesselReg::SeansVesselReg()// : mInvertedTransform(false)
40 {
41  mt_auto_lts = true;
42  mt_ltsRatio = 80;
44  mt_lambda = 0;
45  mt_sigma = 1.0;
46  mt_doOnlyLinear = false;
47  mt_sampleRatio = 1;
50  mt_verbose = false;
51  mt_maximumDurationSeconds = 1E6; // Random high number
52  margin = 40;
53 }
54 
56 {
57 }
58 
64 bool SeansVesselReg::initialize(DataPtr source, DataPtr target, QString logPath)
65 {
66  if (mt_verbose)
67  {
68  reporter()->sendDebug("SOURCE: " + source->getUid());
69  reporter()->sendDebug("TARGET: " + target->getUid());
70  }
71 
72  m_logPath = logPath;
73  mLastRun = this->createContext(source, target);
74  return mLastRun != NULL;
75 }
76 
83 {
84  if (mt_verbose)
85  {
86  std::cout << "stop Threshold:" << mt_distanceDeltaStopThreshold << endl;
87  std::cout << "sigma:" << mt_sigma << endl;
88  std::cout << "lts Ratio:" << mt_ltsRatio << endl;
89  std::cout << "linear flag:" << mt_doOnlyLinear << endl;
90  std::cout << "sample flag:" << mt_sampleRatio << endl;
91  std::cout << "single Point Threshold:" << mt_singlePointThreshold << endl;
92  }
93  QTime start = QTime::currentTime();
94 
95  ContextPtr context = mLastRun;
96 
97  if (!context)
98  return false;
99 
100  if (mt_auto_lts)
101  {
102  context = this->linearRefineAllLTS(context);
103  report(QString("Auto LTS: Found best match using %1\%.").arg(context->mLtsRatio));
104  }
105  else
106  {
107  this->linearRefine(context);
108  }
109 
110  // add a nonlinear step to the end:
111  if (!mt_doOnlyLinear)
112  {
113  this->performOneRegistration(context, false);
114  }
115 
116  printOutResults(m_logPath + "/Vessel_Based_Registration_", context->mConcatenation);
117 
118  if (mt_verbose)
119  std::cout << QString("\n\nV2V Execution time: %1s").arg(start.secsTo(QTime::currentTime())) << endl;
120 
121  mLastRun = context;
122 // mLinearTransformResult = this->getLinearTransform(context->mConcatenation);
123 
124  return true;
125 }
126 
128 {
129  ContextPtr context = mLastRun;
130  if (!context)
131  return false;
132  this->performOneRegistration(context, mt_doOnlyLinear);
133  mLastRun = context;
134  return true;
135 }
136 
138 {
139  return mLastRun && mLastRun->getFixedPoints() && mLastRun->getMovingPoints();
140 }
141 
147 {
148  // all lts values to search along:
149  std::vector<int> lts;
150  lts.push_back(40);
151  lts.push_back(50);
152  lts.push_back(60);
153  lts.push_back(70);
154  lts.push_back(75);
155  lts.push_back(80);
156  lts.push_back(85);
157  lts.push_back(90);
158  lts.push_back(95);
159  lts.push_back(100);
160 
161  std::vector<ContextPtr> paths;
162 
163  // iterate along all paths
164  for (unsigned i=0; i<lts.size(); ++i)
165  {
166  ContextPtr current = this->splitContext(seed);
167  paths.push_back(current);
168  current->mLtsRatio = lts[i];
169  this->linearRefine(current);
170  if (mt_verbose)
171  std::cout << QString("LTS=%1, metric=%2").arg(current->mLtsRatio).arg(current->mMetric) << std::endl;
172  }
173 
174  // search for best path
175  int bestPath = 0;
176  for (unsigned i=0; i<paths.size(); ++i)
177  {
178  if (paths[i]->mMetric < paths[bestPath]->mMetric)
179  bestPath = i;
180  }
181 
182  // return value
183  return paths[bestPath];
184 }
185 
190 {
191  // Perform registrations iteratively until convergence is reached:
192  double previousMetric = 1E6;
193  QDateTime t0 = QDateTime::currentDateTime();
194  for (int iteration = 1; iteration < mt_maximumNumberOfIterations && (t0.msecsTo(QDateTime::currentDateTime()) < mt_maximumDurationSeconds*1000); ++iteration)
195  {
196  this->performOneRegistration(context, true);
197  double difference = context->mMetric - previousMetric;
198 
199  if (mt_verbose)
200  {
201 // std::cout << myNumberOfIterations << " ";
202 // std::cout.flush();
203  std::cout << QString("iteration\t%1\trms:\t%2").arg(iteration).arg(context->mMetric) << std::endl;
204  }
205 
206  // Check for convergence
207  if (fabs(difference) < mt_distanceDeltaStopThreshold)
208  break;
209 
210  previousMetric = context->mMetric;
211  }
212 }
213 
218 {
219  ContextPtr retval = ContextPtr(new Context);
220 
221  retval->mLtsRatio = context->mLtsRatio;
222  retval->mInvertedTransform = context->mInvertedTransform;
223 
224  // constant data: shallow copy
225  retval->mTargetPointLocator = context->mTargetPointLocator;
226  retval->mTargetPoints = context->mTargetPoints;
227 
228  // will be modified: deep copy
229  retval->mSourcePoints = vtkPointsPtr::New();
230  retval->mSourcePoints->DeepCopy(context->mSourcePoints);
231 
232  // return value: initialize only
233  retval->mConcatenation = vtkGeneralTransformPtr::New();;
234 // retval->mTransform = context->mTransform;
235  retval->mMetric = 0;
236 
237  return retval;
238 }
239 
244 {
245  if (!source || !target)
247 
248  vtkPolyDataPtr targetPolyData = this->convertToPolyData(target, "target");
249  vtkPolyDataPtr inputSourcePolyData = this->convertToPolyData(source, "source");
250 
251  if (!targetPolyData || !inputSourcePolyData) // error message has already been emitted in convertToPolyData()
252  return ContextPtr();
253 
254  vtkPolyDataPtr sourcePolyData = this->crop(inputSourcePolyData, targetPolyData, margin);
255 
256  //Make sure we have stuff to work with
257  if (!sourcePolyData->GetNumberOfPoints() || !targetPolyData->GetNumberOfPoints())
258  {
259  std::cerr << "No data after cropping, failed." << std::endl;
260  return ContextPtr();
261  }
262 
263  ContextPtr context = ContextPtr(new Context);
264 
265  context->mConcatenation = vtkGeneralTransformPtr::New();
266  context->mInvertedTransform = false;
267  context->mMetric = -1;
268 
269  // Algorithm requires #source < #target
270  // swap if this is not the case
271  if (sourcePolyData->GetNumberOfPoints() > targetPolyData->GetNumberOfPoints())
272  {
273  //INVERT
274  if (mt_verbose)
275  std::cout << "inverted vessel reg" << std::endl;
276  context->mInvertedTransform = true;
277  std::swap(sourcePolyData, targetPolyData);
278  }
279 
280  vtkIdType numPoints = sourcePolyData->GetNumberOfPoints();
281  if (mt_verbose)
282  {
283  std::cout << "total number of source points:" << numPoints << ", target points: " << targetPolyData->GetNumberOfPoints() << endl;
284  std::cout << "number of source points to be sampled:" << ((int) (numPoints * mt_ltsRatio) / 100) << "\n" << endl;
285  }
286 
287 
288  // Create locator for target points
289  context->mTargetPoints = targetPolyData;
290  context->mTargetPointLocator = vtkCellLocatorPtr::New();
291  context->mTargetPointLocator->SetDataSet(targetPolyData);
292  context->mTargetPointLocator->SetNumberOfCellsPerBucket(1);
293  context->mTargetPointLocator->BuildLocator();
294 
295  //Since we are going to play with the data, we have to make a copy
296  context->mSourcePoints = vtkPointsPtr::New();
297  context->mSourcePoints->DeepCopy(sourcePolyData->GetPoints());
298 
299  context->mLtsRatio = mt_ltsRatio;
300 
301  return context;
302 }
303 
311 {
312  if (!context)
313  context = mLastRun;
314  if (!context)
315  return;
316 
317  if (context->mSortedSourcePoints || context->mSortedTargetPoints)
318  return;
319 
320  // total number of source points:
321  int numPoints = context->mSourcePoints->GetNumberOfPoints();
322  // number of source points used in each iteration (the rest is temporarily rejected from the computation)
323  int nb_points = ((int) (numPoints * context->mLtsRatio) / 100);
324 // std::cout << QString("onestep %1/%2").arg(nb_points).arg(numPoints) << std::endl;
325 
326  // - closestPoint is used so that the internal state of LandmarkTransform remains
327  // correct whenever the iteration process is stopped (hence its source
328  // and landmark points might be used in a vtkThinPlateSplineTransform).
329  vtkPointsPtr closestPoint = vtkPointsPtr::New();
330  closestPoint->SetNumberOfPoints(numPoints);
331 
332  // Fill points with the closest points to each vertex in input
333  vtkFloatArrayPtr residuals = vtkFloatArrayPtr::New();
334  residuals->SetNumberOfValues(numPoints);
335 
336  vtkIdListPtr IdList = vtkIdListPtr::New();
337  IdList->SetNumberOfIds(numPoints);
338  double total_distance = 0;
339  double distanceSquared = 0;
340  //Find closest points to all source points
341  for (int i = 0; i < numPoints; ++i)
342  {
343  //Check the distance to neighbouring points (neighbours should be matched to nearby points)
344  vtkIdType cell_id;
345  int sub_id;
346  double outPoint[3];
347  context->mTargetPointLocator->FindClosestPoint(context->mSourcePoints->GetPoint(i), outPoint, cell_id, sub_id, distanceSquared);
348  closestPoint->SetPoint(i, outPoint);
349  if ((boost::math::isnan)(distanceSquared))
350  {
351  std::cout << "nan found during findClosestPoint!" << std::endl;
352  {
353  context->mMetric = 1E6;
354  return;
355  }
356  }
357  residuals->InsertValue(i, distanceSquared);
358  IdList->InsertId(i, i);
359  total_distance += sqrt(distanceSquared);
360  }
361 
362  // quality of the current iteration
363  context->mMetric = total_distance / numPoints;
364 
365  vtkSortDataArrayPtr sort = vtkSortDataArrayPtr::New();
366  sort->Sort(residuals, IdList);
367  context->mSortedSourcePoints = this->createSortedPoints(IdList, context->mSourcePoints, nb_points);
368  context->mSortedTargetPoints = this->createSortedPoints(IdList, closestPoint, nb_points);
369 }
370 
383 {
384  // compute distances if not already done.
385  this->computeDistances(context);
386 
387  if (!context->mSortedSourcePoints || !context->mSortedTargetPoints)
388  return;
389 
390  if (linear)
391  {
392  context->mTransform = linearRegistration(context->mSortedSourcePoints, context->mSortedTargetPoints);
393  }
394  else
395  {
396  context->mTransform = nonLinearRegistration(context->mSortedSourcePoints, context->mSortedTargetPoints);
397  }
398 
399  // Transform the source points with the transform found during this iteration,
400  // in order to use an updated guess for the next iteration
401  context->mSourcePoints = this->transformPoints(context->mSourcePoints, context->mTransform);
402 
403  // clear sorting data - enables us to call iteratively without fuss.
404  context->mSortedSourcePoints = vtkPointsPtr();
405  context->mSortedTargetPoints = vtkPointsPtr();
406 
407  // add transform from this iteration to the total
408  context->mConcatenation->Concatenate(context->mTransform);
409 
410  this->computeDistances(context);
411 }
412 
417 vtkPointsPtr SeansVesselReg::createSortedPoints(vtkIdListPtr sortedIDList, vtkPointsPtr unsortedPoints, int numPoints)
418 {
419  vtkPointsPtr retval = vtkPointsPtr::New();
420  retval->SetNumberOfPoints(numPoints);
421 
422  double temp_point[3];
423 
424  for (int i = 0; i < numPoints; ++i)
425  {
426  vtkIdType index = sortedIDList->GetId(i);
427  unsortedPoints->GetPoint(index, temp_point); // source points to use in tps
428  retval->SetPoint(i, temp_point);
429  }
430 
431  return retval;
432 }
433 
438 {
439 // QTime start = QTime::currentTime();
440 
441  int N = input->GetNumberOfPoints();
442  vtkPointsPtr retval = vtkPointsPtr::New();
443  retval->SetNumberOfPoints(N);
444 
445  //Transform ALL source points
446  double temp[3];
447  for (int i = 0; i < N; ++i)
448  {
449  transform->InternalTransformPoint(input->GetPoint(i), temp);
450  retval->SetPoint(i, temp);
451  }
452 
453 // std::cout << QString("\nV2V Transform time: %1s").arg(start.secsTo(QTime::currentTime())) << endl;
454  return retval;
455 }
456 
458  vtkPointsPtr sortedTargetPoints)
459 {
460 // std::cout << "linear" << std::endl;
461  //Build landmark transform
462  vtkLandmarkTransformPtr lmt = vtkLandmarkTransformPtr::New();
463  lmt->SetSourceLandmarks(sortedSourcePoints);
464  lmt->SetTargetLandmarks(sortedTargetPoints);
465  lmt->SetModeToRigidBody();
466  lmt->Modified();
467  lmt->Update();
468 
469  return lmt;
470 }
471 
473  vtkPointsPtr sortedTargetPoints)
474 {
475 // std::cout << "nonlinear" << std::endl;
476  vtkPolyDataPtr tpsSourcePolyData = this->convertToPolyData(sortedSourcePoints);
477  vtkPolyDataPtr tpsTargetPolyData = this->convertToPolyData(sortedTargetPoints);
478 
479  vtkMaskPointsPtr mask1 = vtkMaskPointsPtr::New();
480  mask1->SetInputData(tpsSourcePolyData);
481  mask1->SetOnRatio(mt_sampleRatio);
482  mask1->Update();
483  vtkMaskPointsPtr mask2 = vtkMaskPointsPtr::New();
484  mask2->SetInputData(tpsTargetPolyData);
485  mask2->SetOnRatio(mt_sampleRatio);
486  mask2->Update();
487 
488  // Build the thin plate spline transform
489  vtkThinPlateSplineTransformPtr tps = vtkThinPlateSplineTransformPtr::New();
490  tps->SetSourceLandmarks(mask1->GetOutput()->GetPoints());
491  tps->SetTargetLandmarks(mask2->GetOutput()->GetPoints());
492  tps->SetBasisToR();
493  tps->SetSigma(mt_sigma);
494  QTime start = QTime::currentTime();
495  tps->Update();
496  std::cout << QString("\nV2V NL time: %1s").arg(start.secsTo(QTime::currentTime())) << endl;
497 
498  return tps;
499 }
500 
502 {
503  if (!data)
504  {
505  CX_LOG_ERROR(QString("Can't execute: %1 is not set").arg(id));
506  return vtkPolyDataPtr();
507  }
508 
509  // ImagePtr image = boost::dynamic_pointer_cast<Image>(data);
510  MeshPtr mesh = boost::dynamic_pointer_cast<Mesh>(data);
511 
512 // if (image)
513 // {
514 // //Grab the information from the files of target then
515 // //filter out points not fit for the threshold
516 // return this->extractPolyData(image, mt_singlePointThreshold, 0);
517 // }
518 
519  if (!mesh)
520  {
521  CX_LOG_ERROR(QString("Can't execute: %1=%2 is not a mesh").arg(id, data->getName()));
522  return vtkPolyDataPtr();
523  }
524 
525  vtkPolyDataPtr poly = mesh->getTransformedPolyDataCopy(mesh->get_rMd());
526 
527  if (!poly || !poly->GetNumberOfPoints())
528  {
529  CX_LOG_ERROR(QString("Can't execute: %1=%2 is empty").arg(id, data->getName()));
530  return vtkPolyDataPtr();
531  }
532 
533  return poly;
534 }
535 
537 {
538  if (mInvertedTransform)
539  {
540  return mTargetPoints;
541  }
542  else
543  {
544  return SeansVesselReg::convertToPolyData(mSourcePoints);
545  }
546 }
547 
549 {
550  if (mInvertedTransform)
551  {
552  return SeansVesselReg::convertToPolyData(mSourcePoints);
553  }
554  else
555  {
556  return mTargetPoints;
557  }
558 }
559 
561 {
562  vtkPolyDataPtr retval = vtkPolyDataPtr::New();
563  // draw lines
564  retval->Allocate();
565  vtkPointsPtr verts = vtkPointsPtr::New();
566  for (int i = 0; i < this->mSortedSourcePoints->GetNumberOfPoints(); ++i)
567  {
568  verts->InsertNextPoint(this->mSortedSourcePoints->GetPoint(i));
569  verts->InsertNextPoint(this->mSortedTargetPoints->GetPoint(i));
570 
571  vtkIdType connectivity[2];
572  connectivity[0] = 2 * i;
573  connectivity[1] = 2 * i + 1;
574  retval->InsertNextCell(VTK_LINE, 2, connectivity);
575  }
576  retval->SetPoints(verts);
577  return retval;
578 }
579 
580 
582 {
583  vtkCellArrayPtr cellArray = vtkCellArrayPtr::New();
584  int N = input->GetNumberOfPoints();
585 
586  for (int i=0; i<N ; ++i)
587  {
588  cellArray->InsertNextCell(1);
589  cellArray->InsertCellPoint(i);
590  }
591 
592  vtkPolyDataPtr retval = vtkPolyDataPtr::New();
593  retval->SetPoints(input);
594  retval->SetVerts(cellArray);
595 
596  return retval;
597 }
598 
600 {
601  for (int q = 0; q < points->GetNumberOfPoints(); ++q)
602  {
603  Vector3D p(points->GetPoint(q));
604  std::cout << q << "\t" << p[0] << " " << p[1] << " " << p[2] << " " << std::endl;
605  }
606 }
607 
609 {
610  print(data->GetPoints());
611 }
612 
618 {
619  // use the bounding box of target to clip the source data.
620  fixed->GetPoints()->ComputeBounds();
621  double* targetBounds = fixed->GetPoints()->GetBounds();
622  targetBounds[0] -= margin;
623  targetBounds[1] += margin;
624  targetBounds[2] -= margin;
625  targetBounds[3] += margin;
626  targetBounds[4] -= margin;
627  targetBounds[5] += margin;
628 
629  // clip the source data with a box
630  vtkPlanesPtr box = vtkPlanesPtr::New();
631  box->SetBounds(targetBounds);
632  if (mt_verbose)
633  std::cout << "bb: " << DoubleBoundingBox3D(targetBounds) << std::endl;
634  vtkClipPolyDataPtr clipper = vtkClipPolyDataPtr::New();
635  clipper->SetInputData(input);
636  clipper->SetClipFunction(box);
637  clipper->SetInsideOut(true);
638  clipper->Update();
639 
640  int oldSource = input->GetPoints()->GetNumberOfPoints();
641 
642  int clippedSource = 0;
643 
644  if (clipper->GetOutput()->GetPoints())
645  {
646  clippedSource = clipper->GetOutput()->GetPoints()->GetNumberOfPoints();
647  }
648 
649  if (clippedSource < oldSource)
650  {
651  double ratio = double(oldSource - clippedSource) / double(oldSource);
652  if (mt_verbose)
653  std::cout << "Removed " << ratio * 100 << "%" << " of the source data. Outside the target data bounds." << std::endl;
654  }
655 
656  return clipper->GetOutput();
657 }
658 
660 {
661  if (!context)
662  context = mLastRun;
663 
664  Transform3D retval = this->getLinearTransform(context->mConcatenation);
665  if (context->mInvertedTransform)
666  retval = retval.inv();
667 
668  return retval;
669 }
670 
672 {
673  if (!context)
674  context = mLastRun;
675  if (!context)
676  return -1;
677  return context->mMetric;
678 }
679 
681 {
682  if (!context)
683  context = mLastRun;
684  return context->mLtsRatio;
685 }
686 
689 Transform3D SeansVesselReg::getLinearTransform(vtkGeneralTransformPtr concatenation)
690 {
691  Transform3D M_acc = Transform3D::Identity(); // accumulated linear transform
692 
693  for (int i = 0; i < concatenation->GetNumberOfConcatenatedTransforms(); ++i)
694  {
695  vtkLandmarkTransformPtr LM_i = (vtkLandmarkTransform*) concatenation->GetConcatenatedTransform(i);
696  if (!LM_i.GetPointer()) // skip non-landmark transforms
697  continue;
698  Transform3D M_i = Transform3D(LM_i->GetMatrix());
699  M_acc = M_acc * M_i;
700  }
701 
702  return M_acc.inverse();
703 }
704 
705 void SeansVesselReg::printOutResults(QString fileNamePrefix, vtkGeneralTransformPtr myConcatenation)
706 {
707 
708  vtkMatrix4x4Ptr l_tempMatrix = vtkMatrix4x4Ptr::New();
709  vtkMatrix4x4Ptr l_resultMatrix = vtkMatrix4x4Ptr::New();
710 
711  if (mt_doOnlyLinear)
712  l_tempMatrix->DeepCopy(((vtkLandmarkTransform*) myConcatenation->GetConcatenatedTransform(0))->GetMatrix());
713 
714  l_resultMatrix->Identity();
715  for (int i = 1; i < myConcatenation->GetNumberOfConcatenatedTransforms(); ++i)
716  {
717  vtkMatrix4x4::Multiply4x4(l_tempMatrix,
718  ((vtkLandmarkTransform*) myConcatenation->GetConcatenatedTransform(i))->GetMatrix(), l_resultMatrix);
719  l_tempMatrix->DeepCopy(l_resultMatrix);
720 
721  }
722  if (mt_verbose)
723  std::cout << "Filenameprefix: " << fileNamePrefix << std::endl;
724 
725  //std::string logsFolder = string_cast(cx::stateService()->getPatientData()->getActivePatientFolder())+"/Logs/";
726  //std::string logsFolder = "~/Patients/Logs/";
727  std::string nonLinearFile = fileNamePrefix.toStdString();
728  nonLinearFile += "--NonLinear";
729  nonLinearFile += ".txt";
730 
731  std::string linearFile = fileNamePrefix.toStdString();
732  linearFile += "--Linear";
733  linearFile += ".txt";
734 
735  if (mt_verbose)
736  std::cout << "Writing Results to " << nonLinearFile << " and " << linearFile << std::endl;
737 
738  if (!mt_doOnlyLinear)
739  {
740  ofstream file_out(nonLinearFile.c_str());
741 
742  //Non-linear Warped Transform
743  HackTPSTransform* l_theWarpTransform = ((HackTPSTransform*) myConcatenation->GetConcatenatedTransform(0));
744 
745  // Write the header
746  file_out << "MNI Transform File\n" << std::endl;
747  //file_out<<"SeansWarpyTransforms: source"<<std::endl;
748  //file_out<<"SeansWarpyTransforms: target\n"<<std::endl;
749  file_out << "Transform_Type = Thin_Plate_Spline_Transform;" << std::endl;
750  file_out << "Invert_Flag = True;" << std::endl;
751  file_out << "Number_Dimensions = 3;" << std::endl;
752  int n = l_theWarpTransform->GetSourceLandmarks()->GetNumberOfPoints();
753 
754  const double* const * theWarpMatrix = l_theWarpTransform->GetWarpMatrix();
755  double point[3];
756  file_out << "Points = " << std::endl;
757  for (int i = 0; i < n; i++)
758  {
759  l_theWarpTransform->GetSourceLandmarks()->GetPoint(i, point);
760  file_out << point[0] << " " << point[1] << " " << point[2];
761  if (i == n - 1)
762  file_out << ";" << std::endl;
763  else
764  file_out << std::endl;
765  }
766 
767  file_out << "Displacements = " << std::endl;
768  for (int i = 0; i < n + 4; i++)
769  {
770  file_out << theWarpMatrix[i][0] << " " << theWarpMatrix[i][1] << " " << theWarpMatrix[i][2];
771  if (i == n + 3)
772  file_out << ";" << std::endl;
773  else
774  file_out << std::endl;
775  }
776 
777  file_out.close();
778 
779  }
780 
781  ofstream file_out2(linearFile.c_str());
782 
783  //Linear Transform
784  file_out2 << "MNI Transform File\n" << std::endl;
785  //file_out2<<"SeansLinearTransforms: source"<<std::endl;
786  //file_out2<<"SeansLinearTransforms: target\n"<<std::endl;
787  file_out2 << "Transform_Type = Linear;" << std::endl;
788  file_out2 << "Linear_Transform = ";
789  file_out2 << std::endl;
790 
791  for (int i = 0; i < 3; i++)
792  {
793  for (int j = 0; j < 4; j++)
794  {
795 
796  file_out2 << l_resultMatrix->GetElement(i, j);
797  if (j != 3)
798  file_out2 << " ";
799  }
800  if (i == 2)
801  file_out2 << ";" << std::endl;
802  else
803  file_out2 << std::endl;
804  }
805  file_out2.close();
806 }
807 
813 vtkPolyDataPtr SeansVesselReg::extractPolyData(ImagePtr image, int p_neighborhoodFilterThreshold,
814  double p_BoundingBox[6])
815 {
816  vtkPolyDataPtr p_thePolyData = vtkPolyDataPtr::New();
817 
818  int l_dimensions[3];
819  image->getBaseVtkImageData()->GetDimensions(l_dimensions);
820  Vector3D spacing(image->getBaseVtkImageData()->GetSpacing());
821 
822  //set the transform
823  vtkTransformPtr l_dataTransform = vtkTransformPtr::New();
824  l_dataTransform->SetMatrix(image->get_rMd().getVtkMatrix());
825 
826  int l_startPosX, l_startPosY, l_startPosZ; //Beginings of neighborhood offsets
827  int l_stopPosX, l_stopPosY, l_stopPosZ; //Ends of neighborhood offsets
828  int i, j, k, ii, jj, kk, l_counts = 0; //Counter variables
829 
830  bool l_isNeighborFound; //Boolean for loop breakout
831 
832  //dimensions, fast if we first deref it in to variables
833  int l_dimX = l_dimensions[0];
834  int l_dimY = l_dimensions[1];
835  int l_dimZ = l_dimensions[2];
836 
837  double l_tempPoint[3];
838  //Offsets variables
839  int l_offsetI = 0;
840  int l_offsetJ = 0;
841  int l_offsetK = 0;
842  int l_kkjjOffset = 0;
843  int l_kkjjiiOffset = 0;
844  int l_offsetConstIJ = l_dimX * l_dimY;
845 
846  vtkDataArrayPtr l_allTheData = image->getBaseVtkImageData()->GetPointData()->GetScalars();
847  vtkPointsPtr l_dataPoints = vtkPointsPtr::New();
848  vtkCellArrayPtr l_dataCellArray = vtkCellArrayPtr::New();
849 
850  //Loop through the entire volume and precalculate the offsets for each
851  //point along with the points neighborhood offset
852  for (k = 0; k < l_dimZ; ++k)
853  {
854  //the start and stop offsets for the Z neighborhood (a gap in a cube)
855  l_startPosZ = k - p_neighborhoodFilterThreshold;
856  if (l_startPosZ < 0)
857  l_startPosZ = 0;
858  else
859  l_startPosZ *= l_offsetConstIJ;
860 
861  l_stopPosZ = k + p_neighborhoodFilterThreshold;
862  if (l_stopPosZ >= l_dimZ)
863  l_stopPosZ = (l_dimZ - 1) * l_offsetConstIJ;
864  else
865  l_stopPosZ *= l_offsetConstIJ;
866 
867  l_offsetK = k * l_offsetConstIJ;
868 
869  for (j = 0; j < l_dimY; ++j)
870  {
871  //the start and stop offsets for the Y neighborhood (a hole through a cube)
872  l_startPosY = j - p_neighborhoodFilterThreshold;
873  if (l_startPosY < 0)
874  l_startPosY = 0;
875  else
876  l_startPosY *= l_dimX;
877 
878  l_stopPosY = j + p_neighborhoodFilterThreshold;
879  if (l_stopPosY >= l_dimY)
880  l_stopPosY = (l_dimY - 1) * l_dimX;
881  else
882  l_stopPosY *= l_dimX;
883 
884  l_offsetJ = l_offsetK + (j * l_dimX);
885 
886  for (i = 0; i < l_dimX; ++i)
887  {
888  //the start and stop offsets for the X neighborhood (a voxel)
889  l_startPosX = i - p_neighborhoodFilterThreshold;
890  if (l_startPosX < 0)
891  l_startPosX = 0;
892 
893  l_stopPosX = i + p_neighborhoodFilterThreshold;
894  if (l_stopPosX >= l_dimX)
895  l_stopPosX = l_dimX - 1;
896 
897  l_offsetI = l_offsetJ + i;
898  l_isNeighborFound = 0;
899 
900  //See if this voxel contain a vessel center, if so do some more processing
901  if (*(l_allTheData->GetTuple(l_offsetI)))
902  {
903  // added by CA: use spacing when creating point. TODO: check with Ingrid if any other data are affected.
904  l_tempPoint[0] = spacing[0] * i;
905  l_tempPoint[1] = spacing[1] * j;
906  l_tempPoint[2] = spacing[2] * k;
907  l_dataTransform->TransformPoint(l_tempPoint, l_tempPoint);
908 
909  //Do stuff if there is no bounding box, or if there is one check if the
910  //point is in the bounding box
911  if (!p_BoundingBox || (p_BoundingBox[0] < l_tempPoint[0] && p_BoundingBox[1] > l_tempPoint[0]
912  && p_BoundingBox[2] < l_tempPoint[1] && p_BoundingBox[3] > l_tempPoint[1] && p_BoundingBox[4]
913  < l_tempPoint[2] && p_BoundingBox[5] > l_tempPoint[2]))
914  {
915  //Loop through the neigbors of the point and see if they are vessel centers
916  //if one of them is it then write the point down
917  for (kk = l_startPosZ; kk <= l_stopPosZ && !l_isNeighborFound; kk += l_offsetConstIJ)
918  {
919  for (jj = l_startPosY; jj <= l_stopPosY && !l_isNeighborFound; jj += l_dimX)
920  {
921  l_kkjjOffset = kk + jj;
922 
923  for (ii = l_startPosX; ii <= l_stopPosX && !l_isNeighborFound; ++ii)
924  {
925  l_kkjjiiOffset = ii + l_kkjjOffset;
926 
927  if (l_offsetI != l_kkjjiiOffset && //ignore if it is the current point
928  *(l_allTheData->GetTuple(l_kkjjiiOffset))) //check if vessel center
929  {
930  l_isNeighborFound = 1;
931  l_dataPoints->InsertNextPoint(l_tempPoint);
932  l_dataCellArray->InsertNextCell(1);
933  l_dataCellArray->InsertCellPoint(l_counts++);
934  }
935  }
936  }
937  }
938  }
939  }
940  }
941  }
942  }
943 
944  p_thePolyData->SetPoints(l_dataPoints);
945  p_thePolyData->SetVerts(l_dataCellArray);
946 
947  return p_thePolyData;
948 }
949 
955 {
956  // characterize the input perturbation in angle-axis form:
957  Vector3D t_delta = linearTransform.matrix().block<3, 1>(0, 3);
958  Eigen::AngleAxisd angleAxis = Eigen::AngleAxisd(linearTransform.matrix().block<3, 3>(0, 0));
959  double angle = angleAxis.angle();
960 
961  QString qualityText = QString("|t_delta|=%1mm, angle=%2*")
962  .arg(t_delta.length(), 6, 'f', 2)
963  .arg(angle / M_PI * 180.0, 6, 'f', 2);
964 
965  if (t_delta.length() > 20 || fabs(angle) > 10 / 180.0 * M_PI)
966  {
967  reportWarning(qualityText);
968  QString text = QString(
969  "The registration matrix' angle-axis representation shows a large shift. Retry registration.");
970  reportWarning(text);
971  }
972  else
973  {
974  report(qualityText);
975  }
976 }
977 
979 {
980  if (!mLastRun)
981  return vtkPolyDataPtr();
982  this->computeDistances();
983  return mLastRun->getDifferenceLines();
984 }
985 
987 {
988  if (!mLastRun)
989  {
990  reportWarning("ICP not initialized");
991  return;
992  }
993 
994  if(!mLastRun->getMovingPoints())
995  {
996  reportWarning("Moving points not set.");
997  }
998 
999  if(!mLastRun->getFixedPoints())
1000  {
1001  reportWarning("Fixed points not set.");
1002  }
1003 }
1004 
1005 
1006 } //namespace cx
vtkSmartPointer< class vtkMatrix4x4 > vtkMatrix4x4Ptr
Definition: cxMathBase.h:37
vtkPolyDataPtr getDifferenceLines()
Lines connecting the moving and fixed data, according to LTS.
vtkPolyDataPtr getMovingPoints()
the moving data (one of target or source, depending on inversion)
vtkAbstractTransformPtr nonLinearRegistration(vtkPointsPtr sortedSourcePoints, vtkPointsPtr sortedTargetPoints)
DoubleBoundingBox3D transform(const Transform3D &m, const DoubleBoundingBox3D &bb)
void print(vtkPointsPtr points)
vtkSmartPointer< class vtkPlanes > vtkPlanesPtr
A mesh data set.
Definition: cxMesh.h:45
double mt_distanceDeltaStopThreshold
ReporterPtr reporter()
Definition: cxReporter.cpp:36
vtkSmartPointer< class vtkLandmarkTransform > vtkLandmarkTransformPtr
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
void computeDistances(ContextPtr context=ContextPtr())
Compute distances between the two datasets.
vtkSmartPointer< class vtkCellArray > vtkCellArrayPtr
vtkPoints * GetSourceLandmarks()
ContextPtr mLastRun
result from last run of execute()
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
vtkSmartPointer< class vtkFloatArray > vtkFloatArrayPtr
vtkSmartPointer< vtkPoints > vtkPointsPtr
double mt_maximumDurationSeconds
vtkSmartPointer< class vtkThinPlateSplineTransform > vtkThinPlateSplineTransformPtr
vtkSmartPointer< class vtkIdList > vtkIdListPtr
vtkSmartPointer< class vtkClipPolyData > vtkClipPolyDataPtr
vtkSmartPointer< class vtkTransform > vtkTransformPtr
Definition: cxMathBase.h:41
void printOutResults(QString fileNamePrefix, vtkGeneralTransformPtr myConcatenation)
boost::shared_ptr< class Data > DataPtr
vtkAbstractTransformPtr linearRegistration(vtkPointsPtr sortedSourcePoints, vtkPointsPtr sortedTargetPoints)
vtkSmartPointer< class vtkMaskPoints > vtkMaskPointsPtr
vtkPolyDataPtr getDifferenceLines()
Lines connecting the moving and fixed data, according to LTS.
double getResultLtsRatio(ContextPtr context=ContextPtr())
void reportWarning(QString msg)
Definition: cxLogger.cpp:70
ContextPtr linearRefineAllLTS(ContextPtr context)
#define CX_LOG_ERROR
Definition: cxLogger.h:99
SeansVesselReg::ContextPtr splitContext(ContextPtr context)
double getResultMetric(ContextPtr context=ContextPtr())
vtkSmartPointer< class vtkSortDataArray > vtkSortDataArrayPtr
vtkSmartPointer< class vtkAbstractTransform > vtkAbstractTransformPtr
ContextPtr createContext(DataPtr source, DataPtr target)
vtkPointsPtr createSortedPoints(vtkIdListPtr sortedIDList, vtkPointsPtr unsortedPoints, int numPoints)
Representation of a floating-point bounding box in 3D. The data are stored as {xmin,xmax,ymin,ymax,zmin,zmax}, in order to simplify communication with vtk.
void notifyPreRegistrationWarnings()
static vtkPolyDataPtr extractPolyData(ImagePtr image, int p_neighborhoodFilterThreshold, double p_BoundingBox[6])
void linearRefine(ContextPtr context)
vtkSmartPointer< vtkPolyData > vtkPolyDataPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
void report(QString msg)
Definition: cxLogger.cpp:69
vtkSmartPointer< class vtkDataArray > vtkDataArrayPtr
vtkSmartPointer< class vtkGeneralTransform > vtkGeneralTransformPtr
void checkQuality(Transform3D linearTransform)
boost::shared_ptr< Context > ContextPtr
vtkPointsPtr transformPoints(vtkPointsPtr input, vtkAbstractTransformPtr transform)
boost::shared_ptr< class Mesh > MeshPtr
static vtkPolyDataPtr convertToPolyData(vtkPointsPtr input)
vtkPolyDataPtr crop(vtkPolyDataPtr input, vtkPolyDataPtr fixed, double margin)
Transform3D getLinearResult(ContextPtr context=ContextPtr())
bool isValid() const
bool initialize(DataPtr source, DataPtr target, QString logPath)
vtkPolyDataPtr getFixedPoints()
the fixed data (one of target or source, depending on inversion)
#define M_PI
Namespace for all CustusX production code.