11 template <
unsigned int NDimensions>
12 LSWTransformAgregator <NDimensions>::
15 m_EstimationBarycenter.Fill(0);
18 template <
unsigned int NDimensions>
23 return m_EstimationBarycenter;
26 template <
unsigned int NDimensions>
48 this->lswEstimateTranslationsToTranslation();
53 this->lswEstimateTranslationsToRigid();
58 this->lswEstimateTranslationsToAnisotropSim();
63 this->lswEstimateTranslationsToAffine();
76 this->lswEstimateAnyToAffine();
86 this->lswEstimateAnyToAffine();
93 throw itk::ExceptionObject(__FILE__, __LINE__,
"Specific LSW agregation not handled yet...",ITK_LOCATION);
97 template <
unsigned int NDimensions>
100 lswEstimateTranslationsToTranslation()
104 std::vector <PointType> originPoints(nbPts);
105 std::vector <PointType> transformedPoints(nbPts);
107 for (
unsigned int i = 0;i < nbPts;++i)
111 PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
112 originPoints[i] = tmpOrig;
113 transformedPoints[i] = tmpDisp;
116 typename BaseOutputTransformType::Pointer resultTransform;
117 anima::computeTranslationLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>(originPoints,transformedPoints,this->
GetInputWeights(),resultTransform);
121 template <
unsigned int NDimensions>
124 lswEstimateTranslationsToRigid()
129 throw itk::ExceptionObject(__FILE__, __LINE__,
"Dimension not supported for quaternions",ITK_LOCATION);
131 std::vector <PointType> originPoints(nbPts);
132 std::vector <PointType> transformedPoints(nbPts);
134 for (
unsigned int i = 0;i < nbPts;++i)
138 PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
139 originPoints[i] = tmpOrig;
140 transformedPoints[i] = tmpDisp;
143 typename BaseOutputTransformType::Pointer resultTransform;
144 anima::computeRigidLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>(originPoints,transformedPoints,this->
GetInputWeights(),resultTransform);
148 template <
unsigned int NDimensions>
151 lswEstimateTranslationsToAnisotropSim()
157 std::cerr <<
"Dimension not supported for quaternions" << std::endl;
161 std::vector <PointType> originPoints(nbPts);
162 std::vector <PointType> transformedPoints(nbPts);
166 for (
unsigned int i = 0; i < nbPts; ++i)
171 PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
172 originPoints[i] = tmpOrig;
173 transformedPoints[i] = currTrsf->TransformPoint(tmpDisp);
177 typename BaseOutputTransformType::Pointer resultTransform;
179 itk::Matrix<ScalarType, NDimensions, NDimensions> emptyMatrix;
182 vnl_matrix <ScalarType> covPcaOriginPoints(NDimensions, NDimensions, 0);
189 itk::Point <ScalarType, NDimensions> unweightedBarX;
190 vnl_matrix <ScalarType> covOriginPoints(NDimensions, NDimensions, 0);
191 for (
unsigned int i = 0; i < nbPts; ++i)
193 for (
unsigned int j = 0; j < NDimensions; ++j)
194 unweightedBarX[j] += originPoints[i][j] / nbPts;
196 for (
unsigned int i = 0; i < nbPts; ++i)
198 for (
unsigned int j = 0; j < NDimensions; ++j)
200 for (
unsigned int k = 0; k < NDimensions; ++k)
201 covOriginPoints(j, k) += (originPoints[i][j] - unweightedBarX[j])*(originPoints[i][k] - unweightedBarX[k]);
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);
209 covPcaOriginPoints = covPcaOriginPoints.transpose();
210 if (vnl_determinant(covPcaOriginPoints) < 0)
211 covPcaOriginPoints *= -1.0;
214 m_EstimationBarycenter = anima::computeAnisotropSimLSWFromTranslations<InternalScalarType, ScalarType, NDimensions>(originPoints, transformedPoints, this->
GetInputWeights(), resultTransform, covPcaOriginPoints);
218 template <
unsigned int NDimensions>
221 lswEstimateTranslationsToAffine()
225 std::vector <PointType> originPoints(nbPts);
226 std::vector <PointType> transformedPoints(nbPts);
228 for (
unsigned int i = 0;i < nbPts;++i)
232 PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
233 originPoints[i] = tmpOrig;
234 transformedPoints[i] = tmpDisp;
237 typename BaseOutputTransformType::Pointer resultTransform;
238 m_EstimationBarycenter=anima::computeAffineLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>(originPoints,transformedPoints,this->
GetInputWeights(),resultTransform);
242 template <
unsigned int NDimensions>
245 lswEstimateAnyToAffine()
247 typedef itk::MatrixOffsetTransformBase <InternalScalarType, NDimensions, NDimensions> BaseMatrixTransformType;
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;
258 for (
unsigned int i = 0;i < nbPts;++i)
262 BaseMatrixTransformType *tmpTrsf = (BaseMatrixTransformType *)this->
GetInputTransform(i);
263 affinePart = tmpTrsf->GetMatrix();
264 offsetPart = tmpTrsf->GetOffset();
266 for (
unsigned int j = 0;j < NDimensions;++j)
268 tmpMatrix(j,NDimensions) = offsetPart[j];
269 for (
unsigned int k = 0;k < NDimensions;++k)
270 tmpMatrix(j,k) = affinePart(j,k);
274 if (std::isnan(logTransformations[i](0,0)))
276 logTransformations[i].fill(0);
282 LogRigidTransformType *tmpTrsf = (LogRigidTransformType *)this->
GetInputTransform(i);
283 logTransformations[i] = tmpTrsf->GetLogTransform();
287 typename BaseOutputTransformType::Pointer resultTransform;
288 anima::computeLogEuclideanAverage<InternalScalarType,ScalarType,NDimensions>(logTransformations,this->
GetInputWeights(),resultTransform);
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...