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 }