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.AlgebraException;
19  import com.irurueta.algebra.Matrix;
20  import com.irurueta.geometry.Plane;
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.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 projective transformation refiner using plane correspondences.
35   * This class takes into account an initial estimation, inlier line matches
36   * and their residuals to find a solution that minimizes error of inliers in
37   * LMSE 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 PlaneCorrespondenceProjectiveTransformation3DRefiner extends
43          ProjectiveTransformation3DRefiner<Plane, Plane> {
44  
45      /**
46       * Plane to be reused when computing residuals.
47       */
48      private final Plane residualTestPlane = new Plane();
49  
50      /**
51       * Constructor.
52       */
53      public PlaneCorrespondenceProjectiveTransformation3DRefiner() {
54      }
55  
56      /**
57       * Constructor.
58       *
59       * @param initialEstimation           initial estimation to be set.
60       * @param keepCovariance              true if covariance of estimation must be kept after
61       *                                    refinement, false otherwise.
62       * @param inliers                     set indicating which of the provided matches are inliers.
63       * @param residuals                   residuals for matched samples.
64       * @param numInliers                  number of inliers on initial estimation.
65       * @param samples1                    1st set of paired samples.
66       * @param samples2                    2nd set of paired samples.
67       * @param refinementStandardDeviation standard deviation used for
68       *                                    Levenberg-Marquardt fitting.
69       */
70      public PlaneCorrespondenceProjectiveTransformation3DRefiner(
71              final ProjectiveTransformation3D initialEstimation, final boolean keepCovariance, final BitSet inliers,
72              final double[] residuals, final int numInliers, final List<Plane> samples1, final List<Plane> samples2,
73              final double refinementStandardDeviation) {
74          super(initialEstimation, keepCovariance, inliers, residuals, numInliers, samples1, samples2,
75                  refinementStandardDeviation);
76      }
77  
78      /**
79       * Constructor.
80       *
81       * @param initialEstimation           initial estimation to be set.
82       * @param keepCovariance              true if covariance of estimation must be kept after
83       *                                    refinement, false otherwise.
84       * @param inliersData                 inlier data, typically obtained from a robust
85       *                                    estimator.
86       * @param samples1                    1st set of paired samples.
87       * @param samples2                    2nd set of paired samples.
88       * @param refinementStandardDeviation standard deviation used for
89       *                                    Levenberg-Marquardt fitting.
90       */
91      public PlaneCorrespondenceProjectiveTransformation3DRefiner(
92              final ProjectiveTransformation3D initialEstimation, final boolean keepCovariance,
93              final InliersData inliersData, final List<Plane> samples1, final List<Plane> samples2,
94              final double refinementStandardDeviation) {
95          super(initialEstimation, keepCovariance, inliersData, samples1, samples2, refinementStandardDeviation);
96      }
97  
98      /**
99       * Refines provided initial estimation.
100      * This method always sets a value into provided result instance regardless
101      * of the fact that error has actually improved in LMSE terms or not.
102      *
103      * @param result instance where refined estimation will be stored.
104      * @return true if result improves (error decreases) in LMSE terms respect
105      * to initial estimation, false if no improvement has been achieved.
106      * @throws NotReadyException if not enough input data has been provided.
107      * @throws LockedException   if estimator is locked because refinement is
108      *                           already in progress.
109      * @throws RefinerException  if refinement fails for some reason (e.g. unable
110      *                           to converge to a result).
111      */
112     @Override
113     public boolean refine(final ProjectiveTransformation3D result) throws NotReadyException, LockedException,
114             RefinerException {
115         if (isLocked()) {
116             throw new LockedException();
117         }
118         if (!isReady()) {
119             throw new NotReadyException();
120         }
121 
122         locked = true;
123 
124         if (listener != null) {
125             listener.onRefineStart(this, initialEstimation);
126         }
127 
128         final var initialTotalResidual = totalResidual(initialEstimation);
129 
130         try {
131             final var initParams = new double[
132                     ProjectiveTransformation3D.HOM_COORDS * ProjectiveTransformation3D.HOM_COORDS];
133             // copy values
134             System.arraycopy(initialEstimation.getT().getBuffer(), 0, initParams, 0, initParams.length);
135 
136             // output values to be fitted/optimized will contain residuals
137             final var y = new double[numInliers];
138             // input values will contain 2 sets of 2D points to compute residuals
139             final var nDims = 2 * Plane.PLANE_NUMBER_PARAMS;
140             final var x = new Matrix(numInliers, nDims);
141             final var nSamples = inliers.length();
142             var pos = 0;
143             for (var i = 0; i < nSamples; i++) {
144                 if (inliers.get(i)) {
145                     // sample is inlier
146                     final var inputPlane = samples1.get(i);
147                     final var outputPlane = samples2.get(i);
148                     inputPlane.normalize();
149                     outputPlane.normalize();
150                     x.setElementAt(pos, 0, inputPlane.getA());
151                     x.setElementAt(pos, 1, inputPlane.getB());
152                     x.setElementAt(pos, 2, inputPlane.getC());
153                     x.setElementAt(pos, 3, inputPlane.getD());
154                     x.setElementAt(pos, 4, outputPlane.getA());
155                     x.setElementAt(pos, 5, outputPlane.getB());
156                     x.setElementAt(pos, 6, outputPlane.getC());
157                     x.setElementAt(pos, 7, outputPlane.getD());
158 
159                     y[pos] = residuals[i];
160                     pos++;
161                 }
162             }
163 
164             final var evaluator = new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
165 
166                 private final Plane inputPlane = new Plane();
167 
168                 private final Plane outputPlane = new Plane();
169 
170                 private final ProjectiveTransformation3D transformation = new ProjectiveTransformation3D();
171 
172                 private final GradientEstimator gradientEstimator = new GradientEstimator(params -> {
173                     // copy values
174                     System.arraycopy(params, 0, transformation.getT().getBuffer(), 0, params.length);
175                     return residual(transformation, inputPlane, outputPlane);
176                 });
177 
178                 @Override
179                 public int getNumberOfDimensions() {
180                     return nDims;
181                 }
182 
183                 @Override
184                 public double[] createInitialParametersArray() {
185                     return initParams;
186                 }
187 
188                 @Override
189                 public double evaluate(final int i, final double[] point, final double[] params,
190                                        final double[] derivatives) throws EvaluationException {
191                     inputPlane.setParameters(point[0], point[1], point[2], point[3]);
192                     outputPlane.setParameters(point[4], point[5], point[6], point[7]);
193 
194                     // copy values
195                     System.arraycopy(params, 0, transformation.getT().getBuffer(), 0, params.length);
196 
197                     final var y = residual(transformation, inputPlane, outputPlane);
198                     gradientEstimator.gradient(params, derivatives);
199 
200                     return y;
201                 }
202             };
203 
204             final var fitter = new LevenbergMarquardtMultiDimensionFitter(evaluator, x, y,
205                     getRefinementStandardDeviation());
206 
207             fitter.fit();
208 
209             // obtain estimated params
210             final var params = fitter.getA();
211 
212             // update transformation
213 
214             // copy values for A matrix
215             System.arraycopy(params, 0, result.getT().getBuffer(), 0, params.length);
216 
217             if (keepCovariance) {
218                 // keep covariance
219                 covariance = fitter.getCovar();
220             }
221 
222             final var finalTotalResidual = totalResidual(result);
223             final var errorDecreased = finalTotalResidual < initialTotalResidual;
224 
225             if (listener != null) {
226                 listener.onRefineEnd(this, initialEstimation, result, errorDecreased);
227             }
228 
229             return errorDecreased;
230 
231         } catch (final Exception e) {
232             throw new RefinerException(e);
233         } finally {
234             locked = false;
235         }
236     }
237 
238     /**
239      * Computes the residual between the projective transformation and a pair of
240      * matched planes.
241      *
242      * @param transformation a transformation.
243      * @param inputPlane     input 3D plane.
244      * @param outputPlane    output 3D plane.
245      * @return residual.
246      */
247     private double residual(final ProjectiveTransformation3D transformation, final Plane inputPlane,
248                             final Plane outputPlane) {
249         try {
250             inputPlane.normalize();
251             outputPlane.normalize();
252 
253             transformation.transform(inputPlane, residualTestPlane);
254             return 1.0 - Math.abs(outputPlane.dotProduct(residualTestPlane));
255         } catch (final AlgebraException e) {
256             return 1.0;
257         }
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 inputPlane = samples1.get(i);
274                 final var outputPlane = samples2.get(i);
275                 inputPlane.normalize();
276                 outputPlane.normalize();
277                 result += residual(transformation, inputPlane, outputPlane);
278             }
279         }
280 
281         return result;
282     }
283 }