13 template <
unsigned int NDimensions>
14 LTSWTransformAgregator <NDimensions>::
18 m_StoppingThreshold = 1.0e-2;
19 m_EstimationBarycenter.Fill(0);
22 template <
unsigned int NDimensions>
27 return m_EstimationBarycenter;
30 template <
unsigned int NDimensions>
36 bool returnValue =
false;
49 returnValue = this->ltswEstimateTranslationsToAny();
56 returnValue = this->ltswEstimateAnyToAffine();
60 throw itk::ExceptionObject(__FILE__, __LINE__,
"Specific LTSW agregation not handled yet...",ITK_LOCATION);
64 template <
unsigned int NDimensions>
67 ltswEstimateTranslationsToAny()
72 throw itk::ExceptionObject(__FILE__, __LINE__,
"Dimension not supported",ITK_LOCATION);
74 std::vector <PointType> originPoints(nbPts);
75 std::vector <PointType> transformedPoints(nbPts);
82 for (
unsigned int i = 0; i < nbPts; ++i)
86 PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
87 originPoints[i] = tmpOrig;
89 transformedPoints[i] = currTrsf->TransformPoint(tmpDisp);
91 transformedPoints[i] = tmpDisp;
94 vnl_matrix <ScalarType> covPcaOriginPoints(NDimensions, NDimensions, 0);
97 itk::Matrix<ScalarType, NDimensions, NDimensions> emptyMatrix;
106 itk::Point <ScalarType, NDimensions> unweightedBarX;
107 vnl_matrix <ScalarType> covOriginPoints(NDimensions, NDimensions, 0);
108 for (
unsigned int i = 0; i < nbPts; ++i)
110 for (
unsigned int j = 0; j < NDimensions; ++j)
111 unweightedBarX[j] += originPoints[i][j] / nbPts;
113 for (
unsigned int i = 0; i < nbPts; ++i)
115 for (
unsigned int j = 0; j < NDimensions; ++j)
117 for (
unsigned int k = 0; k < NDimensions; ++k)
118 covOriginPoints(j, k) += (originPoints[i][j] - unweightedBarX[j])*(originPoints[i][k] - unweightedBarX[k]);
121 itk::SymmetricEigenAnalysis < vnl_matrix <ScalarType>, vnl_diag_matrix<ScalarType>, vnl_matrix <ScalarType> > eigenSystem(3);
122 vnl_diag_matrix <double> eValsCov(NDimensions);
123 eigenSystem.SetOrderEigenValues(
true);
124 eigenSystem.ComputeEigenValuesAndVectors(covOriginPoints, eValsCov, covPcaOriginPoints);
126 covPcaOriginPoints = covPcaOriginPoints.transpose();
127 if (vnl_determinant(covPcaOriginPoints) < 0)
128 covPcaOriginPoints *= -1.0;
132 std::vector <PointType> originPointsFiltered = originPoints;
133 std::vector <PointType> transformedPointsFiltered = transformedPoints;
134 std::vector <double> weightsFiltered = weights;
136 std::vector < std::pair <unsigned int, double> > residualErrors;
138 itk::Vector <double,3> tmpDiff;
140 bool continueLoop =
true;
141 unsigned int numMaxIter = 100;
142 unsigned int num_itr = 0;
144 typename BaseOutputTransformType::Pointer resultTransform, resultTransformOld;
146 while(num_itr < numMaxIter)
153 anima::computeTranslationLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>
154 (originPointsFiltered,transformedPointsFiltered,weightsFiltered,resultTransform);
158 anima::computeRigidLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>
159 (originPointsFiltered,transformedPointsFiltered,weightsFiltered,resultTransform);
163 m_EstimationBarycenter = anima::computeAnisotropSimLSWFromTranslations<InternalScalarType, ScalarType, NDimensions>
164 (originPointsFiltered, transformedPointsFiltered, weightsFiltered, resultTransform, covPcaOriginPoints);
168 m_EstimationBarycenter = anima::computeAffineLSWFromTranslations<InternalScalarType,ScalarType,NDimensions>
169 (originPointsFiltered,transformedPointsFiltered,weightsFiltered,resultTransform);
173 throw itk::ExceptionObject(__FILE__, __LINE__,
"Not implemented yet...",ITK_LOCATION);
176 continueLoop = endLTSCondition(resultTransformOld,resultTransform);
181 resultTransformOld = resultTransform;
182 residualErrors.clear();
183 for (
unsigned int i = 0;i < nbPts;++i)
188 tmpOutPoint = resultTransform->TransformPoint(originPoints[i]);
189 tmpDiff = tmpOutPoint - transformedPoints[i];
190 residualErrors.push_back(std::make_pair (i,tmpDiff.GetNorm()));
193 unsigned int numLts = floor(residualErrors.size() * m_LTSCut);
195 std::vector < std::pair <unsigned int, double> >::iterator begIt = residualErrors.begin();
196 std::vector < std::pair <unsigned int, double> >::iterator sortPart = begIt + numLts;
200 originPointsFiltered.resize(numLts);
201 transformedPointsFiltered.resize(numLts);
202 weightsFiltered.resize(numLts);
204 for (
unsigned int i = 0;i < numLts;++i)
206 originPointsFiltered[i] = originPoints[residualErrors[i].first];
207 transformedPointsFiltered[i] = transformedPoints[residualErrors[i].first];
208 weightsFiltered[i] = weights[residualErrors[i].first];
216 template <
unsigned int NDimensions>
219 ltswEstimateAnyToAffine()
222 throw itk::ExceptionObject(__FILE__, __LINE__,
"Agregation from affine transforms to rigid is not supported yet...",ITK_LOCATION);
224 typedef itk::MatrixOffsetTransformBase <InternalScalarType, NDimensions> BaseMatrixTransformType;
230 std::vector < vnl_matrix <InternalScalarType> > logTransformations(nbPts);
231 vnl_matrix <InternalScalarType> tmpMatrix(NDimensions+1,NDimensions+1,0), tmpLogMatrix(NDimensions+1,NDimensions+1,0);
232 tmpMatrix(NDimensions,NDimensions) = 1;
233 typename BaseMatrixTransformType::MatrixType affinePart;
234 itk::Vector <InternalScalarType, NDimensions> offsetPart;
236 for (
unsigned int i = 0;i < nbPts;++i)
240 BaseMatrixTransformType *tmpTrsf = (BaseMatrixTransformType *)this->
GetInputTransform(i);
241 affinePart = tmpTrsf->GetMatrix();
242 offsetPart = tmpTrsf->GetOffset();
244 for (
unsigned int j = 0;j < NDimensions;++j)
246 tmpMatrix(j,NDimensions) = offsetPart[j];
247 for (
unsigned int k = 0;k < NDimensions;++k)
248 tmpMatrix(j,k) = affinePart(j,k);
252 if (!std::isfinite(logTransformations[i](0,0)))
254 logTransformations[i].fill(0);
260 LogRigidTransformType *tmpTrsf = (LogRigidTransformType *)this->
GetInputTransform(i);
261 logTransformations[i] = tmpTrsf->GetLogTransform();
265 std::vector < vnl_matrix <InternalScalarType> > logTransformationsFiltered = logTransformations;
266 std::vector <InternalScalarType> weightsFiltered = weights;
269 std::vector < PointType > originPoints(nbPts);
270 std::vector < PointType > transformedPoints(nbPts);
272 for (
unsigned int i = 0;i < nbPts;++i)
276 PointType tmpDisp = tmpTrsf->TransformPoint(tmpOrig);
277 originPoints[i] = tmpOrig;
278 transformedPoints[i] = tmpDisp;
281 std::vector < std::pair <unsigned int, double> > residualErrors;
283 bool continueLoop =
true;
284 unsigned int numMaxIter = 100;
285 unsigned int num_itr = 0;
287 typename BaseOutputTransformType::Pointer resultTransform, resultTransformOld;
289 while(num_itr < numMaxIter)
293 anima::computeLogEuclideanAverage<InternalScalarType,ScalarType,NDimensions>(logTransformationsFiltered,weightsFiltered,resultTransform);
294 continueLoop = endLTSCondition(resultTransformOld,resultTransform);
299 resultTransformOld = resultTransform;
300 residualErrors.clear();
302 BaseMatrixTransformType *tmpTrsf = (BaseMatrixTransformType *)resultTransform.GetPointer();
304 for (
unsigned int i = 0;i < nbPts;++i)
310 PointType tmpDisp = tmpTrsf->TransformPoint(originPoints[i]);
312 for (
unsigned int j = 0;j < NDimensions;++j)
313 tmpDiff += (transformedPoints[i][j] - tmpDisp[j]) * (transformedPoints[i][j] - tmpDisp[j]);
315 residualErrors.push_back(std::make_pair (i,tmpDiff));
318 unsigned int numLts = floor(residualErrors.size() * m_LTSCut);
320 std::vector < std::pair <unsigned int, double> >::iterator begIt = residualErrors.begin();
321 std::vector < std::pair <unsigned int, double> >::iterator sortPart = begIt + numLts;
325 logTransformationsFiltered.resize(numLts);
326 weightsFiltered.resize(numLts);
328 for (
unsigned int i = 0;i < numLts;++i)
330 logTransformationsFiltered[i] = logTransformations[residualErrors[i].first];
331 weightsFiltered[i] = weights[residualErrors[i].first];
339 template <
unsigned int NDimensions>
347 typename BaseOutputTransformType::ParametersType oldParams = oldTrsf->GetParameters();
348 typename BaseOutputTransformType::ParametersType newParams = newTrsf->GetParameters();
350 for (
unsigned int i = 0;i < newParams.GetSize();++i)
352 double diffParam = fabs(newParams[i] - oldParams[i]);
353 if (diffParam > m_StoppingThreshold)
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...