View Javadoc
1   /*
2    * Copyright (C) 2021 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.intervals.thresholdfactor;
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.INSLooselyCoupledKalmanInitializerConfig;
22  import com.irurueta.navigation.inertial.INSTightlyCoupledKalmanInitializerConfig;
23  import com.irurueta.navigation.inertial.calibration.AccelerometerNoiseRootPsdSource;
24  import com.irurueta.navigation.inertial.calibration.BodyKinematicsSequence;
25  import com.irurueta.navigation.inertial.calibration.CalibrationException;
26  import com.irurueta.navigation.inertial.calibration.GyroscopeBiasUncertaintySource;
27  import com.irurueta.navigation.inertial.calibration.GyroscopeCalibrationSource;
28  import com.irurueta.navigation.inertial.calibration.GyroscopeNoiseRootPsdSource;
29  import com.irurueta.navigation.inertial.calibration.StandardDeviationTimedBodyKinematics;
30  import com.irurueta.navigation.inertial.calibration.TimedBodyKinematics;
31  import com.irurueta.navigation.inertial.calibration.accelerometer.KnownBiasAccelerometerCalibrator;
32  import com.irurueta.navigation.inertial.calibration.generators.GyroscopeMeasurementsGenerator;
33  import com.irurueta.navigation.inertial.calibration.generators.GyroscopeMeasurementsGeneratorListener;
34  import com.irurueta.navigation.inertial.calibration.gyroscope.GyroscopeCalibratorMeasurementOrSequenceType;
35  import com.irurueta.navigation.inertial.calibration.gyroscope.GyroscopeNonLinearCalibrator;
36  import com.irurueta.navigation.inertial.calibration.gyroscope.OrderedBodyKinematicsSequenceGyroscopeCalibrator;
37  import com.irurueta.navigation.inertial.calibration.gyroscope.QualityScoredGyroscopeCalibrator;
38  import com.irurueta.navigation.inertial.calibration.gyroscope.UnknownBiasGyroscopeCalibrator;
39  import com.irurueta.navigation.inertial.calibration.intervals.TriadStaticIntervalDetector;
40  import com.irurueta.units.Acceleration;
41  import com.irurueta.units.AccelerationUnit;
42  import com.irurueta.units.Time;
43  
44  import java.util.ArrayList;
45  import java.util.List;
46  
47  /**
48   * Optimizes the threshold factor for interval detection of gyroscope data based
49   * on results of calibration.
50   * Implementations of this class will attempt to find the best threshold factor
51   * between the provided range of values.
52   * Only gyroscope calibrators based on unknown orientation are supported (in other terms,
53   * calibrators must be {@link GyroscopeNonLinearCalibrator} and must support
54   * {@link GyroscopeCalibratorMeasurementOrSequenceType#BODY_KINEMATICS_SEQUENCE}).
55   */
56  public abstract class GyroscopeIntervalDetectorThresholdFactorOptimizer extends
57          IntervalDetectorThresholdFactorOptimizer<TimedBodyKinematics,
58                  GyroscopeIntervalDetectorThresholdFactorOptimizerDataSource>
59          implements GyroscopeNoiseRootPsdSource, GyroscopeBiasUncertaintySource,
60          GyroscopeCalibrationSource, AccelerometerNoiseRootPsdSource {
61  
62      /**
63       * Default minimum threshold factor.
64       */
65      public static final double DEFAULT_MIN_THRESHOLD_FACTOR = 2.0;
66  
67      /**
68       * Default maximum threshold factor.
69       */
70      public static final double DEFAULT_MAX_THRESHOLD_FACTOR = 10.0;
71  
72      /**
73       * Minimum threshold factor.
74       */
75      protected double minThresholdFactor = DEFAULT_MIN_THRESHOLD_FACTOR;
76  
77      /**
78       * Maximum threshold factor.
79       */
80      protected double maxThresholdFactor = DEFAULT_MAX_THRESHOLD_FACTOR;
81  
82      /**
83       * Gyroscope calibrator.
84       */
85      private GyroscopeNonLinearCalibrator calibrator;
86  
87      /**
88       * A measurement generator for gyroscope calibrators.
89       */
90      private GyroscopeMeasurementsGenerator generator;
91  
92      /**
93       * Generated sequences to be used for gyroscope calibration.
94       */
95      private List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences;
96  
97      /**
98       * Mapper to convert {@link BodyKinematicsSequence} sequences of {@link StandardDeviationTimedBodyKinematics}
99       * into quality scores.
100      */
101     private QualityScoreMapper<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> qualityScoreMapper =
102             new DefaultGyroscopeQualityScoreMapper();
103 
104     /**
105      * Estimated norm of gyroscope noise root PSD (Power Spectral Density)
106      * expressed as (rad * s^-0.5).
107      */
108     private double angularSpeedNoiseRootPsd;
109 
110     /**
111      * Accelerometer base noise level that has been detected during
112      * initialization of the best solution that has been found expressed in
113      * meters per squared second (m/s^2).
114      * This is equal to the standard deviation of the accelerometer measurements
115      * during the initialization phase.
116      */
117     private double baseNoiseLevel;
118 
119     /**
120      * Threshold to determine static/dynamic period changes expressed in
121      * meters per squared second (m/s^2) for the best calibration solution that
122      * has been found.
123      */
124     private double threshold;
125 
126     /**
127      * Estimated covariance matrix for estimated parameters.
128      */
129     private Matrix estimatedCovariance;
130 
131     /**
132      * Estimated angular rate biases for each IMU axis expressed in radians per
133      * second (rad/s).
134      */
135     private double[] estimatedBiases;
136 
137     /**
138      * Estimated gyroscope scale factors and cross-coupling errors.
139      * This is the product of matrix Tg containing cross-coupling errors and Kg
140      * containing scaling factors.
141      * So that:
142      * <pre>
143      *     Mg = [sx    mxy  mxz] = Tg*Kg
144      *          [myx   sy   myz]
145      *          [mzx   mzy  sz ]
146      * </pre>
147      * Where:
148      * <pre>
149      *     Kg = [sx 0   0 ]
150      *          [0  sy  0 ]
151      *          [0  0   sz]
152      * </pre>
153      * and
154      * <pre>
155      *     Tg = [1          -alphaXy    alphaXz ]
156      *          [alphaYx    1           -alphaYz]
157      *          [-alphaZx   alphaZy     1       ]
158      * </pre>
159      * Hence:
160      * <pre>
161      *     Mg = [sx    mxy  mxz] = Tg*Kg =  [sx             -sy * alphaXy   sz * alphaXz ]
162      *          [myx   sy   myz]            [sx * alphaYx   sy              -sz * alphaYz]
163      *          [mzx   mzy  sz ]            [-sx * alphaZx  sy * alphaZy    sz           ]
164      * </pre>
165      * This instance allows any 3x3 matrix. However, typically alphaYx, alphaZx and alphaZy
166      * are considered to be zero if the gyroscope z-axis is assumed to be the same
167      * as the body z-axis. When this is assumed, myx = mzx = mzy = 0 and the Mg matrix
168      * becomes upper diagonal:
169      * <pre>
170      *     Mg = [sx    mxy  mxz]
171      *          [0     sy   myz]
172      *          [0     0    sz ]
173      * </pre>
174      * Values of this matrix are unit-less.
175      */
176     private Matrix estimatedMg;
177 
178     /**
179      * Estimated G-dependent cross-biases introduced on the gyroscope by the
180      * specific forces sensed by the accelerometer.
181      * This instance allows any 3x3 matrix.
182      */
183     private Matrix estimatedGg;
184 
185     /**
186      * Constructor.
187      */
188     protected GyroscopeIntervalDetectorThresholdFactorOptimizer() {
189         initialize();
190     }
191 
192     /**
193      * Constructor.
194      *
195      * @param dataSource instance in charge of retrieving data for this optimizer.
196      */
197     protected GyroscopeIntervalDetectorThresholdFactorOptimizer(
198             final GyroscopeIntervalDetectorThresholdFactorOptimizerDataSource dataSource) {
199         super(dataSource);
200         initialize();
201     }
202 
203     /**
204      * Constructor.
205      *
206      * @param calibrator a gyroscope calibrator to be used to optimize its
207      *                   Mean Square Error (MSE).
208      * @throws IllegalArgumentException if gyroscope calibrator does not use
209      *                                  {@link BodyKinematicsSequence} sequences of
210      *                                  {@link StandardDeviationTimedBodyKinematics}.
211      */
212     protected GyroscopeIntervalDetectorThresholdFactorOptimizer(final GyroscopeNonLinearCalibrator calibrator) {
213         try {
214             setCalibrator(calibrator);
215         } catch (final LockedException ignore) {
216             // never happens
217         }
218         initialize();
219     }
220 
221     /**
222      * Constructor.
223      *
224      * @param dataSource instance in charge of retrieving data for this optimizer.
225      * @param calibrator a gyroscope calibrator to be used to optimize its
226      *                   Mean Square Error (MSE).
227      * @throws IllegalArgumentException if gyroscope calibrator does not use
228      *                                  {@link BodyKinematicsSequence} sequences of
229      *                                  {@link StandardDeviationTimedBodyKinematics}.
230      */
231     protected GyroscopeIntervalDetectorThresholdFactorOptimizer(
232             final GyroscopeIntervalDetectorThresholdFactorOptimizerDataSource dataSource,
233             final GyroscopeNonLinearCalibrator calibrator) {
234         super(dataSource);
235         try {
236             setCalibrator(calibrator);
237         } catch (final LockedException ignore) {
238             // never happens
239         }
240         initialize();
241     }
242 
243     /**
244      * Gets provided gyroscope calibrator to be used to optimize its Mean Square Error (MSE).
245      *
246      * @return gyroscope calibrator to be used to optimize its MSE.
247      */
248     public GyroscopeNonLinearCalibrator getCalibrator() {
249         return calibrator;
250     }
251 
252     /**
253      * Sets gyroscope calibrator to be used to optimize its Mean Square Error (MSE).
254      *
255      * @param calibrator gyroscope calibrator to be use dto optimize its MSE.
256      * @throws LockedException          if optimizer is already running.
257      * @throws IllegalArgumentException if gyroscope calibrator does not use
258      *                                  {@link BodyKinematicsSequence} sequences of
259      *                                  {@link StandardDeviationTimedBodyKinematics}.
260      */
261     public void setCalibrator(final GyroscopeNonLinearCalibrator calibrator) throws LockedException {
262         if (running) {
263             throw new LockedException();
264         }
265 
266         if (calibrator != null && calibrator.getMeasurementOrSequenceType()
267                 != GyroscopeCalibratorMeasurementOrSequenceType.BODY_KINEMATICS_SEQUENCE) {
268             throw new IllegalArgumentException();
269         }
270 
271         this.calibrator = calibrator;
272     }
273 
274     /**
275      * Gets mapper to convert {@link BodyKinematicsSequence} sequences of {@link StandardDeviationTimedBodyKinematics}
276      * into quality scores.
277      *
278      * @return mapper to convert sequences into quality scores.
279      */
280     public QualityScoreMapper<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> getQualityScoreMapper() {
281         return qualityScoreMapper;
282     }
283 
284     /**
285      * Sets mapper to convert {@link BodyKinematicsSequence} sequences of {@link StandardDeviationTimedBodyKinematics}
286      * into quality scores.
287      *
288      * @param qualityScoreMapper mapper to convert sequences into quality scores.
289      * @throws LockedException if optimizer is already running.
290      */
291     public void setQualityScoreMapper(
292             final QualityScoreMapper<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> qualityScoreMapper)
293             throws LockedException {
294         if (running) {
295             throw new LockedException();
296         }
297 
298         this.qualityScoreMapper = qualityScoreMapper;
299     }
300 
301     /**
302      * Gets the minimum threshold factor.
303      *
304      * @return minimum threshold factor.
305      */
306     public double getMinThresholdFactor() {
307         return minThresholdFactor;
308     }
309 
310     /**
311      * Gets the maximum threshold factor.
312      *
313      * @return maximum threshold factor.
314      */
315     public double getMaxThresholdFactor() {
316         return maxThresholdFactor;
317     }
318 
319     /**
320      * Sets a range of threshold factor values to get an optimized
321      * threshold factor value.
322      *
323      * @param minThresholdFactor minimum threshold.
324      * @param maxThresholdFactor maximum threshold.
325      * @throws LockedException          if optimizer is already running.
326      * @throws IllegalArgumentException if either minimum or maximum values are
327      *                                  negative, or if the minimum value is larger
328      *                                  than the maximum one.
329      */
330     public void setThresholdFactorRange(final double minThresholdFactor, final double maxThresholdFactor)
331             throws LockedException {
332         if (running) {
333             throw new LockedException();
334         }
335         if (minThresholdFactor < 0.0 || maxThresholdFactor < 0.0 || minThresholdFactor >= maxThresholdFactor) {
336             throw new IllegalArgumentException();
337         }
338 
339         this.minThresholdFactor = minThresholdFactor;
340         this.maxThresholdFactor = maxThresholdFactor;
341     }
342 
343     /**
344      * Indicates whether this optimizer is ready to start optimization.
345      *
346      * @return true if this optimizer is ready, false otherwise.
347      */
348     @Override
349     public boolean isReady() {
350         return super.isReady() && calibrator != null && qualityScoreMapper != null;
351     }
352 
353     /**
354      * Gets the time interval between input measurements provided to the
355      * {@link #getDataSource()} expressed in seconds (s).
356      *
357      * @return time interval between input measurements.
358      */
359     public double getTimeInterval() {
360         return generator.getTimeInterval();
361     }
362 
363     /**
364      * Sets time interval between input measurements provided to the
365      * {@link #getDataSource()} expressed in seconds (s).
366      *
367      * @param timeInterval time interval between input measurements.
368      * @throws LockedException          if optimizer is already running.
369      * @throws IllegalArgumentException if provided value is negative.
370      */
371     public void setTimeInterval(final double timeInterval) throws LockedException {
372         if (running) {
373             throw new LockedException();
374         }
375 
376         generator.setTimeInterval(timeInterval);
377     }
378 
379     /**
380      * Gets the time interval between input measurements provided to the
381      * {@link #getDataSource()}.
382      *
383      * @return time interval between input measurements.
384      */
385     public Time getTimeIntervalAsTime() {
386         return generator.getTimeIntervalAsTime();
387     }
388 
389     /**
390      * Gets the time interval between input measurements provided to the
391      * {@link #getDataSource()}.
392      *
393      * @param result instance where the time interval will be stored.
394      */
395     public void getTimeIntervalAsTime(final Time result) {
396         generator.getTimeIntervalAsTime(result);
397     }
398 
399     /**
400      * Sets time interval between input measurements provided to the
401      * {@link #getDataSource()}.
402      *
403      * @param timeInterval time interval between input measurements.
404      * @throws LockedException          if optimizer is already running.
405      * @throws IllegalArgumentException if provided value is negative.
406      */
407     public void setTimeInterval(final Time timeInterval) throws LockedException {
408         if (running) {
409             throw new LockedException();
410         }
411 
412         generator.setTimeInterval(timeInterval);
413     }
414 
415     /**
416      * Gets minimum number of input measurements provided to the
417      * {@link #getDataSource()} required in a static interval to be taken
418      * into account.
419      * Smaller static intervals will be discarded.
420      *
421      * @return a minimum number of input measurements required in a static interval
422      * to be taken into account.
423      */
424     public int getMinStaticSamples() {
425         return generator.getMinStaticSamples();
426     }
427 
428     /**
429      * Sets minimum number of input measurements provided to the
430      * {@link #getDataSource()} required in a static interval to be taken
431      * into account.
432      * Smaller static intervals will be discarded.
433      *
434      * @param minStaticSamples a minimum number of input measurements required in
435      *                         a static interval to be taken into account.
436      * @throws LockedException          if optimizer is already running.
437      * @throws IllegalArgumentException if provided value is less than 2.
438      */
439     public void setMinStaticSamples(final int minStaticSamples) throws LockedException {
440         if (running) {
441             throw new LockedException();
442         }
443 
444         generator.setMinStaticSamples(minStaticSamples);
445     }
446 
447     /**
448      * Gets maximum number of input measurements provided to the
449      * {@link #getDataSource()} allowed in dynamic intervals.
450      *
451      * @return maximum number of input measurements allowed in dynamic intervals.
452      */
453     public int getMaxDynamicSamples() {
454         return generator.getMaxDynamicSamples();
455     }
456 
457     /**
458      * Sets maximum number of input measurements provided to the
459      * {@link #getDataSource()} allowed in dynamic intervals.
460      *
461      * @param maxDynamicSamples maximum number of input measurements allowed in
462      *                          dynamic intervals.
463      * @throws LockedException          if optimizer is already running.
464      * @throws IllegalArgumentException if provided value is less than 2.
465      */
466     public void setMaxDynamicSamples(final int maxDynamicSamples) throws LockedException {
467         if (running) {
468             throw new LockedException();
469         }
470 
471         generator.setMaxDynamicSamples(maxDynamicSamples);
472     }
473 
474     /**
475      * Gets length of number of input measurements provided to the
476      * {@link #getDataSource()} to keep within the window being processed
477      * to determine instantaneous accelerometer noise level.
478      *
479      * @return length of input measurements to keep within the window.
480      */
481     public int getWindowSize() {
482         return generator.getWindowSize();
483     }
484 
485     /**
486      * Sets length of number of input measurements provided to the
487      * {@link #getDataSource()} to keep within the window being processed
488      * to determine instantaneous noise level.
489      * Window size must always be larger than the allowed minimum value, which is 2
490      * and must have an odd value.
491      *
492      * @param windowSize length of number of samples to keep within the window.
493      * @throws LockedException          if optimizer is already running.
494      * @throws IllegalArgumentException if provided value is not valid.
495      */
496     public void setWindowSize(final int windowSize) throws LockedException {
497         if (running) {
498             throw new LockedException();
499         }
500 
501         generator.setWindowSize(windowSize);
502     }
503 
504     /**
505      * Gets number of input measurements provided to the
506      * {@link #getDataSource()} to be processed initially while keeping the sensor
507      * static to find the base noise level when the device is static.
508      *
509      * @return number of samples to be processed initially.
510      */
511     public int getInitialStaticSamples() {
512         return generator.getInitialStaticSamples();
513     }
514 
515     /**
516      * Sets number of input parameters provided to the {@link #getDataSource()}
517      * to be processed initially while keeping the sensor static to
518      * find the base noise level when the device is static.
519      *
520      * @param initialStaticSamples number of samples ot be processed initially.
521      * @throws LockedException          if optimizer is already running.
522      * @throws IllegalArgumentException if provided value is less than
523      *                                  {@link TriadStaticIntervalDetector#MINIMUM_INITIAL_STATIC_SAMPLES}.
524      */
525     public void setInitialStaticSamples(final int initialStaticSamples) throws LockedException {
526         if (running) {
527             throw new LockedException();
528         }
529 
530         generator.setInitialStaticSamples(initialStaticSamples);
531     }
532 
533     /**
534      * Gets factor to determine that a sudden movement has occurred during
535      * initialization if instantaneous noise level exceeds accumulated noise
536      * level by this factor amount.
537      * This factor is unit-less.
538      *
539      * @return factor to determine that a sudden movement has occurred.
540      */
541     public double getInstantaneousNoiseLevelFactor() {
542         return generator.getInstantaneousNoiseLevelFactor();
543     }
544 
545     /**
546      * Sets factor to determine that a sudden movement has occurred during
547      * initialization if instantaneous noise level exceeds accumulated noise
548      * level by this factor amount.
549      * This factor is unit-less.
550      *
551      * @param instantaneousNoiseLevelFactor factor to determine that a sudden
552      *                                      movement has occurred during
553      *                                      initialization.
554      * @throws LockedException          if optimizer is already running.
555      * @throws IllegalArgumentException if provided value is zero or negative.
556      */
557     public void setInstantaneousNoiseLevelFactor(final double instantaneousNoiseLevelFactor) throws LockedException {
558         if (running) {
559             throw new LockedException();
560         }
561 
562         generator.setInstantaneousNoiseLevelFactor(instantaneousNoiseLevelFactor);
563     }
564 
565     /**
566      * Gets the overall absolute threshold to determine whether there has been
567      * excessive motion during the whole initialization phase.
568      * This threshold is expressed in meters per squared second (m/s^2).
569      *
570      * @return overall absolute threshold to determine whether there has been
571      * excessive motion.
572      */
573     public double getBaseNoiseLevelAbsoluteThreshold() {
574         return generator.getBaseNoiseLevelAbsoluteThreshold();
575     }
576 
577     /**
578      * Sets the overall absolute threshold to determine whether there has been
579      * excessive motion during the whole initialization phase.
580      * This threshold is expressed in meters per squared second (m/s^2).
581      *
582      * @param baseNoiseLevelAbsoluteThreshold overall absolute threshold to
583      *                                        determine whether there has been
584      *                                        excessive motion.
585      * @throws LockedException          if optimizer is already running.
586      * @throws IllegalArgumentException if provided value is zero or negative.
587      */
588     public void setBaseNoiseLevelAbsoluteThreshold(final double baseNoiseLevelAbsoluteThreshold)
589             throws LockedException {
590         if (running) {
591             throw new LockedException();
592         }
593 
594         generator.setBaseNoiseLevelAbsoluteThreshold(baseNoiseLevelAbsoluteThreshold);
595     }
596 
597     /**
598      * Gets the overall absolute threshold to determine whether there has been
599      * excessive motion during the whole initialization phase.
600      *
601      * @return overall absolute threshold to determine whether there has been
602      * excessive motion.
603      */
604     public Acceleration getBaseNoiseLevelAbsoluteThresholdAsMeasurement() {
605         return generator.getBaseNoiseLevelAbsoluteThresholdAsMeasurement();
606     }
607 
608     /**
609      * Gets the overall absolute threshold to determine whether there has been
610      * excessive motion during the whole initialization phase.
611      *
612      * @param result instance where the result will be stored.
613      */
614     public void getBaseNoiseLevelAbsoluteThresholdAsMeasurement(final Acceleration result) {
615         generator.getBaseNoiseLevelAbsoluteThresholdAsMeasurement(result);
616     }
617 
618     /**
619      * Sets the overall absolute threshold to determine whether there has been
620      * excessive motion during the whole initialization phase.
621      *
622      * @param baseNoiseLevelAbsoluteThreshold overall absolute threshold to
623      *                                        determine whether there has been
624      *                                        excessive motion.
625      * @throws LockedException          if optimizer is already running.
626      * @throws IllegalArgumentException if provided value is zero or negative.
627      */
628     public void setBaseNoiseLevelAbsoluteThreshold(final Acceleration baseNoiseLevelAbsoluteThreshold)
629             throws LockedException {
630         if (running) {
631             throw new LockedException();
632         }
633 
634         generator.setBaseNoiseLevelAbsoluteThreshold(baseNoiseLevelAbsoluteThreshold);
635     }
636 
637     /**
638      * Gets accelerometer base noise level that has been detected during
639      * initialization of the best solution that has been found expressed in
640      * meters per squared second (m/s^2).
641      * This is equal to the standard deviation of the accelerometer measurements
642      * during the initialization phase.
643      *
644      * @return accelerometer base noise level of the best solution that has been
645      * found.
646      */
647     public double getAccelerometerBaseNoiseLevel() {
648         return baseNoiseLevel;
649     }
650 
651     /**
652      * Gets accelerometer base noise level that has been detected during
653      * initialization of the best solution that has been found.
654      * This is equal to the standard deviation of the accelerometer measurements
655      * during the initialization phase.
656      *
657      * @return accelerometer base noise level of the best solution that has been
658      * found.
659      */
660     public Acceleration getAccelerometerBaseNoiseLevelAsMeasurement() {
661         return createMeasurement(baseNoiseLevel, getDefaultUnit());
662     }
663 
664     /**
665      * Gets accelerometer base noise level that has been detected during
666      * initialization of the best solution that has been found.
667      * This is equal to the standard deviation of the accelerometer measurements
668      * during the initialization phase.
669      *
670      * @param result instance where the result will be stored.
671      */
672     public void getAccelerometerBaseNoiseLevelAsMeasurement(final Acceleration result) {
673         result.setValue(baseNoiseLevel);
674         result.setUnit(getDefaultUnit());
675     }
676 
677     /**
678      * Gets gyroscope base noise level PSD (Power Spectral Density)
679      * expressed in (rad^2/s).
680      *
681      * @return gyroscope base noise level PSD.
682      */
683     public double getGyroscopeBaseNoiseLevelPsd() {
684         return angularSpeedNoiseRootPsd * angularSpeedNoiseRootPsd;
685     }
686 
687     /**
688      * Gets gyroscope base noise level root PSD (Power Spectral Density)
689      * expressed in (rad * s^-0.5)
690      *
691      * @return gyroscope base noise level root PSD.
692      */
693     @Override
694     public double getGyroscopeBaseNoiseLevelRootPsd() {
695         return angularSpeedNoiseRootPsd;
696     }
697 
698     /**
699      * Gets accelerometer base noise level PSD (Power Spectral Density)
700      * expressed in (m^2 * s^-3).
701      *
702      * @return accelerometer base noise level PSD.
703      */
704     public double getAccelerometerBaseNoiseLevelPsd() {
705         return baseNoiseLevel * baseNoiseLevel * getTimeInterval();
706     }
707 
708     /**
709      * Gets accelerometer base noise level root PSD (Power Spectral Density)
710      * expressed in (m * s^-1.5).
711      *
712      * @return accelerometer base noise level root PSD.
713      */
714     @Override
715     public double getAccelerometerBaseNoiseLevelRootPsd() {
716         return baseNoiseLevel * Math.sqrt(getTimeInterval());
717     }
718 
719     /**
720      * Gets the threshold to determine static/dynamic period changes expressed in
721      * meters per squared second (m/s^2) for the best calibration solution that
722      * has been found.
723      *
724      * @return threshold to determine static/dynamic period changes for the best
725      * solution.
726      */
727     public double getThreshold() {
728         return threshold;
729     }
730 
731     /**
732      * Gets the threshold to determine static/dynamic period changes for the best
733      * calibration solution that has been found.
734      *
735      * @return threshold to determine static/dynamic period changes for the best
736      * solution.
737      */
738     public Acceleration getThresholdAsMeasurement() {
739         return createMeasurement(threshold, getDefaultUnit());
740     }
741 
742     /**
743      * Get the threshold to determine static/dynamic period changes for the best
744      * calibration solution that has been found.
745      *
746      * @param result instance where the result will be stored.
747      */
748     public void getThresholdAsMeasurement(final Acceleration result) {
749         result.setValue(threshold);
750         result.setUnit(getDefaultUnit());
751     }
752 
753     /**
754      * Gets estimated standard deviation norm of gyroscope bias expressed in
755      * radians per second (rad/s).
756      * This can be used as the initial gyroscope bias uncertainty for
757      * {@link INSLooselyCoupledKalmanInitializerConfig} or {@link INSTightlyCoupledKalmanInitializerConfig}.
758      *
759      * @return estimated standard deviation norm of gyroscope bias or null
760      * if not available.
761      */
762     @Override
763     public Double getEstimatedBiasStandardDeviationNorm() {
764         return estimatedCovariance != null
765                 ? Math.sqrt(getEstimatedBiasXVariance() + getEstimatedBiasYVariance() + getEstimatedBiasZVariance())
766                 : null;
767     }
768 
769     /**
770      * Gets estimated x coordinate variance of gyroscope bias expressed in (rad^2/s^2).
771      *
772      * @return estimated x coordinate variance of gyroscope bias or null if not available.
773      */
774     public Double getEstimatedBiasXVariance() {
775         return estimatedCovariance != null ? estimatedCovariance.getElementAt(0, 0) : null;
776     }
777 
778     /**
779      * Gets estimated y coordinate variance of gyroscope bias expressed in (rad^2/s^2).
780      *
781      * @return estimated y coordinate variance of gyroscope bias or null if not available.
782      */
783     public Double getEstimatedBiasYVariance() {
784         return estimatedCovariance != null ? estimatedCovariance.getElementAt(1, 1) : null;
785     }
786 
787     /**
788      * Gets estimated z coordinate variance of gyroscope bias expressed in (rad^2/s^2).
789      *
790      * @return estimated z coordinate variance of gyroscope bias or null if not available.
791      */
792     public Double getEstimatedBiasZVariance() {
793         return estimatedCovariance != null ? estimatedCovariance.getElementAt(2, 2) : null;
794     }
795 
796     /**
797      * Gets the array containing x,y,z components of estimated gyroscope biases
798      * expressed in radians per second (rad/s).
799      *
800      * @return array containing x,y,z components of estimated gyroscope biases.
801      */
802     @Override
803     public double[] getEstimatedBiases() {
804         return estimatedBiases;
805     }
806 
807     /**
808      * Gets estimated gyroscope scale factors and cross-coupling errors.
809      * This is the product of matrix Tg containing cross-coupling errors and Kg
810      * containing scaling factors.
811      * So that:
812      * <pre>
813      *     Mg = [sx    mxy  mxz] = Tg*Kg
814      *          [myx   sy   myz]
815      *          [mzx   mzy  sz ]
816      * </pre>
817      * Where:
818      * <pre>
819      *     Kg = [sx 0   0 ]
820      *          [0  sy  0 ]
821      *          [0  0   sz]
822      * </pre>
823      * and
824      * <pre>
825      *     Tg = [1          -alphaXy    alphaXz ]
826      *          [alphaYx    1           -alphaYz]
827      *          [-alphaZx   alphaZy     1       ]
828      * </pre>
829      * Hence:
830      * <pre>
831      *     Mg = [sx    mxy  mxz] = Tg*Kg =  [sx             -sy * alphaXy   sz * alphaXz ]
832      *          [myx   sy   myz]            [sx * alphaYx   sy              -sz * alphaYz]
833      *          [mzx   mzy  sz ]            [-sx * alphaZx  sy * alphaZy    sz           ]
834      * </pre>
835      * This instance allows any 3x3 matrix. However, typically alphaYx, alphaZx and alphaZy
836      * are considered to be zero if the gyroscope z-axis is assumed to be the same
837      * as the body z-axis. When this is assumed, myx = mzx = mzy = 0 and the Mg matrix
838      * becomes upper diagonal:
839      * <pre>
840      *     Mg = [sx    mxy  mxz]
841      *          [0     sy   myz]
842      *          [0     0    sz ]
843      * </pre>
844      * Values of this matrix are unit-less.
845      *
846      * @return estimated gyroscope scale factors and cross-coupling errors.
847      */
848     @Override
849     public Matrix getEstimatedMg() {
850         return estimatedMg;
851     }
852 
853     /**
854      * Gets estimated G-dependent cross-biases introduced on the gyroscope by the
855      * specific forces sensed by the accelerometer.
856      *
857      * @return a 3x3 matrix containing g-dependent cross biases.
858      */
859     @Override
860     public Matrix getEstimatedGg() {
861         return estimatedGg;
862     }
863 
864     /**
865      * Evaluates calibration Mean Square Error (MSE) for the provided threshold factor.
866      *
867      * @param thresholdFactor threshold factor to be used for interval detection
868      *                        and measurement generation to be used for
869      *                        calibration.
870      * @return calibration MSE.
871      * @throws LockedException                                   if the generator is busy.
872      * @throws CalibrationException                              if calibration fails.
873      * @throws NotReadyException                                 if the calibrator is not ready.
874      * @throws IntervalDetectorThresholdFactorOptimizerException interval detection failed.
875      */
876     protected double evaluateForThresholdFactor(final double thresholdFactor) throws LockedException,
877             CalibrationException, NotReadyException, IntervalDetectorThresholdFactorOptimizerException {
878         if (sequences == null) {
879             sequences = new ArrayList<>();
880         } else {
881             sequences.clear();
882         }
883 
884         generator.reset();
885         generator.setThresholdFactor(thresholdFactor);
886 
887         var count = dataSource.count();
888         var failed = false;
889         for (var i = 0; i < count; i++) {
890             final var bodyKinematics = dataSource.getAt(i);
891             if (!generator.process(bodyKinematics)) {
892                 failed = true;
893                 break;
894             }
895         }
896 
897         if (failed || generator.getStatus() == TriadStaticIntervalDetector.Status.FAILED) {
898             // interval detection failed
899             return Double.MAX_VALUE;
900         }
901 
902         // check that enough measurements have been obtained
903         if (sequences.size() < calibrator.getMinimumRequiredMeasurementsOrSequences()) {
904             return Double.MAX_VALUE;
905         }
906 
907         // set calibrator measurements
908         if (calibrator.getMeasurementOrSequenceType()
909                 == GyroscopeCalibratorMeasurementOrSequenceType.BODY_KINEMATICS_SEQUENCE) {
910             final var cal = (OrderedBodyKinematicsSequenceGyroscopeCalibrator) this.calibrator;
911 
912             cal.setSequences(sequences);
913         } else {
914             throw new IntervalDetectorThresholdFactorOptimizerException();
915         }
916 
917         if (calibrator.isQualityScoresRequired()) {
918             final var cal = (QualityScoredGyroscopeCalibrator) this.calibrator;
919 
920             final var size = sequences.size();
921             final var qualityScores = new double[size];
922             for (var i = 0; i < size; i++) {
923                 qualityScores[i] = qualityScoreMapper.map(sequences.get(i));
924             }
925             cal.setQualityScores(qualityScores);
926         }
927 
928         calibrator.calibrate();
929 
930         final var mse = calibrator.getEstimatedMse();
931         if (mse < minMse) {
932             keepBestResult(mse, thresholdFactor);
933         }
934         return mse;
935     }
936 
937     /**
938      * Initializes gyroscope measurement generator to convert timed body kinematics
939      * measurements after interval detection into measurements used for
940      * gyroscope calibration.
941      */
942     private void initialize() {
943         final var generatorListener = new GyroscopeMeasurementsGeneratorListener() {
944             @Override
945             public void onInitializationStarted(final GyroscopeMeasurementsGenerator generator) {
946                 // not needed
947             }
948 
949             @Override
950             public void onInitializationCompleted(
951                     final GyroscopeMeasurementsGenerator generator,
952                     final double baseNoiseLevel) {
953                 // not needed
954             }
955 
956             @Override
957             public void onError(
958                     final GyroscopeMeasurementsGenerator generator,
959                     final TriadStaticIntervalDetector.ErrorReason reason) {
960                 // not needed
961             }
962 
963             @Override
964             public void onStaticIntervalDetected(final GyroscopeMeasurementsGenerator generator) {
965                 // not needed
966             }
967 
968             @Override
969             public void onDynamicIntervalDetected(final GyroscopeMeasurementsGenerator generator) {
970                 // not needed
971             }
972 
973             @Override
974             public void onStaticIntervalSkipped(final GyroscopeMeasurementsGenerator generator) {
975                 // not needed
976             }
977 
978             @Override
979             public void onDynamicIntervalSkipped(final GyroscopeMeasurementsGenerator generator) {
980                 // not needed
981             }
982 
983             @Override
984             public void onGeneratedMeasurement(
985                     final GyroscopeMeasurementsGenerator generator,
986                     final BodyKinematicsSequence<StandardDeviationTimedBodyKinematics> sequence) {
987                 sequences.add(sequence);
988             }
989 
990             @Override
991             public void onReset(final GyroscopeMeasurementsGenerator generator) {
992                 // not needed
993             }
994         };
995 
996         generator = new GyroscopeMeasurementsGenerator(generatorListener);
997     }
998 
999     /**
1000      * Keeps the best calibration solution found so far.
1001      *
1002      * @param mse             Estimated Mean Square Error during calibration.
1003      * @param thresholdFactor threshold factor to be kept.
1004      */
1005     private void keepBestResult(final double mse, final double thresholdFactor) {
1006         minMse = mse;
1007         optimalThresholdFactor = thresholdFactor;
1008 
1009         angularSpeedNoiseRootPsd = generator.getGyroscopeBaseNoiseLevelRootPsd();
1010         baseNoiseLevel = generator.getAccelerometerBaseNoiseLevel();
1011         threshold = generator.getThreshold();
1012 
1013         if (estimatedCovariance == null) {
1014             estimatedCovariance = new Matrix(calibrator.getEstimatedCovariance());
1015         } else {
1016             estimatedCovariance.copyFrom(calibrator.getEstimatedCovariance());
1017         }
1018         if (estimatedMg == null) {
1019             estimatedMg = new Matrix(calibrator.getEstimatedMg());
1020         } else {
1021             estimatedMg.copyFrom(calibrator.getEstimatedMg());
1022         }
1023         if (estimatedGg == null) {
1024             estimatedGg = new Matrix(calibrator.getEstimatedGg());
1025         } else {
1026             estimatedGg.copyFrom(calibrator.getEstimatedGg());
1027         }
1028         if (calibrator instanceof UnknownBiasGyroscopeCalibrator unknownBiasGyroscopeCalibrator) {
1029             estimatedBiases = unknownBiasGyroscopeCalibrator.getEstimatedBiases();
1030         } else if (calibrator instanceof KnownBiasAccelerometerCalibrator knownBiasAccelerometerCalibrator) {
1031             estimatedBiases = knownBiasAccelerometerCalibrator.getBias();
1032         }
1033     }
1034 
1035     /**
1036      * Creates an acceleration instance using the provided value and unit.
1037      *
1038      * @param value value of measurement.
1039      * @param unit  unit of value.
1040      * @return created acceleration.
1041      */
1042     private Acceleration createMeasurement(final double value, final AccelerationUnit unit) {
1043         return new Acceleration(value, unit);
1044     }
1045 
1046     /**
1047      * Gets the default unit for acceleration, which is meters per
1048      * squared second (m/s^2).
1049      *
1050      * @return default unit for acceleration.
1051      */
1052     private AccelerationUnit getDefaultUnit() {
1053         return AccelerationUnit.METERS_PER_SQUARED_SECOND;
1054     }
1055 }