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