updating raskter to support tiles compositor. this commit puts in some groundwork code to support tiles's pixel processor

6 files changed, 2742 insertions, 1289 deletions

diff --git a/intern/raskter/CMakeLists.txt b/intern/raskter/CMakeLists.txt
index 3e1368d8eb0..06421f0e71f 100644
--- a/intern/raskter/CMakeLists.txt
+++ b/intern/raskter/CMakeLists.txt
@@ -33,8 +33,11 @@ set(INC_SYS
 
 set(SRC
 	raskter.c
+	raskter_mt.c
+	kdtree.c
 
 	raskter.h
+	kdtree.h
 )
 
 blender_add_lib(bf_intern_raskter "${SRC}" "${INC}" "${INC_SYS}")
diff --git a/intern/raskter/kdtree.c b/intern/raskter/kdtree.c
new file mode 100644
index 00000000000..8b84ed412c1
--- /dev/null
+++ b/intern/raskter/kdtree.c
@@ -0,0 +1,836 @@
+/*
+This file is part of ``kdtree'', a library for working with kd-trees.
+Copyright (C) 2007-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+3. The name of the author may not be used to endorse or promote products
+   derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
+WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
+OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+OF SUCH DAMAGE.
+*/
+/* single nearest neighbor search written by Tamas Nepusz <tamas@cs.rhul.ac.uk> */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "kdtree.h"
+
+#if defined(WIN32) || defined(__WIN32__)
+#include <malloc.h>
+#endif
+
+#ifdef USE_LIST_NODE_ALLOCATOR
+
+#ifndef NO_PTHREADS
+#include <pthread.h>
+#else
+
+#ifndef I_WANT_THREAD_BUGS
+#error "You are compiling with the fast list node allocator, with pthreads disabled! This WILL break if used from multiple threads."
+#endif	/* I want thread bugs */
+
+#endif	/* pthread support */
+#endif	/* use list node allocator */
+
+struct kdhyperrect {
+	int dim;
+	double *min, *max;              /* minimum/maximum coords */
+};
+
+struct kdnode {
+	double *pos;
+	int dir;
+	void *data;
+
+	struct kdnode *left, *right;	/* negative/positive side */
+};
+
+struct res_node {
+	struct kdnode *item;
+	double dist_sq;
+	struct res_node *next;
+};
+
+struct kdtree {
+	int dim;
+	struct kdnode *root;
+	struct kdhyperrect *rect;
+	void (*destr)(void*);
+};
+
+struct kdres {
+	struct kdtree *tree;
+	struct res_node *rlist, *riter;
+	int size;
+};
+
+#define SQ(x)			((x) * (x))
+
+
+static void clear_rec(struct kdnode *node, void (*destr)(void*));
+static int insert_rec(struct kdnode **node, const double *pos, void *data, int dir, int dim);
+static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq);
+static void clear_results(struct kdres *set);
+
+static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max);
+static void hyperrect_free(struct kdhyperrect *rect);
+static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect);
+static void hyperrect_extend(struct kdhyperrect *rect, const double *pos);
+static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos);
+
+#ifdef USE_LIST_NODE_ALLOCATOR
+static struct res_node *alloc_resnode(void);
+static void free_resnode(struct res_node*);
+#else
+#define alloc_resnode()		malloc(sizeof(struct res_node))
+#define free_resnode(n)		free(n)
+#endif
+
+
+
+struct kdtree *kd_create(int k)
+{
+	struct kdtree *tree;
+
+	if(!(tree = malloc(sizeof *tree))) {
+		return 0;
+	}
+
+	tree->dim = k;
+	tree->root = 0;
+	tree->destr = 0;
+	tree->rect = 0;
+
+	return tree;
+}
+
+void kd_free(struct kdtree *tree)
+{
+	if(tree) {
+		kd_clear(tree);
+		free(tree);
+	}
+}
+
+static void clear_rec(struct kdnode *node, void (*destr)(void*))
+{
+	if(!node) return;
+
+	clear_rec(node->left, destr);
+	clear_rec(node->right, destr);
+	
+	if(destr) {
+		destr(node->data);
+	}
+	free(node->pos);
+	free(node);
+}
+
+void kd_clear(struct kdtree *tree)
+{
+	clear_rec(tree->root, tree->destr);
+	tree->root = 0;
+
+	if (tree->rect) {
+		hyperrect_free(tree->rect);
+		tree->rect = 0;
+	}
+}
+
+void kd_data_destructor(struct kdtree *tree, void (*destr)(void*))
+{
+	tree->destr = destr;
+}
+
+
+static int insert_rec(struct kdnode **nptr, const double *pos, void *data, int dir, int dim)
+{
+	int new_dir;
+	struct kdnode *node;
+
+	if(!*nptr) {
+		if(!(node = malloc(sizeof *node))) {
+			return -1;
+		}
+		if(!(node->pos = malloc(dim * sizeof *node->pos))) {
+			free(node);
+			return -1;
+		}
+		memcpy(node->pos, pos, dim * sizeof *node->pos);
+		node->data = data;
+		node->dir = dir;
+		node->left = node->right = 0;
+		*nptr = node;
+		return 0;
+	}
+
+	node = *nptr;
+	new_dir = (node->dir + 1) % dim;
+	if(pos[node->dir] < node->pos[node->dir]) {
+		return insert_rec(&(*nptr)->left, pos, data, new_dir, dim);
+	}
+	return insert_rec(&(*nptr)->right, pos, data, new_dir, dim);
+}
+
+int kd_insert(struct kdtree *tree, const double *pos, void *data)
+{
+	if (insert_rec(&tree->root, pos, data, 0, tree->dim)) {
+		return -1;
+	}
+
+	if (tree->rect == 0) {
+		tree->rect = hyperrect_create(tree->dim, pos, pos);
+	} else {
+		hyperrect_extend(tree->rect, pos);
+	}
+
+	return 0;
+}
+
+int kd_insertf(struct kdtree *tree, const float *pos, void *data)
+{
+	static double sbuf[16];
+	double *bptr, *buf = 0;
+	int res, dim = tree->dim;
+
+	if(dim > 16) {
+#ifndef NO_ALLOCA
+		if(dim <= 256)
+			bptr = buf = alloca(dim * sizeof *bptr);
+		else
+#endif
+			if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
+				return -1;
+			}
+	} else {
+		bptr = buf = sbuf;
+	}
+
+	while(dim-- > 0) {
+		*bptr++ = *pos++;
+	}
+
+	res = kd_insert(tree, buf, data);
+#ifndef NO_ALLOCA
+	if(tree->dim > 256)
+#else
+	if(tree->dim > 16)
+#endif
+		free(buf);
+	return res;
+}
+
+int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data)
+{
+	double buf[3];
+	buf[0] = x;
+	buf[1] = y;
+	buf[2] = z;
+	return kd_insert(tree, buf, data);
+}
+
+int kd_insert3f(struct kdtree *tree, float x, float y, float z, void *data)
+{
+	double buf[3];
+	buf[0] = x;
+	buf[1] = y;
+	buf[2] = z;
+	return kd_insert(tree, buf, data);
+}
+
+static int find_nearest(struct kdnode *node, const double *pos, double range, struct res_node *list, int ordered, int dim)
+{
+	double dist_sq, dx;
+	int i, ret, added_res = 0;
+
+	if(!node) return 0;
+
+	dist_sq = 0;
+	for(i=0; i<dim; i++) {
+		dist_sq += SQ(node->pos[i] - pos[i]);
+	}
+	if(dist_sq <= SQ(range)) {
+		if(rlist_insert(list, node, ordered ? dist_sq : -1.0) == -1) {
+			return -1;
+		}
+		added_res = 1;
+	}
+
+	dx = pos[node->dir] - node->pos[node->dir];
+
+	ret = find_nearest(dx <= 0.0 ? node->left : node->right, pos, range, list, ordered, dim);
+	if(ret >= 0 && fabs(dx) < range) {
+		added_res += ret;
+		ret = find_nearest(dx <= 0.0 ? node->right : node->left, pos, range, list, ordered, dim);
+	}
+	if(ret == -1) {
+		return -1;
+	}
+	added_res += ret;
+
+	return added_res;
+}
+
+#if 0
+static int find_nearest_n(struct kdnode *node, const double *pos, double range, int num, struct rheap *heap, int dim)
+{
+	double dist_sq, dx;
+	int i, ret, added_res = 0;
+
+	if(!node) return 0;
+	
+	/* if the photon is close enough, add it to the result heap */
+	dist_sq = 0;
+	for(i=0; i<dim; i++) {
+		dist_sq += SQ(node->pos[i] - pos[i]);
+	}
+	if(dist_sq <= range_sq) {
+		if(heap->size >= num) {
+			/* get furthest element */
+			struct res_node *maxelem = rheap_get_max(heap);
+
+			/* and check if the new one is closer than that */
+			if(maxelem->dist_sq > dist_sq) {
+				rheap_remove_max(heap);
+
+				if(rheap_insert(heap, node, dist_sq) == -1) {
+					return -1;
+				}
+				added_res = 1;
+
+				range_sq = dist_sq;
+			}
+		} else {
+			if(rheap_insert(heap, node, dist_sq) == -1) {
+				return =1;
+			}
+			added_res = 1;
+		}
+	}
+
+
+	/* find signed distance from the splitting plane */
+	dx = pos[node->dir] - node->pos[node->dir];
+
+	ret = find_nearest_n(dx <= 0.0 ? node->left : node->right, pos, range, num, heap, dim);
+	if(ret >= 0 && fabs(dx) < range) {
+		added_res += ret;
+		ret = find_nearest_n(dx <= 0.0 ? node->right : node->left, pos, range, num, heap, dim);
+	}
+
+}
+#endif
+
+static void kd_nearest_i(struct kdnode *node, const double *pos, struct kdnode **result, double *result_dist_sq, struct kdhyperrect* rect)
+{
+	int dir = node->dir;
+	int i;
+	double dummy, dist_sq;
+	struct kdnode *nearer_subtree, *farther_subtree;
+	double *nearer_hyperrect_coord, *farther_hyperrect_coord;
+
+	/* Decide whether to go left or right in the tree */
+	dummy = pos[dir] - node->pos[dir];
+	if (dummy <= 0) {
+		nearer_subtree = node->left;
+		farther_subtree = node->right;
+		nearer_hyperrect_coord = rect->max + dir;
+		farther_hyperrect_coord = rect->min + dir;
+	} else {
+		nearer_subtree = node->right;
+		farther_subtree = node->left;
+		nearer_hyperrect_coord = rect->min + dir;
+		farther_hyperrect_coord = rect->max + dir;
+	}
+
+	if (nearer_subtree) {
+		/* Slice the hyperrect to get the hyperrect of the nearer subtree */
+		dummy = *nearer_hyperrect_coord;
+		*nearer_hyperrect_coord = node->pos[dir];
+		/* Recurse down into nearer subtree */
+		kd_nearest_i(nearer_subtree, pos, result, result_dist_sq, rect);
+		/* Undo the slice */
+		*nearer_hyperrect_coord = dummy;
+	}
+
+	/* Check the distance of the point at the current node, compare it
+	 * with our best so far */
+	dist_sq = 0;
+	for(i=0; i < rect->dim; i++) {
+		dist_sq += SQ(node->pos[i] - pos[i]);
+	}
+	if (dist_sq < *result_dist_sq) {
+		*result = node;
+		*result_dist_sq = dist_sq;
+	}
+
+	if (farther_subtree) {
+		/* Get the hyperrect of the farther subtree */
+		dummy = *farther_hyperrect_coord;
+		*farther_hyperrect_coord = node->pos[dir];
+		/* Check if we have to recurse down by calculating the closest
+		 * point of the hyperrect and see if it's closer than our
+		 * minimum distance in result_dist_sq. */
+		if (hyperrect_dist_sq(rect, pos) < *result_dist_sq) {
+			/* Recurse down into farther subtree */
+			kd_nearest_i(farther_subtree, pos, result, result_dist_sq, rect);
+		}
+		/* Undo the slice on the hyperrect */
+		*farther_hyperrect_coord = dummy;
+	}
+}
+
+struct kdres *kd_nearest(struct kdtree *kd, const double *pos)
+{
+	struct kdhyperrect *rect;
+	struct kdnode *result;
+	struct kdres *rset;
+	double dist_sq;
+	int i;
+
+	if (!kd) return 0;
+	if (!kd->rect) return 0;
+
+	/* Allocate result set */
+	if(!(rset = malloc(sizeof *rset))) {
+		return 0;
+	}
+	if(!(rset->rlist = alloc_resnode())) {
+		free(rset);
+		return 0;
+	}
+	rset->rlist->next = 0;
+	rset->tree = kd;
+
+	/* Duplicate the bounding hyperrectangle, we will work on the copy */
+	if (!(rect = hyperrect_duplicate(kd->rect))) {
+		kd_res_free(rset);
+		return 0;
+	}
+
+	/* Our first guesstimate is the root node */
+	result = kd->root;
+	dist_sq = 0;
+	for (i = 0; i < kd->dim; i++)
+		dist_sq += SQ(result->pos[i] - pos[i]);
+
+	/* Search for the nearest neighbour recursively */
+	kd_nearest_i(kd->root, pos, &result, &dist_sq, rect);
+
+	/* Free the copy of the hyperrect */
+	hyperrect_free(rect);
+
+	/* Store the result */
+	if (result) {
+		if (rlist_insert(rset->rlist, result, -1.0) == -1) {
+			kd_res_free(rset);
+			return 0;
+		}
+		rset->size = 1;
+		kd_res_rewind(rset);
+		return rset;
+	} else {
+		kd_res_free(rset);
+		return 0;
+	}
+}
+
+struct kdres *kd_nearestf(struct kdtree *tree, const float *pos)
+{
+	static double sbuf[16];
+	double *bptr, *buf = 0;
+	int dim = tree->dim;
+	struct kdres *res;
+
+	if(dim > 16) {
+#ifndef NO_ALLOCA
+		if(dim <= 256)
+			bptr = buf = alloca(dim * sizeof *bptr);
+		else
+#endif
+			if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
+				return 0;
+			}
+	} else {
+		bptr = buf = sbuf;
+	}
+
+	while(dim-- > 0) {
+		*bptr++ = *pos++;
+	}
+
+	res = kd_nearest(tree, buf);
+#ifndef NO_ALLOCA
+	if(tree->dim > 256)
+#else
+	if(tree->dim > 16)
+#endif
+		free(buf);
+	return res;
+}
+
+struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z)
+{
+	double pos[3];
+	pos[0] = x;
+	pos[1] = y;
+	pos[2] = z;
+	return kd_nearest(tree, pos);
+}
+
+struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z)
+{
+	double pos[3];
+	pos[0] = x;
+	pos[1] = y;
+	pos[2] = z;
+	return kd_nearest(tree, pos);
+}
+
+/* ---- nearest N search ---- */
+/*
+static kdres *kd_nearest_n(struct kdtree *kd, const double *pos, int num)
+{
+	int ret;
+	struct kdres *rset;
+
+	if(!(rset = malloc(sizeof *rset))) {
+		return 0;
+	}
+	if(!(rset->rlist = alloc_resnode())) {
+		free(rset);
+		return 0;
+	}
+	rset->rlist->next = 0;
+	rset->tree = kd;
+
+	if((ret = find_nearest_n(kd->root, pos, range, num, rset->rlist, kd->dim)) == -1) {
+		kd_res_free(rset);
+		return 0;
+	}
+	rset->size = ret;
+	kd_res_rewind(rset);
+	return rset;
+}*/
+
+struct kdres *kd_nearest_range(struct kdtree *kd, const double *pos, double range)
+{
+	int ret;
+	struct kdres *rset;
+
+	if(!(rset = malloc(sizeof *rset))) {
+		return 0;
+	}
+	if(!(rset->rlist = alloc_resnode())) {
+		free(rset);
+		return 0;
+	}
+	rset->rlist->next = 0;
+	rset->tree = kd;
+
+	if((ret = find_nearest(kd->root, pos, range, rset->rlist, 0, kd->dim)) == -1) {
+		kd_res_free(rset);
+		return 0;
+	}
+	rset->size = ret;
+	kd_res_rewind(rset);
+	return rset;
+}
+
+struct kdres *kd_nearest_rangef(struct kdtree *kd, const float *pos, float range)
+{
+	static double sbuf[16];
+	double *bptr, *buf = 0;
+	int dim = kd->dim;
+	struct kdres *res;
+
+	if(dim > 16) {
+#ifndef NO_ALLOCA
+		if(dim <= 256)
+			bptr = buf = alloca(dim * sizeof *bptr);
+		else
+#endif
+			if(!(bptr = buf = malloc(dim * sizeof *bptr))) {
+				return 0;
+			}
+	} else {
+		bptr = buf = sbuf;
+	}
+
+	while(dim-- > 0) {
+		*bptr++ = *pos++;
+	}
+
+	res = kd_nearest_range(kd, buf, range);
+#ifndef NO_ALLOCA
+	if(kd->dim > 256)
+#else
+	if(kd->dim > 16)
+#endif
+		free(buf);
+	return res;
+}
+
+struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double z, double range)
+{
+	double buf[3];
+	buf[0] = x;
+	buf[1] = y;
+	buf[2] = z;
+	return kd_nearest_range(tree, buf, range);
+}
+
+struct kdres *kd_nearest_range3f(struct kdtree *tree, float x, float y, float z, float range)
+{
+	double buf[3];
+	buf[0] = x;
+	buf[1] = y;
+	buf[2] = z;
+	return kd_nearest_range(tree, buf, range);
+}
+
+void kd_res_free(struct kdres *rset)
+{
+	clear_results(rset);
+	free_resnode(rset->rlist);
+	free(rset);
+}
+
+int kd_res_size(struct kdres *set)
+{
+	return (set->size);
+}
+
+void kd_res_rewind(struct kdres *rset)
+{
+	rset->riter = rset->rlist->next;
+}
+
+int kd_res_end(struct kdres *rset)
+{
+	return rset->riter == 0;
+}
+
+int kd_res_next(struct kdres *rset)
+{
+	rset->riter = rset->riter->next;
+	return rset->riter != 0;
+}
+
+void *kd_res_item(struct kdres *rset, double *pos)
+{
+	if(rset->riter) {
+		if(pos) {
+			memcpy(pos, rset->riter->item->pos, rset->tree->dim * sizeof *pos);
+		}
+		return rset->riter->item->data;
+	}
+	return 0;
+}
+
+void *kd_res_itemf(struct kdres *rset, float *pos)
+{
+	if(rset->riter) {
+		if(pos) {
+			int i;
+			for(i=0; i<rset->tree->dim; i++) {
+				pos[i] = rset->riter->item->pos[i];
+			}
+		}
+		return rset->riter->item->data;
+	}
+	return 0;
+}
+
+void *kd_res_item3(struct kdres *rset, double *x, double *y, double *z)
+{
+	if(rset->riter) {
+		if(*x) *x = rset->riter->item->pos[0];
+		if(*y) *y = rset->riter->item->pos[1];
+		if(*z) *z = rset->riter->item->pos[2];
+	}
+	return 0;
+}
+
+void *kd_res_item3f(struct kdres *rset, float *x, float *y, float *z)
+{
+	if(rset->riter) {
+		if(*x) *x = rset->riter->item->pos[0];
+		if(*y) *y = rset->riter->item->pos[1];
+		if(*z) *z = rset->riter->item->pos[2];
+	}
+	return 0;
+}
+
+void *kd_res_item_data(struct kdres *set)
+{
+	return kd_res_item(set, 0);
+}
+
+/* ---- hyperrectangle helpers ---- */
+static struct kdhyperrect* hyperrect_create(int dim, const double *min, const double *max)
+{
+	size_t size = dim * sizeof(double);
+	struct kdhyperrect* rect = 0;
+
+	if (!(rect = malloc(sizeof(struct kdhyperrect)))) {
+		return 0;
+	}
+
+	rect->dim = dim;
+	if (!(rect->min = malloc(size))) {
+		free(rect);
+		return 0;
+	}
+	if (!(rect->max = malloc(size))) {
+		free(rect->min);
+		free(rect);
+		return 0;
+	}
+	memcpy(rect->min, min, size);
+	memcpy(rect->max, max, size);
+
+	return rect;
+}
+
+static void hyperrect_free(struct kdhyperrect *rect)
+{
+	free(rect->min);
+	free(rect->max);
+	free(rect);
+}
+
+static struct kdhyperrect* hyperrect_duplicate(const struct kdhyperrect *rect)
+{
+	return hyperrect_create(rect->dim, rect->min, rect->max);
+}
+
+static void hyperrect_extend(struct kdhyperrect *rect, const double *pos)
+{
+	int i;
+
+	for (i=0; i < rect->dim; i++) {
+		if (pos[i] < rect->min[i]) {
+			rect->min[i] = pos[i];
+		}
+		if (pos[i] > rect->max[i]) {
+			rect->max[i] = pos[i];
+		}
+	}
+}
+
+static double hyperrect_dist_sq(struct kdhyperrect *rect, const double *pos)
+{
+	int i;
+	double result = 0;
+
+	for (i=0; i < rect->dim; i++) {
+		if (pos[i] < rect->min[i]) {
+			result += SQ(rect->min[i] - pos[i]);
+		} else if (pos[i] > rect->max[i]) {
+			result += SQ(rect->max[i] - pos[i]);
+		}
+	}
+
+	return result;
+}
+
+/* ---- static helpers ---- */
+
+#ifdef USE_LIST_NODE_ALLOCATOR
+/* special list node allocators. */
+static struct res_node *free_nodes;
+
+#ifndef NO_PTHREADS
+static pthread_mutex_t alloc_mutex = PTHREAD_MUTEX_INITIALIZER;
+#endif
+
+static struct res_node *alloc_resnode(void)
+{
+	struct res_node *node;
+
+#ifndef NO_PTHREADS
+	pthread_mutex_lock(&alloc_mutex);
+#endif
+
+	if(!free_nodes) {
+		node = malloc(sizeof *node);
+	} else {
+		node = free_nodes;
+		free_nodes = free_nodes->next;
+		node->next = 0;
+	}
+
+#ifndef NO_PTHREADS
+	pthread_mutex_unlock(&alloc_mutex);
+#endif
+
+	return node;
+}
+
+static void free_resnode(struct res_node *node)
+{
+#ifndef NO_PTHREADS
+	pthread_mutex_lock(&alloc_mutex);
+#endif
+
+	node->next = free_nodes;
+	free_nodes = node;
+
+#ifndef NO_PTHREADS
+	pthread_mutex_unlock(&alloc_mutex);
+#endif
+}
+#endif	/* list node allocator or not */
+
+
+/* inserts the item. if dist_sq is >= 0, then do an ordered insert */
+/* TODO make the ordering code use heapsort */
+static int rlist_insert(struct res_node *list, struct kdnode *item, double dist_sq)
+{
+	struct res_node *rnode;
+
+	if(!(rnode = alloc_resnode())) {
+		return -1;
+	}
+	rnode->item = item;
+	rnode->dist_sq = dist_sq;
+
+	if(dist_sq >= 0.0) {
+		while(list->next && list->next->dist_sq < dist_sq) {
+			list = list->next;
+		}
+	}
+	rnode->next = list->next;
+	list->next = rnode;
+	return 0;
+}
+
+static void clear_results(struct kdres *rset)
+{
+	struct res_node *tmp, *node = rset->rlist->next;
+
+	while(node) {
+		tmp = node;
+		node = node->next;
+		free_resnode(tmp);
+	}
+
+	rset->rlist->next = 0;
+}
diff --git a/intern/raskter/kdtree.h b/intern/raskter/kdtree.h
new file mode 100644
index 00000000000..da80e422a80
--- /dev/null
+++ b/intern/raskter/kdtree.h
@@ -0,0 +1,129 @@
+/*
+This file is part of ``kdtree'', a library for working with kd-trees.
+Copyright (C) 2007-2011 John Tsiombikas <nuclear@member.fsf.org>
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+3. The name of the author may not be used to endorse or promote products
+   derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
+WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
+OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+OF SUCH DAMAGE.
+*/
+#ifndef _KDTREE_H_
+#define _KDTREE_H_
+#define USE_LIST_NODE_ALLOCATOR
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct kdtree;
+struct kdres;
+
+
+/* create a kd-tree for "k"-dimensional data */
+struct kdtree *kd_create(int k);
+
+/* free the struct kdtree */
+void kd_free(struct kdtree *tree);
+
+/* remove all the elements from the tree */
+void kd_clear(struct kdtree *tree);
+
+/* if called with non-null 2nd argument, the function provided
+ * will be called on data pointers (see kd_insert) when nodes
+ * are to be removed from the tree.
+ */
+void kd_data_destructor(struct kdtree *tree, void (*destr)(void*));
+
+/* insert a node, specifying its position, and optional data */
+int kd_insert(struct kdtree *tree, const double *pos, void *data);
+int kd_insertf(struct kdtree *tree, const float *pos, void *data);
+int kd_insert3(struct kdtree *tree, double x, double y, double z, void *data);
+int kd_insert3f(struct kdtree *tree, float x, float y, float z, void *data);
+
+/* Find the nearest node from a given point.
+ *
+ * This function returns a pointer to a result set with at most one element.
+ */
+struct kdres *kd_nearest(struct kdtree *tree, const double *pos);
+struct kdres *kd_nearestf(struct kdtree *tree, const float *pos);
+struct kdres *kd_nearest3(struct kdtree *tree, double x, double y, double z);
+struct kdres *kd_nearest3f(struct kdtree *tree, float x, float y, float z);
+
+/* Find the N nearest nodes from a given point.
+ *
+ * This function returns a pointer to a result set, with at most N elements,
+ * which can be manipulated with the kd_res_* functions.
+ * The returned pointer can be null as an indication of an error. Otherwise
+ * a valid result set is always returned which may contain 0 or more elements.
+ * The result set must be deallocated with kd_res_free after use.
+ */
+/*
+struct kdres *kd_nearest_n(struct kdtree *tree, const double *pos, int num);
+struct kdres *kd_nearest_nf(struct kdtree *tree, const float *pos, int num);
+struct kdres *kd_nearest_n3(struct kdtree *tree, double x, double y, double z);
+struct kdres *kd_nearest_n3f(struct kdtree *tree, float x, float y, float z);
+*/
+
+/* Find any nearest nodes from a given point within a range.
+ *
+ * This function returns a pointer to a result set, which can be manipulated
+ * by the kd_res_* functions.
+ * The returned pointer can be null as an indication of an error. Otherwise
+ * a valid result set is always returned which may contain 0 or more elements.
+ * The result set must be deallocated with kd_res_free after use.
+ */
+struct kdres *kd_nearest_range(struct kdtree *tree, const double *pos, double range);
+struct kdres *kd_nearest_rangef(struct kdtree *tree, const float *pos, float range);
+struct kdres *kd_nearest_range3(struct kdtree *tree, double x, double y, double z, double range);
+struct kdres *kd_nearest_range3f(struct kdtree *tree, float x, float y, float z, float range);
+
+/* frees a result set returned by kd_nearest_range() */
+void kd_res_free(struct kdres *set);
+
+/* returns the size of the result set (in elements) */
+int kd_res_size(struct kdres *set);
+
+/* rewinds the result set iterator */
+void kd_res_rewind(struct kdres *set);
+
+/* returns non-zero if the set iterator reached the end after the last element */
+int kd_res_end(struct kdres *set);
+
+/* advances the result set iterator, returns non-zero on success, zero if
+ * there are no more elements in the result set.
+ */
+int kd_res_next(struct kdres *set);
+
+/* returns the data pointer (can be null) of the current result set item
+ * and optionally sets its position to the pointers(s) if not null.
+ */
+void *kd_res_item(struct kdres *set, double *pos);
+void *kd_res_itemf(struct kdres *set, float *pos);
+void *kd_res_item3(struct kdres *set, double *x, double *y, double *z);
+void *kd_res_item3f(struct kdres *set, float *x, float *y, float *z);
+
+/* equivalent to kd_res_item(set, 0) */
+void *kd_res_item_data(struct kdres *set);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif	/* _KDTREE_H_ */
diff --git a/intern/raskter/raskter.c b/intern/raskter/raskter.c
index b405fde82e8..825287ffe86 100644
--- a/intern/raskter/raskter.c
+++ b/intern/raskter/raskter.c
@@ -30,36 +30,13 @@
 
 #include <stdlib.h>
 #include "raskter.h"
