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 }