ANIMA  4.0
animaLSWTransformAgregator.hxx
Go to the documentation of this file.
1 #pragma once
2 
6 #include <animaMatrixLogExp.h>
7 
8 namespace anima
9 {
10 
11 template <unsigned int NDimensions>
12 LSWTransformAgregator <NDimensions>::
13 LSWTransformAgregator() : Superclass()
14 {
15  m_EstimationBarycenter.Fill(0);
16 }
17 
18 template <unsigned int NDimensions>
22 {
23  return m_EstimationBarycenter;
24 }
25 
26 template <unsigned int NDimensions>
27 bool
30 {
31  this->SetUpToDate(false);
32 
33  if (this->GetInputWeights().size() != this->GetInputTransforms().size())
34  return false;
35 
36  switch (this->GetInputTransformType())
37  {
40  {
41  if ((this->GetInputWeights().size() != this->GetInputOrigins().size())||
42  (this->GetInputTransforms().size() != this->GetInputOrigins().size()))
43  return false;
44 
45  switch (this->GetOutputTransformType())
46  {
48  this->lswEstimateTranslationsToTranslation();
49  this->SetUpToDate(true);
50  return true;
51 
52  case Superclass::RIGID:
53  this->lswEstimateTranslationsToRigid();
54  this->SetUpToDate(true);
55  return true;
56 
58  this->lswEstimateTranslationsToAnisotropSim();
59  this->SetUpToDate(true);
60  return true;
61 
62  case Superclass::AFFINE:
63  this->lswEstimateTranslationsToAffine();
64  this->SetUpToDate(true);
65  return true;
66 
67  default:
68  return false;
69  }
70  break;
71  }
72 
73  case Superclass::RIGID:
75  {
76  this->lswEstimateAnyToAffine();
77  this->SetUpToDate(true);
78  return true;
79  }
80  break;
81 
82  case Superclass::AFFINE:
83  default:
85  {
86  this->lswEstimateAnyToAffine();
87  this->SetUpToDate(true);
88  return true;
89  }
90  break;
91  }
92 
93  throw itk::ExceptionObject(__FILE__, __LINE__,"Specific LSW agregation not handled yet...",ITK_LOCATION);
94  return false;
95 }
96 
97 template <unsigned int NDimensions>
98 void
100 lswEstimateTranslationsToTranslation()
101 {
102  unsigned int nbPts = this->GetInputOrigins().size();
103 
104  std::vector <PointType> originPoints(nbPts);
105  std::vector <PointType> transformedPoints(nbPts);
106 
107  for (unsigned int i = 0;i < nbPts;++i)
108  {
109  PointType tmpOrig = this->GetInputOrigin(i);
110  BaseInputTransformType * tmpTrsf = this->GetInputTransform(i);
111  PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
112  originPoints[i] = tmpOrig;
113  transformedPoints[i] = tmpDisp;
114  }
115 
116  typename BaseOutputTransformType::Pointer resultTransform;
117  anima::computeTranslationLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>(originPoints,transformedPoints,this->GetInputWeights(),resultTransform);
118  this->SetOutput(resultTransform);
119 }
120 
121 template <unsigned int NDimensions>
122 void
124 lswEstimateTranslationsToRigid()
125 {
126  unsigned int nbPts = this->GetInputOrigins().size();
127 
128  if (NDimensions > 3)
129  throw itk::ExceptionObject(__FILE__, __LINE__,"Dimension not supported for quaternions",ITK_LOCATION);
130 
131  std::vector <PointType> originPoints(nbPts);
132  std::vector <PointType> transformedPoints(nbPts);
133 
134  for (unsigned int i = 0;i < nbPts;++i)
135  {
136  PointType tmpOrig = this->GetInputOrigin(i);
137  BaseInputTransformType * tmpTrsf = this->GetInputTransform(i);
138  PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
139  originPoints[i] = tmpOrig;
140  transformedPoints[i] = tmpDisp;
141  }
142 
143  typename BaseOutputTransformType::Pointer resultTransform;
144  anima::computeRigidLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>(originPoints,transformedPoints,this->GetInputWeights(),resultTransform);
145  this->SetOutput(resultTransform);
146 }
147 
148 template <unsigned int NDimensions>
149 void
151 lswEstimateTranslationsToAnisotropSim()
152 {
153  unsigned int nbPts = this->GetInputOrigins().size();
154 
155  if (NDimensions > 3)
156  {
157  std::cerr << "Dimension not supported for quaternions" << std::endl;
158  return;
159  }
160 
161  std::vector <PointType> originPoints(nbPts);
162  std::vector <PointType> transformedPoints(nbPts);
163 
165 
166  for (unsigned int i = 0; i < nbPts; ++i)
167  {
168  PointType tmpOrig = this->GetInputOrigin(i);
169  BaseInputTransformType * tmpTrsf = this->GetInputTransform(i);
170 
171  PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
172  originPoints[i] = tmpOrig;
173  transformedPoints[i] = currTrsf->TransformPoint(tmpDisp);
174 
175  }
176 
177  typename BaseOutputTransformType::Pointer resultTransform;
178 
179  itk::Matrix<ScalarType, NDimensions, NDimensions> emptyMatrix;
180  emptyMatrix.Fill(0);
181 
182  vnl_matrix <ScalarType> covPcaOriginPoints(NDimensions, NDimensions, 0);
183  if (this->GetOrthogonalDirectionMatrix() != emptyMatrix)
184  {
185  covPcaOriginPoints = this->GetOrthogonalDirectionMatrix().GetVnlMatrix().as_matrix();
186  }
187  else
188  {
189  itk::Point <ScalarType, NDimensions> unweightedBarX;
190  vnl_matrix <ScalarType> covOriginPoints(NDimensions, NDimensions, 0);
191  for (unsigned int i = 0; i < nbPts; ++i)
192  {
193  for (unsigned int j = 0; j < NDimensions; ++j)
194  unweightedBarX[j] += originPoints[i][j] / nbPts;
195  }
196  for (unsigned int i = 0; i < nbPts; ++i)
197  {
198  for (unsigned int j = 0; j < NDimensions; ++j)
199  {
200  for (unsigned int k = 0; k < NDimensions; ++k)
201  covOriginPoints(j, k) += (originPoints[i][j] - unweightedBarX[j])*(originPoints[i][k] - unweightedBarX[k]);
202  }
203  }
204  itk::SymmetricEigenAnalysis < vnl_matrix <ScalarType>, vnl_diag_matrix<ScalarType>, vnl_matrix <ScalarType> > eigenSystem(3);
205  vnl_diag_matrix <double> eValsCov(NDimensions);
206  eigenSystem.SetOrderEigenValues(true);
207  eigenSystem.ComputeEigenValuesAndVectors(covOriginPoints, eValsCov, covPcaOriginPoints);
208  /* return eigen vectors in row !!!!!!! */
209  covPcaOriginPoints = covPcaOriginPoints.transpose();
210  if (vnl_determinant(covPcaOriginPoints) < 0)
211  covPcaOriginPoints *= -1.0;
212  }
213 
214  m_EstimationBarycenter = anima::computeAnisotropSimLSWFromTranslations<InternalScalarType, ScalarType, NDimensions>(originPoints, transformedPoints, this->GetInputWeights(), resultTransform, covPcaOriginPoints);
215  this->SetOutput(resultTransform);
216 }
217 
218 template <unsigned int NDimensions>
219 void
221 lswEstimateTranslationsToAffine()
222 {
223  unsigned int nbPts = this->GetInputOrigins().size();
224 
225  std::vector <PointType> originPoints(nbPts);
226  std::vector <PointType> transformedPoints(nbPts);
227 
228  for (unsigned int i = 0;i < nbPts;++i)
229  {
230  PointType tmpOrig = this->GetInputOrigin(i);
231  BaseInputTransformType * tmpTrsf = this->GetInputTransform(i);
232  PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
233  originPoints[i] = tmpOrig;
234  transformedPoints[i] = tmpDisp;
235  }
236 
237  typename BaseOutputTransformType::Pointer resultTransform;
238  m_EstimationBarycenter=anima::computeAffineLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>(originPoints,transformedPoints,this->GetInputWeights(),resultTransform);
239  this->SetOutput(resultTransform);
240 }
241 
242 template <unsigned int NDimensions>
243 void
245 lswEstimateAnyToAffine()
246 {
247  typedef itk::MatrixOffsetTransformBase <InternalScalarType, NDimensions, NDimensions> BaseMatrixTransformType;
248  typedef anima::LogRigid3DTransform <InternalScalarType> LogRigidTransformType;
249 
250  unsigned int nbPts = this->GetInputTransforms().size();
251 
252  std::vector < vnl_matrix <InternalScalarType> > logTransformations(nbPts);
253  vnl_matrix <InternalScalarType> tmpMatrix(NDimensions+1,NDimensions+1,0);
254  tmpMatrix(NDimensions,NDimensions) = 1;
255  typename BaseMatrixTransformType::MatrixType affinePart;
256  itk::Vector <InternalScalarType,NDimensions> offsetPart;
257 
258  for (unsigned int i = 0;i < nbPts;++i)
259  {
261  {
262  BaseMatrixTransformType *tmpTrsf = (BaseMatrixTransformType *)this->GetInputTransform(i);
263  affinePart = tmpTrsf->GetMatrix();
264  offsetPart = tmpTrsf->GetOffset();
265 
266  for (unsigned int j = 0;j < NDimensions;++j)
267  {
268  tmpMatrix(j,NDimensions) = offsetPart[j];
269  for (unsigned int k = 0;k < NDimensions;++k)
270  tmpMatrix(j,k) = affinePart(j,k);
271  }
272 
273  logTransformations[i] = anima::GetLogarithm(tmpMatrix);
274  if (std::isnan(logTransformations[i](0,0)))
275  {
276  logTransformations[i].fill(0);
277  this->SetInputWeight(i,0);
278  }
279  }
280  else
281  {
282  LogRigidTransformType *tmpTrsf = (LogRigidTransformType *)this->GetInputTransform(i);
283  logTransformations[i] = tmpTrsf->GetLogTransform();
284  }
285  }
286 
287  typename BaseOutputTransformType::Pointer resultTransform;
288  anima::computeLogEuclideanAverage<InternalScalarType,ScalarType,NDimensions>(logTransformations,this->GetInputWeights(),resultTransform);
289  this->SetOutput(resultTransform);
290 }
291 
292 } // end of namespace anima
itk::Point< InternalScalarType, NDimensions > PointType
BaseInputTransformType * GetInputTransform(unsigned int i)
std::vector< BaseInputTransformPointer > & GetInputTransforms()
BaseInputTransformPointer & GetCurrentLinearTransform()
PointType GetEstimationBarycenter() ITK_OVERRIDE
std::vector< InternalScalarType > & GetInputWeights()
PointType & GetInputOrigin(unsigned int i)
void SetInputWeight(unsigned int i, double w)
vnl_matrix< T > GetLogarithm(const vnl_matrix< T > &m, const double precision=0.00000000001, const int numApprox=1)
Computation of the matrix logarithm. Algo: inverse scaling and squaring, variant proposed by Cheng et...
void SetOutput(BaseOutputTransformType *output)
std::vector< PointType > & GetInputOrigins()
itk::Transform< InternalScalarType, NDimensions, NDimensions > BaseInputTransformType
virtual bool Update() ITK_OVERRIDE