compositor - EWA filter was blurring too much by default, this caused the displace...
[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 /**
263  * \note \a sampler at the moment is either 'COM_PS_NEAREST' or not, other values won't matter.
264  */
265 void MemoryBuffer::readEWA(float result[4], float fx, float fy, float dx, float dy, PixelSampler sampler)
266 {
267         const int width = this->getWidth(), height = this->getHeight();
268         
269         // scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values,
270         // scaling by aspect ratio alone does the opposite, so try something in between instead...
271         const float ff2 = width, ff = sqrtf(ff2), q = height / ff;
272         const float Ux = dx * ff, Vx = dx * q, Uy = dy * ff, Vy = dy * q;
273         float A = Vx * Vx + Vy * Vy;
274         float B = -2.f * (Ux * Vx + Uy * Vy);
275         float C = Ux * Ux + Uy * Uy;
276         float F = A * C - B * B * 0.25f;
277         float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d;
278         int u, v, u1, u2, v1, v2;
279         // The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C,
280         // so the ellipse always covers at least some texels. But since the filter is now always larger,
281         // it also means that everywhere else it's also more blurry then ideally should be the case.
282         // So instead here the ellipse radii are modified instead whenever either is too low.
283         // Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off,
284         // and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on
285         // (minimum values: const float rmin = intpol ? 1.f : 0.5f;)
286         const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 0.765625f) / ff2;
287         imp2radangle(A, B, C, F, &a, &b, &th, &ecc);
288         if ((b2 = b * b) < rmin) {
289                 if ((a2 = a * a) < rmin) {
290                         B = 0.f;
291                         A = C = rmin;
292                         F = A * C;
293                 }
294                 else {
295                         b2 = rmin;
296                         radangle2imp(a2, b2, th, &A, &B, &C, &F);
297                 }
298         }
299
300         ue = ff * sqrtf(C);
301         ve = ff * sqrtf(A);
302         d = (float)(EWA_MAXIDX + 1) / (F * ff2);
303         A *= d;
304         B *= d;
305         C *= d;
306
307         U0 = fx;
308         V0 = fy;
309         u1 = (int)(floorf(U0 - ue));
310         u2 = (int)(ceilf(U0 + ue));
311         v1 = (int)(floorf(V0 - ve));
312         v2 = (int)(ceilf(V0 + ve));
313         U0 -= 0.5f;
314         V0 -= 0.5f;
315         DDQ = 2.f * A;
316         U = u1 - U0;
317         ac1 = A * (2.f * U + 1.f);
318         ac2 = A * U * U;
319         BU = B * U;
320
321         d = result[0] = result[1] = result[2] = result[3] = 0.f;
322         for (v = v1; v <= v2; ++v) {
323                 const float V = v - V0;
324                 float DQ = ac1 + B * V;
325                 float Q = (C * V + BU) * V + ac2;
326                 for (u = u1; u <= u2; ++u) {
327                         if (Q < (float)(EWA_MAXIDX + 1)) {
328                                 float tc[4];
329                                 const float wt = EWA_WTS[(Q < 0.f) ? 0 : (unsigned int)Q];
330                                 read(tc, clipuv(u, width), clipuv(v, height));
331                                 madd_v3_v3fl(result, tc, wt);
332                                 result[3] += result[3] ? tc[3] * wt : 0.f;
333                                 d += wt;
334                         }
335                         Q += DQ;
336                         DQ += DDQ;
337                 }
338         }
339         
340         // d should hopefully never be zero anymore
341         d = 1.f / d;
342         result[0] *= d;
343         result[1] *= d;
344         result[2] *= d;
345         // clipping can be ignored if alpha used, texr->ta already includes filtered edge
346         result[3] = result[3] ? result[3] * d : 1.f;
347 }