+//#define __PLX__FAKE_AA__
+//#define __PLX_KD_TREE__
+#ifdef __PLX_KD_TREE__
+#include "kdtree.h"
+#endif
+
 
-#define __PLX__FAKE_AA__
-
-/* from BLI_utildefines.h */
-#define MIN2(x, y)      ( (x) < (y) ? (x) : (y) )
-#define MAX2(x, y)      ( (x) > (y) ? (x) : (y) )
-#define ABS(a)          ( (a) < 0 ? (-(a)) : (a) )
-
-struct e_status {
-	int x;
-	int ybeg;
-	int xshift;
-	int xdir;
-	int drift;
-	int drift_inc;
-	int drift_dec;
-	int num;
-	struct e_status *e_next;
-};
-
-struct r_buffer_stats {
-	float *buf;
-	int sizex;
-	int sizey;
-};
-
-struct r_fill_context {
-	struct e_status *all_edges, *possible_edges;
-	struct r_buffer_stats rb;
-};
 
 /*
  * Sort all the edges of the input polygon by Y, then by X, of the "first" vertex encountered.
@@ -69,102 +46,113 @@ struct r_fill_context {
  * just the poly. Since the DEM code could end up being coupled with this, we'll keep it separate
  * for now.
  */
-static void preprocess_all_edges(struct r_fill_context *ctx, struct poly_vert *verts, int num_verts, struct e_status *open_edge)
-{
-	int i;
-	int xbeg;
-	int ybeg;
-	int xend;
-	int yend;
-	int dx;
-	int dy;
-	int temp_pos;
-	int xdist;
-	struct e_status *e_new;
-	struct e_status *next_edge;
-	struct e_status **next_edge_ref;
-	struct poly_vert *v;
-	/* set up pointers */
-	v = verts;
-	ctx->all_edges = NULL;
-	/* loop all verts */
-	for (i = 0; i < num_verts; i++) {
-		/* determine beginnings and endings of edges, linking last vertex to first vertex */
-		xbeg = v[i].x;
-		ybeg = v[i].y;
-		if (i) {
-			/* we're not at the last vert, so end of the edge is the previous vertex */
-			xend = v[i - 1].x;
-			yend = v[i - 1].y;
-		}
-		else {
-			/* we're at the first vertex, so the "end" of this edge is the last vertex */
-			xend = v[num_verts - 1].x;
-			yend = v[num_verts - 1].y;
-		}
-		/* make sure our edges are facing the correct direction */
-		if (ybeg > yend) {
-			/* flip the Xs */
-			temp_pos = xbeg;
-			xbeg = xend;
-			xend = temp_pos;
-			/* flip the Ys */
-			temp_pos = ybeg;
-			ybeg = yend;
-			yend = temp_pos;
-		}
-
-		/* calculate y delta */
-		dy = yend - ybeg;
-		/* dont draw horizontal lines directly, they are scanned as part of the edges they connect, so skip em. :) */
-		if (dy) {
-			/* create the edge and determine it's slope (for incremental line drawing) */
-			e_new = open_edge++;
-
-			/* calculate x delta */
-			dx = xend - xbeg;
-			if (dx > 0) {
-				e_new->xdir = 1;
-				xdist = dx;
-			}
-			else {
-				e_new->xdir = -1;
-				xdist = -dx;
-			}
-
-			e_new->x = xbeg;
-			e_new->ybeg = ybeg;
-			e_new->num = dy;
-			e_new->drift_dec = dy;
-
-			/* calculate deltas for incremental drawing */
-			if (dx >= 0) {
-				e_new->drift = 0;
-			}
-			else {
-				e_new->drift = -dy + 1;
-			}
-			if (dy >= xdist) {
-				e_new->drift_inc = xdist;
-				e_new->xshift = 0;
-			}
-			else {
-				e_new->drift_inc = xdist % dy;
-				e_new->xshift = (xdist / dy) * e_new->xdir;
-			}
-			next_edge_ref = &ctx->all_edges;
-			/* link in all the edges, in sorted order */
-			for (;; ) {
-				next_edge = *next_edge_ref;
-				if (!next_edge || (next_edge->ybeg > ybeg) || ((next_edge->ybeg == ybeg) && (next_edge->x >= xbeg))) {
-					e_new->e_next = next_edge;
-					*next_edge_ref = e_new;
-					break;
-				}
-				next_edge_ref = &next_edge->e_next;
-			}
-		}
-	}
+static void preprocess_all_edges(struct r_fill_context *ctx, struct poly_vert *verts, int num_verts, struct e_status *open_edge) {
+    int i;
+    int xbeg;
+    int ybeg;
+    int xend;
+    int yend;
+    int dx;
+    int dy;
+    int temp_pos;
+    int xdist;
+    struct e_status *e_new;
+    struct e_status *next_edge;
+    struct e_status **next_edge_ref;
+    struct poly_vert *v;
+    /* set up pointers */
+    v = verts;
+    ctx->all_edges = NULL;
+    /* initialize some boundaries */
+    ctx->rb.xmax = v[0].x;
+    ctx->rb.xmin = v[0].x;
+    ctx->rb.ymax = v[0].y;
+    ctx->rb.ymin = v[0].y;
+    /* loop all verts */
+    for(i = 0; i < num_verts; i++) {
+        /* determine beginnings and endings of edges, linking last vertex to first vertex */
+        xbeg = v[i].x;
+        ybeg = v[i].y;
+        /* keep track of our x and y bounds */
+        if(xbeg >= ctx->rb.xmax) {
+            ctx->rb.xmax = xbeg;
+        } else if(xbeg <= ctx->rb.xmin) {
+            ctx->rb.xmin = xbeg;
+        }
+        if(ybeg >= ctx->rb.ymax) {
+            ctx->rb.ymax= ybeg;
+        } else if(ybeg <= ctx->rb.ymin) {
+            ctx->rb.ymin=ybeg;
+        }
+        if(i) {
+            /* we're not at the last vert, so end of the edge is the previous vertex */
+            xend = v[i - 1].x;
+            yend = v[i - 1].y;
+        } else {
+            /* we're at the first vertex, so the "end" of this edge is the last vertex */
+            xend = v[num_verts - 1].x;
+            yend = v[num_verts - 1].y;
+        }
+        /* make sure our edges are facing the correct direction */
+        if(ybeg > yend) {
+            /* flip the Xs */
+            temp_pos = xbeg;
+            xbeg = xend;
+            xend = temp_pos;
+            /* flip the Ys */
+            temp_pos = ybeg;
+            ybeg = yend;
+            yend = temp_pos;
+        }
+
+        /* calculate y delta */
+        dy = yend - ybeg;
+        /* dont draw horizontal lines directly, they are scanned as part of the edges they connect, so skip em. :) */
+        if(dy) {
+            /* create the edge and determine it's slope (for incremental line drawing) */
+            e_new = open_edge++;
+
+            /* calculate x delta */
+            dx = xend - xbeg;
+            if(dx > 0) {
+                e_new->xdir = 1;
+                xdist = dx;
+            } else {
+                e_new->xdir = -1;
+                xdist = -dx;
+            }
+
+            e_new->x = xbeg;
+            e_new->ybeg = ybeg;
+            e_new->num = dy;
+            e_new->drift_dec = dy;
+
+            /* calculate deltas for incremental drawing */
+            if(dx >= 0) {
+                e_new->drift = 0;
+            } else {
+                e_new->drift = -dy + 1;
+            }
+            if(dy >= xdist) {
+                e_new->drift_inc = xdist;
+                e_new->xshift = 0;
+            } else {
+                e_new->drift_inc = xdist % dy;
+                e_new->xshift = (xdist / dy) * e_new->xdir;
+            }
+            next_edge_ref = &ctx->all_edges;
+            /* link in all the edges, in sorted order */
+            for(;;) {
+                next_edge = *next_edge_ref;
+                if(!next_edge || (next_edge->ybeg > ybeg) || ((next_edge->ybeg == ybeg) && (next_edge->x >= xbeg))) {
+                    e_new->e_next = next_edge;
+                    *next_edge_ref = e_new;
+                    break;
+                }
+                next_edge_ref = &next_edge->e_next;
+            }
+        }
+    }
 }
 
 /*
@@ -172,281 +160,275 @@ static void preprocess_all_edges(struct r_fill_context *ctx, struct poly_vert *v
  * for speed, but waiting on final design choices for curve-data before eliminating data the DEM code will need
  * if it ends up being coupled with this function.
  */
-static int rast_scan_fill(struct r_fill_context *ctx, struct poly_vert *verts, int num_verts, float intensity)
-{
-	int x_curr;                 /* current pixel position in X */
-	int y_curr;                 /* current scan line being drawn */
-	int yp;                     /* y-pixel's position in frame buffer */
-	int swixd = 0;              /* whether or not edges switched position in X */
-	float *cpxl;                /* pixel pointers... */
-	float *mpxl;
-	float *spxl;
-	struct e_status *e_curr;    /* edge pointers... */
-	struct e_status *e_temp;
-	struct e_status *edgbuf;
-	struct e_status **edgec;
-
-
-	/*
-	 * If the number of verts specified to render as a polygon is less than 3,
-	 * return immediately. Obviously we cant render a poly with sides < 3. The
-	 * return for this we set to 1, simply so it can be distinguished from the
-	 * next place we could return.
-	 * which is a failure to allocate memory.
-	 */
-	if (num_verts < 3) {
-		return(1);
-	}
-
-	/*
-	 * Try to allocate an edge buffer in memory. needs to be the size of the edge tracking data
-	 * multiplied by the number of edges, which is always equal to the number of verts in
-	 * a 2D polygon. Here we return 0 to indicate a memory allocation failure, as opposed to a 1 for
-	 * the preceeding error, which was a rasterization request on a 2D poly with less than
-	 * 3 sides.
-	 */
-	if ((edgbuf = (struct e_status *)(malloc(sizeof(struct e_status) * num_verts))) == NULL) {
-		return(0);
-	}
-
-	/*
-	 * Do some preprocessing on all edges. This constructs a table structure in memory of all
-	 * the edge properties and can "flip" some edges so sorting works correctly.
-	 */
-	preprocess_all_edges(ctx, verts, num_verts, edgbuf);
-
-	/* can happen with a zero area mask */
-	if (ctx->all_edges == NULL) {
-		free(edgbuf);
-		return(1);
-	}
-
-	/*
-	 * Set the pointer for tracking the edges currently in processing to NULL to make sure
-	 * we don't get some crazy value after initialization.
-	 */
-	ctx->possible_edges = NULL;
-
-	/*
-	 * Loop through all scan lines to be drawn. Since we sorted by Y values during
-	 * preprocess_all_edges(), we can already exact values for the lowest and
-	 * highest Y values we could possibly need by induction. The preprocessing sorted
-	 * out edges by Y position, we can cycle the current edge being processed once
-	 * it runs out of Y pixels. When we have no more edges, meaning the current edge
-	 * is NULL after setting the "current" edge to be the previous current edge's
-	 * "next" edge in the Y sorted edge connection chain, we can stop looping Y values,
-	 * since we can't possibly have more scan lines if we ran out of edges. :)
-	 *
-	 * TODO: This clips Y to the frame buffer, which should be done in the preprocessor, but for now is done here.
-	 *       Will get changed once DEM code gets in.
-	 */
-	for (y_curr = ctx->all_edges->ybeg; (ctx->all_edges || ctx->possible_edges); y_curr++) {
-
-		/*
-		 * Link any edges that start on the current scan line into the list of
-		 * edges currently needed to draw at least this, if not several, scan lines.
-		 */
-
-		/*
-		 * Set the current edge to the beginning of the list of edges to be rasterized
-		 * into this scan line.
-		 *
-		 * We could have lots of edge here, so iterate over all the edges needed. The
-		 * preprocess_all_edges() function sorted edges by X within each chunk of Y sorting
-		 * so we safely cycle edges to thier own "next" edges in order.
-		 *
-		 * At each iteration, make sure we still have a non-NULL edge.
-		 */
-		for (edgec = &ctx->possible_edges; ctx->all_edges && (ctx->all_edges->ybeg == y_curr); ) {
-			x_curr = ctx->all_edges->x;                  /* Set current X position. */
-			for (;; ) {                                  /* Start looping edges. Will break when edges run out. */
-				e_curr = *edgec;                         /* Set up a current edge pointer. */
-				if (!e_curr || (e_curr->x >= x_curr)) {  /* If we have an no edge, or we need to skip some X-span, */
-					e_temp = ctx->all_edges->e_next;     /* set a temp "next" edge to test. */
-					*edgec = ctx->all_edges;             /* Add this edge to the list to be scanned. */
-					ctx->all_edges->e_next = e_curr;     /* Set up the next edge. */
-					edgec = &ctx->all_edges->e_next;     /* Set our list to the next edge's location in memory. */
-					ctx->all_edges = e_temp;             /* Skip the NULL or bad X edge, set pointer to next edge. */
-					break;                               /* Stop looping edges (since we ran out or hit empty X span. */
-				}
-				else {
-					edgec = &e_curr->e_next;             /* Set the pointer to the edge list the "next" edge. */
-				}
-			}
-		}
-
-		/*
-		 * Determine the current scan line's offset in the pixel buffer based on its Y position.
-		 * Basically we just multiply the current scan line's Y value by the number of pixels in each line.
-		 */
-		yp = y_curr * ctx->rb.sizex;
-		/*
-		 * Set a "scan line pointer" in memory. The location of the buffer plus the row offset.
-		 */
-		spxl = ctx->rb.buf + (yp);
-		/*
-		 * Set up the current edge to the first (in X) edge. The edges which could possibly be in this
-		 * list were determined in the preceeding edge loop above. They were already sorted in X by the
-		 * initial processing function.
-		 *
-		 * At each iteration, test for a NULL edge. Since we'll keep cycling edge's to their own "next" edge
-		 * we will eventually hit a NULL when the list runs out.
-		 */
-		for (e_curr = ctx->possible_edges; e_curr; e_curr = e_curr->e_next) {
-			/*
-			 * Calculate a span of pixels to fill on the current scan line.
-			 *
-			 * Set the current pixel pointer by adding the X offset to the scan line's start offset.
-			 * Cycle the current edge the next edge.
-			 * Set the max X value to draw to be one less than the next edge's first pixel. This way we are
-			 * sure not to ever get into a situation where we have overdraw. (drawing the same pixel more than
-			 * one time because it's on a vertex connecting two edges)
-			 *
-			 * Then blast through all the pixels in the span, advancing the pointer and setting the color to white.
-			 *
-			 * TODO: Here we clip to the scan line, this is not efficient, and should be done in the preprocessor,
-			 *       but for now it is done here until the DEM code comes in.
-			 */
-
-			/* set up xmin and xmax bounds on this scan line */
-			cpxl = spxl + MAX2(e_curr->x, 0);
-			e_curr = e_curr->e_next;
-			mpxl = spxl + MIN2(e_curr->x, ctx->rb.sizex) - 1;
-
-			if ((y_curr >= 0) && (y_curr < ctx->rb.sizey)) {
-				/* draw the pixels. */
-				for(; cpxl <= mpxl; *cpxl++ += intensity);
-			}
-		}
-
-		/*
-		 * Loop through all edges of polygon that could be hit by this scan line,
-		 * and figure out their x-intersections with the next scan line.
-		 *
-		 * Either A.) we wont have any more edges to test, or B.) we just add on the
-		 * slope delta computed in preprocessing step. Since this draws non-antialiased
-		 * polygons, we dont have fractional positions, so we only move in x-direction
-		 * when needed to get all the way to the next pixel over...
-		 */
-		for (edgec = &ctx->possible_edges; (e_curr = *edgec); ) {
-			if (!(--(e_curr->num))) {
-				*edgec = e_curr->e_next;
-			}
-			else {
-				e_curr->x += e_curr->xshift;
-				if ((e_curr->drift += e_curr->drift_inc) > 0) {
-					e_curr->x += e_curr->xdir;
-					e_curr->drift -= e_curr->drift_dec;
-				}
-				edgec = &e_curr->e_next;
-			}
-		}
-		/*
-		 * It's possible that some edges may have crossed during the last step, so we'll be sure
-		 * that we ALWAYS intersect scan lines in order by shuffling if needed to make all edges
-		 * sorted by x-intersection coordinate. We'll always scan through at least once to see if
-		 * edges crossed, and if so, we set the 'swixd' flag. If 'swixd' gets set on the initial
-		 * pass, then we know we need to sort by x, so then cycle through edges again and perform
-		 * the sort.-
-		 */
-		if (ctx->possible_edges) {
-			for (edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
-				/* if the current edge hits scan line at greater X than the next edge, we need to exchange the edges */
-				if (e_curr->x > e_curr->e_next->x) {
-					*edgec = e_curr->e_next;
-					/* exchange the pointers */
-					e_temp = e_curr->e_next->e_next;
-					e_curr->e_next->e_next = e_curr;
-					e_curr->e_next = e_temp;
-					/* set flag that we had at least one switch */
-					swixd = 1;
-				}
-			}
-			/* if we did have a switch, look for more (there will more if there was one) */
-			for (;; ) {
-				/* reset exchange flag so it's only set if we encounter another one */
-				swixd = 0;
-				for (edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
-					/* again, if current edge hits scan line at higher X than next edge, exchange the edges and set flag */
-					if (e_curr->x > e_curr->e_next->x) {
-						*edgec = e_curr->e_next;
-						/* exchange the pointers */
-						e_temp = e_curr->e_next->e_next;
-						e_curr->e_next->e_next = e_curr;
-						e_curr->e_next = e_temp;
-						/* flip the exchanged flag */
-						swixd = 1;
-					}
-				}
-				/* if we had no exchanges, we're done reshuffling the pointers */
-				if (!swixd) {
-					break;
-				}
-			}
-		}
-	}
-
-	free(edgbuf);
-	return 1;
+static int rast_scan_fill(struct r_fill_context *ctx, struct poly_vert *verts, int num_verts, float intensity) {
+    int x_curr;                 /* current pixel position in X */
+    int y_curr;                 /* current scan line being drawn */
+    int yp;                     /* y-pixel's position in frame buffer */
+    int swixd = 0;              /* whether or not edges switched position in X */
+    float *cpxl;                /* pixel pointers... */
+    float *mpxl;
+    float *spxl;
+    struct e_status *e_curr;    /* edge pointers... */
+    struct e_status *e_temp;
+    struct e_status *edgbuf;
+    struct e_status **edgec;
+
+
+    /*
+     * If the number of verts specified to render as a polygon is less than 3,
+     * return immediately. Obviously we cant render a poly with sides < 3. The
+     * return for this we set to 1, simply so it can be distinguished from the
+     * next place we could return, /home/guest/blender-svn/soc-2011-tomato/intern/raskter/raskter.
+     * which is a failure to allocate memory.
+     */
+    if(num_verts < 3) {
+        return(1);
+    }
+
+    /*
+     * Try to allocate an edge buffer in memory. needs to be the size of the edge tracking data
+     * multiplied by the number of edges, which is always equal to the number of verts in
+     * a 2D polygon. Here we return 0 to indicate a memory allocation failure, as opposed to a 1 for
+     * the preceeding error, which was a rasterization request on a 2D poly with less than
+     * 3 sides.
+     */
+    if((edgbuf = (struct e_status *)(malloc(sizeof(struct e_status) * num_verts))) == NULL) {
+        return(0);
+    }
+
+    /*
+     * Do some preprocessing on all edges. This constructs a table structure in memory of all
+     * the edge properties and can "flip" some edges so sorting works correctly.
+     */
+    preprocess_all_edges(ctx, verts, num_verts, edgbuf);
+
+    /* can happen with a zero area mask */
+    if (ctx->all_edges == NULL) {
+        free(edgbuf);
+        return(1);
+    }
+    /*
+     * Set the pointer for tracking the edges currently in processing to NULL to make sure
+     * we don't get some crazy value after initialization.
+     */
+    ctx->possible_edges = NULL;
+
+    /*
+     * Loop through all scan lines to be drawn. Since we sorted by Y values during
+     * preprocess_all_edges(), we can already exact values for the lowest and
+     * highest Y values we could possibly need by induction. The preprocessing sorted
+     * out edges by Y position, we can cycle the current edge being processed once
+     * it runs out of Y pixels. When we have no more edges, meaning the current edge
+     * is NULL after setting the "current" edge to be the previous current edge's
+     * "next" edge in the Y sorted edge connection chain, we can stop looping Y values,
+     * since we can't possibly have more scan lines if we ran out of edges. :)
+     *
+     * TODO: This clips Y to the frame buffer, which should be done in the preprocessor, but for now is done here.
+     *       Will get changed once DEM code gets in.
+     */
+    for(y_curr = ctx->all_edges->ybeg; (ctx->all_edges || ctx->possible_edges); y_curr++) {
+
+        /*
+         * Link any edges that start on the current scan line into the list of
+         * edges currently needed to draw at least this, if not several, scan lines.
+         */
+
+        /*
+         * Set the current edge to the beginning of the list of edges to be rasterized
+         * into this scan line.
+         *
+         * We could have lots of edge here, so iterate over all the edges needed. The
+         * preprocess_all_edges() function sorted edges by X within each chunk of Y sorting
+         * so we safely cycle edges to thier own "next" edges in order.
+         *
+         * At each iteration, make sure we still have a non-NULL edge.
+         */
+        for(edgec = &ctx->possible_edges; ctx->all_edges && (ctx->all_edges->ybeg == y_curr);) {
+            x_curr = ctx->all_edges->x;                  /* Set current X position. */
+            for(;;) {                                    /* Start looping edges. Will break when edges run out. */
+                e_curr = *edgec;                         /* Set up a current edge pointer. */
+                if(!e_curr || (e_curr->x >= x_curr)) {   /* If we have an no edge, or we need to skip some X-span, */
+                    e_temp = ctx->all_edges->e_next;     /* set a temp "next" edge to test. */
+                    *edgec = ctx->all_edges;             /* Add this edge to the list to be scanned. */
+                    ctx->all_edges->e_next = e_curr;     /* Set up the next edge. */
+                    edgec = &ctx->all_edges->e_next;     /* Set our list to the next edge's location in memory. */
+                    ctx->all_edges = e_temp;             /* Skip the NULL or bad X edge, set pointer to next edge. */
+                    break;                               /* Stop looping edges (since we ran out or hit empty X span. */
+                } else {
+                    edgec = &e_curr->e_next;             /* Set the pointer to the edge list the "next" edge. */
+                }
+            }
+        }
+
+        /*
+         * Determine the current scan line's offset in the pixel buffer based on its Y position.
+         * Basically we just multiply the current scan line's Y value by the number of pixels in each line.
+         */
+        yp = y_curr * ctx->rb.sizex;
+        /*
+         * Set a "scan line pointer" in memory. The location of the buffer plus the row offset.
+         */
+        spxl = ctx->rb.buf + (yp);
+        /*
+         * Set up the current edge to the first (in X) edge. The edges which could possibly be in this
+         * list were determined in the preceeding edge loop above. They were already sorted in X by the
+         * initial processing function.
+         *
+         * At each iteration, test for a NULL edge. Since we'll keep cycling edge's to their own "next" edge
+         * we will eventually hit a NULL when the list runs out.
+         */
+        for(e_curr = ctx->possible_edges; e_curr; e_curr = e_curr->e_next) {
+            /*
+             * Calculate a span of pixels to fill on the current scan line.
+             *
+             * Set the current pixel pointer by adding the X offset to the scan line's start offset.
+             * Cycle the current edge the next edge.
+             * Set the max X value to draw to be one less than the next edge's first pixel. This way we are
+             * sure not to ever get into a situation where we have overdraw. (drawing the same pixel more than
+             * one time because it's on a vertex connecting two edges)
+             *
+             * Then blast through all the pixels in the span, advancing the pointer and setting the color to white.
+             *
+             * TODO: Here we clip to the scan line, this is not efficient, and should be done in the preprocessor,
+             *       but for now it is done here until the DEM code comes in.
+             */
+
+            /* set up xmin and xmax bounds on this scan line */
+            cpxl = spxl + MAX2(e_curr->x, 0);
+            e_curr = e_curr->e_next;
+            mpxl = spxl + MIN2(e_curr->x, ctx->rb.sizex) - 1;
+
+            if((y_curr >= 0) && (y_curr < ctx->rb.sizey)) {
+                /* draw the pixels. */
+                for(; cpxl <= mpxl; *cpxl++ += intensity);
+            }
+        }
+
+        /*
+         * Loop through all edges of polygon that could be hit by this scan line,
+         * and figure out their x-intersections with the next scan line.
+         *
+         * Either A.) we wont have any more edges to test, or B.) we just add on the
+         * slope delta computed in preprocessing step. Since this draws non-antialiased
+         * polygons, we dont have fractional positions, so we only move in x-direction
+         * when needed to get all the way to the next pixel over...
+         */
+        for(edgec = &ctx->possible_edges; (e_curr = *edgec);) {
+            if(!(--(e_curr->num))) {
+                *edgec = e_curr->e_next;
+            } else {
+                e_curr->x += e_curr->xshift;
+                if((e_curr->drift += e_curr->drift_inc) > 0) {
+                    e_curr->x += e_curr->xdir;
+                    e_curr->drift -= e_curr->drift_dec;
+                }
+                edgec = &e_curr->e_next;
+            }
+        }
+        /*
+         * It's possible that some edges may have crossed during the last step, so we'll be sure
+         * that we ALWAYS intersect scan lines in order by shuffling if needed to make all edges
+         * sorted by x-intersection coordinate. We'll always scan through at least once to see if
+         * edges crossed, and if so, we set the 'swixd' flag. If 'swixd' gets set on the initial
+         * pass, then we know we need to sort by x, so then cycle through edges again and perform
+         * the sort.-
+         */
+        if(ctx->possible_edges) {
+            for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                /* if the current edge hits scan line at greater X than the next edge, we need to exchange the edges */
+                if(e_curr->x > e_curr->e_next->x) {
+                    *edgec = e_curr->e_next;
+                    /* exchange the pointers */
+                    e_temp = e_curr->e_next->e_next;
+                    e_curr->e_next->e_next = e_curr;
+                    e_curr->e_next = e_temp;
+                    /* set flag that we had at least one switch */
+                    swixd = 1;
+                }
+            }
+            /* if we did have a switch, look for more (there will more if there was one) */
+            for(;;) {
+                /* reset exchange flag so it's only set if we encounter another one */
+                swixd = 0;
+                for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                    /* again, if current edge hits scan line at higher X than next edge, exchange the edges and set flag */
+                    if(e_curr->x > e_curr->e_next->x) {
+                        *edgec = e_curr->e_next;
+                        /* exchange the pointers */
+                        e_temp = e_curr->e_next->e_next;
+                        e_curr->e_next->e_next = e_curr;
+                        e_curr->e_next = e_temp;
+                        /* flip the exchanged flag */
+                        swixd = 1;
+                    }
+                }
+                /* if we had no exchanges, we're done reshuffling the pointers */
+                if(!swixd) {
+                    break;
+                }
+            }
+        }
+    }
+
+    free(edgbuf);
+    return 1;
 }
 
