14 #include <boost/math/special_functions/fpclassify.hpp>
16 #include "vtkClipPolyData.h"
17 #include "vtkPlanes.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"
28 #include "vtkSortDataArray.h"
29 #include "vtkMaskPoints.h"
30 #include "vtkPointData.h"
31 #include "vtkLandmarkTransform.h"
32 #include "vtkFloatArray.h"
68 reporter()->sendDebug(
"SOURCE: " + source->getUid());
69 reporter()->sendDebug(
"TARGET: " + target->getUid());
87 std::cout <<
"sigma:" <<
mt_sigma << endl;
93 QTime start = QTime::currentTime();
103 report(QString(
"Auto LTS: Found best match using %1\%.").arg(context->mLtsRatio));
119 std::cout << QString(
"\n\nV2V Execution time: %1s").arg(start.secsTo(QTime::currentTime())) << endl;
149 std::vector<int> lts;
161 std::vector<ContextPtr> paths;
164 for (
unsigned i=0; i<lts.size(); ++i)
167 paths.push_back(current);
168 current->mLtsRatio = lts[i];
171 std::cout << QString(
"LTS=%1, metric=%2").arg(current->mLtsRatio).arg(current->mMetric) << std::endl;
176 for (
unsigned i=0; i<paths.size(); ++i)
178 if (paths[i]->mMetric < paths[bestPath]->mMetric)
183 return paths[bestPath];
192 double previousMetric = 1E6;
193 QDateTime t0 = QDateTime::currentDateTime();
197 double difference = context->mMetric - previousMetric;
203 std::cout << QString(
"iteration\t%1\trms:\t%2").arg(iteration).arg(context->mMetric) << std::endl;
210 previousMetric = context->mMetric;
221 retval->mLtsRatio = context->mLtsRatio;
222 retval->mInvertedTransform = context->mInvertedTransform;
225 retval->mTargetPointLocator = context->mTargetPointLocator;
226 retval->mTargetPoints = context->mTargetPoints;
229 retval->mSourcePoints = vtkPointsPtr::New();
230 retval->mSourcePoints->DeepCopy(context->mSourcePoints);
233 retval->mConcatenation = vtkGeneralTransformPtr::New();;
245 if (!source || !target)
251 if (!targetPolyData || !inputSourcePolyData)
257 if (!sourcePolyData->GetNumberOfPoints() || !targetPolyData->GetNumberOfPoints())
259 std::cerr <<
"No data after cropping, failed." << std::endl;
265 context->mConcatenation = vtkGeneralTransformPtr::New();
266 context->mInvertedTransform =
false;
267 context->mMetric = -1;
271 if (sourcePolyData->GetNumberOfPoints() > targetPolyData->GetNumberOfPoints())
275 std::cout <<
"inverted vessel reg" << std::endl;
276 context->mInvertedTransform =
true;
277 std::swap(sourcePolyData, targetPolyData);
280 vtkIdType numPoints = sourcePolyData->GetNumberOfPoints();
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;
289 context->mTargetPoints = targetPolyData;
290 context->mTargetPointLocator = vtkCellLocatorPtr::New();
291 context->mTargetPointLocator->SetDataSet(targetPolyData);
292 context->mTargetPointLocator->SetNumberOfCellsPerBucket(1);
293 context->mTargetPointLocator->BuildLocator();
296 context->mSourcePoints = vtkPointsPtr::New();
297 context->mSourcePoints->DeepCopy(sourcePolyData->GetPoints());
317 if (context->mSortedSourcePoints || context->mSortedTargetPoints)
321 int numPoints = context->mSourcePoints->GetNumberOfPoints();
323 int nb_points = ((int) (numPoints * context->mLtsRatio) / 100);
330 closestPoint->SetNumberOfPoints(numPoints);
334 residuals->SetNumberOfValues(numPoints);
337 IdList->SetNumberOfIds(numPoints);
338 double total_distance = 0;
339 double distanceSquared = 0;
341 for (
int i = 0; i < numPoints; ++i)
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))
351 std::cout <<
"nan found during findClosestPoint!" << std::endl;
353 context->mMetric = 1E6;
357 residuals->InsertValue(i, distanceSquared);
358 IdList->InsertId(i, i);
359 total_distance += sqrt(distanceSquared);
363 context->mMetric = total_distance / numPoints;
366 sort->Sort(residuals, IdList);
367 context->mSortedSourcePoints = this->
createSortedPoints(IdList, context->mSourcePoints, nb_points);
368 context->mSortedTargetPoints = this->
createSortedPoints(IdList, closestPoint, nb_points);
387 if (!context->mSortedSourcePoints || !context->mSortedTargetPoints)
392 context->mTransform =
linearRegistration(context->mSortedSourcePoints, context->mSortedTargetPoints);
396 context->mTransform =
nonLinearRegistration(context->mSortedSourcePoints, context->mSortedTargetPoints);
401 context->mSourcePoints = this->
transformPoints(context->mSourcePoints, context->mTransform);
408 context->mConcatenation->Concatenate(context->mTransform);
420 retval->SetNumberOfPoints(numPoints);
422 double temp_point[3];
424 for (
int i = 0; i < numPoints; ++i)
426 vtkIdType index = sortedIDList->GetId(i);
427 unsortedPoints->GetPoint(index, temp_point);
428 retval->SetPoint(i, temp_point);
441 int N = input->GetNumberOfPoints();
443 retval->SetNumberOfPoints(N);
447 for (
int i = 0; i < N; ++i)
449 transform->InternalTransformPoint(input->GetPoint(i), temp);
450 retval->SetPoint(i, temp);
463 lmt->SetSourceLandmarks(sortedSourcePoints);
464 lmt->SetTargetLandmarks(sortedTargetPoints);
465 lmt->SetModeToRigidBody();
480 mask1->SetInputData(tpsSourcePolyData);
484 mask2->SetInputData(tpsTargetPolyData);
490 tps->SetSourceLandmarks(mask1->GetOutput()->GetPoints());
491 tps->SetTargetLandmarks(mask2->GetOutput()->GetPoints());
494 QTime start = QTime::currentTime();
496 std::cout << QString(
"\nV2V NL time: %1s").arg(start.secsTo(QTime::currentTime())) << endl;
505 CX_LOG_ERROR(QString(
"Can't execute: %1 is not set").arg(
id));
510 MeshPtr mesh = boost::dynamic_pointer_cast<
Mesh>(data);
521 CX_LOG_ERROR(QString(
"Can't execute: %1=%2 is not a mesh").arg(
id, data->getName()));
525 vtkPolyDataPtr poly = mesh->getTransformedPolyData(mesh->get_rMd());
527 if (!poly || !poly->GetNumberOfPoints())
529 CX_LOG_ERROR(QString(
"Can't execute: %1=%2 is empty").arg(
id, data->getName()));
550 if (mInvertedTransform)
556 return mTargetPoints;
566 for (
int i = 0; i < this->mSortedSourcePoints->GetNumberOfPoints(); ++i)
568 verts->InsertNextPoint(this->mSortedSourcePoints->GetPoint(i));
569 verts->InsertNextPoint(this->mSortedTargetPoints->GetPoint(i));
571 vtkIdType connectivity[2];
572 connectivity[0] = 2 * i;
573 connectivity[1] = 2 * i + 1;
574 retval->InsertNextCell(VTK_LINE, 2, connectivity);
576 retval->SetPoints(verts);
584 int N = input->GetNumberOfPoints();
586 for (
int i=0; i<N ; ++i)
588 cellArray->InsertNextCell(1);
589 cellArray->InsertCellPoint(i);
593 retval->SetPoints(input);
594 retval->SetVerts(cellArray);
601 for (
int q = 0; q < points->GetNumberOfPoints(); ++q)
604 std::cout << q <<
"\t" << p[0] <<
" " << p[1] <<
" " << p[2] <<
" " << std::endl;
610 print(data->GetPoints());
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;
631 box->SetBounds(targetBounds);
635 clipper->SetInputData(input);
636 clipper->SetClipFunction(box);
637 clipper->SetInsideOut(
true);
640 int oldSource = input->GetPoints()->GetNumberOfPoints();
642 int clippedSource = 0;
644 if (clipper->GetOutput()->GetPoints())
646 clippedSource = clipper->GetOutput()->GetPoints()->GetNumberOfPoints();
649 if (clippedSource < oldSource)
651 double ratio = double(oldSource - clippedSource) / double(oldSource);
653 std::cout <<
"Removed " << ratio * 100 <<
"%" <<
" of the source data. Outside the target data bounds." << std::endl;
656 return clipper->GetOutput();
664 Transform3D retval = this->getLinearTransform(context->mConcatenation);
665 if (context->mInvertedTransform)
666 retval = retval.inv();
677 return context->mMetric;
684 return context->mLtsRatio;
693 for (
int i = 0; i < concatenation->GetNumberOfConcatenatedTransforms(); ++i)
696 if (!LM_i.GetPointer())
702 return M_acc.inverse();
712 l_tempMatrix->DeepCopy(((vtkLandmarkTransform*) myConcatenation->GetConcatenatedTransform(0))->GetMatrix());
714 l_resultMatrix->Identity();
715 for (
int i = 1; i < myConcatenation->GetNumberOfConcatenatedTransforms(); ++i)
717 vtkMatrix4x4::Multiply4x4(l_tempMatrix,
718 ((vtkLandmarkTransform*) myConcatenation->GetConcatenatedTransform(i))->GetMatrix(), l_resultMatrix);
719 l_tempMatrix->DeepCopy(l_resultMatrix);
723 std::cout <<
"Filenameprefix: " << fileNamePrefix << std::endl;
727 std::string nonLinearFile = fileNamePrefix.toStdString();
728 nonLinearFile +=
"--NonLinear";
729 nonLinearFile +=
".txt";
731 std::string linearFile = fileNamePrefix.toStdString();
732 linearFile +=
"--Linear";
733 linearFile +=
".txt";
736 std::cout <<
"Writing Results to " << nonLinearFile <<
" and " << linearFile << std::endl;
740 ofstream file_out(nonLinearFile.c_str());
746 file_out <<
"MNI Transform File\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;
754 const double*
const * theWarpMatrix = l_theWarpTransform->GetWarpMatrix();
756 file_out <<
"Points = " << std::endl;
757 for (
int i = 0; i < n; i++)
759 l_theWarpTransform->GetSourceLandmarks()->GetPoint(i, point);
760 file_out << point[0] <<
" " << point[1] <<
" " << point[2];
762 file_out <<
";" << std::endl;
764 file_out << std::endl;
767 file_out <<
"Displacements = " << std::endl;
768 for (
int i = 0; i < n + 4; i++)
770 file_out << theWarpMatrix[i][0] <<
" " << theWarpMatrix[i][1] <<
" " << theWarpMatrix[i][2];
772 file_out <<
";" << std::endl;
774 file_out << std::endl;
781 ofstream file_out2(linearFile.c_str());
784 file_out2 <<
"MNI Transform File\n" << std::endl;
787 file_out2 <<
"Transform_Type = Linear;" << std::endl;
788 file_out2 <<
"Linear_Transform = ";
789 file_out2 << std::endl;
791 for (
int i = 0; i < 3; i++)
793 for (
int j = 0; j < 4; j++)
796 file_out2 << l_resultMatrix->GetElement(i, j);
801 file_out2 <<
";" << std::endl;
803 file_out2 << std::endl;
814 double p_BoundingBox[6])
819 image->getBaseVtkImageData()->GetDimensions(l_dimensions);
820 Vector3D spacing(image->getBaseVtkImageData()->GetSpacing());
824 l_dataTransform->SetMatrix(image->get_rMd().getVtkMatrix());
826 int l_startPosX, l_startPosY, l_startPosZ;
827 int l_stopPosX, l_stopPosY, l_stopPosZ;
828 int i, j, k, ii, jj, kk, l_counts = 0;
830 bool l_isNeighborFound;
833 int l_dimX = l_dimensions[0];
834 int l_dimY = l_dimensions[1];
835 int l_dimZ = l_dimensions[2];
837 double l_tempPoint[3];
842 int l_kkjjOffset = 0;
843 int l_kkjjiiOffset = 0;
844 int l_offsetConstIJ = l_dimX * l_dimY;
846 vtkDataArrayPtr l_allTheData = image->getBaseVtkImageData()->GetPointData()->GetScalars();
852 for (k = 0; k < l_dimZ; ++k)
855 l_startPosZ = k - p_neighborhoodFilterThreshold;
859 l_startPosZ *= l_offsetConstIJ;
861 l_stopPosZ = k + p_neighborhoodFilterThreshold;
862 if (l_stopPosZ >= l_dimZ)
863 l_stopPosZ = (l_dimZ - 1) * l_offsetConstIJ;
865 l_stopPosZ *= l_offsetConstIJ;
867 l_offsetK = k * l_offsetConstIJ;
869 for (j = 0; j < l_dimY; ++j)
872 l_startPosY = j - p_neighborhoodFilterThreshold;
876 l_startPosY *= l_dimX;
878 l_stopPosY = j + p_neighborhoodFilterThreshold;
879 if (l_stopPosY >= l_dimY)
880 l_stopPosY = (l_dimY - 1) * l_dimX;
882 l_stopPosY *= l_dimX;
884 l_offsetJ = l_offsetK + (j * l_dimX);
886 for (i = 0; i < l_dimX; ++i)
889 l_startPosX = i - p_neighborhoodFilterThreshold;
893 l_stopPosX = i + p_neighborhoodFilterThreshold;
894 if (l_stopPosX >= l_dimX)
895 l_stopPosX = l_dimX - 1;
897 l_offsetI = l_offsetJ + i;
898 l_isNeighborFound = 0;
901 if (*(l_allTheData->GetTuple(l_offsetI)))
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);
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]))
917 for (kk = l_startPosZ; kk <= l_stopPosZ && !l_isNeighborFound; kk += l_offsetConstIJ)
919 for (jj = l_startPosY; jj <= l_stopPosY && !l_isNeighborFound; jj += l_dimX)
921 l_kkjjOffset = kk + jj;
923 for (ii = l_startPosX; ii <= l_stopPosX && !l_isNeighborFound; ++ii)
925 l_kkjjiiOffset = ii + l_kkjjOffset;
927 if (l_offsetI != l_kkjjiiOffset &&
928 *(l_allTheData->GetTuple(l_kkjjiiOffset)))
930 l_isNeighborFound = 1;
931 l_dataPoints->InsertNextPoint(l_tempPoint);
932 l_dataCellArray->InsertNextCell(1);
933 l_dataCellArray->InsertCellPoint(l_counts++);
944 p_thePolyData->SetPoints(l_dataPoints);
945 p_thePolyData->SetVerts(l_dataCellArray);
947 return p_thePolyData;
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();
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);
965 if (t_delta.length() > 20 || fabs(angle) > 10 / 180.0 *
M_PI)
968 QString text = QString(
969 "The registration matrix' angle-axis representation shows a large shift. Retry registration.");
983 return mLastRun->getDifferenceLines();
vtkSmartPointer< class vtkMatrix4x4 > vtkMatrix4x4Ptr
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
double mt_distanceDeltaStopThreshold
vtkSmartPointer< class vtkLandmarkTransform > vtkLandmarkTransformPtr
int mt_maximumNumberOfIterations
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
ContextPtr mLastRun
result from last run of execute()
boost::shared_ptr< class Image > ImagePtr
vtkSmartPointer< class vtkFloatArray > vtkFloatArrayPtr
bool mInvertedTransform
the calculated registration goes from target to source instead of source to target ...
double mt_maximumDurationSeconds
vtkSmartPointer< class vtkThinPlateSplineTransform > vtkThinPlateSplineTransformPtr
vtkPolyDataPtr mTargetPoints
input: target data
vtkSmartPointer< class vtkIdList > vtkIdListPtr
vtkSmartPointer< class vtkClipPolyData > vtkClipPolyDataPtr
vtkSmartPointer< class vtkTransform > vtkTransformPtr
int mt_singlePointThreshold
void printOutResults(QString fileNamePrefix, vtkGeneralTransformPtr myConcatenation)
boost::shared_ptr< class Data > DataPtr
vtkAbstractTransformPtr linearRegistration(vtkPointsPtr sortedSourcePoints, vtkPointsPtr sortedTargetPoints)
vtkSmartPointer< class vtkPolyData > vtkPolyDataPtr
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)
ContextPtr linearRefineAllLTS(ContextPtr context)
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)
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
bool performOneRegistration()
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 initialize(DataPtr source, DataPtr target, QString logPath)
vtkPolyDataPtr getFixedPoints()
the fixed data (one of target or source, depending on inversion)
vtkSmartPointer< class vtkPoints > vtkPointsPtr
vtkPointsPtr mSourcePoints
input: current source data, modified according to last iteration