OSDN Git Service

Winter cleanup: suffix our own types with _t
[android-x86/hardware-intel-libsensors.git] / compass-calibration.c
1 /*
2  * Copyright (C) 2014 Intel Corporation.
3  */
4
5 #include <errno.h>
6 #include <fcntl.h>
7 #include <math.h>
8 #include <hardware/sensors.h>
9 #include <sys/stat.h>
10 #include <string.h>
11 #include <utils/Log.h>
12 #include "calibration.h"
13 #include "matrix-ops.h"
14 #include "description.h"
15
16 #ifdef DBG_RAW_DATA
17 #define MAX_RAW_DATA_COUNT 2000
18 #define RAW_DATA_FULL_PATH "/data/raw_compass_data_full_%d.txt"
19 #define RAW_DATA_SELECTED_PATH "/data/raw_compass_data_selected_%d.txt"
20 static FILE *raw_data = NULL;
21 static FILE *raw_data_selected = NULL;
22 static int raw_data_count = 0;
23 int file_no = 0;
24 #endif
25
26 /* compass defines */
27 #define COMPASS_CALIBRATION_PATH "/data/compass.conf"
28 #define EPSILON 0.000000001
29
30 #define MAGNETIC_LOW 960 /* 31 micro tesla squared */
31 #define CAL_STEPS 5
32
33 /* We'll have multiple calibration levels
34 *  so that we can provide an estimation as fast as possible
35 */
36 static const float min_diffs[CAL_STEPS] =  {0.2,  0.25, 0.4, 0.6, 1.0 };
37 static const float max_sqr_errs[CAL_STEPS] = {10.0, 10.0, 8.0, 5.0, 3.5 };
38 static const unsigned int lookback_counts[CAL_STEPS] = {2, 3, 4, 5, 6 };
39
40 /* reset calibration algorithm */
41 static void reset_sample (struct compass_cal_t* data)
42 {
43     int i,j;
44     data->sample_count = 0;
45     for (i = 0; i < MAGN_DS_SIZE; i++)
46         for (j=0; j < 3; j++)
47             data->sample[i][j] = 0;
48
49     data->average[0] = data->average[1] = data->average[2] = 0;
50 }
51
52 static double calc_square_err (struct compass_cal_t* data)
53 {
54     double err = 0;
55     double raw[3][1], result[3][1], mat_diff[3][1];
56     int i;
57     float stdev[3] = {0,0,0};
58
59     for (i = 0; i < MAGN_DS_SIZE; i++) {
60         raw[0][0] = data->sample[i][0];
61         raw[1][0] = data->sample[i][1];
62         raw[2][0] = data->sample[i][2];
63
64         stdev[0] += (raw[0][0] - data->average[0]) * (raw[0][0] - data->average[0]);
65         stdev[1] += (raw[1][0] - data->average[1]) * (raw[1][0] - data->average[1]);
66         stdev[2] += (raw[2][0] - data->average[2]) * (raw[2][0] - data->average[2]);
67
68         substract (3, 1, raw, data->offset, mat_diff);
69         multiply(3, 3, 1, data->w_invert, mat_diff, result);
70
71         double diff = sqrt(result[0][0] * result[0][0] + result[1][0] * result[1][0]
72              + result[2][0] * result[2][0]) - data->bfield;
73
74         err += diff * diff;
75     }
76
77     stdev[0] = sqrt(stdev[0] / MAGN_DS_SIZE);
78     stdev[1] = sqrt(stdev[1] / MAGN_DS_SIZE);
79     stdev[2] = sqrt(stdev[2] / MAGN_DS_SIZE);
80
81     /*
82      * A sanity check - if we have too little variation for an axis
83      * it's best to reject the calibration than risking a wrong calibration.
84      */
85     if (stdev[0] <= 1 || stdev[1] <= 1 || stdev[2] <= 1)
86         return max_sqr_errs[0];
87
88     err /= MAGN_DS_SIZE;
89     return err;
90 }
91
92 /* Given an real symmetric 3x3 matrix A, compute the eigenvalues */
93 static void compute_eigenvalues(double mat[3][3], double* eig1, double* eig2, double* eig3)
94 {
95     double p = mat[0][1] * mat[0][1] + mat[0][2] * mat[0][2] + mat[1][2] * mat[1][2];
96
97     if (p < EPSILON) {
98         *eig1 = mat[0][0];
99         *eig2 = mat[1][1];
100         *eig3 = mat[2][2];
101         return;
102     }
103
104     double q = (mat[0][0] + mat[1][1] + mat[2][2]) / 3;
105     double temp1 = mat[0][0] - q;
106     double temp2 = mat[1][1] - q;
107     double temp3 = mat[2][2] - q;
108
109     p = temp1 * temp1 + temp2 * temp2 + temp3 * temp3 + 2 * p;
110     p = sqrt(p / 6);
111
112     double mat2[3][3];
113     assign(3, 3, mat, mat2);
114     mat2[0][0] -= q;
115     mat2[1][1] -= q;
116     mat2[2][2] -= q;
117     multiply_scalar_inplace(3, 3, mat2, 1/p);
118
119     double r = (mat2[0][0] * mat2[1][1] * mat2[2][2] + mat2[0][1] * mat2[1][2] * mat2[2][0]
120         + mat2[0][2] * mat2[1][0] * mat2[2][1] - mat2[0][2] * mat2[1][1] * mat2[2][0]
121         - mat2[0][0] * mat2[1][2] * mat2[2][1] - mat2[0][1] * mat2[1][0] * mat2[2][2]) / 2;
122
123     double phi;
124     if (r <= -1.0)
125         phi = M_PI/3;
126     else if (r >= 1.0)
127         phi = 0;
128     else
129         phi = acos(r) / 3;
130
131     *eig3 = q + 2 * p * cos(phi);
132     *eig1 = q + 2 * p * cos(phi + 2 * M_PI / 3);
133     *eig2 = 3 * q - *eig1 - *eig3;
134 }
135
136 static void calc_evector(double mat[3][3], double eig, double vec[3][1])
137 {
138     double h[3][3];
139     double x_tmp[2][2];
140     assign(3, 3, mat, h);
141     h[0][0] -= eig;
142     h[1][1] -= eig;
143     h[2][2] -= eig;
144
145     double x[2][2];
146     x[0][0] = h[1][1];
147     x[0][1] = h[1][2];
148     x[1][0] = h[2][1];
149     x[1][1] = h[2][2];
150     invert(2, x, x_tmp);
151     assign(2, 2, x_tmp, x);
152
153     double temp1 = x[0][0] * (-h[1][0]) + x[0][1] * (-h[2][0]);
154     double temp2 = x[1][0] * (-h[1][0]) + x[1][1] * (-h[2][0]);
155     double norm = sqrt(1 + temp1 * temp1 + temp2 * temp2);
156
157     vec[0][0] = 1.0 / norm;
158     vec[1][0] = temp1 / norm;
159     vec[2][0] = temp2 / norm;
160 }
161
162 static int ellipsoid_fit (mat_input_t m, double offset[3][1], double w_invert[3][3], double* bfield)
163 {
164     int i;
165     double h[MAGN_DS_SIZE][9];
166     double w[MAGN_DS_SIZE][1];
167     double h_trans[9][MAGN_DS_SIZE];
168     double p_temp1[9][9];
169     double p_temp2[9][MAGN_DS_SIZE];
170     double temp1[3][3], temp[3][3];
171     double temp1_inv[3][3];
172     double temp2[3][1];
173     double result[9][9];
174     double p[9][1];
175     double a[3][3], sqrt_evals[3][3], evecs[3][3], evecs_trans[3][3];
176     double evec1[3][1], evec2[3][1], evec3[3][1];
177
178     for (i = 0; i < MAGN_DS_SIZE; i++) {
179         w[i][0] = m[i][0] * m[i][0];
180         h[i][0] = m[i][0];
181         h[i][1] = m[i][1];
182         h[i][2] = m[i][2];
183         h[i][3] = -1 * m[i][0] * m[i][1];
184         h[i][4] = -1 * m[i][0] * m[i][2];
185         h[i][5] = -1 * m[i][1] * m[i][2];
186         h[i][6] = -1 * m[i][1] * m[i][1];
187         h[i][7] = -1 * m[i][2] * m[i][2];
188         h[i][8] = 1;
189     }
190     transpose (MAGN_DS_SIZE, 9, h, h_trans);
191     multiply (9, MAGN_DS_SIZE, 9, h_trans, h, result);
192     invert (9, result, p_temp1);
193     multiply (9, 9, MAGN_DS_SIZE, p_temp1, h_trans, p_temp2);
194     multiply (9, MAGN_DS_SIZE, 1, p_temp2, w, p);
195
196     temp1[0][0] = 2;
197     temp1[0][1] = p[3][0];
198     temp1[0][2] = p[4][0];
199     temp1[1][0] = p[3][0];
200     temp1[1][1] = 2 * p[6][0];
201     temp1[1][2] = p[5][0];
202     temp1[2][0] = p[4][0];
203     temp1[2][1] = p[5][0];
204     temp1[2][2] = 2 * p[7][0];
205
206     temp2[0][0] = p[0][0];
207     temp2[1][0] = p[1][0];
208     temp2[2][0] = p[2][0];
209
210     invert(3, temp1, temp1_inv);
211     multiply(3, 3, 1, temp1_inv, temp2, offset);
212     double off_x = offset[0][0];
213     double off_y = offset[1][0];
214     double off_z = offset[2][0];
215
216
217     a[0][0] = 1.0 / (p[8][0] + off_x * off_x + p[6][0] * off_y * off_y
218             + p[7][0] * off_z * off_z + p[3][0] * off_x * off_y
219             + p[4][0] * off_x * off_z + p[5][0] * off_y * off_z);
220
221     a[0][1] = p[3][0] * a[0][0] / 2;
222     a[0][2] = p[4][0] * a[0][0] / 2;
223     a[1][2] = p[5][0] * a[0][0] / 2;
224     a[1][1] = p[6][0] * a[0][0];
225     a[2][2] = p[7][0] * a[0][0];
226     a[2][1] = a[1][2];
227     a[1][0] = a[0][1];
228     a[2][0] = a[0][2];
229
230     double eig1 = 0, eig2 = 0, eig3 = 0;
231     compute_eigenvalues(a, &eig1, &eig2, &eig3);
232
233     if (eig1 <=0 || eig2 <= 0 || eig3 <= 0)
234         return 0;
235
236     sqrt_evals[0][0] = sqrt(eig1);
237     sqrt_evals[1][0] = 0;
238     sqrt_evals[2][0] = 0;
239     sqrt_evals[0][1] = 0;
240     sqrt_evals[1][1] = sqrt(eig2);
241     sqrt_evals[2][1] = 0;
242     sqrt_evals[0][2] = 0;
243     sqrt_evals[1][2] = 0;
244     sqrt_evals[2][2] = sqrt(eig3);
245
246     calc_evector(a, eig1, evec1);
247     calc_evector(a, eig2, evec2);
248     calc_evector(a, eig3, evec3);
249
250     evecs[0][0] = evec1[0][0];
251     evecs[1][0] = evec1[1][0];
252     evecs[2][0] = evec1[2][0];
253     evecs[0][1] = evec2[0][0];
254     evecs[1][1] = evec2[1][0];
255     evecs[2][1] = evec2[2][0];
256     evecs[0][2] = evec3[0][0];
257     evecs[1][2] = evec3[1][0];
258     evecs[2][2] = evec3[2][0];
259
260     multiply (3, 3, 3, evecs, sqrt_evals, temp1);
261     transpose(3, 3, evecs, evecs_trans);
262     multiply (3, 3, 3, temp1, evecs_trans, temp);
263     transpose (3, 3, temp, w_invert);
264     *bfield = pow(sqrt(1/eig1) * sqrt(1/eig2) * sqrt(1/eig3), 1.0/3.0);
265
266     if (*bfield < 0)
267         return 0;
268
269     multiply_scalar_inplace(3, 3, w_invert, *bfield);
270
271     return 1;
272 }
273
274 static void compass_cal_init (FILE* data_file, struct sensor_info_t* info)
275 {
276
277 #ifdef DBG_RAW_DATA
278     if (raw_data) {
279         fclose(raw_data);
280         raw_data = NULL;
281     }
282
283     if (raw_data_selected) {
284         fclose(raw_data_selected);
285         raw_data_selected = NULL;
286     }
287
288     char path[64];
289     snprintf(path, 64, RAW_DATA_FULL_PATH, file_no);
290     raw_data = fopen(path,"w+");
291     snprintf(path, 64, RAW_DATA_SELECTED_PATH, file_no);
292     raw_data_selected = fopen(path,"w+");
293     file_no++;
294     raw_data_count = 0;
295 #endif
296
297     struct compass_cal_t* cal_data = (struct compass_cal_t*) info->cal_data;
298     int cal_steps = (info->max_cal_level && info->max_cal_level <= CAL_STEPS) ?
299         info->max_cal_level : CAL_STEPS;
300     if (cal_data == NULL)
301         return;
302
303     int data_count = 14;
304     reset_sample(cal_data);
305
306     if (!info->cal_level && data_file != NULL) {
307        int ret = fscanf(data_file, "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
308             &info->cal_level, &cal_data->offset[0][0], &cal_data->offset[1][0], &cal_data->offset[2][0],
309             &cal_data->w_invert[0][0], &cal_data->w_invert[0][1], &cal_data->w_invert[0][2],
310             &cal_data->w_invert[1][0], &cal_data->w_invert[1][1], &cal_data->w_invert[1][2],
311             &cal_data->w_invert[2][0], &cal_data->w_invert[2][1], &cal_data->w_invert[2][2],
312             &cal_data->bfield);
313
314         if (ret != data_count || info->cal_level >= cal_steps) {
315             info->cal_level = 0;
316         }
317     }
318
319
320     if (info->cal_level) {
321         ALOGV("CompassCalibration: load old data, caldata: %f %f %f %f %f %f %f %f %f %f %f %f %f",
322             cal_data->offset[0][0], cal_data->offset[1][0], cal_data->offset[2][0],
323             cal_data->w_invert[0][0], cal_data->w_invert[0][1], cal_data->w_invert[0][2], cal_data->w_invert[1][0],
324             cal_data->w_invert[1][1], cal_data->w_invert[1][2], cal_data->w_invert[2][0], cal_data->w_invert[2][1],
325             cal_data->w_invert[2][2], cal_data->bfield);
326
327     } else {
328         cal_data->offset[0][0] = 0;
329         cal_data->offset[1][0] = 0;
330         cal_data->offset[2][0] = 0;
331
332         cal_data->w_invert[0][0] = 1;
333         cal_data->w_invert[1][0] = 0;
334         cal_data->w_invert[2][0] = 0;
335         cal_data->w_invert[0][1] = 0;
336         cal_data->w_invert[1][1] = 1;
337         cal_data->w_invert[2][1] = 0;
338         cal_data->w_invert[0][2] = 0;
339         cal_data->w_invert[1][2] = 0;
340         cal_data->w_invert[2][2] = 1;
341
342         cal_data->bfield = 0;
343     }
344
345 }
346
347 static void compass_store_result(FILE* data_file, struct sensor_info_t* info)
348 {
349     struct compass_cal_t* cal_data = (struct compass_cal_t*) info->cal_data;
350
351     if (data_file == NULL || cal_data == NULL)
352         return;
353
354     int ret = fprintf(data_file, "%d %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
355         info->cal_level, cal_data->offset[0][0], cal_data->offset[1][0], cal_data->offset[2][0],
356         cal_data->w_invert[0][0], cal_data->w_invert[0][1], cal_data->w_invert[0][2],
357         cal_data->w_invert[1][0], cal_data->w_invert[1][1], cal_data->w_invert[1][2],
358         cal_data->w_invert[2][0], cal_data->w_invert[2][1], cal_data->w_invert[2][2],
359         cal_data->bfield);
360
361     if (ret < 0)
362         ALOGE ("compass calibration - store data failed!");
363 }
364
365 static int compass_collect (struct sensors_event_t* event, struct sensor_info_t* info)
366 {
367     float data[3] = {event->magnetic.x, event->magnetic.y, event->magnetic.z};
368     unsigned int index,j;
369     unsigned int lookback_count;
370     float min_diff;
371
372     struct compass_cal_t* cal_data = (struct compass_cal_t*) info->cal_data;
373
374     if (cal_data == NULL)
375         return -1;
376
377     /* Discard the point if not valid */
378     if (data[0] == 0 || data[1] == 0 || data[2] == 0)
379         return -1;
380
381 #ifdef DBG_RAW_DATA
382     if (raw_data && raw_data_count < MAX_RAW_DATA_COUNT) {
383         fprintf(raw_data, "%f %f %f\n", (double)data[0], (double)data[1],
384                  (double)data[2]);
385         raw_data_count++;
386     }
387
388     if (raw_data && raw_data_count >= MAX_RAW_DATA_COUNT) {
389         fclose(raw_data);
390         raw_data = NULL;
391     }
392 #endif
393
394     lookback_count = lookback_counts[info->cal_level];
395     min_diff = min_diffs[info->cal_level];
396
397     /*
398      * For the current point to be accepted, each x/y/z value must be different
399      * enough to the last several collected points.
400      */
401     if (cal_data->sample_count > 0 && cal_data->sample_count < MAGN_DS_SIZE) {
402         unsigned int lookback = lookback_count < cal_data->sample_count ? lookback_count :
403                         cal_data->sample_count;
404         for (index = 0; index < lookback; index++){
405             for (j = 0; j < 3; j++) {
406                 if (fabsf(data[j] - cal_data->sample[cal_data->sample_count-1-index][j]) < min_diff) {
407                     ALOGV("CompassCalibration:point reject: [%f,%f,%f], selected_count=%d",
408                        data[0], data[1], data[2], cal_data->sample_count);
409                        return 0;
410                 }
411             }
412         }
413     }
414
415     if (cal_data->sample_count < MAGN_DS_SIZE) {
416         memcpy(cal_data->sample[cal_data->sample_count], data, sizeof(float) * 3);
417         cal_data->sample_count++;
418         cal_data->average[0] += data[0];
419         cal_data->average[1] += data[1];
420         cal_data->average[2] += data[2];
421         ALOGV("CompassCalibration:point collected [%f,%f,%f], selected_count=%d",
422             (double)data[0], (double)data[1], (double)data[2], cal_data->sample_count);
423 #ifdef DBG_RAW_DATA
424         if (raw_data_selected) {
425             fprintf(raw_data_selected, "%f %f %f\n", (double)data[0], (double)data[1], (double)data[2]);
426         }
427 #endif
428     }
429     return 1;
430 }
431
432 static void scale_event (struct sensors_event_t* event)
433 {
434     float sqr_norm = 0;
435     float sanity_norm = 0;
436     float scale = 1;
437
438     sqr_norm = (event->magnetic.x * event->magnetic.x +
439                 event->magnetic.y * event->magnetic.y +
440                 event->magnetic.z * event->magnetic.z);
441
442     if (sqr_norm < MAGNETIC_LOW)
443         sanity_norm = MAGNETIC_LOW;
444
445     if (sanity_norm && sqr_norm) {
446         scale = sanity_norm / sqr_norm;
447         scale = sqrt(scale);
448         event->magnetic.x = event->magnetic.x * scale;
449         event->magnetic.y = event->magnetic.y * scale;
450         event->magnetic.z = event->magnetic.z * scale;
451
452     }
453 }
454
455 static void compass_compute_cal (struct sensors_event_t* event, struct sensor_info_t* info)
456 {
457     struct compass_cal_t* cal_data = (struct compass_cal_t*) info->cal_data;
458     double result[3][1], raw[3][1], diff[3][1];
459
460     if (!info->cal_level || cal_data == NULL)
461         return;
462
463     raw[0][0] = event->magnetic.x;
464     raw[1][0] = event->magnetic.y;
465     raw[2][0] = event->magnetic.z;
466
467     substract(3, 1, raw, cal_data->offset, diff);
468     multiply (3, 3, 1, cal_data->w_invert, diff, result);
469
470     event->magnetic.x = event->data[0] = result[0][0];
471     event->magnetic.y = event->data[1] = result[1][0];
472     event->magnetic.z = event->data[2] = result[2][0];
473
474     scale_event(event);
475 }
476
477
478 static int compass_ready (struct sensor_info_t* info)
479 {
480     mat_input_t mat;
481     int i;
482     float max_sqr_err;
483
484     struct compass_cal_t* cal_data = (struct compass_cal_t*) info->cal_data;
485     struct compass_cal_t new_cal_data;
486
487     /*
488      *  Some sensors take unrealistically long to calibrate at higher levels.
489      *  We'll use a max_cal_level if we have such a property setup, or go with
490      *  the default settings if not.
491      */
492     int cal_steps = (info->max_cal_level && info->max_cal_level <= CAL_STEPS) ?
493         info->max_cal_level : CAL_STEPS;
494
495     if (cal_data->sample_count < MAGN_DS_SIZE)
496         return info->cal_level;
497
498     max_sqr_err = max_sqr_errs[info->cal_level];
499
500     /* Enough points have been collected, do the ellipsoid calibration */
501
502     /* Compute average per axis */
503     cal_data->average[0] /= MAGN_DS_SIZE;
504     cal_data->average[1] /= MAGN_DS_SIZE;
505     cal_data->average[2] /= MAGN_DS_SIZE;
506
507     for (i = 0; i < MAGN_DS_SIZE; i++) {
508        mat[i][0] = cal_data->sample[i][0];
509        mat[i][1] = cal_data->sample[i][1];
510        mat[i][2] = cal_data->sample[i][2];
511     }
512
513     /* Check if result is good. The sample data must remain the same */
514     new_cal_data = *cal_data;
515
516     if (ellipsoid_fit(mat, new_cal_data.offset, new_cal_data.w_invert, &new_cal_data.bfield)) {
517         double new_err = calc_square_err (&new_cal_data);
518         ALOGI("new err is %f, max sqr err id %f", new_err,max_sqr_err);
519         if (new_err < max_sqr_err) {
520             double err = calc_square_err(cal_data);
521             if (new_err < err) {
522                 /* new cal data is better, so we switch to the new */
523                 memcpy(cal_data->offset, new_cal_data.offset, sizeof(cal_data->offset));
524                 memcpy(cal_data->w_invert, new_cal_data.w_invert, sizeof(cal_data->w_invert));
525                 cal_data->bfield = new_cal_data.bfield;
526                 if (info->cal_level < (cal_steps - 1))
527                     info->cal_level++;
528                 ALOGV("CompassCalibration: ready check success, caldata: %f %f %f %f %f %f %f %f %f %f %f %f %f, err %f",
529                     cal_data->offset[0][0], cal_data->offset[1][0], cal_data->offset[2][0], cal_data->w_invert[0][0],
530                     cal_data->w_invert[0][1], cal_data->w_invert[0][2], cal_data->w_invert[1][0], cal_data->w_invert[1][1],
531                     cal_data->w_invert[1][2], cal_data->w_invert[2][0], cal_data->w_invert[2][1], cal_data->w_invert[2][2],
532                     cal_data->bfield, new_err);
533             }
534         }
535     }
536     reset_sample(cal_data);
537     return info->cal_level;
538 }
539
540
541 void calibrate_compass (struct sensors_event_t* event, struct sensor_info_t* info)
542 {
543     int cal_level;
544
545     /* Calibration is continuous */
546     compass_collect (event, info);
547
548     cal_level = compass_ready(info);
549
550     switch (cal_level) {
551
552         case 0:
553             scale_event(event);
554             event->magnetic.status = SENSOR_STATUS_UNRELIABLE;
555             break;
556
557         case 1:
558             compass_compute_cal (event, info);
559             event->magnetic.status = SENSOR_STATUS_ACCURACY_LOW;
560             break;
561
562         case 2:
563             compass_compute_cal (event, info);
564             event->magnetic.status = SENSOR_STATUS_ACCURACY_MEDIUM;
565             break;
566
567         default:
568             compass_compute_cal (event, info);
569             event->magnetic.status = SENSOR_STATUS_ACCURACY_HIGH;
570             break;
571     }
572 }
573
574 void compass_read_data (struct sensor_info_t* info)
575 {
576     FILE* data_file = fopen (COMPASS_CALIBRATION_PATH, "r");
577
578     compass_cal_init(data_file, info);
579     if (data_file)
580         fclose(data_file);
581 }
582
583 void compass_store_data (struct sensor_info_t* info)
584 {
585     FILE* data_file = fopen (COMPASS_CALIBRATION_PATH, "w");
586
587     compass_store_result(data_file, info);
588     if (data_file)
589         fclose(data_file);
590
591 }