Merging r58196 through r58265 from trunk into soc-2013-depsgraph_mt
[blender.git] / source / blender / blenlib / intern / math_interp.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) 2012 by Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Sergey Sharybin
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  *
27  */
28
29 #include <math.h>
30
31 #include "BLI_math.h"
32
33 /**************************************************************************
34  *                            INTERPOLATIONS
35  *
36  * Reference and docs:
37  * http://wiki.blender.org/index.php/User:Damiles#Interpolations_Algorithms
38  ***************************************************************************/
39
40 /* BICUBIC Interpolation functions
41  *  More info: http://wiki.blender.org/index.php/User:Damiles#Bicubic_pixel_interpolation
42  * function assumes out to be zero'ed, only does RGBA */
43
44 static float P(float k)
45 {
46         float p1, p2, p3, p4;
47         p1 = max_ff(k + 2.0f, 0.0f);
48         p2 = max_ff(k + 1.0f, 0.0f);
49         p3 = max_ff(k, 0.0f);
50         p4 = max_ff(k - 1.0f, 0.0f);
51         return (float)(1.0f / 6.0f) * (p1 * p1 * p1 - 4.0f * p2 * p2 * p2 + 6.0f * p3 * p3 * p3 - 4.0f * p4 * p4 * p4);
52 }
53
54
55 #if 0
56 /* older, slower function, works the same as above */
57 static float P(float k)
58 {
59         return (float)(1.0f / 6.0f) * (pow(MAX2(k + 2.0f, 0), 3.0f) - 4.0f * pow(MAX2(k + 1.0f, 0), 3.0f) + 6.0f * pow(MAX2(k, 0), 3.0f) - 4.0f * pow(MAX2(k - 1.0f, 0), 3.0f));
60 }
61 #endif
62
63 static void vector_from_float(const float *data, float vector[4], int components)
64 {
65         if (components == 1) {
66                 vector[0] = data[0];
67         }
68         else if (components == 3) {
69                 copy_v3_v3(vector, data);
70         }
71         else {
72                 copy_v4_v4(vector, data);
73         }
74 }
75
76 static void vector_from_byte(const unsigned char *data, float vector[4], int components)
77 {
78         if (components == 1) {
79                 vector[0] = data[0];
80         }
81         else if (components == 3) {
82                 vector[0] = data[0];
83                 vector[1] = data[1];
84                 vector[2] = data[2];
85         }
86         else {
87                 vector[0] = data[0];
88                 vector[1] = data[1];
89                 vector[2] = data[2];
90                 vector[3] = data[3];
91         }
92 }
93
94 /* BICUBIC INTERPOLATION */
95 BLI_INLINE void bicubic_interpolation(const unsigned char *byte_buffer, const float *float_buffer,
96                                       unsigned char *byte_output, float *float_output, int width, int height,
97                                       int components, float u, float v)
98 {
99         int i, j, n, m, x1, y1;
100         float a, b, w, wx, wy[4], out[4];
101
102         /* sample area entirely outside image? */
103         if (ceil(u) < 0 || floor(u) > width - 1 || ceil(v) < 0 || floor(v) > height - 1) {
104                 return;
105         }
106
107         i = (int)floor(u);
108         j = (int)floor(v);
109         a = u - i;
110         b = v - j;
111
112         zero_v4(out);
113
114 /* Optimized and not so easy to read */
115
116         /* avoid calling multiple times */
117         wy[0] = P(b - (-1));
118         wy[1] = P(b -  0);
119         wy[2] = P(b -  1);
120         wy[3] = P(b -  2);
121
122         for (n = -1; n <= 2; n++) {
123                 x1 = i + n;
124                 CLAMP(x1, 0, width - 1);
125                 wx = P(n - a);
126                 for (m = -1; m <= 2; m++) {
127                         float data[4];
128
129                         y1 = j + m;
130                         CLAMP(y1, 0, height - 1);
131                         /* normally we could do this */
132                         /* w = P(n-a) * P(b-m); */
133                         /* except that would call P() 16 times per pixel therefor pow() 64 times, better precalc these */
134                         w = wx * wy[m + 1];
135
136                         if (float_output) {
137                                 const float *float_data = float_buffer + width * y1 * components + components * x1;
138
139                                 vector_from_float(float_data, data, components);
140                         }
141                         else {
142                                 const unsigned char *byte_data = byte_buffer + width * y1 * components + components * x1;
143
144                                 vector_from_byte(byte_data, data, components);
145                         }
146
147                         if (components == 1) {
148                                 out[0] += data[0] * w;
149                         }
150                         else if (components == 3) {
151                                 out[0] += data[0] * w;
152                                 out[1] += data[1] * w;
153                                 out[2] += data[2] * w;
154                         }
155                         else {
156                                 out[0] += data[0] * w;
157                                 out[1] += data[1] * w;
158                                 out[2] += data[2] * w;
159                                 out[3] += data[3] * w;
160                         }
161                 }
162         }
163
164 /* Done with optimized part */
165
166 #if 0
167         /* older, slower function, works the same as above */
168         for (n = -1; n <= 2; n++) {
169                 for (m = -1; m <= 2; m++) {
170                         x1 = i + n;
171                         y1 = j + m;
172                         if (x1 > 0 && x1 < width && y1 > 0 && y1 < height) {
173                                 float data[4];
174
175                                 if (float_output) {
176                                         const float *float_data = float_buffer + width * y1 * components + components * x1;
177
178                                         vector_from_float(float_data, data, components);
179                                 }
180                                 else {
181                                         const unsigned char *byte_data = byte_buffer + width * y1 * components + components * x1;
182
183                                         vector_from_byte(byte_data, data, components);
184                                 }
185
186                                 if (components == 1) {
187                                         out[0] += data[0] * P(n - a) * P(b - m);
188                                 }
189                                 else if (components == 3) {
190                                         out[0] += data[0] * P(n - a) * P(b - m);
191                                         out[1] += data[1] * P(n - a) * P(b - m);
192                                         out[2] += data[2] * P(n - a) * P(b - m);
193                                 }
194                                 else {
195                                         out[0] += data[0] * P(n - a) * P(b - m);
196                                         out[1] += data[1] * P(n - a) * P(b - m);
197                                         out[2] += data[2] * P(n - a) * P(b - m);
198                                         out[3] += data[3] * P(n - a) * P(b - m);
199                                 }
200                         }
201                 }
202         }
203 #endif
204
205         if (float_output) {
206                 if (components == 1) {
207                         float_output[0] = out[0];
208                 }
209                 else if (components == 3) {
210                         copy_v3_v3(float_output, out);
211                 }
212                 else {
213                         copy_v4_v4(float_output, out);
214                 }
215         }
216         else {
217                 if (components == 1) {
218                         byte_output[0] = out[0] + 0.5f;
219                 }
220                 else if (components == 3) {
221                         byte_output[0] = out[0] + 0.5f;
222                         byte_output[1] = out[1] + 0.5f;
223                         byte_output[2] = out[2] + 0.5f;
224                 }
225                 else {
226                         byte_output[0] = out[0] + 0.5f;
227                         byte_output[1] = out[1] + 0.5f;
228                         byte_output[2] = out[2] + 0.5f;
229                         byte_output[3] = out[3] + 0.5f;
230                 }
231         }
232 }
233
234 void BLI_bicubic_interpolation_fl(const float *buffer, float *output, int width, int height,
235                                   int components, float u, float v)
236 {
237         bicubic_interpolation(NULL, buffer, NULL, output, width, height, components, u, v);
238 }
239
240 void BLI_bicubic_interpolation_char(const unsigned char *buffer, unsigned char *output, int width, int height,
241                                     int components, float u, float v)
242 {
243         bicubic_interpolation(buffer, NULL, output, NULL, width, height, components, u, v);
244 }
245
246 /* BILINEAR INTERPOLATION */
247 BLI_INLINE void bilinear_interpolation(const unsigned char *byte_buffer, const float *float_buffer,
248                                        unsigned char *byte_output, float *float_output, int width, int height,
249                                        int components, float u, float v)
250 {
251         float a, b;
252         float a_b, ma_b, a_mb, ma_mb;
253         int y1, y2, x1, x2;
254
255         /* ImBuf in must have a valid rect or rect_float, assume this is already checked */
256
257         x1 = (int)floor(u);
258         x2 = (int)ceil(u);
259         y1 = (int)floor(v);
260         y2 = (int)ceil(v);
261
262         /* sample area entirely outside image? */
263         if (x2 < 0 || x1 > width - 1 || y2 < 0 || y1 > height - 1) {
264                 return;
265         }
266
267         if (float_output) {
268                 const float *row1, *row2, *row3, *row4;
269                 float empty[4] = {0.0f, 0.0f, 0.0f, 0.0f};
270
271                 /* sample including outside of edges of image */
272                 if (x1 < 0 || y1 < 0) row1 = empty;
273                 else row1 = float_buffer + width * y1 * components + components * x1;
274
275                 if (x1 < 0 || y2 > height - 1) row2 = empty;
276                 else row2 = float_buffer + width * y2 * components + components * x1;
277
278                 if (x2 > width - 1 || y1 < 0) row3 = empty;
279                 else row3 = float_buffer + width * y1 * components + components * x2;
280
281                 if (x2 > width - 1 || y2 > height - 1) row4 = empty;
282                 else row4 = float_buffer + width * y2 * components + components * x2;
283
284                 a = u - floorf(u);
285                 b = v - floorf(v);
286                 a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
287
288                 if (components == 1) {
289                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
290                 }
291                 else if (components == 3) {
292                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
293                         float_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
294                         float_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
295                 }
296                 else {
297                         float_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
298                         float_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
299                         float_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
300                         float_output[3] = ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3];
301                 }
302         }
303         else {
304                 const unsigned char *row1, *row2, *row3, *row4;
305                 unsigned char empty[4] = {0, 0, 0, 0};
306
307                 /* sample including outside of edges of image */
308                 if (x1 < 0 || y1 < 0) row1 = empty;
309                 else row1 = byte_buffer + width * y1 * components + components * x1;
310
311                 if (x1 < 0 || y2 > height - 1) row2 = empty;
312                 else row2 = byte_buffer + width * y2 * components + components * x1;
313
314                 if (x2 > width - 1 || y1 < 0) row3 = empty;
315                 else row3 = byte_buffer + width * y1 * components + components * x2;
316
317                 if (x2 > width - 1 || y2 > height - 1) row4 = empty;
318                 else row4 = byte_buffer + width * y2 * components + components * x2;
319
320                 a = u - floorf(u);
321                 b = v - floorf(v);
322                 a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
323
324                 if (components == 1) {
325                         byte_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f;
326                 }
327                 else if (components == 3) {
328                         byte_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f;
329                         byte_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1] + 0.5f;
330                         byte_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2] + 0.5f;
331                 }
332                 else {
333                         byte_output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0] + 0.5f;
334                         byte_output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1] + 0.5f;
335                         byte_output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2] + 0.5f;
336                         byte_output[3] = ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3] + 0.5f;
337                 }
338         }
339 }
340
341 void BLI_bilinear_interpolation_fl(const float *buffer, float *output, int width, int height,
342                                    int components, float u, float v)
343 {
344         bilinear_interpolation(NULL, buffer, NULL, output, width, height, components, u, v);
345 }
346
347 void BLI_bilinear_interpolation_char(const unsigned char *buffer, unsigned char *output, int width, int height,
348                                      int components, float u, float v)
349 {
350         bilinear_interpolation(buffer, NULL, output, NULL, width, height, components, u, v);
351 }