fix for bokeh blur using uninitialized memory - it would cause some tiles not to...
[blender.git] / source / blender / compositor / intern / COM_MemoryBuffer.cpp
1 /*
2  * Copyright 2011, Blender Foundation.
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  * Contributor: 
19  *              Jeroen Bakker 
20  *              Monique Dewanchand
21  */
22
23 #include "COM_MemoryBuffer.h"
24 #include "MEM_guardedalloc.h"
25 //#include "BKE_global.h"
26
27 unsigned int MemoryBuffer::determineBufferSize()
28 {
29         return getWidth() * getHeight();
30 }
31
32 int MemoryBuffer::getWidth() const
33 {
34         return this->m_rect.xmax - this->m_rect.xmin;
35 }
36 int MemoryBuffer::getHeight() const
37 {
38         return this->m_rect.ymax - this->m_rect.ymin;
39 }
40
41 MemoryBuffer::MemoryBuffer(MemoryProxy *memoryProxy, unsigned int chunkNumber, rcti *rect)
42 {
43         BLI_rcti_init(&this->m_rect, rect->xmin, rect->xmax, rect->ymin, rect->ymax);
44         this->m_memoryProxy = memoryProxy;
45         this->m_chunkNumber = chunkNumber;
46         this->m_buffer = (float *)MEM_mallocN(sizeof(float) * determineBufferSize() * COM_NUMBER_OF_CHANNELS, "COM_MemoryBuffer");
47         this->m_state = COM_MB_ALLOCATED;
48         this->m_datatype = COM_DT_COLOR;
49         this->m_chunkWidth = this->m_rect.xmax - this->m_rect.xmin;
50 }
51
52 MemoryBuffer::MemoryBuffer(MemoryProxy *memoryProxy, rcti *rect)
53 {
54         BLI_rcti_init(&this->m_rect, rect->xmin, rect->xmax, rect->ymin, rect->ymax);
55         this->m_memoryProxy = memoryProxy;
56         this->m_chunkNumber = -1;
57         this->m_buffer = (float *)MEM_mallocN(sizeof(float) * determineBufferSize() * COM_NUMBER_OF_CHANNELS, "COM_MemoryBuffer");
58         this->m_state = COM_MB_TEMPORARILY;
59         this->m_datatype = COM_DT_COLOR;
60         this->m_chunkWidth = this->m_rect.xmax - this->m_rect.xmin;
61 }
62 MemoryBuffer *MemoryBuffer::duplicate()
63 {
64         MemoryBuffer *result = new MemoryBuffer(this->m_memoryProxy, &this->m_rect);
65         memcpy(result->m_buffer, this->m_buffer, this->determineBufferSize() * COM_NUMBER_OF_CHANNELS * sizeof(float));
66         return result;
67 }
68 void MemoryBuffer::clear()
69 {
70         memset(this->m_buffer, 0, this->determineBufferSize() * COM_NUMBER_OF_CHANNELS * sizeof(float));
71 }
72
73 float *MemoryBuffer::convertToValueBuffer()
74 {
75         const unsigned int size = this->determineBufferSize();
76         unsigned int i;
77
78         float *result = new float[size];
79
80         const float *fp_src = this->m_buffer;
81         float       *fp_dst = result;
82
83         for (i = 0; i < size; i++, fp_dst++, fp_src += COM_NUMBER_OF_CHANNELS) {
84                 *fp_dst = *fp_src;
85         }
86
87         return result;
88 }
89
90 float MemoryBuffer::getMaximumValue()
91 {
92         float result = this->m_buffer[0];
93         const unsigned int size = this->determineBufferSize();
94         unsigned int i;
95
96         const float *fp_src = this->m_buffer;
97
98         for (i = 0; i < size; i++, fp_src += COM_NUMBER_OF_CHANNELS) {
99                 float value = *fp_src;
100                 if (value > result) {
101                         result = value;
102                 }
103         }
104
105         return result;
106 }
107
108 float MemoryBuffer::getMaximumValue(rcti *rect)
109 {
110         rcti rect_clamp;
111
112         /* first clamp the rect by the bounds or we get un-initialized values */
113         BLI_rcti_isect(rect, &this->m_rect, &rect_clamp);
114
115         if (!BLI_rcti_is_empty(&rect_clamp)) {
116                 MemoryBuffer *temp = new MemoryBuffer(NULL, &rect_clamp);
117                 temp->copyContentFrom(this);
118                 float result = temp->getMaximumValue();
119                 delete temp;
120                 return result;
121         }
122         else {
123                 BLI_assert(0);
124                 return 0.0f;
125         }
126 }
127
128 MemoryBuffer::~MemoryBuffer()
129 {
130         if (this->m_buffer) {
131                 MEM_freeN(this->m_buffer);
132                 this->m_buffer = NULL;
133         }
134 }
135
136 void MemoryBuffer::copyContentFrom(MemoryBuffer *otherBuffer)
137 {
138         if (!otherBuffer) {
139                 BLI_assert(0);
140                 return;
141         }
142         unsigned int otherY;
143         unsigned int minX = max(this->m_rect.xmin, otherBuffer->m_rect.xmin);
144         unsigned int maxX = min(this->m_rect.xmax, otherBuffer->m_rect.xmax);
145         unsigned int minY = max(this->m_rect.ymin, otherBuffer->m_rect.ymin);
146         unsigned int maxY = min(this->m_rect.ymax, otherBuffer->m_rect.ymax);
147         int offset;
148         int otherOffset;
149
150
151         for (otherY = minY; otherY < maxY; otherY++) {
152                 otherOffset = ((otherY - otherBuffer->m_rect.ymin) * otherBuffer->m_chunkWidth + minX - otherBuffer->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
153                 offset = ((otherY - this->m_rect.ymin) * this->m_chunkWidth + minX - this->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
154                 memcpy(&this->m_buffer[offset], &otherBuffer->m_buffer[otherOffset], (maxX - minX) * COM_NUMBER_OF_CHANNELS * sizeof(float));
155         }
156 }
157
158 void MemoryBuffer::writePixel(int x, int y, const float color[4])
159 {
160         if (x >= this->m_rect.xmin && x < this->m_rect.xmax &&
161             y >= this->m_rect.ymin && y < this->m_rect.ymax)
162         {
163                 const int offset = (this->m_chunkWidth * (y-this->m_rect.ymin) + x-this->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
164                 copy_v4_v4(&this->m_buffer[offset], color);
165         }
166 }
167
168 void MemoryBuffer::addPixel(int x, int y, const float color[4])
169 {
170         if (x >= this->m_rect.xmin && x < this->m_rect.xmax &&
171             y >= this->m_rect.ymin && y < this->m_rect.ymax)
172         {
173                 const int offset = (this->m_chunkWidth * (y-this->m_rect.ymin) + x-this->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
174                 add_v4_v4(&this->m_buffer[offset], color);
175         }
176 }
177
178
179 // table of (exp(ar) - exp(a)) / (1 - exp(a)) for r in range [0, 1] and a = -2
180 // used instead of actual gaussian, otherwise at high texture magnifications circular artifacts are visible
181 #define EWA_MAXIDX 255
182 static float EWA_WTS[EWA_MAXIDX + 1] = {
183         1.f, 0.990965f, 0.982f, 0.973105f, 0.96428f, 0.955524f, 0.946836f, 0.938216f, 0.929664f,
184         0.921178f, 0.912759f, 0.904405f, 0.896117f, 0.887893f, 0.879734f, 0.871638f, 0.863605f,
185         0.855636f, 0.847728f, 0.839883f, 0.832098f, 0.824375f, 0.816712f, 0.809108f, 0.801564f,
186         0.794079f, 0.786653f, 0.779284f, 0.771974f, 0.76472f, 0.757523f, 0.750382f, 0.743297f,
187         0.736267f, 0.729292f, 0.722372f, 0.715505f, 0.708693f, 0.701933f, 0.695227f, 0.688572f,
188         0.68197f, 0.67542f, 0.66892f, 0.662471f, 0.656073f, 0.649725f, 0.643426f, 0.637176f,
189         0.630976f, 0.624824f, 0.618719f, 0.612663f, 0.606654f, 0.600691f, 0.594776f, 0.588906f,
190         0.583083f, 0.577305f, 0.571572f, 0.565883f, 0.56024f, 0.55464f, 0.549084f, 0.543572f,
191         0.538102f, 0.532676f, 0.527291f, 0.521949f, 0.516649f, 0.511389f, 0.506171f, 0.500994f,
192         0.495857f, 0.490761f, 0.485704f, 0.480687f, 0.475709f, 0.470769f, 0.465869f, 0.461006f,
193         0.456182f, 0.451395f, 0.446646f, 0.441934f, 0.437258f, 0.432619f, 0.428017f, 0.42345f,
194         0.418919f, 0.414424f, 0.409963f, 0.405538f, 0.401147f, 0.39679f, 0.392467f, 0.388178f,
195         0.383923f, 0.379701f, 0.375511f, 0.371355f, 0.367231f, 0.363139f, 0.359079f, 0.355051f,
196         0.351055f, 0.347089f, 0.343155f, 0.339251f, 0.335378f, 0.331535f, 0.327722f, 0.323939f,
197         0.320186f, 0.316461f, 0.312766f, 0.3091f, 0.305462f, 0.301853f, 0.298272f, 0.294719f,
198         0.291194f, 0.287696f, 0.284226f, 0.280782f, 0.277366f, 0.273976f, 0.270613f, 0.267276f,
199         0.263965f, 0.26068f, 0.257421f, 0.254187f, 0.250979f, 0.247795f, 0.244636f, 0.241502f,
200         0.238393f, 0.235308f, 0.232246f, 0.229209f, 0.226196f, 0.223206f, 0.220239f, 0.217296f,
201         0.214375f, 0.211478f, 0.208603f, 0.20575f, 0.20292f, 0.200112f, 0.197326f, 0.194562f,
202         0.191819f, 0.189097f, 0.186397f, 0.183718f, 0.18106f, 0.178423f, 0.175806f, 0.17321f,
203         0.170634f, 0.168078f, 0.165542f, 0.163026f, 0.16053f, 0.158053f, 0.155595f, 0.153157f,
204         0.150738f, 0.148337f, 0.145955f, 0.143592f, 0.141248f, 0.138921f, 0.136613f, 0.134323f,
205         0.132051f, 0.129797f, 0.12756f, 0.125341f, 0.123139f, 0.120954f, 0.118786f, 0.116635f,
206         0.114501f, 0.112384f, 0.110283f, 0.108199f, 0.106131f, 0.104079f, 0.102043f, 0.100023f,
207         0.0980186f, 0.09603f, 0.094057f, 0.0920994f, 0.0901571f, 0.08823f, 0.0863179f, 0.0844208f,
208         0.0825384f, 0.0806708f, 0.0788178f, 0.0769792f, 0.0751551f, 0.0733451f, 0.0715493f, 0.0697676f,
209         0.0679997f, 0.0662457f, 0.0645054f, 0.0627786f, 0.0610654f, 0.0593655f, 0.0576789f, 0.0560055f,
210         0.0543452f, 0.0526979f, 0.0510634f, 0.0494416f, 0.0478326f, 0.0462361f, 0.0446521f, 0.0430805f,
211         0.0415211f, 0.039974f, 0.0384389f, 0.0369158f, 0.0354046f, 0.0339052f, 0.0324175f, 0.0309415f,
212         0.029477f, 0.0280239f, 0.0265822f, 0.0251517f, 0.0237324f, 0.0223242f, 0.020927f, 0.0195408f,
213         0.0181653f, 0.0168006f, 0.0154466f, 0.0141031f, 0.0127701f, 0.0114476f, 0.0101354f, 0.00883339f,
214         0.00754159f, 0.00625989f, 0.00498819f, 0.00372644f, 0.00247454f, 0.00123242f, 0.f
215 };
216
217 static void radangle2imp(float a2, float b2, float th, float *A, float *B, float *C, float *F)
218 {
219         float ct2 = cosf(th);
220         const float st2 = 1.f - ct2 * ct2;    // <- sin(th)^2
221         ct2 *= ct2;
222         *A = a2 * st2 + b2 * ct2;
223         *B = (b2 - a2) * sinf(2.f * th);
224         *C = a2 * ct2 + b2 * st2;
225         *F = a2 * b2;
226 }
227
228 // all tests here are done to make sure possible overflows are hopefully minimized
229 static void imp2radangle(float A, float B, float C, float F, float *a, float *b, float *th, float *ecc)
230 {
231         if (F <= 1e-5f) {   // use arbitrary major radius, zero minor, infinite eccentricity
232                 *a = sqrtf(A > C ? A : C);
233                 *b = 0.f;
234                 *ecc = 1e10f;
235                 *th = 0.5f * (atan2f(B, A - C) + (float)M_PI);
236         }
237         else {
238                 const float AmC = A - C, ApC = A + C, F2 = F * 2.f;
239                 const float r = sqrtf(AmC * AmC + B * B);
240                 float d = ApC - r;
241                 *a = (d <= 0.f) ? sqrtf(A > C ? A : C) : sqrtf(F2 / d);
242                 d = ApC + r;
243                 if (d <= 0.f) {
244                         *b = 0.f;
245                         *ecc = 1e10f;
246                 }
247                 else {
248                         *b = sqrtf(F2 / d);
249                         *ecc = *a / *b;
250                 }
251                 // incr theta by 0.5*pi (angle of major axis)
252                 *th = 0.5f * (atan2f(B, AmC) + (float)M_PI);
253         }
254 }
255
256 float clipuv(float x, float limit)
257 {
258         x = (x < 0) ? 0 : ((x >= limit) ? (limit - 1) : x);
259         return x;
260 }
261
262 void MemoryBuffer::readEWA(float result[4], float fx, float fy, float dx, float dy)
263 {
264         const int width = this->getWidth(), height = this->getHeight();
265         
266         // scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values,
267         // scaling by aspect ratio alone does the opposite, so try something in between instead...
268         const float ff2 = width, ff = sqrtf(ff2), q = height / ff;
269         const float Ux = dx * ff, Vx = dx * q, Uy = dy * ff, Vy = dy * q;
270         float A = Vx * Vx + Vy * Vy;
271         float B = -2.f * (Ux * Vx + Uy * Vy);
272         float C = Ux * Ux + Uy * Uy;
273         float F = A * C - B * B * 0.25f;
274         float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d;
275         int u, v, u1, u2, v1, v2;
276         // The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C,
277         // so the ellipse always covers at least some texels. But since the filter is now always larger,
278         // it also means that everywhere else it's also more blurry then ideally should be the case.
279         // So instead here the ellipse radii are modified instead whenever either is too low.
280         // Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off,
281         // and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on
282         // (minimum values: const float rmin = intpol ? 1.f : 0.5f;)
283         const float rmin = 1.5625f / ff2;
284         imp2radangle(A, B, C, F, &a, &b, &th, &ecc);
285         if ((b2 = b * b) < rmin) {
286                 if ((a2 = a * a) < rmin) {
287                         B = 0.f;
288                         A = C = rmin;
289                         F = A * C;
290                 }
291                 else {
292                         b2 = rmin;
293                         radangle2imp(a2, b2, th, &A, &B, &C, &F);
294                 }
295         }
296
297         ue = ff * sqrtf(C);
298         ve = ff * sqrtf(A);
299         d = (float)(EWA_MAXIDX + 1) / (F * ff2);
300         A *= d;
301         B *= d;
302         C *= d;
303
304         U0 = fx;
305         V0 = fy;
306         u1 = (int)(floorf(U0 - ue));
307         u2 = (int)(ceilf(U0 + ue));
308         v1 = (int)(floorf(V0 - ve));
309         v2 = (int)(ceilf(V0 + ve));
310         U0 -= 0.5f;
311         V0 -= 0.5f;
312         DDQ = 2.f * A;
313         U = u1 - U0;
314         ac1 = A * (2.f * U + 1.f);
315         ac2 = A * U * U;
316         BU = B * U;
317
318         d = result[0] = result[1] = result[2] = result[3] = 0.f;
319         for (v = v1; v <= v2; ++v) {
320                 const float V = v - V0;
321                 float DQ = ac1 + B * V;
322                 float Q = (C * V + BU) * V + ac2;
323                 for (u = u1; u <= u2; ++u) {
324                         if (Q < (float)(EWA_MAXIDX + 1)) {
325                                 float tc[4];
326                                 const float wt = EWA_WTS[(Q < 0.f) ? 0 : (unsigned int)Q];
327                                 read(tc, clipuv(u, width), clipuv(v, height));
328                                 madd_v3_v3fl(result, tc, wt);
329                                 result[3] += result[3] ? tc[3] * wt : 0.f;
330                                 d += wt;
331                         }
332                         Q += DQ;
333                         DQ += DDQ;
334                 }
335         }
336         
337         // d should hopefully never be zero anymore
338         d = 1.f / d;
339         result[0] *= d;
340         result[1] *= d;
341         result[2] *= d;
342         // clipping can be ignored if alpha used, texr->ta already includes filtered edge
343         result[3] = result[3] ? result[3] * d : 1.f;
344 }