7f7429bf2e61e78d07705d75d04999fa45467006
[blender.git] / source / blender / compositor / operations / COM_GlareFogGlowOperation.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_GlareFogGlowOperation.h"
24 #include "MEM_guardedalloc.h"
25
26 /*
27  *  2D Fast Hartley Transform, used for convolution
28  */
29
30 typedef float fREAL;
31
32 // returns next highest power of 2 of x, as well it's log2 in L2
33 static unsigned int nextPow2(unsigned int x, unsigned int *L2)
34 {
35         unsigned int pw, x_notpow2 = x & (x - 1);
36         *L2 = 0;
37         while (x >>= 1) ++(*L2);
38         pw = 1 << (*L2);
39         if (x_notpow2) { (*L2)++;  pw <<= 1; }
40         return pw;
41 }
42
43 //------------------------------------------------------------------------------
44
45 // from FXT library by Joerg Arndt, faster in order bitreversal
46 // use: r = revbin_upd(r, h) where h = N>>1
47 static unsigned int revbin_upd(unsigned int r, unsigned int h)
48 {
49         while (!((r ^= h) & h)) h >>= 1;
50         return r;
51 }
52 //------------------------------------------------------------------------------
53 static void FHT(fREAL *data, unsigned int M, unsigned int inverse)
54 {
55         double tt, fc, dc, fs, ds, a = M_PI;
56         fREAL t1, t2;
57         int n2, bd, bl, istep, k, len = 1 << M, n = 1;
58
59         int i, j = 0;
60         unsigned int Nh = len >> 1;
61         for (i = 1; i < (len - 1); ++i) {
62                 j = revbin_upd(j, Nh);
63                 if (j > i) {
64                         t1 = data[i];
65                         data[i] = data[j];
66                         data[j] = t1;
67                 }
68         }
69
70         do {
71                 fREAL *data_n = &data[n];
72
73                 istep = n << 1;
74                 for (k = 0; k < len; k += istep) {
75                         t1 = data_n[k];
76                         data_n[k] = data[k] - t1;
77                         data[k] += t1;
78                 }
79
80                 n2 = n >> 1;
81                 if (n > 2) {
82                         fc = dc = cos(a);
83                         fs = ds = sqrt(1.0 - fc * fc); //sin(a);
84                         bd = n - 2;
85                         for (bl = 1; bl < n2; bl++) {
86                                 fREAL *data_nbd = &data_n[bd];
87                                 fREAL *data_bd = &data[bd];
88                                 for (k = bl; k < len; k += istep) {
89                                         t1 = fc * data_n[k] + fs * data_nbd[k];
90                                         t2 = fs * data_n[k] - fc * data_nbd[k];
91                                         data_n[k] = data[k] - t1;
92                                         data_nbd[k] = data_bd[k] - t2;
93                                         data[k] += t1;
94                                         data_bd[k] += t2;
95                                 }
96                                 tt = fc * dc - fs * ds;
97                                 fs = fs * dc + fc * ds;
98                                 fc = tt;
99                                 bd -= 2;
100                         }
101                 }
102
103                 if (n > 1) {
104                         for (k = n2; k < len; k += istep) {
105                                 t1 = data_n[k];
106                                 data_n[k] = data[k] - t1;
107                                 data[k] += t1;
108                         }
109                 }
110
111                 n = istep;
112                 a *= 0.5;
113         } while (n < len);
114
115         if (inverse) {
116                 fREAL sc = (fREAL)1 / (fREAL)len;
117                 for (k = 0; k < len; ++k)
118                         data[k] *= sc;
119         }
120 }
121 //------------------------------------------------------------------------------
122 /* 2D Fast Hartley Transform, Mx/My -> log2 of width/height,
123  * nzp -> the row where zero pad data starts,
124  * inverse -> see above */
125 static void FHT2D(fREAL *data, unsigned int Mx, unsigned int My,
126                   unsigned int nzp, unsigned int inverse)
127 {
128         unsigned int i, j, Nx, Ny, maxy;
129         fREAL t;
130
131         Nx = 1 << Mx;
132         Ny = 1 << My;
133
134         // rows (forward transform skips 0 pad data)
135         maxy = inverse ? Ny : nzp;
136         for (j = 0; j < maxy; ++j)
137                 FHT(&data[Nx * j], Mx, inverse);
138
139         // transpose data
140         if (Nx == Ny) {  // square
141                 for (j = 0; j < Ny; ++j)
142                         for (i = j + 1; i < Nx; ++i) {
143                                 unsigned int op = i + (j << Mx), np = j + (i << My);
144                                 t = data[op], data[op] = data[np], data[np] = t;
145                         }
146         }
147         else {  // rectangular
148                 unsigned int k, Nym = Ny - 1, stm = 1 << (Mx + My);
149                 for (i = 0; stm > 0; i++) {
150                         #define PRED(k) (((k & Nym) << Mx) + (k >> My))
151                         for (j = PRED(i); j > i; j = PRED(j)) ;
152                         if (j < i) continue;
153                         for (k = i, j = PRED(i); j != i; k = j, j = PRED(j), stm--) {
154                                 t = data[j], data[j] = data[k], data[k] = t;
155                         }
156                         #undef PRED
157                         stm--;
158                 }
159         }
160         // swap Mx/My & Nx/Ny
161         i = Nx, Nx = Ny, Ny = i;
162         i = Mx, Mx = My, My = i;
163
164         // now columns == transposed rows
165         for (j = 0; j < Ny; ++j)
166                 FHT(&data[Nx * j], Mx, inverse);
167
168         // finalize
169         for (j = 0; j <= (Ny >> 1); j++) {
170                 unsigned int jm = (Ny - j) & (Ny - 1);
171                 unsigned int ji = j << Mx;
172                 unsigned int jmi = jm << Mx;
173                 for (i = 0; i <= (Nx >> 1); i++) {
174                         unsigned int im = (Nx - i) & (Nx - 1);
175                         fREAL A = data[ji + i];
176                         fREAL B = data[jmi + i];
177                         fREAL C = data[ji + im];
178                         fREAL D = data[jmi + im];
179                         fREAL E = (fREAL)0.5 * ((A + D) - (B + C));
180                         data[ji + i] = A - E;
181                         data[jmi + i] = B + E;
182                         data[ji + im] = C + E;
183                         data[jmi + im] = D - E;
184                 }
185         }
186
187 }
188
189 //------------------------------------------------------------------------------
190
191 /* 2D convolution calc, d1 *= d2, M/N - > log2 of width/height */
192 static void fht_convolve(fREAL *d1, fREAL *d2, unsigned int M, unsigned int N)
193 {
194         fREAL a, b;
195         unsigned int i, j, k, L, mj, mL;
196         unsigned int m = 1 << M, n = 1 << N;
197         unsigned int m2 = 1 << (M - 1), n2 = 1 << (N - 1);
198         unsigned int mn2 = m << (N - 1);
199
200         d1[0] *= d2[0];
201         d1[mn2] *= d2[mn2];
202         d1[m2] *= d2[m2];
203         d1[m2 + mn2] *= d2[m2 + mn2];
204         for (i = 1; i < m2; i++) {
205                 k = m - i;
206                 a = d1[i] * d2[i] - d1[k] * d2[k];
207                 b = d1[k] * d2[i] + d1[i] * d2[k];
208                 d1[i] = (b + a) * (fREAL)0.5;
209                 d1[k] = (b - a) * (fREAL)0.5;
210                 a = d1[i + mn2] * d2[i + mn2] - d1[k + mn2] * d2[k + mn2];
211                 b = d1[k + mn2] * d2[i + mn2] + d1[i + mn2] * d2[k + mn2];
212                 d1[i + mn2] = (b + a) * (fREAL)0.5;
213                 d1[k + mn2] = (b - a) * (fREAL)0.5;
214         }
215         for (j = 1; j < n2; j++) {
216                 L = n - j;
217                 mj = j << M;
218                 mL = L << M;
219                 a = d1[mj] * d2[mj] - d1[mL] * d2[mL];
220                 b = d1[mL] * d2[mj] + d1[mj] * d2[mL];
221                 d1[mj] = (b + a) * (fREAL)0.5;
222                 d1[mL] = (b - a) * (fREAL)0.5;
223                 a = d1[m2 + mj] * d2[m2 + mj] - d1[m2 + mL] * d2[m2 + mL];
224                 b = d1[m2 + mL] * d2[m2 + mj] + d1[m2 + mj] * d2[m2 + mL];
225                 d1[m2 + mj] = (b + a) * (fREAL)0.5;
226                 d1[m2 + mL] = (b - a) * (fREAL)0.5;
227         }
228         for (i = 1; i < m2; i++) {
229                 k = m - i;
230                 for (j = 1; j < n2; j++) {
231                         L = n - j;
232                         mj = j << M;
233                         mL = L << M;
234                         a = d1[i + mj] * d2[i + mj] - d1[k + mL] * d2[k + mL];
235                         b = d1[k + mL] * d2[i + mj] + d1[i + mj] * d2[k + mL];
236                         d1[i + mj] = (b + a) * (fREAL)0.5;
237                         d1[k + mL] = (b - a) * (fREAL)0.5;
238                         a = d1[i + mL] * d2[i + mL] - d1[k + mj] * d2[k + mj];
239                         b = d1[k + mj] * d2[i + mL] + d1[i + mL] * d2[k + mj];
240                         d1[i + mL] = (b + a) * (fREAL)0.5;
241                         d1[k + mj] = (b - a) * (fREAL)0.5;
242                 }
243         }
244 }
245 //------------------------------------------------------------------------------
246
247 void convolve(float *dst, MemoryBuffer *in1, MemoryBuffer *in2)
248 {
249         fREAL *data1, *data2, *fp;
250         unsigned int w2, h2, hw, hh, log2_w, log2_h;
251         fRGB wt, *colp;
252         int x, y, ch;
253         int xbl, ybl, nxb, nyb, xbsz, ybsz;
254         int in2done = FALSE;
255         const unsigned int kernelWidth = in2->getWidth();
256         const unsigned int kernelHeight = in2->getHeight();
257         const unsigned int imageWidth = in1->getWidth();
258         const unsigned int imageHeight = in1->getHeight();
259         float *kernelBuffer = in2->getBuffer();
260         float *imageBuffer = in1->getBuffer();
261
262         MemoryBuffer *rdst = new MemoryBuffer(NULL, in1->getRect());
263         memset(rdst->getBuffer(), 0, rdst->getWidth() * rdst->getHeight() * COM_NUMBER_OF_CHANNELS * sizeof(float));
264
265         // convolution result width & height
266         w2 = 2 * kernelWidth - 1;
267         h2 = 2 * kernelHeight - 1;
268         // FFT pow2 required size & log2
269         w2 = nextPow2(w2, &log2_w);
270         h2 = nextPow2(h2, &log2_h);
271
272         // alloc space
273         data1 = (fREAL *)MEM_callocN(3 * w2 * h2 * sizeof(fREAL), "convolve_fast FHT data1");
274         data2 = (fREAL *)MEM_callocN(w2 * h2 * sizeof(fREAL), "convolve_fast FHT data2");
275
276         // normalize convolutor
277         wt[0] = wt[1] = wt[2] = 0.f;
278         for (y = 0; y < kernelHeight; y++) {
279                 colp = (fRGB *)&kernelBuffer[y * kernelWidth * COM_NUMBER_OF_CHANNELS];
280                 for (x = 0; x < kernelWidth; x++)
281                         add_v3_v3(wt, colp[x]);
282         }
283         if (wt[0] != 0.f) wt[0] = 1.f / wt[0];
284         if (wt[1] != 0.f) wt[1] = 1.f / wt[1];
285         if (wt[2] != 0.f) wt[2] = 1.f / wt[2];
286         for (y = 0; y < kernelHeight; y++) {
287                 colp = (fRGB *)&kernelBuffer[y * kernelWidth * COM_NUMBER_OF_CHANNELS];
288                 for (x = 0; x < kernelWidth; x++)
289                         mul_v3_v3(colp[x], wt);
290         }
291
292         // copy image data, unpacking interleaved RGBA into separate channels
293         // only need to calc data1 once
294
295         // block add-overlap
296         hw = kernelWidth >> 1;
297         hh = kernelHeight >> 1;
298         xbsz = (w2 + 1) - kernelWidth;
299         ybsz = (h2 + 1) - kernelHeight;
300         nxb = imageWidth / xbsz;
301         if (imageWidth % xbsz) nxb++;
302         nyb = imageHeight / ybsz;
303         if (imageHeight % ybsz) nyb++;
304         for (ybl = 0; ybl < nyb; ybl++) {
305                 for (xbl = 0; xbl < nxb; xbl++) {
306
307                         // each channel one by one
308                         for (ch = 0; ch < 3; ch++) {
309                                 fREAL *data1ch = &data1[ch * w2 * h2];
310
311                                 // only need to calc fht data from in2 once, can re-use for every block
312                                 if (!in2done) {
313                                         // in2, channel ch -> data1
314                                         for (y = 0; y < kernelHeight; y++) {
315                                                 fp = &data1ch[y * w2];
316                                                 colp = (fRGB *)&kernelBuffer[y * kernelWidth * COM_NUMBER_OF_CHANNELS];
317                                                 for (x = 0; x < kernelWidth; x++)
318                                                         fp[x] = colp[x][ch];
319                                         }
320                                 }
321
322                                 // in1, channel ch -> data2
323                                 memset(data2, 0, w2 * h2 * sizeof(fREAL));
324                                 for (y = 0; y < ybsz; y++) {
325                                         int yy = ybl * ybsz + y;
326                                         if (yy >= imageHeight) continue;
327                                         fp = &data2[y * w2];
328                                         colp = (fRGB *)&imageBuffer[yy * imageWidth * COM_NUMBER_OF_CHANNELS];
329                                         for (x = 0; x < xbsz; x++) {
330                                                 int xx = xbl * xbsz + x;
331                                                 if (xx >= imageWidth) continue;
332                                                 fp[x] = colp[xx][ch];
333                                         }
334                                 }
335
336                                 // forward FHT
337                                 // zero pad data start is different for each == height+1
338                                 if (!in2done) FHT2D(data1ch, log2_w, log2_h, kernelHeight + 1, 0);
339                                 FHT2D(data2, log2_w, log2_h, kernelHeight + 1, 0);
340
341                                 // FHT2D transposed data, row/col now swapped
342                                 // convolve & inverse FHT
343                                 fht_convolve(data2, data1ch, log2_h, log2_w);
344                                 FHT2D(data2, log2_h, log2_w, 0, 1);
345                                 // data again transposed, so in order again
346
347                                 // overlap-add result
348                                 for (y = 0; y < (int)h2; y++) {
349                                         const int yy = ybl * ybsz + y - hh;
350                                         if ((yy < 0) || (yy >= imageHeight)) continue;
351                                         fp = &data2[y * w2];
352                                         colp = (fRGB *)&rdst->getBuffer()[yy * imageWidth * COM_NUMBER_OF_CHANNELS];
353                                         for (x = 0; x < (int)w2; x++) {
354                                                 const int xx = xbl * xbsz + x - hw;
355                                                 if ((xx < 0) || (xx >= imageWidth)) continue;
356                                                 colp[xx][ch] += fp[x];
357                                         }
358                                 }
359
360                         }
361                         in2done = TRUE;
362                 }
363         }
364
365         MEM_freeN(data2);
366         MEM_freeN(data1);
367         memcpy(dst, rdst->getBuffer(), sizeof(float) * imageWidth * imageHeight * COM_NUMBER_OF_CHANNELS);
368         delete(rdst);
369 }
370
371 void GlareFogGlowOperation::generateGlare(float *data, MemoryBuffer *inputTile, NodeGlare *settings)
372 {
373         int x, y;
374         float scale, u, v, r, w, d;
375         fRGB fcol;
376         MemoryBuffer *ckrn;
377         unsigned int sz = 1 << settings->size;
378         const float cs_r = 1.f, cs_g = 1.f, cs_b = 1.f;
379
380         // temp. src image
381         // make the convolution kernel
382         rcti kernelRect;
383         BLI_rcti_init(&kernelRect, 0, sz, 0, sz);
384         ckrn = new MemoryBuffer(NULL, &kernelRect);
385
386         scale = 0.25f * sqrtf((float)(sz * sz));
387
388         for (y = 0; y < sz; ++y) {
389                 v = 2.f * (y / (float)sz) - 1.0f;
390                 for (x = 0; x < sz; ++x) {
391                         u = 2.f * (x / (float)sz) - 1.0f;
392                         r = (u * u + v * v) * scale;
393                         d = -sqrtf(sqrtf(sqrtf(r))) * 9.0f;
394                         fcol[0] = expf(d * cs_r), fcol[1] = expf(d * cs_g), fcol[2] = expf(d * cs_b);
395                         // linear window good enough here, visual result counts, not scientific analysis
396                         //w = (1.f-fabs(u))*(1.f-fabs(v));
397                         // actually, Hanning window is ok, cos^2 for some reason is slower
398                         w = (0.5f + 0.5f * cos((double)u * M_PI)) * (0.5f + 0.5f * cos((double)v * M_PI));
399                         mul_v3_fl(fcol, w);
400                         ckrn->writePixel(x, y, fcol);
401                 }
402         }
403
404         convolve(data, inputTile, ckrn);
405         delete ckrn;
406 }