View Javadoc
1   /*
2    * Copyright (C) 2020 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.navigation.inertial.calibration.magnetometer;
17  
18  import com.irurueta.navigation.LockedException;
19  import com.irurueta.navigation.NotReadyException;
20  import com.irurueta.navigation.inertial.calibration.CalibrationException;
21  import com.irurueta.navigation.inertial.calibration.StandardDeviationFrameBodyMagneticFluxDensity;
22  import com.irurueta.numerical.robust.LMedSRobustEstimator;
23  import com.irurueta.numerical.robust.LMedSRobustEstimatorListener;
24  import com.irurueta.numerical.robust.RobustEstimator;
25  import com.irurueta.numerical.robust.RobustEstimatorException;
26  import com.irurueta.numerical.robust.RobustEstimatorMethod;
27  
28  import java.io.IOException;
29  import java.util.List;
30  
31  /**
32   * Robustly estimates magnetometer hard-iron biases, soft-iron cross
33   * couplings and scaling factors using LMedS algorithm.
34   * <p>
35   * To use this calibrator at least 4 measurements at different known
36   * frames must be provided. In other words, magnetometer samples must
37   * be obtained at 4 different positions or orientations.
38   * Notice that frame velocities are ignored by this calibrator.
39   * <p>
40   * Measured magnetic flux density is assumed to follow the model shown below:
41   * <pre>
42   *     mBmeas = bm + (I + Mm) * mBtrue + w
43   * </pre>
44   * Where:
45   * - mBmeas is the measured magnetic flux density. This is a 3x1 vector.
46   * - bm is magnetometer hard-iron bias. Ideally, on a perfect magnetometer,
47   * this should be a 3x1 zero vector.
48   * - I is the 3x3 identity matrix.
49   * - Mm is the 3x3 soft-iron matrix containing cross-couplings and scaling
50   * factors. Ideally, on a perfect magnetometer, this should be a 3x3 zero
51   * matrix.
52   * - mBtrue is ground-truth magnetic flux density. This is a 3x1 vector.
53   * - w is measurement noise. This is a 3x1 vector.
54   */
55  public class LMedSRobustKnownFrameMagnetometerCalibrator extends RobustKnownFrameMagnetometerCalibrator {
56  
57      /**
58       * Default value to be used for stop threshold. Stop threshold can be used to
59       * avoid keeping the algorithm unnecessarily iterating in case that best
60       * estimated threshold using median of residuals is not small enough. Once a
61       * solution is found that generates a threshold below this value, the
62       * algorithm will stop.
63       * The stop threshold can be used to prevent the LMedS algorithm iterating
64       * too many times in cases where samples have a very similar accuracy.
65       * For instance, in cases where proportion of outliers is very small (close
66       * to 0%), and samples are very accurate (i.e. 1e-6), the algorithm would
67       * iterate for a long time trying to find the best solution when indeed
68       * there is no need to do that if a reasonable threshold has already been
69       * reached.
70       * Because of this behaviour the stop threshold can be set to a value much
71       * lower than the one typically used in RANSAC, and yet the algorithm could
72       * still produce even smaller thresholds in estimated results.
73       */
74      public static final double DEFAULT_STOP_THRESHOLD = 500e-9;
75  
76      /**
77       * Minimum allowed stop threshold value.
78       */
79      public static final double MIN_STOP_THRESHOLD = 0.0;
80  
81      /**
82       * Threshold to be used to keep the algorithm iterating in case that best
83       * estimated threshold using median of residuals is not small enough. Once
84       * a solution is found that generates a threshold below this value, the
85       * algorithm will stop.
86       * The stop threshold can be used to prevent the LMedS algorithm iterating
87       * too many times in cases where samples have a very similar accuracy.
88       * For instance, in cases where proportion of outliers is very small (close
89       * to 0%), and samples are very accurate (i.e. 1e-6), the algorithm would
90       * iterate for a long time trying to find the best solution when indeed
91       * there is no need to do that if a reasonable threshold has already been
92       * reached.
93       * Because of this behaviour the stop threshold can be set to a value much
94       * lower than the one typically used in RANSAC, and yet the algorithm could
95       * still produce even smaller thresholds in estimated results.
96       */
97      private double stopThreshold = DEFAULT_STOP_THRESHOLD;
98  
99      /**
100      * Constructor.
101      */
102     public LMedSRobustKnownFrameMagnetometerCalibrator() {
103         super();
104     }
105 
106     /**
107      * Constructor.
108      *
109      * @param listener listener to be notified of events such as when estimation
110      *                 starts, ends or its progress significantly changes.
111      */
112     public LMedSRobustKnownFrameMagnetometerCalibrator(final RobustKnownFrameMagnetometerCalibratorListener listener) {
113         super(listener);
114     }
115 
116     /**
117      * Constructor.
118      *
119      * @param measurements list of body magnetic flux density measurements with standard
120      *                     deviations taken at different frames (positions and
121      *                     orientations).
122      */
123     public LMedSRobustKnownFrameMagnetometerCalibrator(
124             final List<StandardDeviationFrameBodyMagneticFluxDensity> measurements) {
125         super(measurements);
126     }
127 
128     /**
129      * Constructor.
130      *
131      * @param measurements list of body magnetic flux density measurements with standard
132      *                     deviations taken at different frames (positions and
133      *                     orientations).
134      * @param listener     listener to handle events raised by this calibrator.
135      */
136     public LMedSRobustKnownFrameMagnetometerCalibrator(
137             final List<StandardDeviationFrameBodyMagneticFluxDensity> measurements,
138             final RobustKnownFrameMagnetometerCalibratorListener listener) {
139         super(measurements, listener);
140     }
141 
142     /**
143      * Constructor.
144      *
145      * @param commonAxisUsed indicates whether z-axis is assumed to be common
146      *                       for the accelerometer, gyroscope and magnetometer.
147      */
148     public LMedSRobustKnownFrameMagnetometerCalibrator(final boolean commonAxisUsed) {
149         super(commonAxisUsed);
150     }
151 
152     /**
153      * Constructor.
154      *
155      * @param commonAxisUsed indicates whether z-axis is assumed to be common
156      *                       for the accelerometer, gyroscope and magnetometer.
157      * @param listener       listener to handle events raised by this calibrator.
158      */
159     public LMedSRobustKnownFrameMagnetometerCalibrator(
160             final boolean commonAxisUsed, final RobustKnownFrameMagnetometerCalibratorListener listener) {
161         super(commonAxisUsed, listener);
162     }
163 
164     /**
165      * Constructor.
166      *
167      * @param measurements   list of body magnetic flux density measurements with standard
168      *                       deviations taken at different frames (positions and
169      *                       orientations).
170      * @param commonAxisUsed indicates whether z-axis is assumed to be common
171      *                       for the accelerometer, gyroscope and magnetometer.
172      */
173     public LMedSRobustKnownFrameMagnetometerCalibrator(
174             final List<StandardDeviationFrameBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed) {
175         super(measurements, commonAxisUsed);
176     }
177 
178     /**
179      * Constructor.
180      *
181      * @param measurements   list of body magnetic flux density measurements with standard
182      *                       deviations taken at different frames (positions and
183      *                       orientations).
184      * @param commonAxisUsed indicates whether z-axis is assumed to be common
185      *                       for the accelerometer, gyroscope and magnetometer.
186      * @param listener       listener to handle events raised by this calibrator.
187      */
188     public LMedSRobustKnownFrameMagnetometerCalibrator(
189             final List<StandardDeviationFrameBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
190             final RobustKnownFrameMagnetometerCalibratorListener listener) {
191         super(measurements, commonAxisUsed, listener);
192     }
193 
194     /**
195      * Returns threshold to be used to keep the algorithm iterating in case that
196      * best estimated threshold using median of residuals is not small enough.
197      * Once a solution is found that generates a threshold below this value, the
198      * algorithm will stop.
199      * The stop threshold can be used to prevent the LMedS algorithm to iterate
200      * too many times in cases where samples have a very similar accuracy.
201      * For instance, in cases where proportion of outliers is very small (close
202      * to 0%), and samples are very accurate (i.e. 1e-6), the algorithm would
203      * iterate for a long time trying to find the best solution when indeed
204      * there is no need to do that if a reasonable threshold has already been
205      * reached.
206      * Because of this behaviour the stop threshold can be set to a value much
207      * lower than the one typically used in RANSAC, and yet the algorithm could
208      * still produce even smaller thresholds in estimated results.
209      *
210      * @return stop threshold to stop the algorithm prematurely when a certain
211      * accuracy has been reached.
212      */
213     public double getStopThreshold() {
214         return stopThreshold;
215     }
216 
217     /**
218      * Sets threshold to be used to keep the algorithm iterating in case that
219      * best estimated threshold using median of residuals is not small enough.
220      * Once a solution is found that generates a threshold below this value,
221      * the algorithm will stop.
222      * The stop threshold can be used to prevent the LMedS algorithm to iterate
223      * too many times in cases where samples have a very similar accuracy.
224      * For instance, in cases where proportion of outliers is very small (close
225      * to 0%), and samples are very accurate (i.e. 1e-6), the algorithm would
226      * iterate for a long time trying to find the best solution when indeed
227      * there is no need to do that if a reasonable threshold has already been
228      * reached.
229      * Because of this behaviour the stop threshold can be set to a value much
230      * lower than the one typically used in RANSAC, and yet the algorithm could
231      * still produce even smaller thresholds in estimated results.
232      *
233      * @param stopThreshold stop threshold to stop the algorithm prematurely
234      *                      when a certain accuracy has been reached.
235      * @throws IllegalArgumentException if provided value is zero or negative.
236      * @throws LockedException          if calibrator is currently running.
237      */
238     public void setStopThreshold(final double stopThreshold) throws LockedException {
239         if (running) {
240             throw new LockedException();
241         }
242         if (stopThreshold <= MIN_STOP_THRESHOLD) {
243             throw new IllegalArgumentException();
244         }
245 
246         this.stopThreshold = stopThreshold;
247     }
248 
249     /**
250      * Estimates magnetometer calibration parameters containing hard-iron
251      * bias and soft-iron scale factors and cross-coupling errors.
252      *
253      * @throws LockedException      if calibrator is currently running.
254      * @throws NotReadyException    if calibrator is not ready.
255      * @throws CalibrationException if estimation fails for numerical reasons.
256      */
257     @SuppressWarnings("DuplicatedCode")
258     @Override
259     public void calibrate() throws LockedException, NotReadyException, CalibrationException {
260         if (running) {
261             throw new LockedException();
262         }
263         if (!isReady()) {
264             throw new NotReadyException();
265         }
266 
267         final var innerEstimator = new LMedSRobustEstimator<>(new LMedSRobustEstimatorListener<PreliminaryResult>() {
268             @Override
269             public int getTotalSamples() {
270                 return measurements.size();
271             }
272 
273             @Override
274             public int getSubsetSize() {
275                 return preliminarySubsetSize;
276             }
277 
278             @Override
279             public void estimatePreliminarSolutions(
280                     final int[] samplesIndices, final List<PreliminaryResult> solutions) {
281                 computePreliminarySolutions(samplesIndices, solutions);
282             }
283 
284             @Override
285             public double computeResidual(final PreliminaryResult currentEstimation, final int i) {
286                 return computeError(measurements.get(i), currentEstimation);
287             }
288 
289             @Override
290             public boolean isReady() {
291                 return LMedSRobustKnownFrameMagnetometerCalibrator.this.isReady();
292             }
293 
294             @Override
295             public void onEstimateStart(final RobustEstimator<PreliminaryResult> estimator) {
296                 // no action needed
297             }
298 
299             @Override
300             public void onEstimateEnd(final RobustEstimator<PreliminaryResult> estimator) {
301                 // no action needed
302             }
303 
304             @Override
305             public void onEstimateNextIteration(
306                     final RobustEstimator<PreliminaryResult> estimator, final int iteration) {
307                 if (listener != null) {
308                     listener.onCalibrateNextIteration(
309                             LMedSRobustKnownFrameMagnetometerCalibrator.this, iteration);
310                 }
311             }
312 
313             @Override
314             public void onEstimateProgressChange(
315                     final RobustEstimator<PreliminaryResult> estimator, final float progress) {
316                 if (listener != null) {
317                     listener.onCalibrateProgressChange(
318                             LMedSRobustKnownFrameMagnetometerCalibrator.this, progress);
319                 }
320             }
321         });
322 
323         try {
324             running = true;
325 
326             if (listener != null) {
327                 listener.onCalibrateStart(this);
328             }
329 
330             inliersData = null;
331 
332             setupWmmEstimator();
333 
334             innerEstimator.setConfidence(confidence);
335             innerEstimator.setMaxIterations(maxIterations);
336             innerEstimator.setProgressDelta(progressDelta);
337             innerEstimator.setStopThreshold(stopThreshold);
338             final var preliminaryResult = innerEstimator.estimate();
339             inliersData = innerEstimator.getInliersData();
340 
341             attemptRefine(preliminaryResult);
342 
343             if (listener != null) {
344                 listener.onCalibrateEnd(this);
345             }
346 
347         } catch (final com.irurueta.numerical.LockedException e) {
348             throw new LockedException(e);
349         } catch (final com.irurueta.numerical.NotReadyException e) {
350             throw new NotReadyException(e);
351         } catch (final RobustEstimatorException | IOException e) {
352             throw new CalibrationException(e);
353         } finally {
354             running = false;
355         }
356     }
357 
358     /**
359      * Returns method being used for robust estimation.
360      *
361      * @return method being used for robust estimation.
362      */
363     @Override
364     public RobustEstimatorMethod getMethod() {
365         return RobustEstimatorMethod.LMEDS;
366     }
367 
368     /**
369      * Indicates whether this calibrator requires quality scores for each
370      * measurement or not.
371      *
372      * @return true if quality scores are required, false otherwise.
373      */
374     @Override
375     public boolean isQualityScoresRequired() {
376         return false;
377     }
378 }