-int PLX_raskterize(float (*base_verts)[2], int num_base_verts,
-                   float *buf, int buf_x, int buf_y, int do_mask_AA)
-{
-	int subdiv_AA = (do_mask_AA != 0) ? 8 : 0;
-	int i;                                   /* i: Loop counter. */
-	int sAx;
-	int sAy;
-	struct poly_vert *ply;                   /* ply: Pointer to a list of integer buffer-space vertex coordinates. */
-	struct r_fill_context ctx = {0};
-	const float buf_x_f = (float)(buf_x);
-	const float buf_y_f = (float)(buf_y);
-	float div_offset = (1.0f / (float)(subdiv_AA));
-	float div_offset_static = 0.5f * (float)(subdiv_AA) * div_offset;
-	/*
-	 * Allocate enough memory for our poly_vert list. It'll be the size of the poly_vert
-	 * data structure multiplied by the number of base_verts.
-	 *
-	 * In the event of a failure to allocate the memory, return 0, so this error can
-	 * be distinguished as a memory allocation error.
-	 */
-	if ((ply = (struct poly_vert *)(malloc(sizeof(struct poly_vert) * num_base_verts))) == NULL) {
-		return(0);
-	}
-
-	ctx.rb.buf = buf;                            /* Set the output buffer pointer. */
-	ctx.rb.sizex = buf_x;                        /* Set the output buffer size in X. (width) */
-	ctx.rb.sizey = buf_y;                        /* Set the output buffer size in Y. (height) */
-	/*
-	 * Loop over all verts passed in to be rasterized. Each vertex's X and Y coordinates are
-	 * then converted from normalized screen space (0.0 <= POS <= 1.0) to integer coordinates
-	 * in the buffer-space coordinates passed in inside buf_x and buf_y.
-	 *
-	 * It's worth noting that this function ONLY outputs fully white pixels in a mask. Every pixel
-	 * drawn will be 1.0f in value, there is no anti-aliasing.
-	 */
-
-	if (!subdiv_AA) {
-		for (i = 0; i < num_base_verts; i++) {                     /* Loop over all base_verts. */
-			ply[i].x = (int)((base_verts[i][0] * buf_x_f) + 0.5f); /* Range expand normalized X to integer buffer-space X. */
-			ply[i].y = (int)((base_verts[i][1] * buf_y_f) + 0.5f); /* Range expand normalized Y to integer buffer-space Y. */
-		}
-
-		i = rast_scan_fill(&ctx, ply, num_base_verts, 1.0f);  /* Call our rasterizer, passing in the integer coords for each vert. */
-	}
-	else {
-		for (sAx = 0; sAx < subdiv_AA; sAx++) {
-			for (sAy = 0; sAy < subdiv_AA; sAy++) {
-				for (i = 0; i < num_base_verts; i++) {
-					ply[i].x = (int)((base_verts[i][0] * buf_x_f) + 0.5f - div_offset_static + (div_offset * (float)(sAx)));
-					ply[i].y = (int)((base_verts[i][1] * buf_y_f) + 0.5f - div_offset_static + (div_offset * (float)(sAy)));
-				}
-				i = rast_scan_fill(&ctx, ply, num_base_verts, (1.0f / (float)(subdiv_AA * subdiv_AA)));
-			}
-		}
-	}
-	free(ply);                                      /* Free the memory allocated for the integer coordinate table. */
-	return(i);                                      /* Return the value returned by the rasterizer. */
+int PLX_raskterize(float(*base_verts)[2], int num_base_verts,
+                   float *buf, int buf_x, int buf_y, int do_mask_AA) {
+    int subdiv_AA = (do_mask_AA != 0)? 0:0;
+    int i;                                   /* i: Loop counter. */
+    int sAx;
+    int sAy;
+    struct poly_vert *ply;                   /* ply: Pointer to a list of integer buffer-space vertex coordinates. */
+    struct r_fill_context ctx = {0};
+    const float buf_x_f = (float)(buf_x);
+    const float buf_y_f = (float)(buf_y);
+    float div_offset=(1.0f / (float)(subdiv_AA));
+    float div_offset_static = 0.5f * (float)(subdiv_AA) * div_offset;
+    /*
+     * Allocate enough memory for our poly_vert list. It'll be the size of the poly_vert
+     * data structure multiplied by the number of base_verts.
+     *
+     * In the event of a failure to allocate the memory, return 0, so this error can
+     * be distinguished as a memory allocation error.
+     */
+    if((ply = (struct poly_vert *)(malloc(sizeof(struct poly_vert) * num_base_verts))) == NULL) {
+        return(0);
+    }
+
+    ctx.rb.buf = buf;                            /* Set the output buffer pointer. */
+    ctx.rb.sizex = buf_x;                        /* Set the output buffer size in X. (width) */
+    ctx.rb.sizey = buf_y;                        /* Set the output buffer size in Y. (height) */
+    /*
+     * Loop over all verts passed in to be rasterized. Each vertex's X and Y coordinates are
+     * then converted from normalized screen space (0.0 <= POS <= 1.0) to integer coordinates
+     * in the buffer-space coordinates passed in inside buf_x and buf_y.
+     *
+     * It's worth noting that this function ONLY outputs fully white pixels in a mask. Every pixel
+     * drawn will be 1.0f in value, there is no anti-aliasing.
+     */
+
+    if(!subdiv_AA) {
+        for(i = 0; i < num_base_verts; i++) {                           /* Loop over all base_verts. */
+            ply[i].x = (int)((base_verts[i][0] * buf_x_f) + 0.5f);       /* Range expand normalized X to integer buffer-space X. */
+            ply[i].y = (int)((base_verts[i][1] * buf_y_f) + 0.5f); /* Range expand normalized Y to integer buffer-space Y. */
+        }
+
+        i = rast_scan_fill(&ctx, ply, num_base_verts,1.0f);  /* Call our rasterizer, passing in the integer coords for each vert. */
+    } else {
+        for(sAx=0; sAx < subdiv_AA; sAx++) {
+            for(sAy=0; sAy < subdiv_AA; sAy++) {
+                for(i=0; i < num_base_verts; i++) {
+                    ply[i].x = (int)((base_verts[i][0]*buf_x_f)+0.5f - div_offset_static + (div_offset*(float)(sAx)));
+                    ply[i].y = (int)((base_verts[i][1]*buf_y_f)+0.5f - div_offset_static + (div_offset*(float)(sAy)));
+                }
+                i = rast_scan_fill(&ctx, ply, num_base_verts,(1.0f / (float)(subdiv_AA*subdiv_AA)));
+            }
+        }
+    }
+    free(ply);                                      /* Free the memory allocated for the integer coordinate table. */
+    return(i);                                      /* Return the value returned by the rasterizer. */
 }
 
 /*
@@ -455,911 +437,1087 @@ int PLX_raskterize(float (*base_verts)[2], int num_base_verts,
  * if it ends up being coupled with this function.
  */
 static int rast_scan_feather(struct r_fill_context *ctx,
-                             float (*base_verts_f)[2], int num_base_verts,
-                             struct poly_vert *feather_verts, float(*feather_verts_f)[2], int num_feather_verts)
-{
-	int x_curr;                 /* current pixel position in X */
-	int y_curr;                 /* current scan line being drawn */
-	int yp;                     /* y-pixel's position in frame buffer */
-	int swixd = 0;              /* whether or not edges switched position in X */
-	float *cpxl;                /* pixel pointers... */
-	float *mpxl;
-	float *spxl;
-	struct e_status *e_curr;    /* edge pointers... */
-	struct e_status *e_temp;
-	struct e_status *edgbuf;
-	struct e_status **edgec;
-
-	/* from dem */
-	int a;                          // a = temporary pixel index buffer loop counter
-	float fsz;                        // size of the frame
-	unsigned int rsl;               // long used for finding fast 1.0/sqrt
-	float rsf;                      // float used for finding fast 1.0/sqrt
-	const float rsopf = 1.5f;       // constant float used for finding fast 1.0/sqrt
-
-	//unsigned int gradientFillOffset;
-	float t;
-	float ud;                // ud = unscaled edge distance
-	float dmin;              // dmin = minimun edge distance
-	float odist;                    // odist = current outer edge distance
-	float idist;                    // idist = current inner edge distance
-	float dx;                         // dx = X-delta (used for distance proportion calculation)
-	float dy;                         // dy = Y-delta (used for distance proportion calculation)
-	float xpxw = (1.0f / (float)(ctx->rb.sizex));  // xpxw = normalized pixel width
-	float ypxh = (1.0f / (float)(ctx->rb.sizey));  // ypxh = normalized pixel height
-
-	/*
-	 * If the number of verts specified to render as a polygon is less than 3,
-	 * return immediately. Obviously we cant render a poly with sides < 3. The
-	 * return for this we set to 1, simply so it can be distinguished from the
-	 * next place we could return,
-	 * which is a failure to allocate memory.
-	 */
-	if (num_feather_verts < 3) {
-		return(1);
-	}
-
-	/*
-	 * Try to allocate an edge buffer in memory. needs to be the size of the edge tracking data
-	 * multiplied by the number of edges, which is always equal to the number of verts in
-	 * a 2D polygon. Here we return 0 to indicate a memory allocation failure, as opposed to a 1 for
-	 * the preceeding error, which was a rasterization request on a 2D poly with less than
-	 * 3 sides.
-	 */
-	if ((edgbuf = (struct e_status *)(malloc(sizeof(struct e_status) * num_feather_verts))) == NULL) {
-		return(0);
-	}
-
-	/*
-	 * Do some preprocessing on all edges. This constructs a table structure in memory of all
-	 * the edge properties and can "flip" some edges so sorting works correctly.
-	 */
-	preprocess_all_edges(ctx, feather_verts, num_feather_verts, edgbuf);
-
-	/* can happen with a zero area mask */
-	if (ctx->all_edges == NULL) {
-		free(edgbuf);
-		return(1);
-	}
-
-	/*
-	 * Set the pointer for tracking the edges currently in processing to NULL to make sure
-	 * we don't get some crazy value after initialization.
-	 */
-	ctx->possible_edges = NULL;
-
-	/*
-	 * Loop through all scan lines to be drawn. Since we sorted by Y values during
-	 * preprocess_all_edges(), we can already exact values for the lowest and
-	 * highest Y values we could possibly need by induction. The preprocessing sorted
-	 * out edges by Y position, we can cycle the current edge being processed once
-	 * it runs out of Y pixels. When we have no more edges, meaning the current edge
-	 * is NULL after setting the "current" edge to be the previous current edge's
-	 * "next" edge in the Y sorted edge connection chain, we can stop looping Y values,
-	 * since we can't possibly have more scan lines if we ran out of edges. :)
-	 *
-	 * TODO: This clips Y to the frame buffer, which should be done in the preprocessor, but for now is done here.
-	 *       Will get changed once DEM code gets in.
-	 */
-	for (y_curr = ctx->all_edges->ybeg; (ctx->all_edges || ctx->possible_edges); y_curr++) {
-
-		/*
-		 * Link any edges that start on the current scan line into the list of
-		 * edges currently needed to draw at least this, if not several, scan lines.
-		 */
-
-		/*
-		 * Set the current edge to the beginning of the list of edges to be rasterized
-		 * into this scan line.
-		 *
-		 * We could have lots of edge here, so iterate over all the edges needed. The
-		 * preprocess_all_edges() function sorted edges by X within each chunk of Y sorting
-		 * so we safely cycle edges to thier own "next" edges in order.
-		 *
-		 * At each iteration, make sure we still have a non-NULL edge.
-		 */
-		for (edgec = &ctx->possible_edges; ctx->all_edges && (ctx->all_edges->ybeg == y_curr); ) {
-			x_curr = ctx->all_edges->x;                  /* Set current X position. */
-			for (;; ) {                                  /* Start looping edges. Will break when edges run out. */
-				e_curr = *edgec;                         /* Set up a current edge pointer. */
-				if (!e_curr || (e_curr->x >= x_curr)) {  /* If we have an no edge, or we need to skip some X-span, */
-					e_temp = ctx->all_edges->e_next;     /* set a temp "next" edge to test. */
-					*edgec = ctx->all_edges;             /* Add this edge to the list to be scanned. */
-					ctx->all_edges->e_next = e_curr;     /* Set up the next edge. */
-					edgec = &ctx->all_edges->e_next;     /* Set our list to the next edge's location in memory. */
-					ctx->all_edges = e_temp;             /* Skip the NULL or bad X edge, set pointer to next edge. */
-					break;                               /* Stop looping edges (since we ran out or hit empty X span. */
-				}
-				else {
-					edgec = &e_curr->e_next;             /* Set the pointer to the edge list the "next" edge. */
-				}
-			}
-		}
-
-		/*
-		 * Determine the current scan line's offset in the pixel buffer based on its Y position.
-		 * Basically we just multiply the current scan line's Y value by the number of pixels in each line.
-		 */
-		yp = y_curr * ctx->rb.sizex;
-		/*
-		 * Set a "scan line pointer" in memory. The location of the buffer plus the row offset.
-		 */
-		spxl = ctx->rb.buf + (yp);
-		/*
-		 * Set up the current edge to the first (in X) edge. The edges which could possibly be in this
-		 * list were determined in the preceeding edge loop above. They were already sorted in X by the
-		 * initial processing function.
-		 *
-		 * At each iteration, test for a NULL edge. Since we'll keep cycling edge's to their own "next" edge
-		 * we will eventually hit a NULL when the list runs out.
-		 */
-		for (e_curr = ctx->possible_edges; e_curr; e_curr = e_curr->e_next) {
-			/*
-			 * Calculate a span of pixels to fill on the current scan line.
-			 *
-			 * Set the current pixel pointer by adding the X offset to the scan line's start offset.
-			 * Cycle the current edge the next edge.
-			 * Set the max X value to draw to be one less than the next edge's first pixel. This way we are
-			 * sure not to ever get into a situation where we have overdraw. (drawing the same pixel more than
-			 * one time because it's on a vertex connecting two edges)
-			 *
-			 * Then blast through all the pixels in the span, advancing the pointer and setting the color to white.
-			 *
-			 * TODO: Here we clip to the scan line, this is not efficient, and should be done in the preprocessor,
-			 *       but for now it is done here until the DEM code comes in.
-			 */
-
-			/* set up xmin and xmax bounds on this scan line */
-			cpxl = spxl + MAX2(e_curr->x, 0);
-			e_curr = e_curr->e_next;
-			mpxl = spxl + MIN2(e_curr->x, ctx->rb.sizex) - 1;
-
-			if ((y_curr >= 0) && (y_curr < ctx->rb.sizey)) {
-				t = ((float)((cpxl - spxl) % ctx->rb.sizex) + 0.5f) * xpxw;
-				fsz = ((float)(y_curr) + 0.5f) * ypxh;
-				/* draw the pixels. */
-				for (; cpxl <= mpxl; cpxl++, t += xpxw) {
-					//do feather check
-					// first check that pixel isn't already full, and only operate if it is not
-					if (*cpxl < 0.9999f) {
-
-						dmin = 2.0f;                        // reset min distance to edge pixel
-						for (a = 0; a < num_feather_verts; a++) { // loop through all outer edge buffer pixels
-							dy = t - feather_verts_f[a][0];          // set dx to gradient pixel column - outer edge pixel row
-							dx = fsz - feather_verts_f[a][1];        // set dy to gradient pixel row - outer edge pixel column
-							ud = dx * dx + dy * dy;               // compute sum of squares
-							if (ud < dmin) {                      // if our new sum of squares is less than the current minimum
-								dmin = ud;                        // set a new minimum equal to the new lower value
-							}
-						}
-						odist = dmin;                    // cast outer min to a float
-						rsf = odist * 0.5f;                       //
-						rsl = *(unsigned int *)&odist;            // use some peculiar properties of the way bits are stored
-						rsl = 0x5f3759df - (rsl >> 1);            // in floats vs. unsigned ints to compute an approximate
-						odist = *(float *)&rsl;                   // reciprocal square root
-						odist = odist * (rsopf - (rsf * odist * odist));  // -- ** this line can be iterated for more accuracy ** --
-						odist = odist * (rsopf - (rsf * odist * odist));
-						dmin = 2.0f;                        // reset min distance to edge pixel
-						for (a = 0; a < num_base_verts; a++) {    // loop through all inside edge pixels
-							dy = t - base_verts_f[a][0];             // compute delta in Y from gradient pixel to inside edge pixel
-							dx = fsz - base_verts_f[a][1];           // compute delta in X from gradient pixel to inside edge pixel
-							ud = dx * dx + dy * dy;   // compute sum of squares
-							if (ud < dmin) {          // if our new sum of squares is less than the current minimum we've found
-								dmin = ud;            // set a new minimum equal to the new lower value
-							}
-						}
-						idist = dmin;                    // cast inner min to a float
-						rsf = idist * 0.5f;                       //
-						rsl = *(unsigned int *)&idist;            //
-						rsl = 0x5f3759df - (rsl >> 1);            // see notes above
-						idist = *(float *)&rsl;                   //
-						idist = idist * (rsopf - (rsf * idist * idist));  //
-						idist = idist * (rsopf - (rsf * idist * idist));
-						/*
-						 * Note once again that since we are using reciprocals of distance values our
-						 * proportion is already the correct intensity, and does not need to be
-						 * subracted from 1.0 like it would have if we used real distances.
-						 */
-
-						/* set intensity, do the += so overlapping gradients are additive */
-						*cpxl = (idist / (idist + odist));
-					}
-				}
-			}
-		}
-
-		/*
-		 * Loop through all edges of polygon that could be hit by this scan line,
-		 * and figure out their x-intersections with the next scan line.
-		 *
-		 * Either A.) we wont have any more edges to test, or B.) we just add on the
-		 * slope delta computed in preprocessing step. Since this draws non-antialiased
-		 * polygons, we dont have fractional positions, so we only move in x-direction
-		 * when needed to get all the way to the next pixel over...
-		 */
-		for (edgec = &ctx->possible_edges; (e_curr = *edgec); ) {
-			if (!(--(e_curr->num))) {
-				*edgec = e_curr->e_next;
-			}
-			else {
-				e_curr->x += e_curr->xshift;
-				if ((e_curr->drift += e_curr->drift_inc) > 0) {
-					e_curr->x += e_curr->xdir;
-					e_curr->drift -= e_curr->drift_dec;
-				}
-				edgec = &e_curr->e_next;
-			}
-		}
-		/*
-		 * It's possible that some edges may have crossed during the last step, so we'll be sure
-		 * that we ALWAYS intersect scan lines in order by shuffling if needed to make all edges
-		 * sorted by x-intersection coordinate. We'll always scan through at least once to see if
-		 * edges crossed, and if so, we set the 'swixd' flag. If 'swixd' gets set on the initial
-		 * pass, then we know we need to sort by x, so then cycle through edges again and perform
-		 * the sort.-
-		 */
-		if (ctx->possible_edges) {
-			for (edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
-				/* if the current edge hits scan line at greater X than the next edge, we need to exchange the edges */
-				if (e_curr->x > e_curr->e_next->x) {
-					*edgec = e_curr->e_next;
-					/* exchange the pointers */
-					e_temp = e_curr->e_next->e_next;
-					e_curr->e_next->e_next = e_curr;
-					e_curr->e_next = e_temp;
-					/* set flag that we had at least one switch */
-					swixd = 1;
-				}
-			}
-			/* if we did have a switch, look for more (there will more if there was one) */
-			for (;; ) {
-				/* reset exchange flag so it's only set if we encounter another one */
-				swixd = 0;
-				for (edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
-					/* again, if current edge hits scan line at higher X than next edge,
-					 * exchange the edges and set flag */
-					if (e_curr->x > e_curr->e_next->x) {
-						*edgec = e_curr->e_next;
-						/* exchange the pointers */
-						e_temp = e_curr->e_next->e_next;
-						e_curr->e_next->e_next = e_curr;
-						e_curr->e_next = e_temp;
-						/* flip the exchanged flag */
-						swixd = 1;
-					}
-				}
-				/* if we had no exchanges, we're done reshuffling the pointers */
-				if (!swixd) {
-					break;
-				}
-			}
-		}
-	}
-
-	free(edgbuf);
-	return 1;
+                             float(*base_verts_f)[2], int num_base_verts,
+                             struct poly_vert *feather_verts, float(*feather_verts_f)[2], int num_feather_verts) {
+    int x_curr;                 /* current pixel position in X */
+    int y_curr;                 /* current scan line being drawn */
+    int yp;                     /* y-pixel's position in frame buffer */
+    int swixd = 0;              /* whether or not edges switched position in X */
+    float *cpxl;                /* pixel pointers... */
+    float *mpxl;
+    float *spxl;
+    struct e_status *e_curr;    /* edge pointers... */
+    struct e_status *e_temp;
+    struct e_status *edgbuf;
+    struct e_status **edgec;
+
+    /* from dem */
+    int a;                          // a = temporary pixel index buffer loop counter
+    float fsz;                        // size of the frame
+    unsigned int rsl;               // long used for finding fast 1.0/sqrt
+    float rsf;                      // float used for finding fast 1.0/sqrt
+    const float rsopf = 1.5f;       // constant float used for finding fast 1.0/sqrt
+
+    //unsigned int gradientFillOffset;
+    float t;
+    float ud;                // ud = unscaled edge distance
+    float dmin;              // dmin = minimun edge distance
+    float odist;                    // odist = current outer edge distance
+    float idist;                    // idist = current inner edge distance
+    float dx;                         // dx = X-delta (used for distance proportion calculation)
+    float dy;                         // dy = Y-delta (used for distance proportion calculation)
+    float xpxw = (1.0f / (float)(ctx->rb.sizex));  // xpxw = normalized pixel width
+    float ypxh = (1.0f / (float)(ctx->rb.sizey));  // ypxh = normalized pixel height
+#ifdef __PLX_KD_TREE__
+    void *res_kdi;
+    void *res_kdo;
+    float clup[2];
+#endif
+
+    /*
+     * If the number of verts specified to render as a polygon is less than 3,
+     * return immediately. Obviously we cant render a poly with sides < 3. The
+     * return for this we set to 1, simply so it can be distinguished from the
+     * next place we could return,
+     * which is a failure to allocate memory.
+     */
+    if(num_feather_verts < 3) {
+        return(1);
+    }
+
+    /*
+     * Try to allocate an edge buffer in memory. needs to be the size of the edge tracking data
+     * multiplied by the number of edges, which is always equal to the number of verts in
+     * a 2D polygon. Here we return 0 to indicate a memory allocation failure, as opposed to a 1 for
+     * the preceeding error, which was a rasterization request on a 2D poly with less than
+     * 3 sides.
+     */
+    if((edgbuf = (struct e_status *)(malloc(sizeof(struct e_status) * num_feather_verts))) == NULL) {
+        return(0);
+    }
+
+    /*
+     * Do some preprocessing on all edges. This constructs a table structure in memory of all
+     * the edge properties and can "flip" some edges so sorting works correctly.
+     */
+    preprocess_all_edges(ctx, feather_verts, num_feather_verts, edgbuf);
+
+    /* can happen with a zero area mask */
+    if (ctx->all_edges == NULL) {
+        free(edgbuf);
+        return(1);
+    }
+
+    /*
+     * Set the pointer for tracking the edges currently in processing to NULL to make sure
+     * we don't get some crazy value after initialization.
+     */
+    ctx->possible_edges = NULL;
+
+    /*
+     * Loop through all scan lines to be drawn. Since we sorted by Y values during
+     * preprocess_all_edges(), we can already exact values for the lowest and
+     * highest Y values we could possibly need by induction. The preprocessing sorted
+     * out edges by Y position, we can cycle the current edge being processed once
+     * it runs out of Y pixels. When we have no more edges, meaning the current edge
+     * is NULL after setting the "current" edge to be the previous current edge's
+     * "next" edge in the Y sorted edge connection chain, we can stop looping Y values,
+     * since we can't possibly have more scan lines if we ran out of edges. :)
+     *
+     * TODO: This clips Y to the frame buffer, which should be done in the preprocessor, but for now is done here.
+     *       Will get changed once DEM code gets in.
+     */
+    for(y_curr = ctx->all_edges->ybeg; (ctx->all_edges || ctx->possible_edges); y_curr++) {
+
+        /*
+         * Link any edges that start on the current scan line into the list of
+         * edges currently needed to draw at least this, if not several, scan lines.
+         */
+
+        /*
+         * Set the current edge to the beginning of the list of edges to be rasterized
+         * into this scan line.
+         *
+         * We could have lots of edge here, so iterate over all the edges needed. The
+         * preprocess_all_edges() function sorted edges by X within each chunk of Y sorting
+         * so we safely cycle edges to thier own "next" edges in order.
+         *
+         * At each iteration, make sure we still have a non-NULL edge.
+         */
+        for(edgec = &ctx->possible_edges; ctx->all_edges && (ctx->all_edges->ybeg == y_curr);) {
+            x_curr = ctx->all_edges->x;                  /* Set current X position. */
+            for(;;) {                                    /* Start looping edges. Will break when edges run out. */
+                e_curr = *edgec;                         /* Set up a current edge pointer. */
+                if(!e_curr || (e_curr->x >= x_curr)) {   /* If we have an no edge, or we need to skip some X-span, */
+                    e_temp = ctx->all_edges->e_next;     /* set a temp "next" edge to test. */
+                    *edgec = ctx->all_edges;             /* Add this edge to the list to be scanned. */
+                    ctx->all_edges->e_next = e_curr;     /* Set up the next edge. */
+                    edgec = &ctx->all_edges->e_next;     /* Set our list to the next edge's location in memory. */
+                    ctx->all_edges = e_temp;             /* Skip the NULL or bad X edge, set pointer to next edge. */
+                    break;                               /* Stop looping edges (since we ran out or hit empty X span. */
+                } else {
+                    edgec = &e_curr->e_next;             /* Set the pointer to the edge list the "next" edge. */
+                }
+            }
+        }
+
+        /*
+         * Determine the current scan line's offset in the pixel buffer based on its Y position.
+         * Basically we just multiply the current scan line's Y value by the number of pixels in each line.
+         */
+        yp = y_curr * ctx->rb.sizex;
+        /*
+         * Set a "scan line pointer" in memory. The location of the buffer plus the row offset.
+         */
+        spxl = ctx->rb.buf + (yp);
+        /*
+         * Set up the current edge to the first (in X) edge. The edges which could possibly be in this
+         * list were determined in the preceeding edge loop above. They were already sorted in X by the
+         * initial processing function.
+         *
+         * At each iteration, test for a NULL edge. Since we'll keep cycling edge's to their own "next" edge
+         * we will eventually hit a NULL when the list runs out.
+         */
+        for(e_curr = ctx->possible_edges; e_curr; e_curr = e_curr->e_next) {
+            /*
+             * Calculate a span of pixels to fill on the current scan line.
+             *
+             * Set the current pixel pointer by adding the X offset to the scan line's start offset.
+             * Cycle the current edge the next edge.
+             * Set the max X value to draw to be one less than the next edge's first pixel. This way we are
+             * sure not to ever get into a situation where we have overdraw. (drawing the same pixel more than
+             * one time because it's on a vertex connecting two edges)
+             *
+             * Then blast through all the pixels in the span, advancing the pointer and setting the color to white.
+             *
+             * TODO: Here we clip to the scan line, this is not efficient, and should be done in the preprocessor,
+             *       but for now it is done here until the DEM code comes in.
+             */
+
+            /* set up xmin and xmax bounds on this scan line */
+            cpxl = spxl + MAX2(e_curr->x, 0);
+            e_curr = e_curr->e_next;
+            mpxl = spxl + MIN2(e_curr->x, ctx->rb.sizex) - 1;
+
+            if((y_curr >= 0) && (y_curr < ctx->rb.sizey)) {
+                t = ((float)((cpxl - spxl) % ctx->rb.sizex) + 0.5f) * xpxw;
+                fsz = ((float)(y_curr) + 0.5f) * ypxh;
+                /* draw the pixels. */
+                for(; cpxl <= mpxl; cpxl++, t += xpxw) {
+                    //do feather check
+                    // first check that pixel isn't already full, and only operate if it is not
+                    if(*cpxl < 0.9999f) {
+#ifndef __PLX_KD_TREE__
+                        dmin = 2.0f;                        // reset min distance to edge pixel
+                        for(a = 0; a < num_feather_verts; a++) {  // loop through all outer edge buffer pixels
+                            dx = t - feather_verts_f[a][0];          // set dx to gradient pixel column - outer edge pixel row
+                            dy = fsz - feather_verts_f[a][1];        // set dy to gradient pixel row - outer edge pixel column
+                            ud = dx * dx + dy * dy;               // compute sum of squares
+                            if(ud < dmin) {                       // if our new sum of squares is less than the current minimum
+                                dmin = ud;                        // set a new minimum equal to the new lower value
+                            }
+                        }
+                        odist = dmin;                    // cast outer min to a float
+                        rsf = odist * 0.5f;                       //
+                        rsl = *(unsigned int *)&odist;            // use some peculiar properties of the way bits are stored
+                        rsl = 0x5f3759df - (rsl >> 1);            // in floats vs. unsigned ints to compute an approximate
+                        odist = *(float *)&rsl;                   // reciprocal square root
+                        odist = odist * (rsopf - (rsf * odist * odist));  // -- ** this line can be iterated for more accuracy ** --
+                        odist = odist * (rsopf - (rsf * odist * odist));
+                        dmin = 2.0f;                        // reset min distance to edge pixel
+                        for(a = 0; a < num_base_verts; a++) {     // loop through all inside edge pixels
+                            dx = t - base_verts_f[a][0];             // compute delta in Y from gradient pixel to inside edge pixel
+                            dy = fsz - base_verts_f[a][1];           // compute delta in X from gradient pixel to inside edge pixel
+                            ud = dx * dx + dy * dy;   // compute sum of squares
+                            if(ud < dmin) {           // if our new sum of squares is less than the current minimum we've found
+                                dmin = ud;            // set a new minimum equal to the new lower value
+                            }
+                        }
+                        idist = dmin;                    // cast inner min to a float
+                        rsf = idist * 0.5f;                       //
+                        rsl = *(unsigned int *)&idist;            //
+                        rsl = 0x5f3759df - (rsl >> 1);            // see notes above
+                        idist = *(float *)&rsl;                   //
+                        idist = idist * (rsopf - (rsf * idist * idist));  //
+                        idist = idist * (rsopf - (rsf * idist * idist));
+                        /*
+                         * Note once again that since we are using reciprocals of distance values our
+                         * proportion is already the correct intensity, and does not need to be
+                         * subracted from 1.0 like it would have if we used real distances.
+                         */
+#else
+                        clup[0]=t;
+                        clup[1]=fsz;
+                        res_kdi=kd_nearestf(ctx->kdi,clup);
+                        res_kdo=kd_nearestf(ctx->kdo,clup);
+                        kd_res_itemf(res_kdi,clup);
+                        dx=t-clup[0];
+                        dy=fsz-clup[1];
+                        idist=dx*dx+dy*dy;
+                        rsf = idist * 0.5f;                       //
+                        rsl = *(unsigned int *)&idist;            //
+                        rsl = 0x5f3759df - (rsl >> 1);            // see notes above
+                        idist = *(float *)&rsl;                   //
+                        idist = idist * (rsopf - (rsf * idist * idist));  //
+                        idist = idist * (rsopf - (rsf * idist * idist));
+                        kd_res_itemf(res_kdo,clup);
+                        dx=t-clup[0];
+                        dy=fsz-clup[1];
+                        odist=dx*dx+dy*dy;
+                        rsf = odist * 0.5f;                       //
+                        rsl = *(unsigned int *)&odist;            // use some peculiar properties of the way bits are stored
+                        rsl = 0x5f3759df - (rsl >> 1);            // in floats vs. unsigned ints to compute an approximate
+                        odist = *(float *)&rsl;                   // reciprocal square root
+                        odist = odist * (rsopf - (rsf * odist * odist));  // -- ** this line can be iterated for more accuracy ** --
+                        odist = odist * (rsopf - (rsf * odist * odist));
+
+#endif
+                        /* set intensity, do the += so overlapping gradients are additive */
+                        *cpxl = (idist / (idist+odist));
+                    }
+                }
+            }
+        }
+
+        /*
+         * Loop through all edges of polygon that could be hit by this scan line,
+         * and figure out their x-intersections with the next scan line.
+         *
+         * Either A.) we wont have any more edges to test, or B.) we just add on the
+         * slope delta computed in preprocessing step. Since this draws non-antialiased
+         * polygons, we dont have fractional positions, so we only move in x-direction
+         * when needed to get all the way to the next pixel over...
+         */
+        for(edgec = &ctx->possible_edges; (e_curr = *edgec);) {
+            if(!(--(e_curr->num))) {
+                *edgec = e_curr->e_next;
+            } else {
+                e_curr->x += e_curr->xshift;
+                if((e_curr->drift += e_curr->drift_inc) > 0) {
+                    e_curr->x += e_curr->xdir;
+                    e_curr->drift -= e_curr->drift_dec;
+                }
+                edgec = &e_curr->e_next;
+            }
+        }
+        /*
+         * It's possible that some edges may have crossed during the last step, so we'll be sure
+         * that we ALWAYS intersect scan lines in order by shuffling if needed to make all edges
+         * sorted by x-intersection coordinate. We'll always scan through at least once to see if
+         * edges crossed, and if so, we set the 'swixd' flag. If 'swixd' gets set on the initial
+         * pass, then we know we need to sort by x, so then cycle through edges again and perform
+         * the sort.-
+         */
+        if(ctx->possible_edges) {
+            for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                /* if the current edge hits scan line at greater X than the next edge, we need to exchange the edges */
+                if(e_curr->x > e_curr->e_next->x) {
+                    *edgec = e_curr->e_next;
+                    /* exchange the pointers */
+                    e_temp = e_curr->e_next->e_next;
+                    e_curr->e_next->e_next = e_curr;
+                    e_curr->e_next = e_temp;
+                    /* set flag that we had at least one switch */
+                    swixd = 1;
+                }
+            }
+            /* if we did have a switch, look for more (there will more if there was one) */
+            for(;;) {
+                /* reset exchange flag so it's only set if we encounter another one */
+                swixd = 0;
+                for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                    /* again, if current edge hits scan line at higher X than next edge,
+                     * exchange the edges and set flag */
+                    if(e_curr->x > e_curr->e_next->x) {
+                        *edgec = e_curr->e_next;
+                        /* exchange the pointers */
+                        e_temp = e_curr->e_next->e_next;
+                        e_curr->e_next->e_next = e_curr;
+                        e_curr->e_next = e_temp;
+                        /* flip the exchanged flag */
+                        swixd = 1;
+                    }
+                }
+                /* if we had no exchanges, we're done reshuffling the pointers */
+                if(!swixd) {
+                    break;
+                }
+            }
+        }
+    }
+
+    free(edgbuf);
+    return 1;
 }
 
