Solved issue with distorted compositor results in some cases
[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 = MAX2(k + 2.0f, 0);
48         p2 = MAX2(k + 1.0f, 0);
49         p3 = MAX2(k, 0);
50         p4 = MAX2(k - 1.0f, 0);
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 /* BICUBIC INTERPOLATION */
64 void BLI_bicubic_interpolation(const float *buffer, float *output, int width, int height, int components, float u, float v)
65 {
66         int i, j, n, m, x1, y1;
67         float a, b, w, wx, wy[4], out[4];
68         const float *data;
69
70         /* sample area entirely outside image? */
71         if (ceil(u) < 0 || floor(u) > width - 1 || ceil(v) < 0 || floor(v) > height - 1) {
72                 return;
73         }
74
75         i = (int)floor(u);
76         j = (int)floor(v);
77         a = u - i;
78         b = v - j;
79
80         zero_v4(out);
81
82 /* Optimized and not so easy to read */
83
84         /* avoid calling multiple times */
85         wy[0] = P(b - (-1));
86         wy[1] = P(b -  0);
87         wy[2] = P(b -  1);
88         wy[3] = P(b -  2);
89
90         for (n = -1; n <= 2; n++) {
91                 x1 = i + n;
92                 CLAMP(x1, 0, width - 1);
93                 wx = P(n - a);
94                 for (m = -1; m <= 2; m++) {
95                         y1 = j + m;
96                         CLAMP(y1, 0, height - 1);
97                         /* normally we could do this */
98                         /* w = P(n-a) * P(b-m); */
99                         /* except that would call P() 16 times per pixel therefor pow() 64 times, better precalc these */
100                         w = wx * wy[m + 1];
101
102                         data = buffer + width * y1 * 4 + 4 * x1;
103
104                         if (components == 1) {
105                                 out[0] += data[0] * w;
106                         }
107                         else if (components == 2) {
108                                 out[0] += data[0] * w;
109                                 out[1] += data[1] * w;
110                         }
111                         else if (components == 3) {
112                                 out[0] += data[0] * w;
113                                 out[1] += data[1] * w;
114                                 out[2] += data[2] * w;
115                         }
116                         else {
117                                 out[0] += data[0] * w;
118                                 out[1] += data[1] * w;
119                                 out[2] += data[2] * w;
120                                 out[3] += data[3] * w;
121                         }
122                 }
123         }
124
125 /* Done with optimized part */
126
127 #if 0
128         /* older, slower function, works the same as above */
129         for (n = -1; n <= 2; n++) {
130                 for (m = -1; m <= 2; m++) {
131                         x1 = i + n;
132                         y1 = j + m;
133                         if (x1 > 0 && x1 < width && y1 > 0 && y1 < height) {
134                                 data = in->rect_float + width * y1 * 4 + 4 * x1;
135
136                                 if (components == 1) {
137                                         out[0] += data[0] * P(n - a) * P(b - m);
138                                 }
139                                 else if (components == 2) {
140                                         out[0] += data[0] * P(n - a) * P(b - m);
141                                         out[1] += data[1] * P(n - a) * P(b - m);
142                                 }
143                                 else if (components == 3) {
144                                         out[0] += data[0] * P(n - a) * P(b - m);
145                                         out[1] += data[1] * P(n - a) * P(b - m);
146                                         out[2] += data[2] * P(n - a) * P(b - m);
147                                 }
148                                 else {
149                                         out[0] += data[0] * P(n - a) * P(b - m);
150                                         out[1] += data[1] * P(n - a) * P(b - m);
151                                         out[2] += data[2] * P(n - a) * P(b - m);
152                                         out[3] += data[3] * P(n - a) * P(b - m);
153                                 }
154                         }
155                 }
156         }
157 #endif
158
159         if (components == 1) {
160                 output[0] = out[0];
161         }
162         else if (components == 2) {
163                 output[0] = out[0];
164                 output[1] = out[1];
165         }
166         else if (components == 3) {
167                 output[0] = out[0];
168                 output[1] = out[1];
169                 output[2] = out[2];
170         }
171         else {
172                 output[0] = out[0];
173                 output[1] = out[1];
174                 output[2] = out[2];
175                 output[3] = out[3];
176         }
177 }
178
179 /* BILINEAR INTERPOLATION */
180 void BLI_bilinear_interpolation(const float *buffer, float *output, int width, int height, int components, float u, float v)
181 {
182         const float *row1, *row2, *row3, *row4;
183         float a, b;
184         float a_b, ma_b, a_mb, ma_mb;
185         float empty[4] = {0.0f, 0.0f, 0.0f, 0.0f};
186         int y1, y2, x1, x2;
187
188         /* ImBuf in must have a valid rect or rect_float, assume this is already checked */
189
190         x1 = (int)floor(u);
191         x2 = (int)ceil(u);
192         y1 = (int)floor(v);
193         y2 = (int)ceil(v);
194
195         /* sample area entirely outside image? */
196         if (x2 < 0 || x1 > width - 1 || y2 < 0 || y1 > height - 1) {
197                 return;
198         }
199
200         /* sample including outside of edges of image */
201         if (x1 < 0 || y1 < 0) row1 = empty;
202         else row1 = buffer + width * y1 * 4 + 4 * x1;
203
204         if (x1 < 0 || y2 > height - 1) row2 = empty;
205         else row2 = buffer + width * y2 * 4 + 4 * x1;
206
207         if (x2 > width - 1 || y1 < 0) row3 = empty;
208         else row3 = buffer + width * y1 * 4 + 4 * x2;
209
210         if (x2 > width - 1 || y2 > height - 1) row4 = empty;
211         else row4 = buffer + width * y2 * 4 + 4 * x2;
212
213         a = u - floorf(u);
214         b = v - floorf(v);
215         a_b = a * b; ma_b = (1.0f - a) * b; a_mb = a * (1.0f - b); ma_mb = (1.0f - a) * (1.0f - b);
216
217         if (components == 1) {
218                 output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
219         }
220         else if (components == 2) {
221                 output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
222                 output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
223         }
224         else if (components == 3) {
225                 output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
226                 output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
227                 output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
228         }
229         else {
230                 output[0] = ma_mb * row1[0] + a_mb * row3[0] + ma_b * row2[0] + a_b * row4[0];
231                 output[1] = ma_mb * row1[1] + a_mb * row3[1] + ma_b * row2[1] + a_b * row4[1];
232                 output[2] = ma_mb * row1[2] + a_mb * row3[2] + ma_b * row2[2] + a_b * row4[2];
233                 output[3] = ma_mb * row1[3] + a_mb * row3[3] + ma_b * row2[3] + a_b * row4[3];
234         }
235 }