ANIMA  4.0
animaLinearTransformEstimationTools.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <animaMatrixLogExp.h>
5 
6 #include <vnl/vnl_inverse.h>
7 #include <itkVector.h>
8 
9 #include <itkMatrixOffsetTransformBase.h>
10 #include <itkTranslationTransform.h>
11 #include <itkContinuousIndex.h>
12 
13 namespace anima
14 {
15 
16 template <class TInput, class TScalarType, unsigned int NDimensions>
17 void computeTranslationLSWFromTranslations(std::vector < itk::Point<TInput,NDimensions> > &inputOrigins,
18  std::vector < itk::Point<TInput,NDimensions> > &inputTransformed,
19  std::vector <TInput> &weights,
20  typename itk::AffineTransform<TScalarType,NDimensions>::Pointer &resultTransform)
21 {
22  unsigned int nbPts = inputOrigins.size();
23  itk::Point <TInput,NDimensions> barX, barY;
24  // See Pennec PhD
25  vnl_matrix <TInput> AMatrix(4,4,0);
26  vnl_matrix <TInput> tmpMatrix(4,4);
27 
28  for (unsigned int j = 0;j < NDimensions;++j)
29  {
30  barX[j] = 0;
31  barY[j] = 0;
32  }
33 
34  double sumWeights = 0;
35  for (unsigned int i = 0;i < nbPts;++i)
36  {
37  for (unsigned int j = 0;j < NDimensions;++j)
38  {
39  barX[j] += weights[i] * inputOrigins[i][j];
40  barY[j] += weights[i] * inputTransformed[i][j];
41  }
42 
43  sumWeights += weights[i];
44  }
45 
46  for (unsigned int j = 0;j < NDimensions;++j)
47  {
48  barX[j] /= sumWeights;
49  barY[j] /= sumWeights;
50  }
51 
52  itk::Vector <TScalarType,NDimensions> translationPart;
53  for (unsigned int i = 0;i < NDimensions;++i)
54  translationPart[i] = barY[i] - barX[i];
55 
56  resultTransform = itk::AffineTransform<TScalarType,NDimensions>::New();
57 
58  resultTransform->SetIdentity();
59  resultTransform->SetOffset(translationPart);
60 }
61 
62 template <class TInput, class TScalarType, unsigned int NDimensions>
63 void computeRigidLSWFromTranslations(std::vector < itk::Point<TInput,NDimensions> > &inputOrigins,
64  std::vector < itk::Point<TInput,NDimensions> > &inputTransformed,
65  std::vector <TInput> &weights,
66  typename itk::AffineTransform<TScalarType,NDimensions>::Pointer &resultTransform)
67 {
68  unsigned int nbPts = inputOrigins.size();
69  itk::Point <TInput, NDimensions> barX, barY;
70  // See Pennec PhD
71  vnl_matrix <TInput> AMatrix(4,4,0);
72  vnl_matrix <TInput> tmpMatrix(4,4);
73 
74  for (unsigned int j = 0;j < NDimensions;++j)
75  {
76  barX[j] = 0;
77  barY[j] = 0;
78  }
79 
80  double sumWeights = 0;
81  for (unsigned int i = 0;i < nbPts;++i)
82  {
83  for (unsigned int j = 0;j < NDimensions;++j)
84  {
85  barX[j] += weights[i] * inputOrigins[i][j];
86  barY[j] += weights[i] * inputTransformed[i][j];
87  }
88 
89  sumWeights += weights[i];
90  }
91 
92  for (unsigned int j = 0;j < NDimensions;++j)
93  {
94  barX[j] /= sumWeights;
95  barY[j] /= sumWeights;
96  }
97 
98  itk::Vector <TInput,NDimensions> xVector, yVector;
99  for (unsigned int i = 0;i < nbPts;++i)
100  {
101  xVector = inputOrigins[i] - barX;
102  yVector = inputTransformed[i] - barY;
103  anima::pairingToQuaternion(xVector,yVector,tmpMatrix);
104  AMatrix += weights[i]*tmpMatrix.transpose()*tmpMatrix;
105  }
106 
107  itk::SymmetricEigenAnalysis < vnl_matrix <TInput>, vnl_diag_matrix<TInput>, vnl_matrix <TInput> > eigenSystem(4);
108  vnl_matrix <double> eVec(4,4);
109  vnl_diag_matrix <double> eVals(4);
110 
111  eigenSystem.SetOrderEigenValues(true);
112  eigenSystem.ComputeEigenValuesAndVectors(AMatrix, eVals, eVec);
113 
114  vnl_matrix <TScalarType> rotationMatrix = anima::computeRotationFromQuaternion<TInput,TScalarType>(eVec.get_row(0));
115 
116  itk::Vector <TScalarType,NDimensions> translationPart;
117  for (unsigned int i = 0;i < NDimensions;++i)
118  {
119  translationPart[i] = barY[i];
120  for (unsigned int j = 0;j < NDimensions;++j)
121  translationPart[i] -= rotationMatrix(i,j)*barX[j];
122  }
123 
124  resultTransform = itk::AffineTransform<TScalarType,NDimensions>::New();
125 
126  resultTransform->SetMatrix(rotationMatrix);
127  resultTransform->SetOffset(translationPart);
128 }
129 
130 template <class TInput, class TScalarType, unsigned int NDimensions>
131 itk::Point <TInput, NDimensions> computeAnisotropSimLSWFromTranslations(std::vector < itk::Point<TInput, NDimensions> > &inputOrigins,
132  std::vector < itk::Point<TInput, NDimensions> > &inputTransformed,
133  std::vector <TInput> &weights,
134  typename itk::AffineTransform<TScalarType, NDimensions>::Pointer &resultTransform,
135  vnl_matrix <double> &UMatrix)
136 {
137  unsigned int nbPts = inputOrigins.size();
138  itk::Point <TInput, NDimensions> barX, barY;
139  // See Pennec PhD
140  vnl_matrix <TInput> AMatrix(4, 4, 0);
141  vnl_matrix <TInput> covInputOrigins(NDimensions, NDimensions, 0);
142  vnl_matrix <TInput> tmpMatrix(4, 4);
143 
144  for (unsigned int j = 0; j < NDimensions; ++j)
145  {
146  barX[j] = 0;
147  barY[j] = 0;
148  }
149 
150  double sumWeights = 0;
151  for (unsigned int i = 0; i < nbPts; ++i)
152  {
153  for (unsigned int j = 0; j < NDimensions; ++j)
154  {
155  barX[j] += weights[i] * inputOrigins[i][j];
156  barY[j] += weights[i] * inputTransformed[i][j];
157  }
158  sumWeights += weights[i];
159  }
160 
161  for (unsigned int j = 0; j < NDimensions; ++j)
162  {
163  barX[j] /= sumWeights;
164  barY[j] /= sumWeights;
165  }
166 
167  vnl_vector_fixed <TInput, NDimensions> xVector, xiVector, yVector;
168  unsigned int iter = 0;
169  double eps = 1.0e-8;
170  double diffRot = eps + 1;
171  double diffScal = eps + 1;
172  vnl_diag_matrix <double> scal(NDimensions);
173  scal.fill_diagonal(1);
174  vnl_diag_matrix <double> oldscal(NDimensions);
175  vnl_vector <double> q, oldq;
176  std::vector<vnl_matrix<double>> sumdScal(NDimensions);
177  for (unsigned int j = 0; j < NDimensions; ++j)
178  sumdScal[j].set_size(NDimensions + 1, NDimensions + 1);
179  vnl_vector < double > sumx2(3, 0);
180  itk::SymmetricEigenAnalysis < vnl_matrix <TInput>, vnl_diag_matrix<TInput>, vnl_matrix <TInput> > eigenSystem2(4);
181  vnl_matrix <double> eVec(4, 4);
182  vnl_diag_matrix <double> eVals(4);
183  eigenSystem2.SetOrderEigenValues(true);
184 
185  // Loop over rotation and scalings optimization
186  while (((diffRot > eps) || (diffScal > eps)) && (iter < 100))
187  {
188  // optimize rotation
189  AMatrix.fill(0);
190 
191  for (unsigned int i = 0;i < nbPts;++i)
192  {
193  xiVector = scal*UMatrix.transpose()*(inputOrigins[i].GetVnlVector() - barX.GetVnlVector());
194  yVector = inputTransformed[i].GetVnlVector() - barY.GetVnlVector();
195  anima::pairingToQuaternion(xiVector, yVector, tmpMatrix);
196  AMatrix += weights[i] * tmpMatrix.transpose()*tmpMatrix;
197  }
198 
199  eigenSystem2.ComputeEigenValuesAndVectors(AMatrix, eVals, eVec);
200 
201  q = eVec.get_row(0);
202 
203  if (iter > 0)
204  diffRot = (q - oldq).two_norm();
205  oldq = q;
206 
207  // optimize scalings
208  for (unsigned int j = 0;j < NDimensions;++j)
209  sumdScal[j].fill(0);
210  sumx2.fill(0);
211 
212  for (unsigned int i = 0;i < nbPts;++i)
213  {
214  yVector = inputTransformed[i].GetVnlVector() - barY.GetVnlVector();
215  xVector = UMatrix.transpose()*(inputOrigins[i].GetVnlVector() - barX.GetVnlVector());
216  for (unsigned int j = 0;j < NDimensions;++j)
217  {
218  anima::pairingToQuaternionScalingsDerivative(xVector, yVector, tmpMatrix, j);
219  sumdScal[j] += weights[i] * tmpMatrix;
220  sumx2[j] += weights[i] * xVector[j] * xVector[j] ;
221  }
222  }
223 
224  scal.fill(0);
225  for (unsigned int j = 0;j < NDimensions;++j)
226  {
227  for (unsigned int k = 0;k < NDimensions + 1;++k)
228  {
229  for (unsigned int l = 0;l < NDimensions + 1;++l)
230  scal[j] += (1 / sumx2[j]) * q[k] * sumdScal[j](k, l) * q[l];
231  }
232  }
233 
234  if (iter > 0)
235  diffScal = (scal.get_diagonal() - oldscal.get_diagonal()).two_norm();
236  oldscal = scal;
237 
238  iter++;
239  }
240 
241  vnl_matrix <TScalarType> rotationMatrix = anima::computeRotationFromQuaternion<TInput, TScalarType>(q);
242  vnl_matrix <TScalarType> linearPartMatrix = rotationMatrix * scal * UMatrix.transpose();
243 
244  itk::Vector <TScalarType, NDimensions> translationPart;
245  for (unsigned int i = 0;i < NDimensions;++i)
246  {
247  translationPart[i] = barY[i];
248  for (unsigned int j = 0;j < NDimensions;++j)
249  translationPart[i] -= linearPartMatrix(i,j)*barX[j];
250  }
251 
252  resultTransform = itk::AffineTransform <TScalarType, NDimensions>::New();
253 
254  resultTransform->SetMatrix(linearPartMatrix);
255  resultTransform->SetOffset(translationPart);
256 
257  return barX;
258 }
259 
260 template <class TInput, class TScalarType, unsigned int NDimensions>
261 void computeLogEuclideanAverage(std::vector < vnl_matrix <TInput> > &inputTransforms, std::vector <TInput> &weights,
262  typename itk::AffineTransform<TScalarType,NDimensions>::Pointer &resultTransform)
263 {
264  unsigned int nbPts = inputTransforms.size();
265 
266  vnl_matrix <TInput> resultMatrix(NDimensions+1,NDimensions+1,0);
267  vnl_matrix <TInput> tmpMatrix;
268  double sumWeights = 0;
269  for (unsigned int i = 0;i < nbPts;++i)
270  {
271  sumWeights += weights[i];
272  resultMatrix += weights[i] * inputTransforms[i];
273  }
274 
275  resultMatrix /= sumWeights;
276 
277  resultMatrix = anima::GetExponential(resultMatrix);
278 
279  resultTransform = itk::AffineTransform<TScalarType,NDimensions>::New();
280 
281  vnl_matrix <TScalarType> affinePart(NDimensions,NDimensions);
282  itk::Vector <TScalarType,NDimensions> translationPart;
283 
284  for (unsigned int i = 0;i < NDimensions;++i)
285  {
286  translationPart[i] = resultMatrix(i,NDimensions);
287  for (unsigned int j = 0;j < NDimensions;++j)
288  affinePart(i,j) = resultMatrix(i,j);
289  }
290 
291  resultTransform->SetMatrix(affinePart);
292  resultTransform->SetOffset(translationPart);
293 }
294 
295 template <class TInput, class TScalarType, unsigned int NDimensions>
296 itk::Point <TInput, NDimensions> computeAffineLSWFromTranslations(std::vector < itk::Point<TInput,NDimensions> > &inputOrigins,
297  std::vector < itk::Point<TInput,NDimensions> > &inputTransformed,
298  std::vector <TInput> &weights,
299  typename itk::AffineTransform<TScalarType,NDimensions>::Pointer &resultTransform)
300 {
301  unsigned int nbPts = inputOrigins.size();
302  itk::Point <TInput,NDimensions> barX, barY;
303  // See Pennec PhD
304  vnl_matrix <TInput> SigmaXXMatrix(NDimensions,NDimensions,0);
305  vnl_matrix <TInput> SigmaYXMatrix(NDimensions,NDimensions,0);
306 
307  for (unsigned int j = 0;j < NDimensions;++j)
308  {
309  barX[j] = 0;
310  barY[j] = 0;
311  }
312 
313  double sumWeights = 0;
314  for (unsigned int i = 0;i < nbPts;++i)
315  {
316  for (unsigned int j = 0;j < NDimensions;++j)
317  {
318  barX[j] += weights[i] * inputOrigins[i][j];
319  barY[j] += weights[i] * inputTransformed[i][j];
320  }
321 
322  sumWeights += weights[i];
323  }
324 
325  for (unsigned int j = 0;j < NDimensions;++j)
326  {
327  barX[j] /= sumWeights;
328  barY[j] /= sumWeights;
329  }
330 
331  for (unsigned int i = 0;i < nbPts;++i)
332  {
333  for (unsigned int j = 0;j < NDimensions;++j)
334  for (unsigned int k = 0;k < NDimensions;++k)
335  {
336  SigmaXXMatrix(j,k) += weights[i]*(inputOrigins[i][j] - barX[j])*(inputOrigins[i][k] - barX[k]);
337  SigmaYXMatrix(j,k) += weights[i]*(inputTransformed[i][j] - barY[j])*(inputOrigins[i][k] - barX[k]);
338  }
339  }
340 
341  vnl_matrix <TInput> affineMatrix = SigmaYXMatrix * vnl_inverse (SigmaXXMatrix);
342 
343  vnl_matrix <TScalarType> outMatrix(NDimensions,NDimensions);
344  itk::Vector <TScalarType,NDimensions> translationPart;
345  for (unsigned int i = 0;i < NDimensions;++i)
346  {
347  translationPart[i] = barY[i];
348  for (unsigned int j = 0;j < NDimensions;++j)
349  {
350  outMatrix(i,j) = affineMatrix(i,j);
351  translationPart[i] -= affineMatrix(i,j)*barX[j];
352  }
353  }
354 
355  resultTransform = itk::AffineTransform<TScalarType,NDimensions>::New();
356 
357  resultTransform->SetMatrix(outMatrix);
358  resultTransform->SetOffset(translationPart);
359 
360  return barX;
361 }
362 
363 template <class TInput, class TOutput> vnl_matrix <TOutput> computeRotationFromQuaternion(vnl_vector <TInput> eigenVector)
364 {
365  double normVector = 0;
366 
367  for (unsigned int i = 0;i < eigenVector.size();++i)
368  normVector += eigenVector[i]*eigenVector[i];
369 
370  vnl_matrix <TOutput> resVal(3,3,0);
371 
372  resVal(0,0) = (1.0/normVector) * (eigenVector[0]*eigenVector[0] + eigenVector[1]*eigenVector[1] - eigenVector[2]*eigenVector[2] - eigenVector[3]*eigenVector[3]);
373  resVal(0,1) = (2.0/normVector) * (eigenVector[1]*eigenVector[2] - eigenVector[0]*eigenVector[3]);
374  resVal(0,2) = (2.0/normVector) * (eigenVector[1]*eigenVector[3] + eigenVector[0]*eigenVector[2]);
375 
376  resVal(1,0) = (2.0/normVector) * (eigenVector[1]*eigenVector[2] + eigenVector[0]*eigenVector[3]);
377  resVal(1,1) = (1.0/normVector) * (eigenVector[0]*eigenVector[0] - eigenVector[1]*eigenVector[1] + eigenVector[2]*eigenVector[2] - eigenVector[3]*eigenVector[3]);
378  resVal(1,2) = (2.0/normVector) * (eigenVector[2]*eigenVector[3] - eigenVector[0]*eigenVector[1]);
379 
380  resVal(2,0) = (2.0/normVector) * (eigenVector[1]*eigenVector[3] - eigenVector[0]*eigenVector[2]);
381  resVal(2,1) = (2.0/normVector) * (eigenVector[2]*eigenVector[3] + eigenVector[0]*eigenVector[1]);
382  resVal(2,2) = (1.0/normVector) * (eigenVector[0]*eigenVector[0] - eigenVector[1]*eigenVector[1] - eigenVector[2]*eigenVector[2] + eigenVector[3]*eigenVector[3]);
383 
384  return resVal;
385 }
386 
387 template <class TInput, class TOutput, unsigned int NDimensions>
388 void pairingToQuaternion(const vnl_vector_fixed <TInput,NDimensions> &inputPoint, const vnl_vector_fixed <TInput,NDimensions> &inputTransformedPoint,
389  vnl_matrix <TOutput> &outputMatrix)
390 {
391  return pairingToQuaternion(inputPoint,inputTransformedPoint,outputMatrix,NDimensions);
392 }
393 
394 template <class TInput, class TOutput, unsigned int NDimensions>
395 void pairingToQuaternion(const itk::Vector <TInput,NDimensions> &inputPoint, const itk::Vector <TInput,NDimensions> &inputTransformedPoint,
396  vnl_matrix <TOutput> &outputMatrix)
397 {
398  return pairingToQuaternion(inputPoint,inputTransformedPoint,outputMatrix,NDimensions);
399 }
400 
401 template <class PointType, class TOutput>
402 void pairingToQuaternion(const PointType &inputPoint, const PointType &inputTransformedPoint, vnl_matrix <TOutput> &outputMatrix, unsigned int ndim)
403 {
404  outputMatrix.set_size(ndim+1,ndim+1);
405  outputMatrix.fill(0);
406 
407  for (unsigned int i = 0;i < ndim;++i)
408  {
409  switch (i)
410  {
411  case 0:
412  outputMatrix(0,1) = inputPoint[i] - inputTransformedPoint[i];
413  outputMatrix(2,3) = - (inputPoint[i] + inputTransformedPoint[i]);
414  outputMatrix(1,0) = - outputMatrix(0,1);
415  outputMatrix(3,2) = - outputMatrix(2,3);
416  break;
417 
418  case 1:
419  outputMatrix(0,2) = inputPoint[i] - inputTransformedPoint[i];
420  outputMatrix(1,3) = inputPoint[i] + inputTransformedPoint[i];
421  outputMatrix(2,0) = - outputMatrix(0,2);
422  outputMatrix(3,1) = - outputMatrix(1,3);
423  break;
424 
425  case 2:
426  outputMatrix(0,3) = inputPoint[i] - inputTransformedPoint[i];
427  outputMatrix(1,2) = - (inputPoint[i] + inputTransformedPoint[i]);
428  outputMatrix(2,1) = - outputMatrix(1,2);
429  outputMatrix(3,0) = - outputMatrix(0,3);
430  break;
431 
432  default:
433  break;
434  }
435  }
436 }
437 
438 template <class TInput, class TOutput, unsigned int NDimensions>
439 void pairingToQuaternionScalingsDerivative(const vnl_vector_fixed <TInput, NDimensions> &inputPoint, const vnl_vector_fixed <TInput, NDimensions> &inputTransformedPoint,
440  vnl_matrix <TOutput> &outputMatrix, const int &dimScal)
441 {
442  return pairingToQuaternionScalingsDerivative(inputPoint, inputTransformedPoint, outputMatrix, NDimensions, dimScal);
443 }
444 
445 template <class TInput, class TOutput, unsigned int NDimensions>
446 void pairingToQuaternionScalingsDerivative(const itk::Vector <TInput, NDimensions> &inputPoint, const itk::Vector <TInput, NDimensions> &inputTransformedPoint,
447  vnl_matrix <TOutput> &outputMatrix, const int &dimScal)
448 {
449  return pairingToQuaternionScalingsDerivative(inputPoint, inputTransformedPoint, outputMatrix, NDimensions, dimScal);
450 }
451 
452 template <class PointType, class TOutput>
453 void pairingToQuaternionScalingsDerivative(const PointType &inputPoint, const PointType &inputTransformedPoint, vnl_matrix <TOutput> &outputMatrix, unsigned int ndim, const int &dimScal)
454 {
455  outputMatrix.set_size(ndim + 1, ndim + 1);
456  outputMatrix.fill(0);
457 
458  switch (dimScal)
459  {
460  case 0:
461  outputMatrix(0, 2) = -inputTransformedPoint[2];
462  outputMatrix(0, 3) = inputTransformedPoint[1];
463  outputMatrix(1, 2) = inputTransformedPoint[1];
464  outputMatrix(1, 3) = inputTransformedPoint[2];
465  outputMatrix += outputMatrix.transpose();
466  outputMatrix(0, 0) = inputTransformedPoint[0];
467  outputMatrix(1, 1) = inputTransformedPoint[0];
468  outputMatrix(2, 2) = -inputTransformedPoint[0];
469  outputMatrix(3, 3) = -inputTransformedPoint[0];
470  outputMatrix *= inputPoint[0];
471  break;
472 
473  case 1:
474  outputMatrix(0, 1) = inputTransformedPoint[2];
475  outputMatrix(0, 3) = -inputTransformedPoint[0];
476  outputMatrix(1, 2) = inputTransformedPoint[0];
477  outputMatrix(2, 3) = inputTransformedPoint[2];
478  outputMatrix += outputMatrix.transpose();
479  outputMatrix(0, 0) = inputTransformedPoint[1];
480  outputMatrix(1, 1) = -inputTransformedPoint[1];
481  outputMatrix(2, 2) = inputTransformedPoint[1];
482  outputMatrix(3, 3) = -inputTransformedPoint[1];
483  outputMatrix *= inputPoint[1];
484  break;
485 
486  case 2:
487  outputMatrix(0, 1) = -inputTransformedPoint[1];
488  outputMatrix(0, 2) = inputTransformedPoint[0];
489  outputMatrix(1, 3) = inputTransformedPoint[0];
490  outputMatrix(2, 3) = inputTransformedPoint[1];
491  outputMatrix += outputMatrix.transpose();
492  outputMatrix(0, 0) = inputTransformedPoint[2];
493  outputMatrix(1, 1) = -inputTransformedPoint[2];
494  outputMatrix(2, 2) = -inputTransformedPoint[2];
495  outputMatrix(3, 3) = inputTransformedPoint[2];
496  outputMatrix *= inputPoint[2];
497  break;
498 
499  default:
500  break;
501  }
502 
503 }
504 
505 } // end of namespace anima
void pairingToQuaternionScalingsDerivative(const vnl_vector_fixed< TInput, NDimensions > &inputPoint, const vnl_vector_fixed< TInput, NDimensions > &inputTransformedPoint, vnl_matrix< TOutput > &outputMatrix, const int &dimScal)
itk::Point< TInput, NDimensions > computeAffineLSWFromTranslations(std::vector< itk::Point< TInput, NDimensions > > &inputOrigins, std::vector< itk::Point< TInput, NDimensions > > &inputTransformed, std::vector< TInput > &weights, typename itk::AffineTransform< TScalarType, NDimensions >::Pointer &resultTransform)
void computeLogEuclideanAverage(std::vector< vnl_matrix< TInput > > &inputTransforms, std::vector< TInput > &weights, typename itk::AffineTransform< TScalarType, NDimensions >::Pointer &resultTransform)
void computeTranslationLSWFromTranslations(std::vector< itk::Point< TInput, NDimensions > > &inputOrigins, std::vector< itk::Point< TInput, NDimensions > > &inputTransformed, std::vector< TInput > &weights, typename itk::AffineTransform< TScalarType, NDimensions >::Pointer &resultTransform)
void pairingToQuaternion(const vnl_vector_fixed< TInput, NDimensions > &inputPoint, const vnl_vector_fixed< TInput, NDimensions > &inputTransformedPoint, vnl_matrix< TOutput > &outputMatrix)
vnl_matrix< T > GetExponential(const vnl_matrix< T > &m, const int numApprox=1)
Computation of the matrix exponential. Algo: classical scaling and squaring, as in Matlab...
void computeRigidLSWFromTranslations(std::vector< itk::Point< TInput, NDimensions > > &inputOrigins, std::vector< itk::Point< TInput, NDimensions > > &inputTransformed, std::vector< TInput > &weights, typename itk::AffineTransform< TScalarType, NDimensions >::Pointer &resultTransform)
vnl_matrix< TOutput > computeRotationFromQuaternion(vnl_vector< TInput > eigenVector)
itk::Point< TInput, NDimensions > computeAnisotropSimLSWFromTranslations(std::vector< itk::Point< TInput, NDimensions > > &inputOrigins, std::vector< itk::Point< TInput, NDimensions > > &inputTransformed, std::vector< TInput > &weights, typename itk::AffineTransform< TScalarType, NDimensions >::Pointer &resultTransform, vnl_matrix< double > &UMatrix)