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.CoordinatesType;
20  import com.irurueta.geometry.Point3D;
21  import com.irurueta.geometry.ProjectiveTransformation3D;
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.MultiDimensionFunctionEvaluatorListener;
27  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFitter;
28  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFunctionEvaluator;
29  import com.irurueta.numerical.robust.InliersData;
30  
31  import java.util.BitSet;
32  import java.util.List;
33  
34  /**
35   * A 3D projective transformation refiner using point correspondences.
36   * This class takes into account an initial estimation, inlier point matches
37   * and their residuals to find a solution that minimizes error of inliers in
38   * LMSE terms.
39   * Typically, a refiner is used by a robust estimator, however it can also be
40   * useful in some other situations.
41   */
42  @SuppressWarnings("DuplicatedCode")
43  public class PointCorrespondenceProjectiveTransformation3DRefiner extends
44          ProjectiveTransformation3DRefiner<Point3D, Point3D> {
45  
46      /**
47       * Point to be reused when computing residuals.
48       */
49      private final Point3D residualTestPoint = Point3D.create(CoordinatesType.HOMOGENEOUS_COORDINATES);
50  
51      /**
52       * Constructor.
53       */
54      public PointCorrespondenceProjectiveTransformation3DRefiner() {
55      }
56  
57      /**
58       * Constructor.
59       *
60       * @param initialEstimation           initial estimation to be set.
61       * @param keepCovariance              true if covariance of estimation must be kept after
62       *                                    refinement, false otherwise.
63       * @param inliers                     set indicating which of the provided matches are inliers.
64       * @param residuals                   residuals for matched samples.
65       * @param numInliers                  number of inliers on initial estimation.
66       * @param samples1                    1st set of paired samples.
67       * @param samples2                    2nd set of paired samples.
68       * @param refinementStandardDeviation standard deviation used for
69       *                                    Levenberg-Marquardt fitting.
70       */
71      public PointCorrespondenceProjectiveTransformation3DRefiner(
72              final ProjectiveTransformation3D initialEstimation, final boolean keepCovariance, final BitSet inliers,
73              final double[] residuals, final int numInliers, final List<Point3D> samples1, final List<Point3D> samples2,
74              final double refinementStandardDeviation) {
75          super(initialEstimation, keepCovariance, inliers, residuals, numInliers, samples1, samples2,
76                  refinementStandardDeviation);
77      }
78  
79      /**
80       * Constructor.
81       *
82       * @param initialEstimation           initial estimation to be set.
83       * @param keepCovariance              true if covariance of estimation must be kept after
84       *                                    refinement, false otherwise.
85       * @param inliersData                 inlier data, typically obtained from a robust
86       *                                    estimator.
87       * @param samples1                    1st set of paired samples.
88       * @param samples2                    2nd set of paired samples.
89       * @param refinementStandardDeviation standard deviation used for
90       *                                    Levenberg-Marquardt fitting.
91       */
92      public PointCorrespondenceProjectiveTransformation3DRefiner(
93              final ProjectiveTransformation3D initialEstimation, final boolean keepCovariance,
94              final InliersData inliersData, final List<Point3D> samples1, final List<Point3D> samples2,
95              final double refinementStandardDeviation) {
96          super(initialEstimation, keepCovariance, inliersData, samples1, samples2, refinementStandardDeviation);
97      }
98  
99      /**
100      * Refines provided initial estimation.
101      * This method always sets a value into provided result instance regardless
102      * of the fact that error has actually improved in LMSE terms or not.
103      *
104      * @param result instance where refined estimation will be stored.
105      * @return true if result improves (error decreases) in LMSE terms respect
106      * to initial estimation, false if no improvement has been achieved.
107      * @throws NotReadyException if not enough input data has been provided.
108      * @throws LockedException   if estimator is locked because refinement is
109      *                           already in progress.
110      * @throws RefinerException  if refinement fails for some reason (e.g. unable
111      *                           to converge to a result).
112      */
113     @Override
114     public boolean refine(final ProjectiveTransformation3D result) throws NotReadyException, LockedException,
115             RefinerException {
116         if (isLocked()) {
117             throw new LockedException();
118         }
119         if (!isReady()) {
120             throw new NotReadyException();
121         }
122 
123         locked = true;
124 
125         if (listener != null) {
126             listener.onRefineStart(this, initialEstimation);
127         }
128 
129         initialEstimation.normalize();
130 
131         final var initialTotalResidual = totalResidual(initialEstimation);
132 
133         try {
134             final var initParams = new double[ProjectiveTransformation3D.HOM_COORDS
135                     * ProjectiveTransformation3D.HOM_COORDS];
136             // copy values
137             System.arraycopy(initialEstimation.getT().getBuffer(), 0, initParams, 0, initParams.length);
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 ProjectiveTransformation3D transformation = new ProjectiveTransformation3D();
174 
175                 private final GradientEstimator gradientEstimator = new GradientEstimator(params -> {
176                     // copy values
177                     System.arraycopy(params, 0, transformation.getT().getBuffer(), 0, params.length);
178                     return residual(transformation, inputPoint, outputPoint);
179                 });
180 
181                 @Override
182                 public int getNumberOfDimensions() {
183                     return nDims;
184                 }
185 
186                 @Override
187                 public double[] createInitialParametersArray() {
188                     return initParams;
189                 }
190 
191                 @Override
192                 public double evaluate(final int i, final double[] point, final double[] params,
193                                        final double[] derivatives) throws EvaluationException {
194                     inputPoint.setHomogeneousCoordinates(point[0], point[1], point[2], point[3]);
195                     outputPoint.setHomogeneousCoordinates(point[4], point[5], point[6], point[7]);
196 
197                     // copy values
198                     System.arraycopy(params, 0, transformation.getT().getBuffer(), 0, params.length);
199 
200                     final var y = residual(transformation, inputPoint, outputPoint);
201                     gradientEstimator.gradient(params, derivatives);
202 
203                     return y;
204 
205                 }
206             };
207 
208             final var fitter = new LevenbergMarquardtMultiDimensionFitter(evaluator, x, y,
209                     getRefinementStandardDeviation());
210 
211             fitter.fit();
212 
213             // obtain estimated params
214             final var params = fitter.getA();
215 
216             // update transformation
217 
218             // copy values
219             System.arraycopy(params, 0, result.getT().getBuffer(), 0, params.length);
220 
221             if (keepCovariance) {
222                 // keep covariance
223                 covariance = fitter.getCovar();
224             }
225 
226             final var finalTotalResidual = totalResidual(result);
227             final var errorDecreased = finalTotalResidual < initialTotalResidual;
228 
229             if (listener != null) {
230                 listener.onRefineEnd(this, initialEstimation, result, errorDecreased);
231             }
232 
233             return errorDecreased;
234 
235         } catch (final Exception e) {
236             throw new RefinerException(e);
237         } finally {
238             locked = false;
239         }
240     }
241 
242     /**
243      * Computes the residual between the affine transformation and a pair of
244      * matched points.
245      *
246      * @param transformation a transformation.
247      * @param inputPoint     input 2D point.
248      * @param outputPoint    output 2D point.
249      * @return residual.
250      */
251     private double residual(final ProjectiveTransformation3D transformation, final Point3D inputPoint,
252                             final Point3D outputPoint) {
253         inputPoint.normalize();
254         outputPoint.normalize();
255 
256         transformation.transform(inputPoint, residualTestPoint);
257         return residualTestPoint.distanceTo(outputPoint);
258     }
259 
260     /**
261      * Computes total residual among all provided inlier samples.
262      *
263      * @param transformation a transformation.
264      * @return total residual.
265      */
266     private double totalResidual(final ProjectiveTransformation3D transformation) {
267         var result = 0.0;
268 
269         final var nSamples = inliers.length();
270         for (var i = 0; i < nSamples; i++) {
271             if (inliers.get(i)) {
272                 // sample is inlier
273                 final var inputPoint = samples1.get(i);
274                 final var outputPoint = samples2.get(i);
275                 inputPoint.normalize();
276                 outputPoint.normalize();
277                 result += residual(transformation, inputPoint, outputPoint);
278             }
279         }
280 
281         return result;
282     }
283 }