Files
blender/source/blender/nodes/intern/CMP_nodes/CMP_defocus.c
Ton Roosendaal 5c760e4811 Changed node type definitions to use a dynamic list.
This will allow python or plugin defined nodes to work as well.
(And fixes compile issues with MSVC in yesterdays commit for nodes)

Code provided by Nathan L.
Fixes in his code:
- free_nodesystem() was called too late (after guarded alloc was closed)
- free_nodesystem() was freeing nodes that were not malloced even
- free_nodesystem was using free, not freeN :)
- the typedefs needed to be malloced yes, to allow duplicate nodes like
  group but also for dynamic nodes.
2007-03-26 15:07:38 +00:00

823 lines
25 KiB
C

/**
* $Id$
*
* ***** BEGIN GPL LICENSE BLOCK *****
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software Foundation,
* Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*
* The Original Code is Copyright (C) 2006 Blender Foundation.
* All rights reserved.
*
* The Original Code is: all of this file.
*
* Contributor(s): none yet.
*
* ***** END GPL LICENSE BLOCK *****
*/
#include "../CMP_util.h"
/* ************ qdn: Defocus node ****************** */
static bNodeSocketType cmp_node_defocus_in[]= {
{ SOCK_RGBA, 1, "Image", 0.8f, 0.8f, 0.8f, 1.0f, 0.0f, 1.0f},
{ SOCK_VALUE, 1, "Z", 0.8f, 0.8f, 0.8f, 1.0f, 0.0f, 1.0f},
{ -1, 0, "" }
};
static bNodeSocketType cmp_node_defocus_out[]= {
{ SOCK_RGBA, 0, "Image", 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 1.0f},
{ -1, 0, "" }
};
// line coefs for point sampling & scancon. data.
typedef struct BokehCoeffs {
float x0, y0, dx, dy;
float ls_x, ls_y;
float min_x, min_y, max_x, max_y;
} BokehCoeffs;
// returns array of BokehCoeffs
// returns length of array in 'len_bkh',
// radius squared of inscribed disk in 'inradsq', needed in getWeight() test,
// BKH[8] is the data returned for the bokeh shape & bkh_b[4] is it's 2d bound
static void makeBokeh(char bktype, char ro, int* len_bkh, float* inradsq, BokehCoeffs BKH[8], float bkh_b[4])
{
float x0, x1, y0, y1, dx, dy, iDxy, w = ro*M_PI/180.f;
float wi = (360.f/bktype)*M_PI/180.f;
int i, ov, nv;
// bktype must be at least 3 & <= 8
bktype = (bktype<3) ? 3 : ((bktype>8) ? 8 : bktype);
*len_bkh = bktype;
*inradsq = -1.f;
for (i=0; i<(*len_bkh); i++) {
x0 = cos(w);
y0 = sin(w);
w += wi;
x1 = cos(w);
y1 = sin(w);
if ((*inradsq)<0.f) {
// radius squared of inscribed disk
float idx=(x0+x1)*0.5f, idy=(y0+y1)*0.5f;
*inradsq = idx*idx + idy*idy;
}
BKH[i].x0 = x0;
BKH[i].y0 = y0;
dx = x1-x0, dy = y1-y0;
iDxy = 1.f / sqrt(dx*dx + dy*dy);
dx *= iDxy;
dy *= iDxy;
BKH[i].dx = dx;
BKH[i].dy = dy;
}
// precalc scanconversion data
// bokeh bound, not transformed, for scanconvert
bkh_b[0] = bkh_b[2] = 1e10f; // xmin/ymin
bkh_b[1] = bkh_b[3] = -1e10f; // xmax/ymax
ov = (*len_bkh) - 1;
for (nv=0; nv<(*len_bkh); nv++) {
bkh_b[0] = MIN2(bkh_b[0], BKH[nv].x0); // xmin
bkh_b[1] = MAX2(bkh_b[1], BKH[nv].x0); // xmax
bkh_b[2] = MIN2(bkh_b[2], BKH[nv].y0); // ymin
bkh_b[3] = MAX2(bkh_b[3], BKH[nv].y0); // ymax
BKH[nv].min_x = MIN2(BKH[ov].x0, BKH[nv].x0);
BKH[nv].max_x = MAX2(BKH[ov].x0, BKH[nv].x0);
BKH[nv].min_y = MIN2(BKH[ov].y0, BKH[nv].y0);
BKH[nv].max_y = MAX2(BKH[ov].y0, BKH[nv].y0);
dy = BKH[nv].y0 - BKH[ov].y0;
BKH[nv].ls_x = (BKH[nv].x0 - BKH[ov].x0) / ((dy==0.f) ? 1.f : dy);
BKH[nv].ls_y = (BKH[nv].ls_x==0.f) ? 1.f : (1.f/BKH[nv].ls_x);
ov = nv;
}
}
// test if u/v inside shape & returns weight value
static float getWeight(BokehCoeffs* BKH, int len_bkh, float u, float v, float rad, float inradsq)
{
BokehCoeffs* bc = BKH;
float cdist, irad = (rad==0.f) ? 1.f : (1.f/rad);
u *= irad;
v *= irad;
// early out test1: if point outside outer unit disk, it cannot be inside shape
cdist = u*u + v*v;
if (cdist>1.f) return 0.f;
// early out test2: if point inside or on inner disk, point must be inside shape
if (cdist<=inradsq) return 1.f;
while (len_bkh--) {
if ((bc->dy*(u - bc->x0) - bc->dx*(v - bc->y0)) > 0.f) return 0.f;
bc++;
}
return 1.f;
}
// QMC.seq. for sampling, A.Keller, EMS
static float RI_vdC(unsigned int bits, unsigned int r)
{
bits = ( bits << 16) | ( bits >> 16);
bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
bits ^= r;
return (float)((double)bits / 4294967296.0);
}
// single channel IIR gaussian filtering
// much faster than anything else, constant time independent of width
// should extend to multichannel and make this a node, could be useful
static void IIR_gauss(CompBuf* buf, float sigma)
{
double q, q2, sc, cf[4], tsM[9], tsu[3], tsv[3];
float *X, *Y, *W;
int i, x, y, sz;
// single channel only for now
if (buf->type != CB_VAL) return;
// <0.5 not valid, though can have a possibly useful sort of sharpening effect
if (sigma < 0.5) return;
// see "Recursive Gabor Filtering" by Young/VanVliet
// all factors here in double.prec. Required, because for single.prec it seems to blow up if sigma > ~200
if (sigma >= 3.556)
q = 0.9804*(sigma - 3.556) + 2.5091;
else // sigma >= 0.5
q = (0.0561*sigma + 0.5784)*sigma - 0.2568;
q2 = q*q;
sc = (1.1668 + q)*(3.203729649 + (2.21566 + q)*q);
// no gabor filtering here, so no complex multiplies, just the regular coefs.
// all negated here, so as not to have to recalc Triggs/Sdika matrix
cf[1] = q*(5.788961737 + (6.76492 + 3.0*q)*q)/ sc;
cf[2] = -q2*(3.38246 + 3.0*q)/sc;
// 0 & 3 unchanged
cf[3] = q2*q/sc;
cf[0] = 1.0 - cf[1] - cf[2] - cf[3];
// Triggs/Sdika border corrections,
// it seems to work, not entirely sure if it is actually totally correct,
// Besides J.M.Geusebroek's anigauss.c (see http://www.science.uva.nl/~mark),
// found one other implementation by Cristoph Lampert,
// but neither seem to be quite the same, result seems to be ok sofar anyway.
// Extra scale factor here to not have to do it in filter,
// though maybe this had something to with the precision errors
sc = cf[0]/((1.0 + cf[1] - cf[2] + cf[3])*(1.0 - cf[1] - cf[2] - cf[3])*(1.0 + cf[2] + (cf[1] - cf[3])*cf[3]));
tsM[0] = sc*(-cf[3]*cf[1] + 1.0 - cf[3]*cf[3] - cf[2]);
tsM[1] = sc*((cf[3] + cf[1])*(cf[2] + cf[3]*cf[1]));
tsM[2] = sc*(cf[3]*(cf[1] + cf[3]*cf[2]));
tsM[3] = sc*(cf[1] + cf[3]*cf[2]);
tsM[4] = sc*(-(cf[2] - 1.0)*(cf[2] + cf[3]*cf[1]));
tsM[5] = sc*(-(cf[3]*cf[1] + cf[3]*cf[3] + cf[2] - 1.0)*cf[3]);
tsM[6] = sc*(cf[3]*cf[1] + cf[2] + cf[1]*cf[1] - cf[2]*cf[2]);
tsM[7] = sc*(cf[1]*cf[2] + cf[3]*cf[2]*cf[2] - cf[1]*cf[3]*cf[3] - cf[3]*cf[3]*cf[3] - cf[3]*cf[2] + cf[3]);
tsM[8] = sc*(cf[3]*(cf[1] + cf[3]*cf[2]));
#define YVV(L)\
{\
W[0] = cf[0]*X[0] + cf[1]*X[0] + cf[2]*X[0] + cf[3]*X[0];\
W[1] = cf[0]*X[1] + cf[1]*W[0] + cf[2]*X[0] + cf[3]*X[0];\
W[2] = cf[0]*X[2] + cf[1]*W[1] + cf[2]*W[0] + cf[3]*X[0];\
for (i=3; i<L; i++)\
W[i] = cf[0]*X[i] + cf[1]*W[i-1] + cf[2]*W[i-2] + cf[3]*W[i-3];\
tsu[0] = W[L-1] - X[L-1];\
tsu[1] = W[L-2] - X[L-1];\
tsu[2] = W[L-3] - X[L-1];\
tsv[0] = tsM[0]*tsu[0] + tsM[1]*tsu[1] + tsM[2]*tsu[2] + X[L-1];\
tsv[1] = tsM[3]*tsu[0] + tsM[4]*tsu[1] + tsM[5]*tsu[2] + X[L-1];\
tsv[2] = tsM[6]*tsu[0] + tsM[7]*tsu[1] + tsM[8]*tsu[2] + X[L-1];\
Y[L-1] = cf[0]*W[L-1] + cf[1]*tsv[0] + cf[2]*tsv[1] + cf[3]*tsv[2];\
Y[L-2] = cf[0]*W[L-2] + cf[1]*Y[L-1] + cf[2]*tsv[0] + cf[3]*tsv[1];\
Y[L-3] = cf[0]*W[L-3] + cf[1]*Y[L-2] + cf[2]*Y[L-1] + cf[3]*tsv[0];\
for (i=L-4; i>=0; i--)\
Y[i] = cf[0]*W[i] + cf[1]*Y[i+1] + cf[2]*Y[i+2] + cf[3]*Y[i+3];\
}
// intermediate buffers
sz = MAX2(buf->x, buf->y);
Y = MEM_callocN(sz*sizeof(float), "IIR_gauss Y buf");
W = MEM_callocN(sz*sizeof(float), "IIR_gauss W buf");
// H
for (y=0; y<buf->y; y++) {
X = &buf->rect[y*buf->x];
YVV(buf->x);
memcpy(X, Y, sizeof(float)*buf->x);
}
// V
X = MEM_callocN(buf->y*sizeof(float), "IIR_gauss X buf");
for (x=0; x<buf->x; x++) {
for (y=0; y<buf->y; y++)
X[y] = buf->rect[x + y*buf->x];
YVV(buf->y);
for (y=0; y<buf->y; y++)
buf->rect[x + y*buf->x] = Y[y];
}
MEM_freeN(X);
MEM_freeN(W);
MEM_freeN(Y);
#undef YVV
}
static void defocus_blur(CompBuf* new, CompBuf* img, CompBuf* zbuf, float inpval, NodeDefocus* nqd)
{
CompBuf *wts; // weights buffer
CompBuf *crad; // CoC radius buffer
BokehCoeffs BKH[8]; // bokeh shape data, here never > 8 pts.
float bkh_b[4] = {0}; // shape 2D bound
unsigned int p, px, p4, zp, cp, cp4;
float *ctcol, u, v, iZ, ct_crad, bcrad, lwt, wt=0, cR2=0;
float dof_sp, maxfgc, nmaxc, scf, bk_hn_theta=0, inradsq=0;
float cam_fdist=1, cam_invfdist=1, cam_lens=35;
int x, y, sx, sy, len_bkh=0;
float aspect, aperture;
int minsz;
// get some required params from the current scene camera
Object* camob = G.scene->camera;
if (camob && camob->type==OB_CAMERA) {
Camera* cam = (Camera*)camob->data;
cam_lens = cam->lens;
cam_fdist = (cam->YF_dofdist==0.f) ? 1e10f : cam->YF_dofdist;
cam_invfdist = 1.f/cam_fdist;
}
// guess work here.. best match with raytraced result
minsz = MIN2(img->x, img->y);
dof_sp = (float)minsz / (16.f / cam_lens); // <- == aspect * MIN2(img->x, img->y) / tan(0.5f * fov);
// aperture
aspect = (img->x > img->y) ? (img->y / (float)img->x) : (img->x / (float)img->y);
aperture = 0.5f*(cam_lens / (aspect*32.f)) / nqd->fstop;
// if not disk, make bokeh coefficients and other needed data
if (nqd->bktype!=0) {
makeBokeh(nqd->bktype, nqd->rotation, &len_bkh, &inradsq, BKH, bkh_b);
bk_hn_theta = 0.5 * nqd->bktype * sin(2.0 * M_PI / nqd->bktype); // weight factor
}
// accumulated weights
wts = alloc_compbuf(img->x, img->y, CB_VAL, 1);
// CoC radius buffer
crad = alloc_compbuf(img->x, img->y, CB_VAL, 1);
// if 'no_zbuf' flag set (which is always set if input is not an image),
// values are instead interpreted directly as blur radius values
if (nqd->no_zbuf) {
// to prevent *reaaallly* big radius values and impossible calculation times,
// limit the maximum to half the image width or height, whichever is smaller
float maxr = 0.5f*(float)MIN2(img->x, img->y);
for (p=0; p<(unsigned int)(img->x*img->y); p++) {
crad->rect[p] = zbuf ? (zbuf->rect[p]*nqd->scale) : inpval;
// bug #5921, limit minimum
crad->rect[p] = MAX2(1e-5f, crad->rect[p]);
crad->rect[p] = MIN2(crad->rect[p], maxr);
// if maxblur!=0, limit maximum
if (nqd->maxblur != 0.f) crad->rect[p] = MIN2(crad->rect[p], nqd->maxblur);
}
}
else {
// actual zbuffer.
// separate foreground from background CoC's
// then blur background and blend in again with foreground,
// improves the 'blurred foreground overlapping in-focus midground' sharp boundary problem.
// wts buffer here used for blendmask
maxfgc = 0.f; // maximum foreground CoC radius
for (y=0; y<img->y; y++) {
p = y * img->x;
for (x=0; x<img->x; x++) {
px = p + x;
iZ = (zbuf->rect[px]==0.f) ? 0.f : (1.f/zbuf->rect[px]);
crad->rect[px] = 0.5f*(aperture*(dof_sp*(cam_invfdist - iZ) - 1.f));
if (crad->rect[px] <= 0.f) {
wts->rect[px] = 1.f;
crad->rect[px] = -crad->rect[px];
if (crad->rect[px] > maxfgc) maxfgc = crad->rect[px];
}
else crad->rect[px] = wts->rect[px] = 0;
}
}
// fast blur...
IIR_gauss(crad, 2.f*maxfgc);
IIR_gauss(wts, 2.f*maxfgc);
// find new maximum to scale it back to original
// (could skip this, not strictly necessary, in general, difference is quite small, but just in case...)
nmaxc = 0;
for (p=0; p<(img->x*img->y); p++)
if (crad->rect[p] > nmaxc) nmaxc = crad->rect[p];
// rescale factor
scf = (nmaxc==0.f) ? 1.f: (maxfgc / nmaxc);
// and blend...
for (y=0; y<img->y; y++) {
p = y*img->x;
for (x=0; x<img->x; x++) {
px = p + x;
if (zbuf->rect[px]!=0.f) {
iZ = (zbuf->rect[px]==0.f) ? 0.f : (1.f/zbuf->rect[px]);
bcrad = 0.5f*fabs(aperture*(dof_sp*(cam_invfdist - iZ) - 1.f));
// scale crad back to original maximum and blend
crad->rect[px] = bcrad + wts->rect[px]*(scf*crad->rect[px] - bcrad);
if (crad->rect[px] < 0.01f) crad->rect[px] = 0.01f;
// if maxblur!=0, limit maximum
if (nqd->maxblur != 0.f) crad->rect[px] = MIN2(crad->rect[px], nqd->maxblur);
}
else crad->rect[px] = 0.f;
// clear weights for next part
wts->rect[px] = 0.f;
}
}
}
//------------------------------------------------------------------
// main loop
for (y=0; y<img->y; y++) {
// some sort of visual feedback would be nice, or at least this text in the renderwin header
// but for now just print some info in the console every 8 scanlines.
if (((y & 7)==0) || (y==(img->y-1))) {
printf("\rdefocus: Processing Line %d of %d ... ", y+1, img->y);
fflush(stdout);
}
zp = y * img->x;
for (x=0; x<img->x; x++) {
cp = zp + x;
cp4 = cp * img->type;
// Circle of Confusion radius for current pixel
cR2 = ct_crad = crad->rect[cp];
// skip if zero (border render)
if (ct_crad==0.f) {
// related to bug #5921, forgot output image when skipping 0 radius values
new->rect[cp4] = img->rect[cp4];
if (new->type != CB_VAL) {
new->rect[cp4+1] = img->rect[cp4+1];
new->rect[cp4+2] = img->rect[cp4+2];
new->rect[cp4+3] = img->rect[cp4+3];
}
continue;
}
cR2 *= cR2;
// pixel color
ctcol = &img->rect[cp4];
if (!nqd->preview) {
int xs, xe, ys, ye;
float lwt, wtcol[4] = {0}, aacol[4] = {0};
// shape weight
if (nqd->bktype==0) // disk
wt = 1.f/((float)M_PI*cR2);
else
wt = 1.f/(cR2*bk_hn_theta);
// weighted color
wtcol[0] = wt*ctcol[0];
if (new->type != CB_VAL) {
wtcol[1] = wt*ctcol[1];
wtcol[2] = wt*ctcol[2];
wtcol[3] = wt*ctcol[3];
}
// macro for background blur overlap test
// unfortunately, since this is done per pixel,
// it has a very significant negative impact on processing time...
// (eg. aa disk blur without test: 112 sec, vs with test: 176 sec...)
// iff center blur radius > threshold
// and if overlap pixel in focus, do nothing, else add color/weigbt
// (threshold constant is dependant on amount of blur)
#define TESTBG1(c, w) {\
if (ct_crad > nqd->bthresh) {\
if (crad->rect[p] > nqd->bthresh) {\
new->rect[p] += c[0];\
wts->rect[p] += w;\
}\
}\
else {\
new->rect[p] += c[0];\
wts->rect[p] += w;\
}\
}
#define TESTBG4(c, w) {\
if (ct_crad > nqd->bthresh) {\
if (crad->rect[p] > nqd->bthresh) {\
new->rect[p4] += c[0];\
new->rect[p4+1] += c[1];\
new->rect[p4+2] += c[2];\
new->rect[p4+3] += c[3];\
wts->rect[p] += w;\
}\
}\
else {\
new->rect[p4] += c[0];\
new->rect[p4+1] += c[1];\
new->rect[p4+2] += c[2];\
new->rect[p4+3] += c[3];\
wts->rect[p] += w;\
}\
}
if (nqd->bktype == 0) {
// Disk
int _x, i, j, di;
float Dj, T;
// AA pixel
#define AAPIX(a, b) {\
int _ny = b;\
if ((_ny >= 0) && (_ny < new->y)) {\
int _nx = a;\
if ((_nx >=0) && (_nx < new->x)) {\
p = _ny*new->x + _nx;\
if (new->type==CB_VAL) {\
TESTBG1(aacol, lwt);\
}\
else {\
p4 = p * new->type;\
TESTBG4(aacol, lwt);\
}\
}\
}\
}
// circle scanline
#define CSCAN(a, b) {\
int _ny = y + b;\
if ((_ny >= 0) && (_ny < new->y)) {\
xs = x - a + 1;\
if (xs < 0) xs = 0;\
xe = x + a;\
if (xe > new->x) xe = new->x;\
p = _ny*new->x + xs;\
if (new->type==CB_VAL) {\
for (_x=xs; _x<xe; _x++, p++) TESTBG1(wtcol, wt);\
}\
else {\
p4 = p * new->type;\
for (_x=xs; _x<xe; _x++, p++, p4+=new->type) TESTBG4(wtcol, wt);\
}\
}\
}
i = ceil(ct_crad);
j = 0;
T = 0;
while (i > j) {
Dj = sqrt(cR2 - j*j);
Dj -= floor(Dj);
di = 0;
if (Dj > T) { i--; di = 1; }
T = Dj;
aacol[0] = wtcol[0]*Dj;
if (new->type != CB_VAL) {
aacol[1] = wtcol[1]*Dj;
aacol[2] = wtcol[2]*Dj;
aacol[3] = wtcol[3]*Dj;
}
lwt = wt*Dj;
if (i!=j) {
// outer pixels
AAPIX(x+j, y+i);
AAPIX(x+j, y-i);
if (j) {
AAPIX(x-j, y+i); // BL
AAPIX(x-j, y-i); // TL
}
if (di) { // only when i changed, interior of outer section
CSCAN(j, i); // bottom
CSCAN(j, -i); // top
}
}
// lower mid section
AAPIX(x+i, y+j);
if (i) AAPIX(x-i, y+j);
CSCAN(i, j);
// upper mid section
if (j) {
AAPIX(x+i, y-j);
if (i) AAPIX(x-i, y-j);
CSCAN(i, -j);
}
j++;
}
#undef CSCAN
#undef AAPIX
}
else {
// n-agonal
int ov, nv;
float mind, maxd, lwt;
ys = MAX2((int)floor(bkh_b[2]*ct_crad + y), 0);
ye = MIN2((int)ceil(bkh_b[3]*ct_crad + y), new->y - 1);
for (sy=ys; sy<=ye; sy++) {
float fxs = 1e10f, fxe = -1e10f;
float yf = (sy - y)/ct_crad;
int found = 0;
ov = len_bkh - 1;
mind = maxd = 0;
for (nv=0; nv<len_bkh; nv++) {
if ((BKH[nv].max_y >= yf) && (BKH[nv].min_y <= yf)) {
float tx = BKH[ov].x0 + BKH[nv].ls_x*(yf - BKH[ov].y0);
if (tx < fxs) { fxs = tx; mind = BKH[nv].ls_x; }
if (tx > fxe) { fxe = tx; maxd = BKH[nv].ls_x; }
if (++found == 2) break;
}
ov = nv;
}
if (found) {
fxs = fxs*ct_crad + x;
fxe = fxe*ct_crad + x;
xs = (int)floor(fxs), xe = (int)ceil(fxe);
// AA hack for first and last x pixel, near vertical edges only
if (fabs(mind) <= 1.f) {
if ((xs >= 0) && (xs < new->x)) {
lwt = 1.f-(fxs - xs);
aacol[0] = wtcol[0]*lwt;
p = xs + sy*new->x;
if (new->type==CB_VAL) {
lwt *= wt;
TESTBG1(aacol, lwt);
}
else {
p4 = p * new->type;
aacol[1] = wtcol[1]*lwt;
aacol[2] = wtcol[2]*lwt;
aacol[3] = wtcol[3]*lwt;
lwt *= wt;
TESTBG4(aacol, lwt);
}
}
}
if (fabs(maxd) <= 1.f) {
if ((xe >= 0) && (xe < new->x)) {
lwt = 1.f-(xe - fxe);
aacol[0] = wtcol[0]*lwt;
p = xe + sy*new->x;
if (new->type==CB_VAL) {
lwt *= wt;
TESTBG1(aacol, lwt);
}
else {
p4 = p * new->type;
aacol[1] = wtcol[1]*lwt;
aacol[2] = wtcol[2]*lwt;
aacol[3] = wtcol[3]*lwt;
lwt *= wt;
TESTBG4(aacol, lwt);
}
}
}
xs = MAX2(xs+1, 0);
xe = MIN2(xe, new->x);
// remaining interior scanline
p = sy*new->x + xs;
if (new->type==CB_VAL) {
for (sx=xs; sx<xe; sx++, p++) TESTBG1(wtcol, wt);
}
else {
p4 = p * new->type;
for (sx=xs; sx<xe; sx++, p++, p4+=new->type) TESTBG4(wtcol, wt);
}
}
}
// now traverse in opposite direction, y scanlines,
// but this time only draw the near horizontal edges,
// applying same AA hack as above
xs = MAX2((int)floor(bkh_b[0]*ct_crad + x), 0);
xe = MIN2((int)ceil(bkh_b[1]*ct_crad + x), img->x - 1);
for (sx=xs; sx<=xe; sx++) {
float xf = (sx - x)/ct_crad;
float fys = 1e10f, fye = -1e10f;
int found = 0;
ov = len_bkh - 1;
mind = maxd = 0;
for (nv=0; nv<len_bkh; nv++) {
if ((BKH[nv].max_x >= xf) && (BKH[nv].min_x <= xf)) {
float ty = BKH[ov].y0 + BKH[nv].ls_y*(xf - BKH[ov].x0);
if (ty < fys) { fys = ty; mind = BKH[nv].ls_y; }
if (ty > fye) { fye = ty; maxd = BKH[nv].ls_y; }
if (++found == 2) break;
}
ov = nv;
}
if (found) {
fys = fys*ct_crad + y;
fye = fye*ct_crad + y;
// near horizontal edges only, line slope <= 1
if (fabs(mind) <= 1.f) {
int iys = (int)floor(fys);
if ((iys >= 0) && (iys < new->y)) {
lwt = 1.f - (fys - iys);
aacol[0] = wtcol[0]*lwt;
p = sx + iys*new->x;
if (new->type==CB_VAL) {
lwt *= wt;
TESTBG1(aacol, lwt);
}
else {
p4 = p * new->type;
aacol[1] = wtcol[1]*lwt;
aacol[2] = wtcol[2]*lwt;
aacol[3] = wtcol[3]*lwt;
lwt *= wt;
TESTBG4(aacol, lwt);
}
}
}
if (fabs(maxd) <= 1.f) {
int iye = ceil(fye);
if ((iye >= 0) && (iye < new->y)) {
lwt = 1.f - (iye - fye);
aacol[0] = wtcol[0]*lwt;
p = sx + iye*new->x;
if (new->type==CB_VAL) {
lwt *= wt;
TESTBG1(aacol, lwt);
}
else {
p4 = p * new->type;
aacol[1] = wtcol[1]*lwt;
aacol[2] = wtcol[2]*lwt;
aacol[3] = wtcol[3]*lwt;
lwt *= wt;
TESTBG4(aacol, lwt);
}
}
}
}
}
}
#undef TESTBG4
#undef TESTBG1
}
else {
// sampled, simple rejection sampling here, good enough
unsigned int maxsam, s, ui = BLI_rand()*BLI_rand();
float wcor, cpr = BLI_frand();
if (nqd->no_zbuf)
maxsam = nqd->samples; // no zbuffer input, use sample value directly
else {
// depth adaptive sampling hack, the more out of focus, the more samples taken, 16 minimum.
maxsam = (int)(0.5f + nqd->samples*(1.f-(float)exp(-fabs(zbuf->rect[cp] - cam_fdist))));
if (maxsam < 16) maxsam = 16;
}
wcor = 1.f/(float)maxsam;
for (s=0; s<maxsam; ++s) {
u = ct_crad*(2.f*RI_vdC(s, ui) - 1.f);
v = ct_crad*(2.f*(s + cpr)/(float)maxsam - 1.f);
sx = (int)(x + u + 0.5f), sy = (int)(y + v + 0.5f);
if ((sx<0) || (sx >= new->x) || (sy<0) || (sy >= new->y)) continue;
p = sx + sy*new->x;
p4 = p * new->type;
if (nqd->bktype==0) // Disk
lwt = ((u*u + v*v)<=cR2) ? wcor : 0.f;
else // AA not needed here
lwt = wcor * getWeight(BKH, len_bkh, u, v, ct_crad, inradsq);
// prevent background bleeding onto in-focus pixels, user-option
if (ct_crad > nqd->bthresh) { // if center blur > threshold
if (crad->rect[p] > nqd->bthresh) { // if overlap pixel in focus, do nothing, else add color/weigbt
new->rect[p4] += ctcol[0] * lwt;
if (new->type != CB_VAL) {
new->rect[p4+1] += ctcol[1] * lwt;
new->rect[p4+2] += ctcol[2] * lwt;
new->rect[p4+3] += ctcol[3] * lwt;
}
wts->rect[p] += lwt;
}
}
else {
new->rect[p4] += ctcol[0] * lwt;
if (new->type != CB_VAL) {
new->rect[p4+1] += ctcol[1] * lwt;
new->rect[p4+2] += ctcol[2] * lwt;
new->rect[p4+3] += ctcol[3] * lwt;
}
wts->rect[p] += lwt;
}
}
}
}
}
// finally, normalize
for (y=0; y<new->y; y++) {
p = y * new->x;
p4 = p * new->type;
for (x=0; x<new->x; x++) {
float dv = (wts->rect[p]==0.f) ? 1.f : (1.f/wts->rect[p]);
new->rect[p4] *= dv;
if (new->type!=CB_VAL) {
new->rect[p4+1] *= dv;
new->rect[p4+2] *= dv;
new->rect[p4+3] *= dv;
}
p++;
p4 += new->type;
}
}
free_compbuf(crad);
free_compbuf(wts);
printf("Done\n");
}
static void node_composit_exec_defocus(void *data, bNode *node, bNodeStack **in, bNodeStack **out)
{
CompBuf *new, *old, *zbuf_use = NULL, *img = in[0]->data, *zbuf = in[1]->data;
NodeDefocus* nqd = node->storage;
if ((img==NULL) || (out[0]->hasoutput==0)) return;
// if image not valid type or fstop==infinite (128), nothing to do, pass in to out
if (((img->type!=CB_RGBA) && (img->type!=CB_VAL)) || ((nqd->no_zbuf==0) && (nqd->fstop==128.f))) {
new = alloc_compbuf(img->x, img->y, img->type, 0);
new->rect = img->rect;
out[0]->data = new;
return;
}
if (zbuf!=NULL) {
// Zbuf input, check to make sure, single channel, same size
// doesn't have to be actual zbuffer, but must be value type
if ((zbuf->x != img->x) || (zbuf->y != img->y)) {
// could do a scale here instead...
printf("Z input must be same size as image !\n");
return;
}
zbuf_use = typecheck_compbuf(zbuf, CB_VAL);
}
else nqd->no_zbuf = 1; // no zbuffer input
// ok, process
old = img;
if (nqd->gamco) {
// gamma correct, blender func is simplified, fixed value & RGBA only, should make user param
old = dupalloc_compbuf(img);
gamma_correct_compbuf(old, 0);
}
new = alloc_compbuf(old->x, old->y, old->type, 1);
defocus_blur(new, old, zbuf_use, in[1]->vec[0]*nqd->scale, node->storage);
if (nqd->gamco) {
gamma_correct_compbuf(new, 1);
free_compbuf(old);
}
out[0]->data = new;
if (zbuf_use && (zbuf_use != zbuf)) free_compbuf(zbuf_use);
}
static void node_composit_init_defocus(bNode* node)
{
/* qdn: defocus node */
NodeDefocus *nbd = MEM_callocN(sizeof(NodeDefocus), "node defocus data");
nbd->bktype = 0;
nbd->rotation = 0.f;
nbd->preview = 1;
nbd->gamco = 0;
nbd->samples = 16;
nbd->fstop = 128.f;
nbd->maxblur = 0;
nbd->bthresh = 1.f;
nbd->scale = 1.f;
nbd->no_zbuf = 1;
node->storage = nbd;
};
bNodeType cmp_node_defocus = {
/* *next,*prev */ NULL, NULL,
/* type code */ CMP_NODE_DEFOCUS,
/* name */ "Defocus",
/* width+range */ 150, 120, 200,
/* class+opts */ NODE_CLASS_OP_FILTER, NODE_OPTIONS,
/* input sock */ cmp_node_defocus_in,
/* output sock */ cmp_node_defocus_out,
/* storage */ "NodeDefocus",
/* execfunc */ node_composit_exec_defocus,
/* butfunc */ NULL,
/* initfunc */ node_composit_init_defocus
};