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