View Javadoc
1   /*
2    * Copyright (C) 2017 Alberto Irurueta Carro (alberto@irurueta.com)
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *         http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
15   */
16  package com.irurueta.geometry.refiners;
17  
18  import com.irurueta.algebra.Matrix;
19  import com.irurueta.geometry.AffineTransformation3D;
20  import com.irurueta.geometry.CoordinatesType;
21  import com.irurueta.geometry.Point3D;
22  import com.irurueta.geometry.estimators.LockedException;
23  import com.irurueta.geometry.estimators.NotReadyException;
24  import com.irurueta.numerical.EvaluationException;
25  import com.irurueta.numerical.GradientEstimator;
26  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFitter;
27  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFunctionEvaluator;
28  import com.irurueta.numerical.robust.InliersData;
29  
30  import java.util.BitSet;
31  import java.util.List;
32  
33  /**
34   * A 3D affine transformation refiner using point correspondences.
35   * This class takes into account an initial estimation, inlier point matches and
36   * their residuals to find a solution that minimizes error of inliers in LMSE
37   * terms.
38   * Typically, a refiner is used by a robust estimator, however it can also be
39   * useful in some other situations.
40   */
41  @SuppressWarnings("DuplicatedCode")
42  public class PointCorrespondenceAffineTransformation3DRefiner extends AffineTransformation3DRefiner<Point3D, Point3D> {
43  
44      /**
45       * Point to be reused when computing residuals.
46       */
47      private final Point3D residualTestPoint = Point3D.create(CoordinatesType.HOMOGENEOUS_COORDINATES);
48  
49      /**
50       * Constructor.
51       */
52      public PointCorrespondenceAffineTransformation3DRefiner() {
53      }
54  
55      /**
56       * Constructor.
57       *
58       * @param initialEstimation           initial estimation to be set.
59       * @param keepCovariance              true if covariance of estimation must be kept after
60       *                                    refinement, false otherwise.
61       * @param inliers                     set indicating which of the provided matches are inliers.
62       * @param residuals                   residuals for matched samples.
63       * @param numInliers                  number of inliers on initial estimation.
64       * @param samples1                    1st set of paired samples.
65       * @param samples2                    2nd set of paired samples.
66       * @param refinementStandardDeviation standard deviation used for
67       *                                    Levenberg-Marquardt fitting.
68       */
69      public PointCorrespondenceAffineTransformation3DRefiner(
70              final AffineTransformation3D initialEstimation, final boolean keepCovariance, final BitSet inliers,
71              final double[] residuals, final int numInliers, final List<Point3D> samples1, final List<Point3D> samples2,
72              final double refinementStandardDeviation) {
73          super(initialEstimation, keepCovariance, inliers, residuals, numInliers, samples1, samples2,
74                  refinementStandardDeviation);
75      }
76  
77      /**
78       * Constructor.
79       *
80       * @param initialEstimation           initial estimation to be set.
81       * @param keepCovariance              true if covariance of estimation must be kept after
82       *                                    refinement, false otherwise.
83       * @param inliersData                 inlier data, typically obtained from a robust
84       *                                    estimator.
85       * @param samples1                    1st set of paired samples.
86       * @param samples2                    2nd set of paired samples.
87       * @param refinementStandardDeviation standard deviation used for
88       *                                    Levenberg-Marquardt fitting.
89       */
90      public PointCorrespondenceAffineTransformation3DRefiner(
91              final AffineTransformation3D initialEstimation, final boolean keepCovariance, final InliersData inliersData,
92              final List<Point3D> samples1, final List<Point3D> samples2, final double refinementStandardDeviation) {
93          super(initialEstimation, keepCovariance, inliersData, samples1, samples2, refinementStandardDeviation);
94      }
95  
96      /**
97       * Refines provide initial estimation.
98       * This method always sets a value into provided result instance regardless
99       * of the fact that error has actually improved in LMSE terms or not.
100      *
101      * @param result instance where refined estimation will be stored.
102      * @return true if result improves (error decreases) in LMSE terms respect
103      * to initial estimation, false if no improvement has been achieved.
104      * @throws NotReadyException if not enough input data has been provided.
105      * @throws LockedException   if estimator is locked because refinement is
106      *                           already in progress.
107      * @throws RefinerException  if refinement fails for some reason (e.g. unable
108      *                           to converge to a result).
109      */
110     @Override
111     public boolean refine(final AffineTransformation3D result) throws NotReadyException, LockedException,
112             RefinerException {
113         if (isLocked()) {
114             throw new LockedException();
115         }
116         if (!isReady()) {
117             throw new NotReadyException();
118         }
119 
120         locked = true;
121 
122         if (listener != null) {
123             listener.onRefineStart(this, initialEstimation);
124         }
125 
126         final var initialTotalResidual = totalResidual(initialEstimation);
127 
128         try {
129             final var initParams = new double[AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS
130                     + AffineTransformation3D.NUM_TRANSLATION_COORDS];
131             // copy values for A matrix
132             System.arraycopy(initialEstimation.getA().getBuffer(), 0, initParams, 0,
133                     AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS);
134             // copy values for translation
135             System.arraycopy(initialEstimation.getTranslation(), 0, initParams,
136                     AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS,
137                     AffineTransformation3D.NUM_TRANSLATION_COORDS);
138 
139             // output values to be fitted/optimized will contain residuals
140             final var y = new double[numInliers];
141             // input values will contain 2 sets of 2D points to compute residuals
142             final var nDims = 2 * Point3D.POINT3D_HOMOGENEOUS_COORDINATES_LENGTH;
143             final var x = new Matrix(numInliers, nDims);
144             final var nSamples = inliers.length();
145             var pos = 0;
146             for (var i = 0; i < nSamples; i++) {
147                 if (inliers.get(i)) {
148                     // sample is inlier
149                     final var inputPoint = samples1.get(i);
150                     final var outputPoint = samples2.get(i);
151                     inputPoint.normalize();
152                     outputPoint.normalize();
153                     x.setElementAt(pos, 0, inputPoint.getHomX());
154                     x.setElementAt(pos, 1, inputPoint.getHomY());
155                     x.setElementAt(pos, 2, inputPoint.getHomZ());
156                     x.setElementAt(pos, 3, inputPoint.getHomW());
157                     x.setElementAt(pos, 4, outputPoint.getHomX());
158                     x.setElementAt(pos, 5, outputPoint.getHomY());
159                     x.setElementAt(pos, 6, outputPoint.getHomZ());
160                     x.setElementAt(pos, 7, outputPoint.getHomW());
161 
162                     y[pos] = residuals[i];
163                     pos++;
164                 }
165             }
166 
167             final var evaluator = new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
168 
169                 private final Point3D inputPoint = Point3D.create(CoordinatesType.HOMOGENEOUS_COORDINATES);
170 
171                 private final Point3D outputPoint = Point3D.create(CoordinatesType.HOMOGENEOUS_COORDINATES);
172 
173                 private final AffineTransformation3D transformation = new AffineTransformation3D();
174 
175                 private final GradientEstimator gradientEstimator = new GradientEstimator(params -> {
176                     // copy values for A matrix
177                     System.arraycopy(params, 0, transformation.getA().getBuffer(), 0,
178                             AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS);
179                     // copy values for translation
180                     System.arraycopy(params,
181                             AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS,
182                             transformation.getTranslation(), 0, AffineTransformation3D.NUM_TRANSLATION_COORDS);
183 
184                     return residual(transformation, inputPoint, outputPoint);
185 
186                 });
187 
188                 @Override
189                 public int getNumberOfDimensions() {
190                     return nDims;
191                 }
192 
193                 @Override
194                 public double[] createInitialParametersArray() {
195                     return initParams;
196                 }
197 
198                 @Override
199                 public double evaluate(final int i, final double[] point, final double[] params,
200                                        final double[] derivatives) throws EvaluationException {
201                     inputPoint.setHomogeneousCoordinates(point[0], point[1], point[2], point[3]);
202                     outputPoint.setHomogeneousCoordinates(point[4], point[5], point[6], point[7]);
203 
204                     // copy values for A matrix
205                     System.arraycopy(params, 0, transformation.getA().getBuffer(), 0,
206                             AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS);
207                     // copy values for translation
208                     System.arraycopy(params,
209                             AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS,
210                             transformation.getTranslation(), 0, AffineTransformation3D.NUM_TRANSLATION_COORDS);
211 
212                     final var y = residual(transformation, inputPoint, outputPoint);
213                     gradientEstimator.gradient(params, derivatives);
214 
215                     return y;
216                 }
217             };
218 
219             final var fitter = new LevenbergMarquardtMultiDimensionFitter(evaluator, x, y,
220                     getRefinementStandardDeviation());
221 
222             fitter.fit();
223 
224             // obtain estimated params
225             final var params = fitter.getA();
226 
227             // update transformation
228 
229             // copy values for A matrix
230             System.arraycopy(params, 0, result.getA().getBuffer(), 0,
231                     AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS);
232             // copy values for translation
233             System.arraycopy(params, AffineTransformation3D.INHOM_COORDS * AffineTransformation3D.INHOM_COORDS,
234                     result.getTranslation(), 0, AffineTransformation3D.NUM_TRANSLATION_COORDS);
235 
236             if (keepCovariance) {
237                 // keep covariance
238                 covariance = fitter.getCovar();
239             }
240 
241             final var finalTotalResidual = totalResidual(result);
242             final var errorDecreased = finalTotalResidual < initialTotalResidual;
243 
244             if (listener != null) {
245                 listener.onRefineEnd(this, initialEstimation, result, errorDecreased);
246             }
247 
248             return errorDecreased;
249 
250         } catch (final Exception e) {
251             throw new RefinerException(e);
252         } finally {
253             locked = false;
254         }
255     }
256 
257     /**
258      * Computes the residual between the affine transformation and a pair of
259      * matched points.
260      *
261      * @param transformation a transformation.
262      * @param inputPoint     input 3D point.
263      * @param outputPoint    output 3D point.
264      * @return residual.
265      */
266     private double residual(final AffineTransformation3D transformation, final Point3D inputPoint,
267                             final Point3D outputPoint) {
268         inputPoint.normalize();
269         outputPoint.normalize();
270 
271         transformation.transform(inputPoint, residualTestPoint);
272         return residualTestPoint.distanceTo(outputPoint);
273     }
274 
275     /**
276      * Computes total residual among all provided inlier samples.
277      *
278      * @param transformation a transformation.
279      * @return total residual.
280      */
281     private double totalResidual(final AffineTransformation3D transformation) {
282         var result = 0.0;
283 
284         final var nSamples = inliers.length();
285         for (var i = 0; i < nSamples; i++) {
286             if (inliers.get(i)) {
287                 // sample is inlier
288                 final var inputPoint = samples1.get(i);
289                 final var outputPoint = samples2.get(i);
290                 inputPoint.normalize();
291                 outputPoint.normalize();
292                 result += residual(transformation, inputPoint, outputPoint);
293             }
294         }
295 
296         return result;
297     }
298 }