Curve Fitting: offset based fallback to calculate cubics
[blender.git] / extern / curve_fit_nd / intern / curve_fit_inline.h
1 /*
2  * Copyright (c) 2016, Blender Foundation.
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *     * Redistributions of source code must retain the above copyright
8  *       notice, this list of conditions and the following disclaimer.
9  *     * Redistributions in binary form must reproduce the above copyright
10  *       notice, this list of conditions and the following disclaimer in the
11  *       documentation and/or other materials provided with the distribution.
12  *     * Neither the name of the <organization> nor the
13  *       names of its contributors may be used to endorse or promote products
14  *       derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  */
27
28
29 /** \file curve_fit_inline.h
30  *  \ingroup curve_fit
31  */
32
33 /** \name Simple Vector Math Lib
34  * \{ */
35
36 #ifdef _MSC_VER
37 #  define MINLINE static __forceinline
38 #else
39 #  define MINLINE static inline
40 #endif
41
42 MINLINE double sq(const double d)
43 {
44         return d * d;
45 }
46
47 #ifndef _MSC_VER
48 MINLINE double min(const double a, const double b)
49 {
50         return b < a ? b : a;
51 }
52
53 MINLINE double max(const double a, const double b)
54 {
55         return a < b ? b : a;
56 }
57 #endif
58
59 MINLINE void zero_vn(
60         double v0[], const uint dims)
61 {
62         for (uint j = 0; j < dims; j++) {
63                 v0[j] = 0.0;
64         }
65 }
66
67 MINLINE void flip_vn_vnvn(
68         double v_out[], const double v0[], const double v1[], const uint dims)
69 {
70         for (uint j = 0; j < dims; j++) {
71                 v_out[j] = v0[j] + (v0[j] - v1[j]);
72         }
73 }
74
75 MINLINE void copy_vnvn(
76         double v0[], const double v1[], const uint dims)
77 {
78         for (uint j = 0; j < dims; j++) {
79                 v0[j] = v1[j];
80         }
81 }
82
83 MINLINE void copy_vnfl_vndb(
84         float v0[], const double v1[], const uint dims)
85 {
86         for (uint j = 0; j < dims; j++) {
87                 v0[j] = (float)v1[j];
88         }
89 }
90
91 MINLINE void copy_vndb_vnfl(
92         double v0[], const float v1[], const uint dims)
93 {
94         for (uint j = 0; j < dims; j++) {
95                 v0[j] = (double)v1[j];
96         }
97 }
98
99 MINLINE double dot_vnvn(
100         const double v0[], const double v1[], const uint dims)
101 {
102         double d = 0.0;
103         for (uint j = 0; j < dims; j++) {
104                 d += v0[j] * v1[j];
105         }
106         return d;
107 }
108
109 MINLINE void add_vn_vnvn(
110         double v_out[], const double v0[], const double v1[], const uint dims)
111 {
112         for (uint j = 0; j < dims; j++) {
113                 v_out[j] = v0[j] + v1[j];
114         }
115 }
116
117 MINLINE void sub_vn_vnvn(
118         double v_out[], const double v0[], const double v1[], const uint dims)
119 {
120         for (uint j = 0; j < dims; j++) {
121                 v_out[j] = v0[j] - v1[j];
122         }
123 }
124
125 MINLINE void iadd_vnvn(
126         double v0[], const double v1[], const uint dims)
127 {
128         for (uint j = 0; j < dims; j++) {
129                 v0[j] += v1[j];
130         }
131 }
132
133 MINLINE void isub_vnvn(
134         double v0[], const double v1[], const uint dims)
135 {
136         for (uint j = 0; j < dims; j++) {
137                 v0[j] -= v1[j];
138         }
139 }
140
141 MINLINE void madd_vn_vnvn_fl(
142         double v_out[],
143         const double v0[], const double v1[],
144         const double f, const uint dims)
145 {
146         for (uint j = 0; j < dims; j++) {
147                 v_out[j] = v0[j] + v1[j] * f;
148         }
149 }
150
151 MINLINE void msub_vn_vnvn_fl(
152         double v_out[],
153         const double v0[], const double v1[],
154         const double f, const uint dims)
155 {
156         for (uint j = 0; j < dims; j++) {
157                 v_out[j] = v0[j] - v1[j] * f;
158         }
159 }
160
161 MINLINE void miadd_vn_vn_fl(
162         double v_out[], const double v0[], double f, const uint dims)
163 {
164         for (uint j = 0; j < dims; j++) {
165                 v_out[j] += v0[j] * f;
166         }
167 }
168
169 #if 0
170 MINLINE void misub_vn_vn_fl(
171         double v_out[], const double v0[], double f, const uint dims)
172 {
173         for (uint j = 0; j < dims; j++) {
174                 v_out[j] -= v0[j] * f;
175         }
176 }
177 #endif
178
179 MINLINE void mul_vnvn_fl(
180         double v_out[],
181         const double v0[], const double f, const uint dims)
182 {
183         for (uint j = 0; j < dims; j++) {
184                 v_out[j] = v0[j] * f;
185         }
186 }
187
188 MINLINE void imul_vn_fl(double v0[], const double f, const uint dims)
189 {
190         for (uint j = 0; j < dims; j++) {
191                 v0[j] *= f;
192         }
193 }
194
195
196 MINLINE double len_squared_vnvn(
197         const double v0[], const double v1[], const uint dims)
198 {
199         double d = 0.0;
200         for (uint j = 0; j < dims; j++) {
201                 d += sq(v0[j] - v1[j]);
202         }
203         return d;
204 }
205
206 MINLINE double len_squared_vn(
207         const double v0[], const uint dims)
208 {
209         double d = 0.0;
210         for (uint j = 0; j < dims; j++) {
211                 d += sq(v0[j]);
212         }
213         return d;
214 }
215
216 MINLINE double len_vnvn(
217         const double v0[], const double v1[], const uint dims)
218 {
219         return sqrt(len_squared_vnvn(v0, v1, dims));
220 }
221
222 MINLINE double len_vn(
223         const double v0[], const uint dims)
224 {
225         return sqrt(len_squared_vn(v0, dims));
226 }
227
228 /* special case, save us negating a copy, then getting the length */
229 MINLINE double len_squared_negated_vnvn(
230         const double v0[], const double v1[], const uint dims)
231 {
232         double d = 0.0;
233         for (uint j = 0; j < dims; j++) {
234                 d += sq(v0[j] + v1[j]);
235         }
236         return d;
237 }
238
239 MINLINE double len_negated_vnvn(
240         const double v0[], const double v1[], const uint dims)
241 {
242         return sqrt(len_squared_negated_vnvn(v0, v1, dims));
243 }
244
245 MINLINE double normalize_vn(
246         double v0[], const uint dims)
247 {
248         double d = len_squared_vn(v0, dims);
249         if (d != 0.0 && ((d = sqrt(d)) != 0.0)) {
250                 imul_vn_fl(v0, 1.0 / d, dims);
251         }
252         return d;
253 }
254
255 /* v_out = (v0 - v1).normalized() */
256 MINLINE double normalize_vn_vnvn(
257         double v_out[],
258         const double v0[], const double v1[], const uint dims)
259 {
260         double d = 0.0;
261         for (uint j = 0; j < dims; j++) {
262                 double a = v0[j] - v1[j];
263                 d += sq(a);
264                 v_out[j] = a;
265         }
266         if (d != 0.0 && ((d = sqrt(d)) != 0.0)) {
267                 imul_vn_fl(v_out, 1.0 / d, dims);
268         }
269         return d;
270 }
271
272 MINLINE bool is_almost_zero_ex(double val, double eps)
273 {
274         return (-eps < val) && (val < eps);
275 }
276
277 MINLINE bool is_almost_zero(double val)
278 {
279         return is_almost_zero_ex(val, 1e-8);
280 }
281
282 MINLINE bool equals_vnvn(
283                 const double v0[], const double v1[], const uint dims)
284 {
285         for (uint j = 0; j < dims; j++) {
286                 if (v0[j] != v1[j]) {
287                         return false;
288                 }
289         }
290         return true;
291 }
292
293 #if 0
294 MINLINE void project_vn_vnvn(
295         double v_out[], const double p[], const double v_proj[], const uint dims)
296 {
297         const double mul = dot_vnvn(p, v_proj, dims) / dot_vnvn(v_proj, v_proj, dims);
298         mul_vnvn_fl(v_out, v_proj, mul, dims);
299 }
300 #endif
301
302 MINLINE void project_vn_vnvn_normalized(
303         double v_out[], const double p[], const double v_proj[], const uint dims)
304 {
305         const double mul = dot_vnvn(p, v_proj, dims);
306         mul_vnvn_fl(v_out, v_proj, mul, dims);
307 }
308
309 MINLINE void project_plane_vn_vnvn_normalized(
310         double v_out[], const double v[], const double v_plane[], const uint dims)
311 {
312         assert(v != v_out);
313         project_vn_vnvn_normalized(v_out, v, v_plane, dims);
314         sub_vn_vnvn(v_out, v, v_out, dims);
315 }
316
317 /** \} */