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