Initial commit. Not in build system so shouldn't interfere with anything at this point. Will commit modified versions of existing files once build system is tested.

1 files changed, 875 insertions, 0 deletions

diff --git a/source/blender/nodes/intern/CMP_nodes/CMP_defocus.c b/source/blender/nodes/intern/CMP_nodes/CMP_defocus.c
new file mode 100644
index 00000000000..157b84a6d6a
--- /dev/null
+++ b/source/blender/nodes/intern/CMP_nodes/CMP_defocus.c
@@ -0,0 +1,875 @@
+/**
+ * $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);
+}
+
+/* qdn: defocus node */
+static int node_composit_buts_defocus(uiBlock *block, bNodeTree *ntree, bNode *node, rctf *butr)
+{
+   if(block) {
+      NodeDefocus *nqd = node->storage;
+      short dy = butr->ymin + 209;
+      short dx = butr->xmax - butr->xmin; 
+      char* mstr1 = "Bokeh Type%t|Octagon %x8|Heptagon %x7|Hexagon %x6|Pentagon %x5|Square %x4|Triangle %x3|Disk %x0";
+
+      uiDefBut(block, LABEL, B_NOP, "Bokeh Type", butr->xmin, dy, dx, 19, NULL, 0, 0, 0, 0, "");
+      uiDefButC(block, MENU, B_NODE_EXEC+node->nr, mstr1,
+         butr->xmin, dy-19, dx, 19,
+         &nqd->bktype, 0, 0, 0, 0, "Bokeh type");
+      if (nqd->bktype) { /* for some reason rotating a disk doesn't seem to work... ;) */
+         uiDefButC(block, NUM, B_NODE_EXEC+node->nr, "Rotate:",
+            butr->xmin, dy-38, dx, 19,
+            &nqd->rotation, 0, 90, 0, 0, "Bokeh shape rotation offset in degrees");
+      }
+      uiDefButC(block, TOG, B_NODE_EXEC+node->nr, "Gamma Correct",
+         butr->xmin, dy-57, dx, 19,
+         &nqd->gamco, 0, 0, 0, 0, "Enable gamma correction before and after main process");
+      if (nqd->no_zbuf==0) {
+         // only needed for zbuffer input
+         uiDefButF(block, NUM, B_NODE_EXEC+node->nr, "fStop:",
+            butr->xmin, dy-76, dx, 19,
+            &nqd->fstop, 0.5, 128, 10, 0, "Amount of focal blur, 128=infinity=perfect focus, half the value doubles the blur radius");
+      }
+      uiDefButF(block, NUM, B_NODE_EXEC+node->nr, "Maxblur:",
+         butr->xmin, dy-95, dx, 19,
+         &nqd->maxblur, 0, 10000, 1000, 0, "blur limit, maximum CoC radius, 0=no limit");
+      uiDefButF(block, NUM, B_NODE_EXEC+node->nr, "BThreshold:",
+         butr->xmin, dy-114, dx, 19,
+         &nqd->bthresh, 0, 100, 100, 0, "CoC radius threshold, prevents background bleed on in-focus midground, 0=off");
+      uiDefButC(block, TOG, B_NODE_EXEC+node->nr, "Preview",
+         butr->xmin, dy-142, dx, 19,
+         &nqd->preview, 0, 0, 0, 0, "Enable sampling mode, useful for preview when using low samplecounts");
+      if (nqd->preview) {
+         /* only visible when sampling mode enabled */
+         uiDefButS(block, NUM, B_NODE_EXEC+node->nr, "Samples:",
+            butr->xmin, dy-161, dx, 19,
+            &nqd->samples, 16, 256, 0, 0, "Number of samples (16=grainy, higher=less noise)");
+      }
+      uiDefButS(block, TOG, B_NODE_EXEC+node->nr, "No zbuffer",
+         butr->xmin, dy-190, dx, 19,
+         &nqd->no_zbuf, 0, 0, 0, 0, "Enable when using an image as input instead of actual zbuffer (auto enabled if node not image based, eg. time node)");
+      if (nqd->no_zbuf) {
+         uiDefButF(block, NUM, B_NODE_EXEC+node->nr, "Zscale:",
+            butr->xmin, dy-209, dx, 19,
+            &nqd->scale, 0, 1000, 100, 0, "Scales the Z input when not using a zbuffer, controls maximum blur designated by the color white or input value 1");
+      }
+   }
+   return 228;
+}
+
+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 = {
+	/* 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     */ node_composit_buts_defocus,
+                     node_composit_init_defocus
+};
+
+