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