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 }