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