archive/blender-archive
Archived
1
1
Fork 0
Code Pull Requests Packages Activity
This repository has been archived on 2023-10-09. You can view files and clone it. You cannot open issues or pull requests or push a commit.
Files
367cf47dd68f210eeaf9ba6b51a1fdbcf93d209b
blender-archive/source/blender/nodes/intern/CMP_nodes/CMP_defocus.c
Nathan Letwory 19a9cffdf0 doxygen: blender/nodes tagged.
2011-02-27 20:13:22 +00:00

893 lines
27 KiB
C
Raw Blame History

/*
* $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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 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 *****
*/
/** \file blender/nodes/intern/CMP_nodes/CMP_defocus.c
* \ingroup cmpnodes
*/
#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;
float w = MAX2(1e-5f, ro)*M_PI/180.f; // never reported stangely enough, but a zero offset causes missing center line...
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_single(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(bNode *node, CompBuf *new, CompBuf *img, CompBuf *zbuf, float inpval, int no_zbuf)
{
NodeDefocus *nqd = node->storage;
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
float cam_fdist=1, cam_invfdist=1, cam_lens=35;
float dof_sp, maxfgc, bk_hn_theta=0, inradsq=0;
int y, len_bkh=0, ydone=0;
float aspect, aperture;
int minsz;
//float bcrad, nmaxc, scf;
// get some required params from the current scene camera
// (ton) this is wrong, needs fixed
Scene *scene= (Scene*)node->id;
Object* camob = (scene)? scene->camera: NULL;
if (camob && camob->type==OB_CAMERA) {
Camera* cam = (Camera*)camob->data;
cam_lens = cam->lens;
cam_fdist = dof_camera(camob);
if (cam_fdist==0.0) cam_fdist = 1e10f; /* if the dof is 0.0 then set it be be far away */
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 (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);
unsigned int p;
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 {
float wt;
// 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++) {
unsigned int p = y * img->x;
int x;
for (x=0; x<img->x; x++) {
unsigned int px = p + x;
float 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...
// bug #6656 part 1, probably when previous node_composite.c was split into separate files, it was not properly updated
// to include recent cvs commits (well, at least not defocus node), so this part was missing...
wt = aperture*128.f;
IIR_gauss_single(crad, wt);
IIR_gauss_single(wts, wt);
// bug #6656 part 2a, although foreground blur is not based anymore on closest object,
// the rescaling op below was still based on that anyway, and unlike the comment in below code,
// the difference is therefore not always that small at all...
// so for now commented out, not sure if this is going to cause other future problems, lets just wait and see...
/*
// 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++) {
unsigned int p = y*img->x;
int x;
for (x=0; x<img->x; x++) {
unsigned px = p + x;
if (zbuf->rect[px]!=0.f) {
float iZ = (zbuf->rect[px]==0.f) ? 0.f : (1.f/zbuf->rect[px]);
// bug #6656 part 2b, do not rescale
/*
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);
*/
crad->rect[px] = 0.5f*fabs(aperture*(dof_sp*(cam_invfdist - iZ) - 1.f));
// 'bug' #6615, limit minimum radius to 1 pixel, not really a solution, but somewhat mitigates the problem
crad->rect[px] = MAX2(crad->rect[px], 0.5f);
// 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;
}
// esc set by main calling process
if(node->exec & NODE_BREAK)
break;
}
}
//------------------------------------------------------------------
// main loop
#ifndef __APPLE__ /* can crash on Mac, see bug #22856, disabled for now */
#ifdef __INTEL_COMPILER /* icc doesn't like the compound statement -- internal error: 0_1506 */
#pragma omp parallel for private(y) if(!nqd->preview) schedule(guided)
#else
#pragma omp parallel for private(y) if(!nqd->preview && img->y*img->x > 16384) schedule(guided)
#endif
#endif
for (y=0; y<img->y; y++) {
unsigned int p, p4, zp, cp, cp4;
float *ctcol, u, v, ct_crad, cR2=0;
int x, sx, sy;
// 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.
#pragma omp critical
{
if (((ydone & 7)==0) || (ydone==(img->y-1))) {
if(G.background==0) {
printf("\rdefocus: Processing Line %d of %d ... ", ydone+1, img->y);
fflush(stdout);
}
}
ydone++;
}
// esc set by main calling process. don't break because openmp doesn't
// allow it, just continue and do nothing
if(node->exec & NODE_BREAK)
continue;
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};
float wt;
// 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(), lwt;
if (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++) {
unsigned int p = y * new->x;
unsigned int p4 = p * new->type;
int x;
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 *UNUSED(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;
int no_zbuf = nqd->no_zbuf;
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)) || ((no_zbuf==0) && (nqd->fstop==128.f))) {
out[0]->data = pass_on_compbuf(img);
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 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. also depremul and premul afterwards, gamma
// correction can't work with premul alpha
old = dupalloc_compbuf(img);
premul_compbuf(old, 1);
gamma_correct_compbuf(old, 0);
premul_compbuf(old, 0);
}
new = alloc_compbuf(old->x, old->y, old->type, 1);
defocus_blur(node, new, old, zbuf_use, in[1]->vec[0]*nqd->scale, no_zbuf);
if (nqd->gamco) {
premul_compbuf(new, 1);
gamma_correct_compbuf(new, 1);
premul_compbuf(new, 0);
free_compbuf(old);
}
if(node->exec & NODE_BREAK) {
free_compbuf(new);
new= NULL;
}
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;
}
void register_node_type_cmp_defocus(ListBase *lb)
{
static bNodeType ntype;
node_type_base(&ntype, CMP_NODE_DEFOCUS, "Defocus", NODE_CLASS_OP_FILTER, NODE_OPTIONS,
cmp_node_defocus_in, cmp_node_defocus_out);
node_type_size(&ntype, 150, 120, 200);
node_type_init(&ntype, node_composit_init_defocus);
node_type_storage(&ntype, "NodeDefocus", node_free_standard_storage, node_copy_standard_storage);
node_type_exec(&ntype, node_composit_exec_defocus);
nodeRegisterType(lb, &ntype);
}
View Git Blame Copy Permalink