Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkAnisotropicRegistrationCommon.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
14 #include <mitkPointSet.h>
15 #include <vtkPoints.h>
16 
18  const CovarianceMatrix &sigma_X, const CovarianceMatrix &sigma_Y)
19 {
20  WeightMatrix returnValue;
21 
22  WeightMatrix sum = sigma_X + sigma_Y;
23  vnl_svd<double> svd(sum.GetVnlMatrix());
24 
25  WeightMatrix diag;
26  diag.Fill(0.0);
27  diag[0][0] = 1.0 / sqrt(svd.W(0));
28  diag[1][1] = 1.0 / sqrt(svd.W(1));
29  diag[2][2] = 1.0 / sqrt(svd.W(2));
30 
31  WeightMatrix V; // convert vnl matrix to itk matrix...
32  for (unsigned int i = 0; i < 3; ++i)
33  for (unsigned int j = 0; j < 3; ++j)
34  V[i][j] = svd.V()[i][j];
35 
36  // add weighting matrix for point j1 (corresponding to identity transform)
37  returnValue = V * diag * V.GetTranspose();
38 
39  return returnValue;
40 }
41 
43  vtkPoints *dst,
44  const Rotation &rotation,
45  const Translation &translation)
46 {
47 #pragma omp parallel for
48  for (int i = 0; i < src->GetNumberOfPoints(); ++i)
49  {
50  double p_in[3];
51  double p_out[3];
52  src->GetPoint(i, p_in);
53 
54  for (unsigned int j = 0; j < 3; ++j)
55  {
56  p_out[j] = p_in[0] * rotation[j][0] + p_in[1] * rotation[j][1] + p_in[2] * rotation[j][2] + translation[j];
57  }
58 
59  dst->SetPoint(i, p_out);
60  }
61 }
62 
64  MatrixList &dst,
65  const Rotation &rotation)
66 {
67  const vnl_matrix_fixed<double, 3, 3> rotationT = rotation.GetTranspose();
68 
69 #pragma omp parallel for
70  for (int i = 0; i < static_cast<int>(src.size()); ++i)
71  {
72  dst[i] = rotation * src[i] * rotationT;
73  }
74 }
75 
77  const mitk::PointSet *fixedTargets,
78  const Rotation &rotation,
79  const Translation &translation)
80 {
81  double tre = 0.0;
82 
83  for (int i = 0; i < movingTargets->GetSize(); ++i)
84  {
85  mitk::Point3D pm = movingTargets->GetPoint(i);
86  mitk::Point3D ps = fixedTargets->GetPoint(i);
87 
88  // transform point
89  pm = rotation * pm + translation;
90 
91  const double dist =
92  (ps[0] - pm[0]) * (ps[0] - pm[0]) + (ps[1] - pm[1]) * (ps[1] - pm[1]) + (ps[2] - pm[2]) * (ps[2] - pm[2]);
93 
94  tre += dist;
95  }
96 
97  tre /= movingTargets->GetSize();
98 
99  return sqrt(tre);
100 }
static Matrix3D rotation
static double ComputeTargetRegistrationError(const mitk::PointSet *movingTargets, const mitk::PointSet *fixedTargets, const Rotation &rotation, const Translation &translation)
Compute the target registration error between two point sets.
virtual int GetSize(unsigned int t=0) const
returns the current size of the point-list
static void PropagateMatrices(const MatrixList &src, MatrixList &dst, const Rotation &rotation)
Propagate a list of matrices with a rotation matrix.
Data structure which stores a set of points. Superclass of mitk::Mesh.
Definition: mitkPointSet.h:75
static void TransformPoints(vtkPoints *src, vtkPoints *dst, const Rotation &rotation, const Translation &translation)
Transforms a point cloud with a Rotation and Translation.
static WeightMatrix CalculateWeightMatrix(const CovarianceMatrix &sigma_X, const CovarianceMatrix &sigma_Y)
Method that computes a WeightMatrix with two CovarianceMatrices.
PointType GetPoint(PointIdentifier id, int t=0) const
Get the point with ID id in world coordinates.