BLI_math_rotation: properly name the quaternion power function.
[blender.git] / source / blender / blenlib / intern / math_statistics.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2015 by Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  * */
25
26 /** \file blender/blenlib/intern/math_statistics.c
27  *  \ingroup bli
28  */
29
30 #include "MEM_guardedalloc.h"
31
32 #include "BLI_math.h"
33 #include "BLI_task.h"
34 #include "BLI_utildefines.h"
35
36 #include "BLI_strict_flags.h"
37
38 /********************************** Covariance Matrices *********************************/
39
40 typedef struct CovarianceData {
41         const float *cos_vn;
42         const float *center;
43         float *r_covmat;
44         float covfac;
45         int n;
46         int nbr_cos_vn;
47 } CovarianceData;
48
49 static void covariance_m_vn_ex_task_cb(
50         void *__restrict userdata,
51         const int a,
52         const ParallelRangeTLS *__restrict UNUSED(tls))
53 {
54         CovarianceData *data = userdata;
55         const float *cos_vn = data->cos_vn;
56         const float *center = data->center;
57         float *r_covmat = data->r_covmat;
58         const int n = data->n;
59         const int nbr_cos_vn = data->nbr_cos_vn;
60
61         int k;
62
63         /* Covariance matrices are always symetrical, so we can compute only one half of it,
64          * and mirror it to the other half (at the end of the func).
65          *
66          * This allows using a flat loop of n*n with same results as imbricated one over half the matrix:
67          *
68          *     for (i = 0; i < n; i++) {
69          *         for (j = i; j < n; j++) {
70          *             ...
71          *         }
72          *      }
73          */
74         const int i = a / n;
75         const int j = a % n;
76         if (j < i)
77                 return;
78
79         if (center) {
80                 for (k = 0; k < nbr_cos_vn; k++) {
81                         r_covmat[a] += (cos_vn[k * n + i] - center[i]) * (cos_vn[k * n + j] - center[j]);
82                 }
83         }
84         else {
85                 for (k = 0; k < nbr_cos_vn; k++) {
86                         r_covmat[a] += cos_vn[k * n + i] * cos_vn[k * n + j];
87                 }
88         }
89         r_covmat[a] *= data->covfac;
90         if (j != i) {
91                 /* Mirror result to other half... */
92                 r_covmat[j * n + i] = r_covmat[a];
93         }
94 }
95
96 /**
97  * \brief Compute the covariance matrix of given set of nD coordinates.
98  *
99  * \param n: the dimension of the vectors (and hence, of the covariance matrix to compute).
100  * \param cos_vn: the nD points to compute covariance from.
101  * \param nbr_cos_vn: the number of nD coordinates in cos_vn.
102  * \param center: the center (or mean point) of cos_vn. If NULL, it is assumed cos_vn is already centered.
103  * \param use_sample_correction: whether to apply sample correction
104  *                              (i.e. get 'sample varince' instead of 'population variance').
105  * \return r_covmat the computed covariance matrix.
106  */
107 void BLI_covariance_m_vn_ex(
108         const int n, const float *cos_vn, const int nbr_cos_vn, const float *center, const bool use_sample_correction,
109         float *r_covmat)
110 {
111         /* Note about that division: see https://en.wikipedia.org/wiki/Bessel%27s_correction.
112          * In a nutshell, it must be 1 / (n - 1) for 'sample data', and 1 / n for 'population data'...
113          */
114         const float covfac = 1.0f / (float)(use_sample_correction ? nbr_cos_vn - 1 : nbr_cos_vn);
115
116         memset(r_covmat, 0, sizeof(*r_covmat) * (size_t)(n * n));
117
118         CovarianceData data = {
119                 .cos_vn = cos_vn, .center = center, .r_covmat = r_covmat,
120                 .covfac = covfac, .n = n, .nbr_cos_vn = nbr_cos_vn,
121         };
122
123         ParallelRangeSettings settings;
124         BLI_parallel_range_settings_defaults(&settings);
125         settings.use_threading = ((nbr_cos_vn * n * n) >= 10000);
126         BLI_task_parallel_range(
127                     0, n * n,
128                     &data,
129                     covariance_m_vn_ex_task_cb,
130                     &settings);
131 }
132
133 /**
134  * \brief Compute the covariance matrix of given set of 3D coordinates.
135  *
136  * \param cos_v3: the 3D points to compute covariance from.
137  * \param nbr_cos_v3: the number of 3D coordinates in cos_v3.
138  * \return r_covmat the computed covariance matrix.
139  * \return r_center the computed center (mean) of 3D points (may be NULL).
140  */
141 void BLI_covariance_m3_v3n(
142         const float (*cos_v3)[3], const int nbr_cos_v3, const bool use_sample_correction,
143         float r_covmat[3][3], float r_center[3])
144 {
145         float center[3];
146         const float mean_fac = 1.0f / (float)nbr_cos_v3;
147         int i;
148
149         zero_v3(center);
150         for (i = 0; i < nbr_cos_v3; i++) {
151                 /* Applying mean_fac here rather than once at the end reduce compute errors... */
152                 madd_v3_v3fl(center, cos_v3[i], mean_fac);
153         }
154
155         if (r_center) {
156                 copy_v3_v3(r_center, center);
157         }
158
159         BLI_covariance_m_vn_ex(3, (const float *)cos_v3, nbr_cos_v3, center, use_sample_correction, (float *)r_covmat);
160 }