-int PLX_raskterize_feather(float (*base_verts)[2], int num_base_verts, float (*feather_verts)[2], int num_feather_verts,
-						   float *buf, int buf_x, int buf_y) {
-	int i;                            /* i: Loop counter. */
-	struct poly_vert *fe;             /* fe: Pointer to a list of integer buffer-space feather vertex coords. */
-	struct r_fill_context ctx = {0};
-
-	/* for faster multiply */
-	const float buf_x_f = (float)buf_x;
-	const float buf_y_f = (float)buf_y;
-
-	/*
-	 * Allocate enough memory for our poly_vert list. It'll be the size of the poly_vert
-	 * data structure multiplied by the number of verts.
-	 *
-	 * In the event of a failure to allocate the memory, return 0, so this error can
-	 * be distinguished as a memory allocation error.
-	 */
-	if ((fe = (struct poly_vert *)(malloc(sizeof(struct poly_vert) * num_feather_verts))) == NULL) {
-		return(0);
-	}
-
-	/*
-	 * Loop over all verts passed in to be rasterized. Each vertex's X and Y coordinates are
-	 * then converted from normalized screen space (0.0 <= POS <= 1.0) to integer coordinates
-	 * in the buffer-space coordinates passed in inside buf_x and buf_y.
-	 *
-	 * It's worth noting that this function ONLY outputs fully white pixels in a mask. Every pixel
-	 * drawn will be 1.0f in value, there is no anti-aliasing.
-	 */
-	for (i = 0; i < num_feather_verts; i++) {            /* Loop over all verts. */
-		fe[i].x = (int)((feather_verts[i][0] * buf_x_f) + 0.5f);  /* Range expand normalized X to integer buffer-space X. */
-		fe[i].y = (int)((feather_verts[i][1] * buf_y_f) + 0.5f);  /* Range expand normalized Y to integer buffer-space Y. */
-	}
-
-	ctx.rb.buf = buf;                            /* Set the output buffer pointer. */
-	ctx.rb.sizex = buf_x;                        /* Set the output buffer size in X. (width) */
-	ctx.rb.sizey = buf_y;                        /* Set the output buffer size in Y. (height) */
-
-	/* Call our rasterizer, passing in the integer coords for each vert. */
-	i = rast_scan_feather(&ctx, base_verts, num_base_verts, fe, feather_verts, num_feather_verts);
-	free(fe);
-	return i;                                   /* Return the value returned by the rasterizer. */
+int PLX_raskterize_feather(float(*base_verts)[2], int num_base_verts, float(*feather_verts)[2], int num_feather_verts,
+                           float *buf, int buf_x, int buf_y) {
+    //void plx_floatsort(float(*f)[2], unsigned int n, int sortby);
+    int i;                            /* i: Loop counter. */
+    struct poly_vert *fe;             /* fe: Pointer to a list of integer buffer-space feather vertex coords. */
+    struct r_fill_context ctx = {0};
+
+    /* for faster multiply */
+    const float buf_x_f = (float)buf_x;
+    const float buf_y_f = (float)buf_y;
+#ifdef __PLX_KD_TREE__
+    ctx.kdi = kd_create(2);
+    ctx.kdo = kd_create(2);
+#endif
+    /*
+     * Allocate enough memory for our poly_vert list. It'll be the size of the poly_vert
+     * data structure multiplied by the number of verts.
+     *
+     * In the event of a failure to allocate the memory, return 0, so this error can
+     * be distinguished as a memory allocation error.
+     */
+    if((fe = (struct poly_vert *)(malloc(sizeof(struct poly_vert) * num_feather_verts))) == NULL) {
+        return(0);
+    }
+
+    /*
+     * Loop over all verts passed in to be rasterized. Each vertex's X and Y coordinates are
+     * then converted from normalized screen space (0.0 <= POS <= 1.0) to integer coordinates
+     * in the buffer-space coordinates passed in inside buf_x and buf_y.
+     *
+     * It's worth noting that this function ONLY outputs fully white pixels in a mask. Every pixel
+     * drawn will be 1.0f in value, there is no anti-aliasing.
+     */
+    for(i = 0; i < num_feather_verts; i++) {             /* Loop over all verts. */
+        fe[i].x = (int)((feather_verts[i][0] * buf_x_f) + 0.5f);  /* Range expand normalized X to integer buffer-space X. */
+        fe[i].y = (int)((feather_verts[i][1] * buf_y_f) + 0.5f);  /* Range expand normalized Y to integer buffer-space Y. */
+#ifdef __PLX_KD_TREE__
+        kd_insertf(ctx.kdo,feather_verts[i],NULL);
+    }
+    for(i=0;i<num_base_verts;i++){
+        kd_insertf(ctx.kdi,base_verts[i],NULL);
+#endif
+    }
+
+    ctx.rb.buf = buf;                            /* Set the output buffer pointer. */
+    ctx.rb.sizex = buf_x;                        /* Set the output buffer size in X. (width) */
+    ctx.rb.sizey = buf_y;                        /* Set the output buffer size in Y. (height) */
+    /* pre-sort the sets of edge verts on y */
+    //plx_floatsort(base_verts,num_base_verts,0);
+    //plx_floatsort(base_verts,num_base_verts,1);
+    //plx_floatsort(feather_verts,num_feather_verts,0);
+    //plx_floatsort(feather_verts,num_feather_verts,1);
+    /* Call our rasterizer, passing in the integer coords for each vert. */
+    i = rast_scan_feather(&ctx, base_verts, num_base_verts, fe, feather_verts, num_feather_verts);
+    free(fe);
+    return i;                                   /* Return the value returned by the rasterizer. */
 }
 
 #ifndef __PLX__FAKE_AA__
 
-static int get_range_expanded_pixel_coord(float normalized_value, int max_value)
-{
-	return (int)((normalized_value * (float)(max_value)) + 0.5f);
+int get_range_expanded_pixel_coord(float normalized_value, int max_value) {
+    return (int)((normalized_value * (float)(max_value)) + 0.5f);
 }
 
-static float get_pixel_intensity(float *buf, int buf_x, int buf_y, int pos_x, int pos_y)
-{
-	if(pos_x < 0 || pos_x >= buf_x || pos_y < 0 || pos_y >= buf_y) {
-		return 0.0f;
-	}
-	return buf[(pos_y * buf_y) + buf_x];
+float get_pixel_intensity(float *buf, int buf_x, int buf_y, int pos_x, int pos_y) {
+    if(pos_x < 0 || pos_x >= buf_x || pos_y < 0 || pos_y >= buf_y) {
+        return 0.0f;
+    }
+    return buf[(pos_y * buf_x) + pos_x];
 }
 
-static float get_pixel_intensity_bilinear(float *buf, int buf_x, int buf_y, float u, float v)
-{
-	int a;
-	int b;
-	int a_plus_1;
-	int b_plus_1;
-	float prop_u;
-	float prop_v;
-	float inv_prop_u;
-	float inv_prop_v;
-	if(u<0.0f || u>1.0f || v<0.0f || v>1.0f) {
-		return 0.0f;
-	}
-	u = u * (float)(buf_x) - 0.5f;
-	v = v * (float)(buf_y) - 0.5f;
-	a = (int)(u);
-	b = (int)(v);
-	prop_u = u - (float)(a);
-	prop_v = v - (float)(b);
-	inv_prop_u = 1.0f - prop_u;
-	inv_prop_v = 1.0f - prop_v;
-	a_plus_1 = MIN2((buf_x-1),a+1);
-	b_plus_1 = MIN2((buf_y-1),b+1);
-	return (buf[(b * buf_y) + a] * inv_prop_u + buf[(b*buf_y)+(a_plus_1)] * prop_u)*inv_prop_v+(buf[((b_plus_1) * buf_y)+a] * inv_prop_u + buf[((b_plus_1)*buf_y)+(a_plus_1)] * prop_u) * prop_v;
+float get_pixel_intensity_bilinear(float *buf, int buf_x, int buf_y, float u, float v) {
+    int a;
+    int b;
+    int a_plus_1;
+    int b_plus_1;
+    float prop_u;
+    float prop_v;
+    float inv_prop_u;
+    float inv_prop_v;
+    if(u<0.0f || u>1.0f || v<0.0f || v>1.0f) {
+        return 0.0f;
+    }
+    u = u * (float)(buf_x) - 0.5f;
+    v = v * (float)(buf_y) - 0.5f;
+    a = (int)(u);
+    b = (int)(v);
+    prop_u = u - (float)(a);
+    prop_v = v - (float)(b);
+    inv_prop_u = 1.0f - prop_u;
+    inv_prop_v = 1.0f - prop_v;
+    a_plus_1 = MIN2((buf_x-1),a+1);
+    b_plus_1 = MIN2((buf_y-1),b+1);
+    return (buf[(b * buf_x) + a] * inv_prop_u + buf[(b*buf_x)+(a_plus_1)] * prop_u)*inv_prop_v+(buf[((b_plus_1) * buf_x)+a] * inv_prop_u + buf[((b_plus_1)*buf_x)+(a_plus_1)] * prop_u) * prop_v;
 
 }
 
-static void set_pixel_intensity(float *buf, int buf_x, int buf_y, int pos_x, int pos_y, float intensity)
-{
-	if(pos_x < 0 || pos_x >= buf_x || pos_y < 0 || pos_y >= buf_y) {
-		return;
-	}
-	buf[(pos_y * buf_y) + buf_x] = intensity;
+void set_pixel_intensity(float *buf, int buf_x, int buf_y, int pos_x, int pos_y, float intensity) {
+    if(pos_x < 0 || pos_x >= buf_x || pos_y < 0 || pos_y >= buf_y) {
+        return;
+    }
+    buf[(pos_y * buf_x) + pos_x] = intensity;
 }
 #endif
 
-int PLX_antialias_buffer(float *buf, int buf_x, int buf_y)
-{
+int PLX_antialias_buffer(float *buf, int buf_x, int buf_y) {
 #ifdef __PLX__FAKE_AA__
 #ifdef __PLX_GREY_AA__
-	int i=0;
-	int sz = buf_x * buf_y;
-	for(i=0; i<sz; i++) {
-		buf[i] *= 0.5f;
-	}
+    int i=0;
+    int sz = buf_x * buf_y;
+    for(i=0; i<sz; i++) {
+        buf[i] *= 0.5f;
+    }
 #endif
-	(void)buf_x;
-	(void)buf_y;
-	(void)buf;
-	return 1;
+    (void)buf_x;
+    (void)buf_y;
+    (void)buf;
+    return 1;
+#else
+    const float p0 = 1.0f;
+    const float p1 = 1.0f;
+    const float p2 = 1.0f;
+    const float p3 = 1.0f;
+    const float p4 = 1.0f;
+    const float p5 = 1.5f;
+    const float p6 = 2.0f;
+    const float p7 = 2.0f;
+    const float p8 = 2.0f;
+    const float p9 = 2.0f;
+    const float p10 = 4.0f;
+    const float p11 = 8.0f;
+
+    const float edge_threshold = 0.063f;
+    const float edge_threshold_min = 0.0312f;
+    const float quality_subpix = 1.0f;
+
+    float posM_x,posM_y;
+    float posB_x,posB_y;
+    float posN_x,posN_y;
+    float posP_x,posP_y;
+    float offNP_x,offNP_y;
+    float lumaM;
+    float lumaS;
+    float lumaE;
+    float lumaN;
+    float lumaW;
+    float lumaNW;
+    float lumaSE;
+    float lumaNE;
+    float lumaSW;
+    float lumaNS;
+    float lumaWE;
+    float lumaNESE;
+    float lumaNWNE;
+    float lumaNWSW;
+    float lumaSWSE;
+    float lumaNN;
+    float lumaSS;
+    float lumaEndN;
+    float lumaEndP;
+    float lumaMM;
+    float lumaMLTZero;
+    float subpixNWSWNESE;
+    float subpixRcpRange;
+    float subpixNSWE;
+    float maxSM;
+    float minSM;
+    float maxESM;
+    float minESM;
+    float maxWN;
+    float minWN;
+    float rangeMax;
+    float rangeMin;
+    float rangeMaxScaled;
+    float range;
+    float rangeMaxClamped;
+    float edgeHorz;
+    float edgeVert;
+    float edgeHorz1;
+    float edgeVert1;
+    float edgeHorz2;
+    float edgeVert2;
+    float edgeHorz3;
+    float edgeVert3;
+    float edgeHorz4;
+    float edgeVert4;
+    float lengthSign;
+    float subpixA;
+    float subpixB;
+    float subpixC;
+    float subpixD;
+    float subpixE;
+    float subpixF;
+    float subpixG;
+    float subpixH;
+    float gradientN;
+    float gradientS;
+    float gradient;
+    float gradientScaled;
+    float dstN;
+    float dstP;
+    float dst;
+    float spanLength;
+    float spanLengthRcp;
+    float pixelOffset;
+    float pixelOffsetGood;
+    float pixelOffsetSubpix;
+    float opx;
+    float opy;
+    int directionN;
+    int goodSpan;
+    int goodSpanN;
+    int goodSpanP;
+    int horzSpan;
+    int earlyExit;
+    int pairN;
+    int doneN;
+    int doneP;
+    int doneNP;
+    int curr_x=0;
+    int curr_y=0;
+    opx = (1.0f / (float)(buf_x));
+    opy = (1.0f / (float)(buf_y));
+    for(curr_y=0; curr_y < buf_y; curr_y++) {
+        for(curr_x=0; curr_x < buf_x; curr_x++) {
+            posM_x = ((float)(curr_x) + 0.5f) * opx;
+            posM_y = ((float)(curr_y) + 0.5f) * opy;
+//#define __PLX_BILINEAR_INITIAL_SAMPLES__ 1
+#ifdef __PLX_BILINEAR_INITIAL_SAMPLES__
+            lumaM = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x, posM_y);
+            lumaS = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x, posM_y + opy);
+            lumaE = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x + opx, posM_y);
+            lumaN = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x, posM_y - opy);
+            lumaW = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x - opx, posM_y);
 #else
-	/*XXX - TODO: THIS IS NOT FINAL CODE - IT DOES NOT WORK - DO NOT ENABLE IT */
-	const float p0 = 1.0f;
-	const float p1 = 1.0f;
-	const float p2 = 1.0f;
-	const float p3 = 1.0f;
-	const float p4 = 1.0f;
-	const float p5 = 1.5f;
-	const float p6 = 2.0f;
-	const float p7 = 2.0f;
-	const float p8 = 2.0f;
-	const float p9 = 2.0f;
-	const float p10 = 4.0f;
-	const float p11 = 8.0f;
-
-	const float edge_threshold = 0.063f;
-	const float edge_threshold_min = 0.0312f;
-	const float quality_subpix = 1.0f;
-//	int px_x;
-//	int px_y;
-
-	float posM_x,posM_y;
-	float posB_x,posB_y;
-	float posN_x,posN_y;
-	float posP_x,posP_y;
-	float offNP_x,offNP_y;
-	float lumaM;
-	float lumaS;
-	float lumaE;
-	float lumaN;
-	float lumaW;
-	float lumaNW;
-	float lumaSE;
-	float lumaNE;
-	float lumaSW;
-	float lumaNS;
-	float lumaWE;
-	float lumaNESE;
-	float lumaNWNE;
-	float lumaNWSW;
-	float lumaSWSE;
-	float lumaNN;
-	float lumaSS;
-	float lumaEndN;
-	float lumaEndP;
-	float lumaMM;
-	float lumaMLTZero;
-	float subpixNWSWNESE;
-	float subpixRcpRange;
-	float subpixNSWE;
-	float maxSM;
-	float minSM;
-	float maxESM;
-	float minESM;
-	float maxWN;
-	float minWN;
-	float rangeMax;
-	float rangeMin;
-	float rangeMaxScaled;
-	float range;
-	float rangeMaxClamped;
-	float edgeHorz;
-	float edgeVert;
-	float edgeHorz1;
-	float edgeVert1;
-	float edgeHorz2;
-	float edgeVert2;
-	float edgeHorz3;
-	float edgeVert3;
-	float edgeHorz4;
-	float edgeVert4;
-	float lengthSign;
-	float subpixA;
-	float subpixB;
-	float subpixC;
-	float subpixD;
-	float subpixE;
-	float subpixF;
-	float subpixG;
-	float subpixH;
-	float gradientN;
-	float gradientS;
-	float gradient;
-	float gradientScaled;
-	float dstN;
-	float dstP;
-	float dst;
-	float spanLength;
-	float spanLengthRcp;
-	float pixelOffset;
-	float pixelOffsetGood;
-	float pixelOffsetSubpix;
-	int directionN;
-	int goodSpan;
-	int goodSpanN;
-	int goodSpanP;
-	int horzSpan;
-	int earlyExit;
-	int pairN;
-	int doneN;
-	int doneP;
-	int doneNP;
-	int curr_x=0;
-	int curr_y=0;
-	for(curr_y=0; curr_y < buf_y; curr_y++) {
-		for(curr_x=0; curr_x < buf_x; curr_x++) {
-			posM_x = ((float)(curr_x) + 0.5f) * (1.0f/(float)(buf_x));
-			posM_y = ((float)(curr_y) + 0.5f) * (1.0f/(float)(buf_y));
-
-			lumaM = get_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y);
-			lumaS = get_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y - 1);
-			lumaE = get_pixel_intensity(buf, buf_x, buf_y, curr_x + 1, curr_y);
-			lumaN = get_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y + 1);
-			lumaW = get_pixel_intensity(buf, buf_x, buf_y, curr_x - 1, curr_y);
-
-			maxSM = MAX2(lumaS, lumaM);
-			minSM = MIN2(lumaS, lumaM);
-			maxESM = MAX2(lumaE, maxSM);
-			minESM = MIN2(lumaE, minSM);
-			maxWN = MAX2(lumaN, lumaW);
-			minWN = MIN2(lumaN, lumaW);
-			rangeMax = MAX2(maxWN, maxESM);
-			rangeMin = MIN2(minWN, minESM);
-			rangeMaxScaled = rangeMax * edge_threshold;
-			range = rangeMax - rangeMin;
-			rangeMaxClamped = MAX2(edge_threshold_min, rangeMaxScaled);
-
-			earlyExit = range < rangeMaxClamped ? 1:0;
-			if(earlyExit) {
-				set_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y, lumaM);
-			}
-
-			lumaNW = get_pixel_intensity(buf, buf_x, buf_y, curr_x + 1, curr_y - 1);
-			lumaSE = get_pixel_intensity(buf, buf_x, buf_y, curr_x - 1, curr_y + 1);
-			lumaNE = get_pixel_intensity(buf, buf_x, buf_y, curr_x + 1, curr_y + 1);
-			lumaSW = get_pixel_intensity(buf, buf_x, buf_y, curr_x - 1, curr_y - 1);
-
-			lumaNS = lumaN + lumaS;
-			lumaWE = lumaW + lumaE;
-			subpixRcpRange = 1.0f/range;
-			subpixNSWE = lumaNS + lumaWE;
-			edgeHorz1 = (-2.0f * lumaM) + lumaNS;
-			edgeVert1 = (-2.0f * lumaM) + lumaWE;
-
-			lumaNESE = lumaNE + lumaSE;
-			lumaNWNE = lumaNW + lumaNE;
-			edgeHorz2 = (-2.0f * lumaE) + lumaNESE;
-			edgeVert2 = (-2.0f * lumaN) + lumaNWNE;
-
-			lumaNWSW = lumaNW + lumaSW;
-			lumaSWSE = lumaSW + lumaSE;
-			edgeHorz4 = (ABS(edgeHorz1) * 2.0f) + ABS(edgeHorz2);
-			edgeVert4 = (ABS(edgeVert1) * 2.0f) + ABS(edgeVert2);
-			edgeHorz3 = (-2.0f * lumaW) + lumaNWSW;
-			edgeVert3 = (-2.0f * lumaS) + lumaSWSE;
-			edgeHorz = ABS(edgeHorz3) + edgeHorz4;
-			edgeVert = ABS(edgeVert3) + edgeVert4;
-
-			subpixNWSWNESE = lumaNWSW + lumaNESE;
-			lengthSign = 1.0f / (float)(buf_x);
-			horzSpan = edgeHorz >= edgeVert ? 1:0;
-			subpixA = subpixNSWE * 2.0f + subpixNWSWNESE;
-
-			if(!horzSpan) {
-				lumaN = lumaW;
-				lumaS = lumaE;
-			}
-			else {
-				lengthSign = 1.0f / (float)(buf_y);
-			}
-			subpixB = (subpixA * (1.0f/12.0f)) - lumaM;
-
-			gradientN = lumaN - lumaM;
-			gradientS = lumaS - lumaM;
-			lumaNN = lumaN + lumaM;
-			lumaSS = lumaS + lumaM;
-			pairN = (ABS(gradientN)) >= (ABS(gradientS)) ? 1:0;
-			gradient = MAX2(ABS(gradientN), ABS(gradientS));
-			if(pairN) {
-				lengthSign = -lengthSign;
-			}
-			subpixC = MAX2(MIN2(ABS(subpixB) * subpixRcpRange,1.0f),0.0f);
-
-			posB_x = posM_x;
-			posB_y = posM_y;
-			offNP_x = (!horzSpan) ? 0.0f:(1.0f / (float)(buf_x));
-			offNP_y = (horzSpan) ? 0.0f:(1.0f / (float)(buf_y));
-			if(!horzSpan) {
-				posB_x += lengthSign * 0.5f;
-			}
-			else {
-				posB_y += lengthSign * 0.5f;
-			}
-
-			posN_x = posB_x - offNP_x * p0;
-			posN_y = posB_y - offNP_y * p0;
-			posP_x = posB_x + offNP_x * p0;
-			posP_y = posB_y + offNP_y * p0;
-			subpixD = ((-2.0f)*subpixC) + 3.0f;
-			//may need bilinear filtered get_pixel_intensity() here...done
-			lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-			subpixE = subpixC * subpixC;
-			//may need bilinear filtered get_pixel_intensity() here...done
-			lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-
-			if(!pairN) {
-				lumaNN = lumaSS;
-			}
-			gradientScaled = gradient * 1.0f/4.0f;
-			lumaMM =lumaM - lumaNN * 0.5f;
-			subpixF = subpixD * subpixE;
-			lumaMLTZero = lumaMM < 0.0f ? 1:0;
-
-			lumaEndN -= lumaNN * 0.5f;
-			lumaEndP -= lumaNN * 0.5f;
-			doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-			doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-			if(!doneN) {
-				posN_x -= offNP_x * p1;
-				posN_y -= offNP_y * p1;
-			}
-			doneNP = (!doneN) || (!doneP) ? 1:0;
-			if(!doneP) {
-				posP_x += offNP_x * p1;
-				posP_y += offNP_y * p1;
-			}
-
-			if(doneNP) {
-				if(!doneN) {
-					lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posN_x,posN_y);
-				}
-				if(!doneP) {
-					lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x, posP_y);
-				}
-				if(!doneN) {
-					lumaEndN = lumaEndN - lumaNN * 0.5;
-				}
-				if(!doneP) {
-					lumaEndP = lumaEndP - lumaNN * 0.5;
-				}
-				doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-				doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-				if(!doneN) {
-					posN_x -= offNP_x * p2;
-					posN_y -= offNP_y * p2;
-				}
-				doneNP = (!doneN) || (!doneP) ? 1:0;
-				if(!doneP) {
-					posP_x += offNP_x * p2;
-					posP_y += offNP_y * p2;
-				}
-				if(doneNP) {
-					if(!doneN) {
-						lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-					}
-					if(!doneP) {
-						lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-					}
-					if(!doneN) {
-						lumaEndN = lumaEndN - lumaNN * 0.5;
-					}
-					if(!doneP) {
-						lumaEndP = lumaEndP - lumaNN * 0.5;
-					}
-					doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-					doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-					if(!doneN) {
-						posN_x -= offNP_x * p3;
-						posN_y -= offNP_y * p3;
-					}
-					doneNP = (!doneN) || (!doneP) ? 1:0;
-					if(!doneP) {
-						posP_x += offNP_x * p3;
-						posP_y += offNP_y * p3;
-					}
-					if(doneNP) {
-						if(!doneN) {
-							lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-						}
-						if(!doneP) {
-							lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-						}
-						if(!doneN) {
-							lumaEndN = lumaEndN - lumaNN * 0.5;
-						}
-						if(!doneP) {
-							lumaEndP = lumaEndP - lumaNN * 0.5;
-						}
-						doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-						doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-						if(!doneN) {
-							posN_x -= offNP_x * p4;
-							posN_y -= offNP_y * p4;
-						}
-						doneNP = (!doneN) || (!doneP) ? 1:0;
-						if(!doneP) {
-							posP_x += offNP_x * p4;
-							posP_y += offNP_y * p4;
-						}
-						if(doneNP) {
-							if(!doneN) {
-								lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-							}
-							if(!doneP) {
-								lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-							}
-							if(!doneN) {
-								lumaEndN = lumaEndN - lumaNN * 0.5;
-							}
-							if(!doneP) {
-								lumaEndP = lumaEndP - lumaNN * 0.5;
-							}
-							doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-							doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-							if(!doneN) {
-								posN_x -= offNP_x * p5;
-								posN_y -= offNP_y * p5;
-							}
-							doneNP = (!doneN) || (!doneP) ? 1:0;
-							if(!doneP) {
-								posP_x += offNP_x * p5;
-								posP_y += offNP_y * p5;
-							}
-							if(doneNP) {
-								if(!doneN) {
-									lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-								}
-								if(!doneP) {
-									lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-								}
-								if(!doneN) {
-									lumaEndN = lumaEndN - lumaNN * 0.5;
-								}
-								if(!doneP) {
-									lumaEndP = lumaEndP - lumaNN * 0.5;
-								}
-								doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-								doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-								if(!doneN) {
-									posN_x -= offNP_x * p6;
-									posN_y -= offNP_y * p6;
-								}
-								doneNP = (!doneN) || (!doneP) ? 1:0;
-								if(!doneP) {
-									posP_x += offNP_x * p6;
-									posP_y += offNP_y * p6;
-								}
-								if(doneNP) {
-									if(!doneN) {
-										lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-									}
-									if(!doneP) {
-										lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-									}
-									if(!doneN) {
-										lumaEndN = lumaEndN - lumaNN * 0.5;
-									}
-									if(!doneP) {
-										lumaEndP = lumaEndP - lumaNN * 0.5;
-									}
-									doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-									doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-									if(!doneN) {
-										posN_x -= offNP_x * p7;
-										posN_y -= offNP_y * p7;
-									}
-									doneNP = (!doneN) || (!doneP) ? 1:0;
-									if(!doneP) {
-										posP_x += offNP_x * p7;
-										posP_y += offNP_y * p7;
-									}
-									if(doneNP) {
-										if(!doneN) {
-											lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-										}
-										if(!doneP) {
-											lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-										}
-										if(!doneN) {
-											lumaEndN = lumaEndN - lumaNN * 0.5;
-										}
-										if(!doneP) {
-											lumaEndP = lumaEndP - lumaNN * 0.5;
-										}
-										doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-										doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-										if(!doneN) {
-											posN_x -= offNP_x * p8;
-											posN_y -= offNP_y * p8;
-										}
-										doneNP = (!doneN) || (!doneP) ? 1:0;
-										if(!doneP) {
-											posP_x += offNP_x * p8;
-											posP_y += offNP_y * p8;
-										}
-										if(doneNP) {
-											if(!doneN) {
-												lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-											}
-											if(!doneP) {
-												lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-											}
-											if(!doneN) {
-												lumaEndN = lumaEndN - lumaNN * 0.5;
-											}
-											if(!doneP) {
-												lumaEndP = lumaEndP - lumaNN * 0.5;
-											}
-											doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-											doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-											if(!doneN) {
-												posN_x -= offNP_x * p9;
-												posN_y -= offNP_y * p9;
-											}
-											doneNP = (!doneN) || (!doneP) ? 1:0;
-											if(!doneP) {
-												posP_x += offNP_x * p9;
-												posP_y += offNP_y * p9;
-											}
-											if(doneNP) {
-												if(!doneN) {
-													lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-												}
-												if(!doneP) {
-													lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-												}
-												if(!doneN) {
-													lumaEndN = lumaEndN - lumaNN * 0.5;
-												}
-												if(!doneP) {
-													lumaEndP = lumaEndP - lumaNN * 0.5;
-												}
-												doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-												doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-												if(!doneN) {
-													posN_x -= offNP_x * p10;
-													posN_y -= offNP_y * p10;
-												}
-												doneNP = (!doneN) || (!doneP) ? 1:0;
-												if(!doneP) {
-													posP_x += offNP_x * p10;
-													posP_y += offNP_y * p10;
-												}
-												if(doneNP) {
-													if(!doneN) {
-														lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
-													}
-													if(!doneP) {
-														lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
-													}
-													if(!doneN) {
-														lumaEndN = lumaEndN - lumaNN * 0.5;
-													}
-													if(!doneP) {
-														lumaEndP = lumaEndP - lumaNN * 0.5;
-													}
-													doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
-													doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
-													if(!doneN) {
-														posN_x -= offNP_x * p11;
-														posN_y -= offNP_y * p11;
-													}
-													doneNP = (!doneN) || (!doneP) ? 1:0;
-													if(!doneP) {
-														posP_x += offNP_x * p11;
-														posP_y += offNP_y * p11;
-													}
-												}
-											}
-										}
-									}
-								}
-							}
-						}
-					}
-				}
-			}
-			dstN = posM_x - posN_x;
-			dstP = posP_x - posM_x;
-			if(!horzSpan) {
-				dstN = posM_y - posN_y;
-				dstP = posP_y - posM_y;
-			}
-
-			goodSpanN = ((lumaEndN < 0.0f) ? 1:0) != lumaMLTZero ? 1:0;
-			spanLength = (dstP + dstN);
-			goodSpanP = ((lumaEndP < 0.0f) ? 1:0) != lumaMLTZero ? 1:0;
-			spanLengthRcp = 1.0f/spanLength;
-
-			directionN = dstN < dstP ? 1:0;
-			dst = MIN2(dstN, dstP);
-			goodSpan = (directionN==1) ? goodSpanN:goodSpanP;
-			subpixG = subpixF * subpixF;
-			pixelOffset = (dst * (-spanLengthRcp)) + 0.5f;
-			subpixH = subpixG * quality_subpix;
-
-			pixelOffsetGood = (goodSpan==1) ? pixelOffset : 0.0f;
-			pixelOffsetSubpix = MAX2(pixelOffsetGood, subpixH);
-			if(!horzSpan) {
-				posM_x += pixelOffsetSubpix * lengthSign;
-			}
-			else {
-				posM_y += pixelOffsetSubpix * lengthSign;
-			}
-			//may need bilinear filtered get_pixel_intensity() here...
-			set_pixel_intensity(buf,buf_x,buf_y,curr_x,curr_y,get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x,posM_y)* lumaM);
-
-		}
-	}
-	return 1;
+    lumaM = get_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y);
+    lumaS = get_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y + 1);
+    lumaE = get_pixel_intensity(buf, buf_x, buf_y, curr_x + 1, curr_y);
+    lumaN = get_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y - 1);
+    lumaW = get_pixel_intensity(buf, buf_x, buf_y, curr_x - 1, curr_y);
+#endif
+            maxSM = MAX2(lumaS, lumaM);
+            minSM = MIN2(lumaS, lumaM);
+            maxESM = MAX2(lumaE, maxSM);
+            minESM = MIN2(lumaE, minSM);
+            maxWN = MAX2(lumaN, lumaW);
+            minWN = MIN2(lumaN, lumaW);
+            rangeMax = MAX2(maxWN, maxESM);
+            rangeMin = MIN2(minWN, minESM);
+            rangeMaxScaled = rangeMax * edge_threshold;
+            range = rangeMax - rangeMin;
+            rangeMaxClamped = MAX2(edge_threshold_min, rangeMaxScaled);
+
+            earlyExit = range < rangeMaxClamped ? 1:0;
+            if(earlyExit) {
+                set_pixel_intensity(buf, buf_x, buf_y, curr_x, curr_y, lumaM);
+            } else {
+
+                lumaNW = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x - opx, posM_y - opy);
+                lumaSE = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x + opx, posM_y + opy);
+                lumaNE = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x + opx, posM_y - opy);
+                lumaSW = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x - opx, posM_y + opy);
+
+                lumaNS = lumaN + lumaS;
+                lumaWE = lumaW + lumaE;
+                subpixRcpRange = 1.0f/range;
+                subpixNSWE = lumaNS + lumaWE;
+                edgeHorz1 = (-2.0f * lumaM) + lumaNS;
+                edgeVert1 = (-2.0f * lumaM) + lumaWE;
+
+                lumaNESE = lumaNE + lumaSE;
+                lumaNWNE = lumaNW + lumaNE;
+                edgeHorz2 = (-2.0f * lumaE) + lumaNESE;
+                edgeVert2 = (-2.0f * lumaN) + lumaNWNE;
+
+                lumaNWSW = lumaNW + lumaSW;
+                lumaSWSE = lumaSW + lumaSE;
+                edgeHorz4 = (ABS(edgeHorz1) * 2.0f) + ABS(edgeHorz2);
+                edgeVert4 = (ABS(edgeVert1) * 2.0f) + ABS(edgeVert2);
+                edgeHorz3 = (-2.0f * lumaW) + lumaNWSW;
+                edgeVert3 = (-2.0f * lumaS) + lumaSWSE;
+                edgeHorz = ABS(edgeHorz3) + edgeHorz4;
+                edgeVert = ABS(edgeVert3) + edgeVert4;
+
+                subpixNWSWNESE = lumaNWSW + lumaNESE;
+                lengthSign = 1.0f / (float)(buf_x);
+                horzSpan = edgeHorz >= edgeVert ? 1:0;
+                subpixA = subpixNSWE * 2.0f + subpixNWSWNESE;
+
+                if(!horzSpan) {
+                    lumaN = lumaW;
+                    lumaS = lumaE;
+                } else {
+                    lengthSign = 1.0f / (float)(buf_y);
+                }
+                subpixB = (subpixA * (1.0f/12.0f)) - lumaM;
+
+                gradientN = lumaN - lumaM;
+                gradientS = lumaS - lumaM;
+                lumaNN = lumaN + lumaM;
+                lumaSS = lumaS + lumaM;
+                pairN = (ABS(gradientN)) >= (ABS(gradientS)) ? 1:0;
+                gradient = MAX2(ABS(gradientN), ABS(gradientS));
+                if(pairN) {
+                    lengthSign = -lengthSign;
+                }
+                subpixC = MAX2(MIN2(ABS(subpixB) * subpixRcpRange,1.0f),0.0f);
+
+                posB_x = posM_x;
+                posB_y = posM_y;
+                offNP_x = (!horzSpan) ? 0.0f:(1.0f / (float)(buf_x));
+                offNP_y = (horzSpan) ? 0.0f:(1.0f / (float)(buf_y));
+                if(!horzSpan) {
+                    posB_x += lengthSign * 0.5f;
+                } else {
+                    posB_y += lengthSign * 0.5f;
+                }
+
+                posN_x = posB_x - offNP_x * p0;
+                posN_y = posB_y - offNP_y * p0;
+                posP_x = posB_x + offNP_x * p0;
+                posP_y = posB_y + offNP_y * p0;
+                subpixD = ((-2.0f)*subpixC) + 3.0f;
+                lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                subpixE = subpixC * subpixC;
+                lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+
+                if(!pairN) {
+                    lumaNN = lumaSS;
+                }
+                gradientScaled = gradient * 1.0f/4.0f;
+                lumaMM =lumaM - lumaNN * 0.5f;
+                subpixF = subpixD * subpixE;
+                lumaMLTZero = lumaMM < 0.0f ? 1:0;
+
+                lumaEndN -= lumaNN * 0.5f;
+                lumaEndP -= lumaNN * 0.5f;
+                doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                if(!doneN) {
+                    posN_x -= offNP_x * p1;
+                    posN_y -= offNP_y * p1;
+                }
+                doneNP = (!doneN) || (!doneP) ? 1:0;
+                if(!doneP) {
+                    posP_x += offNP_x * p1;
+                    posP_y += offNP_y * p1;
+                }
+
+                if(doneNP) {
+                    if(!doneN) {
+                        lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posN_x,posN_y);
+                    }
+                    if(!doneP) {
+                        lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x, posP_y);
+                    }
+                    if(!doneN) {
+                        lumaEndN = lumaEndN - lumaNN * 0.5;
+                    }
+                    if(!doneP) {
+                        lumaEndP = lumaEndP - lumaNN * 0.5;
+                    }
+                    doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                    doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                    if(!doneN) {
+                        posN_x -= offNP_x * p2;
+                        posN_y -= offNP_y * p2;
+                    }
+                    doneNP = (!doneN) || (!doneP) ? 1:0;
+                    if(!doneP) {
+                        posP_x += offNP_x * p2;
+                        posP_y += offNP_y * p2;
+                    }
+                    if(doneNP) {
+                        if(!doneN) {
+                            lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                        }
+                        if(!doneP) {
+                            lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                        }
+                        if(!doneN) {
+                            lumaEndN = lumaEndN - lumaNN * 0.5;
+                        }
+                        if(!doneP) {
+                            lumaEndP = lumaEndP - lumaNN * 0.5;
+                        }
+                        doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                        doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                        if(!doneN) {
+                            posN_x -= offNP_x * p3;
+                            posN_y -= offNP_y * p3;
+                        }
+                        doneNP = (!doneN) || (!doneP) ? 1:0;
+                        if(!doneP) {
+                            posP_x += offNP_x * p3;
+                            posP_y += offNP_y * p3;
+                        }
+                        if(doneNP) {
+                            if(!doneN) {
+                                lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                            }
+                            if(!doneP) {
+                                lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                            }
+                            if(!doneN) {
+                                lumaEndN = lumaEndN - lumaNN * 0.5;
+                            }
+                            if(!doneP) {
+                                lumaEndP = lumaEndP - lumaNN * 0.5;
+                            }
+                            doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                            doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                            if(!doneN) {
+                                posN_x -= offNP_x * p4;
+                                posN_y -= offNP_y * p4;
+                            }
+                            doneNP = (!doneN) || (!doneP) ? 1:0;
+                            if(!doneP) {
+                                posP_x += offNP_x * p4;
+                                posP_y += offNP_y * p4;
+                            }
+                            if(doneNP) {
+                                if(!doneN) {
+                                    lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                }
+                                if(!doneP) {
+                                    lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                }
+                                if(!doneN) {
+                                    lumaEndN = lumaEndN - lumaNN * 0.5;
+                                }
+                                if(!doneP) {
+                                    lumaEndP = lumaEndP - lumaNN * 0.5;
+                                }
+                                doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                if(!doneN) {
+                                    posN_x -= offNP_x * p5;
+                                    posN_y -= offNP_y * p5;
+                                }
+                                doneNP = (!doneN) || (!doneP) ? 1:0;
+                                if(!doneP) {
+                                    posP_x += offNP_x * p5;
+                                    posP_y += offNP_y * p5;
+                                }
+                                if(doneNP) {
+                                    if(!doneN) {
+                                        lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                    }
+                                    if(!doneP) {
+                                        lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                    }
+                                    if(!doneN) {
+                                        lumaEndN = lumaEndN - lumaNN * 0.5;
+                                    }
+                                    if(!doneP) {
+                                        lumaEndP = lumaEndP - lumaNN * 0.5;
+                                    }
+                                    doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                    doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                    if(!doneN) {
+                                        posN_x -= offNP_x * p6;
+                                        posN_y -= offNP_y * p6;
+                                    }
+                                    doneNP = (!doneN) || (!doneP) ? 1:0;
+                                    if(!doneP) {
+                                        posP_x += offNP_x * p6;
+                                        posP_y += offNP_y * p6;
+                                    }
+                                    if(doneNP) {
+                                        if(!doneN) {
+                                            lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                        }
+                                        if(!doneP) {
+                                            lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                        }
+                                        if(!doneN) {
+                                            lumaEndN = lumaEndN - lumaNN * 0.5;
+                                        }
+                                        if(!doneP) {
+                                            lumaEndP = lumaEndP - lumaNN * 0.5;
+                                        }
+                                        doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                        doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                        if(!doneN) {
+                                            posN_x -= offNP_x * p7;
+                                            posN_y -= offNP_y * p7;
+                                        }
+                                        doneNP = (!doneN) || (!doneP) ? 1:0;
+                                        if(!doneP) {
+                                            posP_x += offNP_x * p7;
+                                            posP_y += offNP_y * p7;
+                                        }
+                                        if(doneNP) {
+                                            if(!doneN) {
+                                                lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                            }
+                                            if(!doneP) {
+                                                lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                            }
+                                            if(!doneN) {
+                                                lumaEndN = lumaEndN - lumaNN * 0.5;
+                                            }
+                                            if(!doneP) {
+                                                lumaEndP = lumaEndP - lumaNN * 0.5;
+                                            }
+                                            doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                            doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                            if(!doneN) {
+                                                posN_x -= offNP_x * p8;
+                                                posN_y -= offNP_y * p8;
+                                            }
+                                            doneNP = (!doneN) || (!doneP) ? 1:0;
+                                            if(!doneP) {
+                                                posP_x += offNP_x * p8;
+                                                posP_y += offNP_y * p8;
+                                            }
+                                            if(doneNP) {
+                                                if(!doneN) {
+                                                    lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                                }
+                                                if(!doneP) {
+                                                    lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                                }
+                                                if(!doneN) {
+                                                    lumaEndN = lumaEndN - lumaNN * 0.5;
+                                                }
+                                                if(!doneP) {
+                                                    lumaEndP = lumaEndP - lumaNN * 0.5;
+                                                }
+                                                doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                                doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                                if(!doneN) {
+                                                    posN_x -= offNP_x * p9;
+                                                    posN_y -= offNP_y * p9;
+                                                }
+                                                doneNP = (!doneN) || (!doneP) ? 1:0;
+                                                if(!doneP) {
+                                                    posP_x += offNP_x * p9;
+                                                    posP_y += offNP_y * p9;
+                                                }
+                                                if(doneNP) {
+                                                    if(!doneN) {
+                                                        lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                                    }
+                                                    if(!doneP) {
+                                                        lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                                    }
+                                                    if(!doneN) {
+                                                        lumaEndN = lumaEndN - lumaNN * 0.5;
+                                                    }
+                                                    if(!doneP) {
+                                                        lumaEndP = lumaEndP - lumaNN * 0.5;
+                                                    }
+                                                    doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                                    doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                                    if(!doneN) {
+                                                        posN_x -= offNP_x * p10;
+                                                        posN_y -= offNP_y * p10;
+                                                    }
+                                                    doneNP = (!doneN) || (!doneP) ? 1:0;
+                                                    if(!doneP) {
+                                                        posP_x += offNP_x * p10;
+                                                        posP_y += offNP_y * p10;
+                                                    }
+                                                    if(doneNP) {
+                                                        if(!doneN) {
+                                                            lumaEndN = get_pixel_intensity_bilinear(buf, buf_x, buf_y,posN_x,posN_y);
+                                                        }
+                                                        if(!doneP) {
+                                                            lumaEndP = get_pixel_intensity_bilinear(buf, buf_x, buf_y, posP_x,posP_y);
+                                                        }
+                                                        if(!doneN) {
+                                                            lumaEndN = lumaEndN - lumaNN * 0.5;
+                                                        }
+                                                        if(!doneP) {
+                                                            lumaEndP = lumaEndP - lumaNN * 0.5;
+                                                        }
+                                                        doneN = (ABS(lumaEndN)) >= gradientScaled ? 1:0;
+                                                        doneP = (ABS(lumaEndP)) >= gradientScaled ? 1:0;
+                                                        if(!doneN) {
+                                                            posN_x -= offNP_x * p11;
+                                                            posN_y -= offNP_y * p11;
+                                                        }
+                                                        doneNP = (!doneN) || (!doneP) ? 1:0;
+                                                        if(!doneP) {
+                                                            posP_x += offNP_x * p11;
+                                                            posP_y += offNP_y * p11;
+                                                        }
+                                                    }
+                                                }
+                                            }
+                                        }
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+                dstN = posM_x - posN_x;
+                dstP = posP_x - posM_x;
+                if(!horzSpan) {
+                    dstN = posM_y - posN_y;
+                    dstP = posP_y - posM_y;
+                }
+
+                goodSpanN = ((lumaEndN < 0.0f) ? 1:0) != lumaMLTZero ? 1:0;
+                spanLength = (dstP + dstN);
+                goodSpanP = ((lumaEndP < 0.0f) ? 1:0) != lumaMLTZero ? 1:0;
+                spanLengthRcp = 1.0f/spanLength;
+
+                directionN = dstN < dstP ? 1:0;
+                dst = MIN2(dstN, dstP);
+                goodSpan = (directionN==1) ? goodSpanN:goodSpanP;
+                subpixG = subpixF * subpixF;
+                pixelOffset = (dst * (-spanLengthRcp)) + 0.5f;
+                subpixH = subpixG * quality_subpix;
+
+                pixelOffsetGood = (goodSpan==1) ? pixelOffset : 0.0f;
+                pixelOffsetSubpix = MAX2(pixelOffsetGood, subpixH);
+                if(!horzSpan) {
+                    posM_x += pixelOffsetSubpix * lengthSign;
+                } else {
+                    posM_y += pixelOffsetSubpix * lengthSign;
+                }
+                //may need bilinear filtered get_pixel_intensity() here...
+                set_pixel_intensity(buf,buf_x,buf_y,curr_x,curr_y,get_pixel_intensity_bilinear(buf, buf_x, buf_y, posM_x,posM_y));
+            }
+        }
+    }
+    return 1;
 
 #endif
 }
 
+#define SWAP_POLYVERT(a,b)	ptemp[0]=(a)[0]; ptemp[1]=(a)[1]; (a)[0]=(b)[0]; (a)[1]=(b)[1]; (b)[0]=ptemp[0]; (b)[1]=ptemp[1];
+#define __PLX_SMALL_COUNT__ 13
+void plx_floatsort(float(*f)[2], unsigned int n, int sortby) {
+    unsigned int i;
+    unsigned int ir;
+    unsigned int j;
+    unsigned int k;
+    unsigned int l=1;
+    unsigned int istack[50];
+    int jstack=0;
+    float a[2];
+    float ptemp[2];
+
+    ir=n;
+    for(;;) {
+        if(ir-l < __PLX_SMALL_COUNT__) {
+            for(j=l+1; j<=ir; j++) {
+                a[1]=f[j][1];
+                a[0]=f[j][0];
+                for(i=j-1; i>=l; i--) {
+                    if(f[i][sortby] <= a[sortby]) {
+                        break;
+                    }
+                    f[i+1][1]=f[i][1];
+                    f[i+1][0]=f[i][0];
+                }
+                f[i+1][1]=a[1];
+                f[i+1][0]=a[0];
+            }
+            if(jstack < 0) {
+                break;
+            }
+            ir=istack[jstack--];
+            l=istack[jstack--];
+        } else {
+            k=(l+ir) >> 1;
+            SWAP_POLYVERT(f[k],f[l+1])
+            if(f[l][sortby] > f[ir][sortby]) {
+                SWAP_POLYVERT(f[l],f[ir])
+            }
+            if(f[l+1][sortby] > f[ir][sortby]) {
+                SWAP_POLYVERT(f[l+1],f[ir])
+            }
+            if(f[l][sortby] > f[l+1][sortby]) {
+                SWAP_POLYVERT(f[l],f[l+1])
+            }
+            i=l+1;
+            j=ir;
+            a[0]=f[l+1][0];
+            a[1]=f[l+1][1];
+            for(;;) {
+                do i++;
+                while(f[i][sortby] < a[sortby]);
+                do j--;
+                while(f[j][sortby] > a[sortby]);
+                if(j < i) {
+                    break;
+                }
+                SWAP_POLYVERT(f[i],f[j])
+            }
+            f[l+1][0]=f[j][0];
+            f[l+1][1]=f[j][1];
+            f[j][0]=a[0];
+            f[j][1]=a[1];
+            jstack+=2;
+            if(jstack > __PLX_SMALL_COUNT__) {
+                return;
+            }
+            if(ir-i+1 >= j-l) {
+                istack[jstack]=ir;
+                istack[jstack-1]=i;
+                ir=j-1;
+            } else {
+                istack[jstack]=j-1;
+                istack[jstack-1]=l;
+                l=i;
+            }
+        }
+    }
+}
+
+int plx_find_lower_bound(float v, float(*a)[2], int num_feather_verts) {
+    int x;
+    int l;
+    int r;
+    l=1;
+    r=num_feather_verts;
+    for(;;) {
+        // interpolation style search
+        //x=l+(v-a[l][1])*(r-l) / (a[r][1]-a[l][1]);
+
+        // binary search
+        x=(l+r) / 2;
+        if(v<a[x][1]) {
+            r=x-1;
+        } else {
+            l=x+1;
+        }
+        if((v>a[x-1][1] && v <= a[x][1]) || l>r) {
+            break;
+        }
+    }
+    if(v>a[x-1][1] && v <= a[x][1]) {
+        return x;
+    } else {
+        return num_feather_verts;
+    }
+}
+
+int plx_find_upper_bound(float v, float(*a)[2], int num_feather_verts) {
+    int x;
+    int l;
+    int r;
+    l=1;
+    r=num_feather_verts;
+    for(;;) {
+        // interpolation style search
+        //x=l+(v-a[l][1])*(r-l) / (a[r][1]-a[l][1]);
+
+        // binary search
+        x=(l+r) / 2;
+        if(v<a[x][1]) {
+            r=x-1;
+        } else {
+            l=x+1;
+        }
+        if((v>=a[x-1][1] && v < a[x][1]) || l>r) {
+            break;
+        }
+    }
+    if(v>=a[x-1][1] && v < a[x][1]) {
+        return x-1;
+    } else {
+        return num_feather_verts;
+    }
+}
diff --git a/intern/raskter/raskter.h b/intern/raskter/raskter.h
index e078b0d26be..676bb56a2ae 100644
--- a/intern/raskter/raskter.h
+++ b/intern/raskter/raskter.h
@@ -27,27 +27,64 @@
 /** \file raskter.h
  *  \ingroup RASKTER
  */
+/* from BLI_utildefines.h */
+#define MIN2(x, y)               ( (x) < (y) ? (x) : (y) )
+#define MAX2(x, y)               ( (x) > (y) ? (x) : (y) )
+#define ABS(a)                   ( (a) < 0 ? (-(a)) : (a) )
 
 struct poly_vert {
-	int x;
-	int y;
+    int x;
+    int y;
 };
 
 struct scan_line {
-	int xstart;
-	int xend;
+    int xstart;
+    int xend;
 };
 
 struct scan_line_batch {
-	int num;
-	int ystart;
-	struct scan_line *slines;
+    int num;
+    int ystart;
+    struct scan_line *slines;
+};
+
+struct e_status {
+    int x;
+    int ybeg;
+    int xshift;
+    int xdir;
+    int drift;
+    int drift_inc;
+    int drift_dec;
+    int num;
+    struct e_status *e_next;
+};
+
+struct r_buffer_stats {
+    float *buf;
+    int sizex;
+    int sizey;
+    int ymin;
+    int ymax;
+    int xmin;
+    int xmax;
+};
+
+struct r_fill_context {
+    struct e_status *all_edges, *possible_edges;
+    struct r_buffer_stats rb;
+    struct scan_line *bounds;
+    void *kdo;  //only used with kd tree
+    void *kdi;  //only used with kd tree
+    int *bound_indexes;
+    int bounds_length;
 };
 
 #ifdef __cplusplus
 extern "C" {
 #endif
 
+static void preprocess_all_edges(struct r_fill_context *ctx, struct poly_vert *verts, int num_verts, struct e_status *open_edge);
 int PLX_raskterize(float (*base_verts)[2], int num_base_verts,
                    float *buf, int buf_x, int buf_y, int do_mask_AA);
 int PLX_raskterize_feather(float (*base_verts)[2], int num_base_verts,
diff --git a/intern/raskter/raskter_mt.c b/intern/raskter/raskter_mt.c
new file mode 100644
index 00000000000..0b8f8c7c673
--- /dev/null
+++ b/intern/raskter/raskter_mt.c
@@ -0,0 +1,290 @@
+/*
+ * ***** 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) 2012 Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Peter Larabell.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+/** \file raskter_mt.c
+ *  \ingroup RASKTER
+ */
+#include <stdlib.h>
+
+#include "raskter.h"
+static int rast_scan_init(struct r_fill_context *ctx, struct poly_vert *verts, int num_verts) {
+    int x_curr;                 /* current pixel position in X */
+    int y_curr;                 /* current scan line being drawn */
+    int yp;                     /* y-pixel's position in frame buffer */
+    int swixd = 0;              /* whether or not edges switched position in X */
+    int i=0;                    /* counter */
+    float *cpxl;                /* pixel pointers... */
+    float *mpxl;
+    float *spxl;
+    struct e_status *e_curr;    /* edge pointers... */
+    struct e_status *e_temp;
+    struct e_status *edgbuf;
+    struct e_status **edgec;
+
+
+    if(num_verts < 3) {
+        return(1);
+    }
+
+    if((edgbuf = (struct e_status *)(malloc(sizeof(struct e_status) * num_verts))) == NULL) {
+        return(0);
+    }
+
+    /* set initial bounds length to 0 */
+    ctx->bounds_length=0;
+
+    /* round 1, count up all the possible spans in the base buffer */
+    preprocess_all_edges(ctx, verts, num_verts, edgbuf);
+
+    /* can happen with a zero area mask */
+    if (ctx->all_edges == NULL) {
+        free(edgbuf);
+        return(1);
+    }
+    ctx->possible_edges = NULL;
+
+    for(y_curr = ctx->all_edges->ybeg; (ctx->all_edges || ctx->possible_edges); y_curr++) {
+
+        for(edgec = &ctx->possible_edges; ctx->all_edges && (ctx->all_edges->ybeg == y_curr);) {
+            x_curr = ctx->all_edges->x;                  /* Set current X position. */
+            for(;;) {                                    /* Start looping edges. Will break when edges run out. */
+                e_curr = *edgec;                         /* Set up a current edge pointer. */
+                if(!e_curr || (e_curr->x >= x_curr)) {   /* If we have an no edge, or we need to skip some X-span, */
+                    e_temp = ctx->all_edges->e_next;     /* set a temp "next" edge to test. */
+                    *edgec = ctx->all_edges;             /* Add this edge to the list to be scanned. */
+                    ctx->all_edges->e_next = e_curr;     /* Set up the next edge. */
+                    edgec = &ctx->all_edges->e_next;     /* Set our list to the next edge's location in memory. */
+                    ctx->all_edges = e_temp;             /* Skip the NULL or bad X edge, set pointer to next edge. */
+                    break;                               /* Stop looping edges (since we ran out or hit empty X span. */
+                } else {
+                    edgec = &e_curr->e_next;             /* Set the pointer to the edge list the "next" edge. */
+                }
+            }
+        }
+
+        yp = y_curr * ctx->rb.sizex;
+        spxl = ctx->rb.buf + (yp);
+
+        for(e_curr = ctx->possible_edges; e_curr; e_curr = e_curr->e_next) {
+
+            /* set up xmin and xmax bounds on this scan line */
+            cpxl = spxl + MAX2(e_curr->x, 0);
+            e_curr = e_curr->e_next;
+            mpxl = spxl + MIN2(e_curr->x, ctx->rb.sizex) - 1;
+
+            if((y_curr >= 0) && (y_curr < ctx->rb.sizey)) {
+                ctx->bounds_length++;
+            }
+        }
+
+        for(edgec = &ctx->possible_edges; (e_curr = *edgec);) {
+            if(!(--(e_curr->num))) {
+                *edgec = e_curr->e_next;
+            } else {
+                e_curr->x += e_curr->xshift;
+                if((e_curr->drift += e_curr->drift_inc) > 0) {
+                    e_curr->x += e_curr->xdir;
+                    e_curr->drift -= e_curr->drift_dec;
+                }
+                edgec = &e_curr->e_next;
+            }
+        }
+        if(ctx->possible_edges) {
+            for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                /* if the current edge hits scan line at greater X than the next edge, we need to exchange the edges */
+                if(e_curr->x > e_curr->e_next->x) {
+                    *edgec = e_curr->e_next;
+                    /* exchange the pointers */
+                    e_temp = e_curr->e_next->e_next;
+                    e_curr->e_next->e_next = e_curr;
+                    e_curr->e_next = e_temp;
+                    /* set flag that we had at least one switch */
+                    swixd = 1;
+                }
+            }
+            /* if we did have a switch, look for more (there will more if there was one) */
+            for(;;) {
+                /* reset exchange flag so it's only set if we encounter another one */
+                swixd = 0;
+                for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                    /* again, if current edge hits scan line at higher X than next edge, exchange the edges and set flag */
+                    if(e_curr->x > e_curr->e_next->x) {
+                        *edgec = e_curr->e_next;
+                        /* exchange the pointers */
+                        e_temp = e_curr->e_next->e_next;
+                        e_curr->e_next->e_next = e_curr;
+                        e_curr->e_next = e_temp;
+                        /* flip the exchanged flag */
+                        swixd = 1;
+                    }
+                }
+                /* if we had no exchanges, we're done reshuffling the pointers */
+                if(!swixd) {
+                    break;
+                }
+            }
+        }
+    }
+
+
+/*initialize index buffer and bounds buffers*/
+    //gets the +1 for dummy at the end
+    if((ctx->bound_indexes = (int *)(malloc(sizeof(int) * ctx->rb.sizey+1)))==NULL) {
+        return(0);
+    }
+    //gets the +1 for dummy at the start
+    if((ctx->bounds = (struct scan_line *)(malloc(sizeof(struct scan_line) * ctx->bounds_length+1)))==NULL){
+        return(0);
+    }
+    //init all the indexes to zero (are they already zeroed from malloc???)
+    for(i=0;i<ctx->rb.sizey+1;i++){
+        ctx->bound_indexes[i]=0;
+    }
+    /* round 2, fill in the full list of bounds, and create indexes to the list... */
+    preprocess_all_edges(ctx, verts, num_verts, edgbuf);
+
+    /* can happen with a zero area mask */
+    if (ctx->all_edges == NULL) {
+        free(edgbuf);
+        return(1);
+    }
+    ctx->possible_edges = NULL;
+
+    /* restart i as a counter for total span placement in buffer */
+    i=1;
+    for(y_curr = ctx->all_edges->ybeg; (ctx->all_edges || ctx->possible_edges); y_curr++) {
+
+        for(edgec = &ctx->possible_edges; ctx->all_edges && (ctx->all_edges->ybeg == y_curr);) {
+            x_curr = ctx->all_edges->x;                  /* Set current X position. */
+            for(;;) {                                    /* Start looping edges. Will break when edges run out. */
+                e_curr = *edgec;                         /* Set up a current edge pointer. */
+                if(!e_curr || (e_curr->x >= x_curr)) {   /* If we have an no edge, or we need to skip some X-span, */
+                    e_temp = ctx->all_edges->e_next;     /* set a temp "next" edge to test. */
+                    *edgec = ctx->all_edges;             /* Add this edge to the list to be scanned. */
+                    ctx->all_edges->e_next = e_curr;     /* Set up the next edge. */
+                    edgec = &ctx->all_edges->e_next;     /* Set our list to the next edge's location in memory. */
+                    ctx->all_edges = e_temp;             /* Skip the NULL or bad X edge, set pointer to next edge. */
+                    break;                               /* Stop looping edges (since we ran out or hit empty X span. */
+                } else {
+                    edgec = &e_curr->e_next;             /* Set the pointer to the edge list the "next" edge. */
+                }
+            }
+        }
+
+        yp = y_curr * ctx->rb.sizex;
+        spxl = ctx->rb.buf + (yp);
+        if((y_curr >=0) && (y_curr < ctx->rb.sizey)){
+            ctx->bound_indexes[y_curr]=i;
+        }
+        for(e_curr = ctx->possible_edges; e_curr; e_curr = e_curr->e_next) {
+
+            /* set up xmin and xmax bounds on this scan line */
+            cpxl = spxl + MAX2(e_curr->x, 0);
+            e_curr = e_curr->e_next;
+            mpxl = spxl + MIN2(e_curr->x, ctx->rb.sizex) - 1;
+
+            if((y_curr >= 0) && (y_curr < ctx->rb.sizey)) {
+                ctx->bounds[i].xstart=cpxl-spxl;
+                ctx->bounds[i].xend=mpxl-spxl;
+                i++;
+            }
+        }
+
+        for(edgec = &ctx->possible_edges; (e_curr = *edgec);) {
+            if(!(--(e_curr->num))) {
+                *edgec = e_curr->e_next;
+            } else {
+                e_curr->x += e_curr->xshift;
+                if((e_curr->drift += e_curr->drift_inc) > 0) {
+                    e_curr->x += e_curr->xdir;
+                    e_curr->drift -= e_curr->drift_dec;
+                }
+                edgec = &e_curr->e_next;
+            }
+        }
+        if(ctx->possible_edges) {
+            for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                /* if the current edge hits scan line at greater X than the next edge, we need to exchange the edges */
+                if(e_curr->x > e_curr->e_next->x) {
+                    *edgec = e_curr->e_next;
+                    /* exchange the pointers */
+                    e_temp = e_curr->e_next->e_next;
+                    e_curr->e_next->e_next = e_curr;
+                    e_curr->e_next = e_temp;
+                    /* set flag that we had at least one switch */
+                    swixd = 1;
+                }
+            }
+            /* if we did have a switch, look for more (there will more if there was one) */
+            for(;;) {
+                /* reset exchange flag so it's only set if we encounter another one */
+                swixd = 0;
+                for(edgec = &ctx->possible_edges; (e_curr = *edgec)->e_next; edgec = &(*edgec)->e_next) {
+                    /* again, if current edge hits scan line at higher X than next edge, exchange the edges and set flag */
+                    if(e_curr->x > e_curr->e_next->x) {
+                        *edgec = e_curr->e_next;
+                        /* exchange the pointers */
+                        e_temp = e_curr->e_next->e_next;
+                        e_curr->e_next->e_next = e_curr;
+                        e_curr->e_next = e_temp;
+                        /* flip the exchanged flag */
+                        swixd = 1;
+                    }
+                }
+                /* if we had no exchanges, we're done reshuffling the pointers */
+                if(!swixd) {
+                    break;
+                }
+            }
+        }
+    }
+
+    free(edgbuf);
+    return 1;
+}
+
+static void init_base_data(float(*base_verts)[2], int num_base_verts,
+                   float *buf, int buf_x, int buf_y) {
+    int i;                                   /* i: Loop counter. */
+    struct poly_vert *ply;                   /* ply: Pointer to a list of integer buffer-space vertex coordinates. */
+    struct r_fill_context ctx = {0};
+    const float buf_x_f = (float)(buf_x);
+    const float buf_y_f = (float)(buf_y);
+    if((ply = (struct poly_vert *)(malloc(sizeof(struct poly_vert) * num_base_verts))) == NULL) {
+        return(0);
+    }
+    ctx.rb.buf = buf;                            /* Set the output buffer pointer. */
+    ctx.rb.sizex = buf_x;                        /* Set the output buffer size in X. (width) */
+    ctx.rb.sizey = buf_y;                        /* Set the output buffer size in Y. (height) */
+    for(i = 0; i < num_base_verts; i++) {                           /* Loop over all base_verts. */
+        ply[i].x = (int)((base_verts[i][0] * buf_x_f) + 0.5f);       /* Range expand normalized X to integer buffer-space X. */
+        ply[i].y = (int)((base_verts[i][1] * buf_y_f) + 0.5f); /* Range expand normalized Y to integer buffer-space Y. */
+    }
+    i = rast_scan_init(&ctx, ply, num_base_verts);  /* Call our rasterizer, passing in the integer coords for each vert. */
+    free(ply);                                      /* Free the memory allocated for the integer coordinate table. */
+    return(i);                                      /* Return the value returned by the rasterizer. */
+}
+