Fix for EWA (elliptical weighted average) sampling in the compositor.

EWA sampling is designed for downsampling images, i.e. scaling down the size of
input image pixels, which happens regularly in compositing. While the standard
sampling methods (linear, cubic) work reasonably well for linear
transformations, they don't yield good results in non-linear cases like
perspective projection or arbitrary displacement. EWA sampling is comparable to
mipmapping, but avoids problems with discontinuities.

To work correctly the EWA algorithm needs partial derivatives of the mapping
functions which convert output pixel coordinates back into the input image
space (2x2 Jacobian matrix). With these derivatives the EWA algorithm
projects ellipses into the input space and accumulates colors over their
area. This calculation was not done correctly in the compositor, only the
derivatives du/dx and dv/dy were calculation, basically this means it only
worked for non-rotated input images.

The patch introduces full derivative calculations du/dx, du/dy, dv/dx, dv/dy for
the 3 nodes which use EWA sampling currently: PlaneTrackWarp, MapUV and
Displace. In addition the calculation of ellipsis area and axis-aligned
bounding boxes has been fixed.

For the MapUV and Displace nodes the derivatives have to be estimated by
evaluating the UV/displacement inputs with 1-pixel offsets, which can still have
problems on discontinuities and sub-pixel variations. These potential problems
can only be alleviated by more radical design changes in the compositor
functions, which are out of scope for now. Basically the values passed to the
UV/Displacement inputs would need to be associated with their 1st order
derivatives, which requires a general approach to derivatives in all nodes.

11 files changed, 291 insertions, 276 deletions

diff --git a/source/blender/compositor/intern/COM_MemoryBuffer.cpp b/source/blender/compositor/intern/COM_MemoryBuffer.cpp
index f10e6696c6a..8a010561adf 100644
--- a/source/blender/compositor/intern/COM_MemoryBuffer.cpp
+++ b/source/blender/compositor/intern/COM_MemoryBuffer.cpp
@@ -214,136 +214,108 @@ static const float EWA_WTS[EWA_MAXIDX + 1] = {
 	0.00754159f, 0.00625989f, 0.00498819f, 0.00372644f, 0.00247454f, 0.00123242f, 0.f
 };
 
-static void radangle2imp(float a2, float b2, float th, float *A, float *B, float *C, float *F)
+static void ellipse_bounds(float A, float B, float C, float F, float &xmax, float &ymax)
 {
-	float ct2 = cosf(th);
-	const float st2 = 1.f - ct2 * ct2;    // <- sin(th)^2
-	ct2 *= ct2;
-	*A = a2 * st2 + b2 * ct2;
-	*B = (b2 - a2) * sinf(2.f * th);
-	*C = a2 * ct2 + b2 * st2;
-	*F = a2 * b2;
-}
-
-// all tests here are done to make sure possible overflows are hopefully minimized
-static void imp2radangle(float A, float B, float C, float F, float *a, float *b, float *th, float *ecc)
-{
-	if (F <= 1e-5f) {   // use arbitrary major radius, zero minor, infinite eccentricity
-		*a = sqrtf(A > C ? A : C);
-		*b = 0.f;
-		*ecc = 1e10f;
-		*th = 0.5f * (atan2f(B, A - C) + (float)M_PI);
+	float denom = 4.0f*A*C - B*B;
+	if (denom > 0.0f && A != 0.0f && C != 0.0f) {
+		xmax = sqrt(F)/(2.0f*A) * (sqrt(F*(4.0f*A - B*B/C)) + B*B*sqrt(F/(C*denom)));
+		ymax = sqrt(F)/(2.0f*C) * (sqrt(F*(4.0f*C - B*B/A)) + B*B*sqrt(F/(A*denom)));
 	}
 	else {
-		const float AmC = A - C, ApC = A + C, F2 = F * 2.f;
-		const float r = sqrtf(AmC * AmC + B * B);
-		float d = ApC - r;
-		*a = (d <= 0.f) ? sqrtf(A > C ? A : C) : sqrtf(F2 / d);
-		d = ApC + r;
-		if (d <= 0.f) {
-			*b = 0.f;
-			*ecc = 1e10f;
-		}
-		else {
-			*b = sqrtf(F2 / d);
-			*ecc = *a / *b;
-		}
-		/* incr theta by 0.5 * pi (angle of major axis) */
-		*th = 0.5f * (atan2f(B, AmC) + (float)M_PI);
+		xmax = 0.0f;
+		ymax = 0.0f;
 	}
 }
 
-static float clipuv(float x, float limit)
+static void ellipse_params(float Ux, float Uy, float Vx, float Vy, float &A, float &B, float &C, float &F, float &umax, float &vmax)
 {
-	x = (x < 0) ? 0 : ((x >= limit) ? (limit - 1) : x);
-	return x;
+	A = Vx*Vx + Vy*Vy;
+	B = -2.0f * (Ux*Vx + Uy*Vy);
+	C = Ux*Ux + Uy*Uy;
+	F = A*C - B*B * 0.25f;
+
+	float factor = (F != 0.0f ? (float)(EWA_MAXIDX+1) / F : 0.0f);
+	A *= factor;
+	B *= factor;
+	C *= factor;
+	F = (float)(EWA_MAXIDX+1);
+
+	ellipse_bounds(A, B, C, sqrtf(F), umax, vmax);
 }
 
 /**
- * \note \a sampler at the moment is either 'COM_PS_NEAREST' or not, other values won't matter.
+ * Filtering method based on
+ * "Creating raster omnimax images from multiple perspective views using the elliptical weighted average filter"
+ * by Ned Greene and Paul S. Heckbert (1986)
  */
-void MemoryBuffer::readEWA(float result[4], float fx, float fy, float dx, float dy, PixelSampler sampler)
+void MemoryBuffer::readEWA(float result[4], const float uv[2], const float derivatives[2][2], PixelSampler sampler)
 {
-	const int width = this->getWidth(), height = this->getHeight();
-	
-	// scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values,
-	// scaling by aspect ratio alone does the opposite, so try something in between instead...
-	const float ff2 = width, ff = sqrtf(ff2), q = height / ff;
-	const float Ux = dx * ff, Vx = dx * q, Uy = dy * ff, Vy = dy * q;
-	float A = Vx * Vx + Vy * Vy;
-	float B = -2.f * (Ux * Vx + Uy * Vy);
-	float C = Ux * Ux + Uy * Uy;
-	float F = A * C - B * B * 0.25f;
-	float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d;
-	int u, v, u1, u2, v1, v2;
-	// The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C,
-	// so the ellipse always covers at least some texels. But since the filter is now always larger,
-	// it also means that everywhere else it's also more blurry then ideally should be the case.
-	// So instead here the ellipse radii are modified instead whenever either is too low.
-	// Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off,
-	// and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on
-	// (minimum values: const float rmin = intpol ? 1.f : 0.5f;)
-
-	/* note: 0.765625f is too sharp, 1.0 will not blur with an exact pixel sample
-	 * useful to avoid blurring when there is no distortion */
-#if 0
-	const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 0.765625f) / ff2;
-#else
-	const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 1.0f     ) / ff2;
-#endif
-	imp2radangle(A, B, C, F, &a, &b, &th, &ecc);
-	if ((b2 = b * b) < rmin) {
-		if ((a2 = a * a) < rmin) {
-			B = 0.f;
-			A = C = rmin;
-			F = A * C;
-		}
-		else {
-			b2 = rmin;
-			radangle2imp(a2, b2, th, &A, &B, &C, &F);
-		}
-	}
+	zero_v4(result);
+	int width = this->getWidth(), height = this->getHeight();
+	if (width == 0 || height == 0)
+		return;
 
-	ue = ff * sqrtf(C);
-	ve = ff * sqrtf(A);
-	d = (float)(EWA_MAXIDX + 1) / (F * ff2);
-	A *= d;
-	B *= d;
-	C *= d;
-
-	U0 = fx;
-	V0 = fy;
-	u1 = (int)(floorf(U0 - ue));
-	u2 = (int)(ceilf(U0 + ue));
-	v1 = (int)(floorf(V0 - ve));
-	v2 = (int)(ceilf(V0 + ve));
-	U0 -= 0.5f;
-	V0 -= 0.5f;
-	DDQ = 2.f * A;
-	U = u1 - U0;
-	ac1 = A * (2.f * U + 1.f);
-	ac2 = A * U * U;
-	BU = B * U;
-
-	d = result[0] = result[1] = result[2] = result[3] = 0.f;
-	for (v = v1; v <= v2; ++v) {
-		const float V = v - V0;
-		float DQ = ac1 + B * V;
-		float Q = (C * V + BU) * V + ac2;
-		for (u = u1; u <= u2; ++u) {
-			if (Q < (float)(EWA_MAXIDX + 1)) {
+	float u = uv[0], v = uv[1];
+	float Ux = derivatives[0][0], Vx = derivatives[1][0], Uy = derivatives[0][1], Vy = derivatives[1][1];
+	float A, B, C, F, ue, ve;
+	ellipse_params(Ux, Uy, Vx, Vy, A, B, C, F, ue, ve);
+
+	/* Note: highly eccentric ellipses can lead to large texture space areas to filter!
+	 * This is limited somewhat by the EWA_WTS size in the loop, but a nicer approach
+	 * could be the one found in
+	 * "High Quality Elliptical Texture Filtering on GPU"
+	 * by Pavlos Mavridis and Georgios Papaioannou
+	 * in which the eccentricity of the ellipse is clamped.
+	 */
+
+	int U0 = (int)u;
+	int V0 = (int)v;
+	/* pixel offset for interpolation */
+	float ufac = u - floor(u), vfac = v - floor(v);
+	/* filter size */
+	int u1 = (int)(u - ue);
+	int u2 = (int)(u + ue);
+	int v1 = (int)(v - ve);
+	int v2 = (int)(v + ve);
+
+	/* sane clamping to avoid unnecessarily huge loops */
+	/* note: if eccentricity gets clamped (see above),
+	 * the ue/ve limits can also be lowered accordingly
+	 */
+	if (U0-u1 > EWA_MAXIDX) u1 = U0 - EWA_MAXIDX;
+	if (u2-U0 > EWA_MAXIDX) u2 = U0 + EWA_MAXIDX;
+	if (V0-v1 > EWA_MAXIDX) v1 = V0 - EWA_MAXIDX;
+	if (v2-V0 > EWA_MAXIDX) v2 = V0 + EWA_MAXIDX;
+
+	float DDQ = 2.f * A;
+	float U = u1 - U0;
+	float ac1 = A * (2.f*U + 1.f);
+	float ac2 = A * U*U;
+	float BU = B * U;
+
+	float sum = 0.0f;
+	for (int v = v1; v <= v2; ++v) {
+		float V = v - V0;
+
+		float DQ = ac1 + B*V;
+		float Q = (C*V + BU)*V + ac2;
+		for (int u = u1; u <= u2; ++u) {
+			if (Q < F) {
 				float tc[4];
-				const float wt = EWA_WTS[(Q < 0.f) ? 0 : (unsigned int)Q];
-				read(tc, clipuv(u, width), clipuv(v, height));
+				const float wt = EWA_WTS[CLAMPIS((int)Q, 0, EWA_MAXIDX)];
+				switch (sampler) {
+					case COM_PS_NEAREST: read(tc, u, v); break;
+					case COM_PS_BILINEAR: readBilinear(tc, (float)u+ufac, (float)v+vfac); break;
+					case COM_PS_BICUBIC: readBilinear(tc, (float)u+ufac, (float)v+vfac); break; /* XXX no readBicubic method yet */
+					default: zero_v4(tc); break;
+				}
 				madd_v4_v4fl(result, tc, wt);
-				d += wt;
+				sum += wt;
 			}
 			Q += DQ;
 			DQ += DDQ;
 		}
 	}
 	
-	// d should hopefully never be zero anymore
-	d = 1.f / d;
-	mul_v4_fl(result, d);
+	mul_v4_fl(result, (sum != 0.0f ? 1.0f / sum : 0.0f));
 }
diff --git a/source/blender/compositor/intern/COM_MemoryBuffer.h b/source/blender/compositor/intern/COM_MemoryBuffer.h
index 548744f6660..f9cde79c07f 100644
--- a/source/blender/compositor/intern/COM_MemoryBuffer.h
+++ b/source/blender/compositor/intern/COM_MemoryBuffer.h
@@ -245,7 +245,7 @@ public:
 		result[3] = color1[3] * mvaluex + color3[3] * valuex;
 	}
 
-	void readEWA(float result[4], float fx, float fy, float dx, float dy, PixelSampler sampler);
+	void readEWA(float result[4], const float uv[2], const float derivatives[2][2], PixelSampler sampler);
 	
 	/**
 	 * @brief is this MemoryBuffer a temporarily buffer (based on an area, not on a chunk)
diff --git a/source/blender/compositor/intern/COM_SocketReader.h b/source/blender/compositor/intern/COM_SocketReader.h
index 2168611d0c0..c996ef5bbeb 100644
--- a/source/blender/compositor/intern/COM_SocketReader.h
+++ b/source/blender/compositor/intern/COM_SocketReader.h
@@ -88,7 +88,7 @@ protected:
 	 * @param dy
 	 * @param inputBuffers chunks that can be read by their ReadBufferOperation.
 	 */
-	virtual void executePixelFiltered(float output[4], float x, float y, float dx, float dy, PixelSampler sampler) {}
+	virtual void executePixelFiltered(float output[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler) {}
 
 public:
 	inline void readSampled(float result[4], float x, float y, PixelSampler sampler) {
@@ -97,7 +97,7 @@ public:
 	inline void read(float result[4], int x, int y, void *chunkData) {
 		executePixel(result, x, y, chunkData);
 	}
-	inline void readFiltered(float result[4], float x, float y, float dx, float dy, PixelSampler sampler) {
+	inline void readFiltered(float result[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler) {
 		executePixelFiltered(result, x, y, dx, dy, sampler);
 	}
 
diff --git a/source/blender/compositor/operations/COM_DisplaceOperation.cpp b/source/blender/compositor/operations/COM_DisplaceOperation.cpp
index 17692bd18ef..96031e7acc4 100644
--- a/source/blender/compositor/operations/COM_DisplaceOperation.cpp
+++ b/source/blender/compositor/operations/COM_DisplaceOperation.cpp
@@ -49,54 +49,92 @@ void DisplaceOperation::initExecution()
 	this->m_height_x4 = this->getHeight() * 4;
 }
 
-
-/* minimum distance (in pixels) a pixel has to be displaced
- * in order to take effect */
-#define DISPLACE_EPSILON    0.01f
-
-void DisplaceOperation::executePixel(float output[4], int x, int y, void *data)
+void DisplaceOperation::executePixelSampled(float output[4], float x, float y, PixelSampler sampler)
 {
-	float inVector[4];
-	float inScale[4];
+	float xy[2] = { x, y };
+	float uv[2], deriv[2][2];
 
-	float p_dx, p_dy;   /* main displacement in pixel space */
-	float d_dx, d_dy;
-	float dxt, dyt;
-	float u, v;
+	pixelTransform(xy, uv, deriv);
 
-	this->m_inputScaleXProgram->readSampled(inScale, x, y, COM_PS_NEAREST);
-	float xs = inScale[0];
-	this->m_inputScaleYProgram->readSampled(inScale, x, y, COM_PS_NEAREST);
-	float ys = inScale[0];
+	/* EWA filtering (without nearest it gets blurry with NO distortion) */
+	this->m_inputColorProgram->readFiltered(output, uv[0], uv[1], deriv[0], deriv[1], COM_PS_BILINEAR);
+}
 
+bool DisplaceOperation::read_displacement(float x, float y, float xscale, float yscale, const float origin[2], float &r_u, float &r_v)
+{
+	float width = m_inputVectorProgram->getWidth();
+	float height = m_inputVectorProgram->getHeight();
+	if (x < 0.0f || x >= width || y < 0.0f || y >= height) {
+		r_u = 0.0f;
+		r_v = 0.0f;
+		return false;
+	}
+	else {
+		float col[4];
+		m_inputVectorProgram->readSampled(col, x, y, COM_PS_BILINEAR);
+		r_u = origin[0] - col[0]*xscale + 0.5f;
+		r_v = origin[1] - col[1]*yscale + 0.5f;
+		return true;
+	}
+}
+
+void DisplaceOperation::pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2])
+{
+	float col[4];
+	float uv[2]; /* temporary variables for derivative estimation */
+	int num;
+
+	m_inputScaleXProgram->readSampled(col, xy[0], xy[1], COM_PS_NEAREST);
+	float xs = col[0];
+	m_inputScaleYProgram->readSampled(col, xy[0], xy[1], COM_PS_NEAREST);
+	float ys = col[0];
 	/* clamp x and y displacement to triple image resolution - 
 	 * to prevent hangs from huge values mistakenly plugged in eg. z buffers */
-	CLAMP(xs, -this->m_width_x4, this->m_width_x4);
-	CLAMP(ys, -this->m_height_x4, this->m_height_x4);
-
-	this->m_inputVectorProgram->readSampled(inVector, x, y, COM_PS_NEAREST);
-	p_dx = inVector[0] * xs;
-	p_dy = inVector[1] * ys;
+	CLAMP(xs, -m_width_x4, m_width_x4);
+	CLAMP(ys, -m_height_x4, m_height_x4);
 
 	/* displaced pixel in uv coords, for image sampling */
-	u = x - p_dx + 0.5f;
-	v = y - p_dy + 0.5f;
-
-	/* calc derivatives */
-	this->m_inputVectorProgram->readSampled(inVector, x + 1, y, COM_PS_NEAREST);
-	d_dx = inVector[0] * xs;
-	this->m_inputVectorProgram->readSampled(inVector, x, y + 1, COM_PS_NEAREST);
-	d_dy = inVector[1] * ys;
-
-	/* clamp derivatives to minimum displacement distance in UV space */
-	dxt = p_dx - d_dx;
-	dyt = p_dy - d_dy;
-
-	dxt = signf(dxt) * max_ff(fabsf(dxt), DISPLACE_EPSILON) / this->getWidth();
-	dyt = signf(dyt) * max_ff(fabsf(dyt), DISPLACE_EPSILON) / this->getHeight();
+	read_displacement(xy[0], xy[1], xs, ys, xy, r_uv[0], r_uv[1]);
+
+	/* Estimate partial derivatives using 1-pixel offsets */
+	const float epsilon[2] = { 1.0f, 1.0f };
+	
+	zero_v2(r_deriv[0]);
+	zero_v2(r_deriv[1]);
+
+	num = 0;
+	if (read_displacement(xy[0] + epsilon[0], xy[1], xs, ys, xy, uv[0], uv[1])) {
+		r_deriv[0][0] += uv[0] - r_uv[0];
+		r_deriv[1][0] += uv[1] - r_uv[1];
+		++num;
+	}
+	if (read_displacement(xy[0] - epsilon[0], xy[1], xs, ys, xy, uv[0], uv[1])) {
+		r_deriv[0][0] += r_uv[0] - uv[0];
+		r_deriv[1][0] += r_uv[1] - uv[1];
+		++num;
+	}
+	if (num > 0) {
+		float numinv = 1.0f / (float)num;
+		r_deriv[0][0] *= numinv;
+		r_deriv[1][0] *= numinv;
+	}
 
-	/* EWA filtering (without nearest it gets blurry with NO distortion) */
-	this->m_inputColorProgram->readFiltered(output, u, v, dxt, dyt, COM_PS_NEAREST);
+	num = 0;
+	if (read_displacement(xy[0], xy[1] + epsilon[1], xs, ys, xy, uv[0], uv[1])) {
+		r_deriv[0][1] += uv[0] - r_uv[0];
+		r_deriv[1][1] += uv[1] - r_uv[1];
+		++num;
+	}
+	if (read_displacement(xy[0], xy[1] - epsilon[1], xs, ys, xy, uv[0], uv[1])) {
+		r_deriv[0][1] += r_uv[0] - uv[0];
+		r_deriv[1][1] += r_uv[1] - uv[1];
+		++num;
+	}
+	if (num > 0) {
+		float numinv = 1.0f / (float)num;
+		r_deriv[0][1] *= numinv;
+		r_deriv[1][1] *= numinv;
+	}
 }
 
 void DisplaceOperation::deinitExecution()
@@ -126,10 +164,10 @@ bool DisplaceOperation::determineDependingAreaOfInterest(rcti *input, ReadBuffer
 
 	/* vector */
 	operation = getInputOperation(1);
-	vectorInput.xmax = input->xmax + 2;
-	vectorInput.xmin = input->xmin;
-	vectorInput.ymax = input->ymax + 2;
-	vectorInput.ymin = input->ymin;
+	vectorInput.xmax = input->xmax + 1;
+	vectorInput.xmin = input->xmin - 1;
+	vectorInput.ymax = input->ymax + 1;
+	vectorInput.ymin = input->ymin - 1;
 	if (operation->determineDependingAreaOfInterest(&vectorInput, readOperation, output)) {
 		return true;
 	}
diff --git a/source/blender/compositor/operations/COM_DisplaceOperation.h b/source/blender/compositor/operations/COM_DisplaceOperation.h
index daf54294aa1..cec7937d9d6 100644
--- a/source/blender/compositor/operations/COM_DisplaceOperation.h
+++ b/source/blender/compositor/operations/COM_DisplaceOperation.h
@@ -48,7 +48,9 @@ public:
 	/**
 	 * the inner loop of this program
 	 */
-	void executePixel(float output[4], int x, int y, void *data);
+	void executePixelSampled(float output[4], float x, float y, PixelSampler sampler);
+	
+	void pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2]);
 	
 	/**
 	 * Initialize the execution
@@ -59,5 +61,8 @@ public:
 	 * Deinitialize the execution
 	 */
 	void deinitExecution();
+
+private:
+	bool read_displacement(float x, float y, float xscale, float yscale, const float origin[2], float &r_u, float &r_v);
 };
 #endif
diff --git a/source/blender/compositor/operations/COM_MapUVOperation.cpp b/source/blender/compositor/operations/COM_MapUVOperation.cpp
index 289b447dec7..292f073548a 100644
--- a/source/blender/compositor/operations/COM_MapUVOperation.cpp
+++ b/source/blender/compositor/operations/COM_MapUVOperation.cpp
@@ -42,46 +42,28 @@ void MapUVOperation::initExecution()
 
 void MapUVOperation::executePixelSampled(float output[4], float x, float y, PixelSampler sampler)
 {
-	float inputUV[4];
-	float uv_a[4], uv_b[4];
-	float u, v;
+	float xy[2] = { x, y };
+	float uv[2], deriv[2][2], alpha;
 
-	float dx, dy;
-	float uv_l, uv_r;
-	float uv_u, uv_d;
-
-	this->m_inputUVProgram->readSampled(inputUV, x, y, sampler);
-	if (inputUV[2] == 0.f) {
+	pixelTransform(xy, uv, deriv, alpha);
+	if (alpha == 0.0f) {
 		zero_v4(output);
 		return;
 	}
-	/* adaptive sampling, red (U) channel */
-	this->m_inputUVProgram->readSampled(uv_a, x - 1, y, COM_PS_NEAREST);
-	this->m_inputUVProgram->readSampled(uv_b, x + 1, y, COM_PS_NEAREST);
-	uv_l = uv_a[2] != 0.f ? fabsf(inputUV[0] - uv_a[0]) : 0.f;
-	uv_r = uv_b[2] != 0.f ? fabsf(inputUV[0] - uv_b[0]) : 0.f;
-
-	dx = 0.5f * (uv_l + uv_r);
-
-	/* adaptive sampling, green (V) channel */
-	this->m_inputUVProgram->readSampled(uv_a, x, y - 1, COM_PS_NEAREST);
-	this->m_inputUVProgram->readSampled(uv_b, x, y + 1, COM_PS_NEAREST);
-	uv_u = uv_a[2] != 0.f ? fabsf(inputUV[1] - uv_a[1]) : 0.f;
-	uv_d = uv_b[2] != 0.f ? fabsf(inputUV[1] - uv_b[1]) : 0.f;
-
-	dy = 0.5f * (uv_u + uv_d);
 
+	/* EWA filtering */
+	this->m_inputColorProgram->readFiltered(output, uv[0], uv[1], deriv[0], deriv[1], COM_PS_BILINEAR);
+	
 	/* UV to alpha threshold */
 	const float threshold = this->m_alpha * 0.05f;
-	float alpha = 1.0f - threshold * (dx + dy);
-	if (alpha < 0.f) alpha = 0.f;
-	else alpha *= inputUV[2];
-
-	/* EWA filtering */
-	u = inputUV[0] * this->m_inputColorProgram->getWidth();
-	v = inputUV[1] * this->m_inputColorProgram->getHeight();
-
-	this->m_inputColorProgram->readFiltered(output, u, v, dx, dy, COM_PS_NEAREST);
+	/* XXX alpha threshold is used to fade out pixels on boundaries with invalid derivatives.
+	 * this calculation is not very well defined, should be looked into if it becomes a problem ...
+	 */
+	float du = len_v2(deriv[0]);
+	float dv = len_v2(deriv[1]);
+	float factor = 1.0f - threshold * (du + dv);
+	if (factor < 0.f) alpha = 0.f;
+	else alpha *= factor;
 
 	/* "premul" */
 	if (alpha < 1.0f) {
@@ -89,6 +71,74 @@ void MapUVOperation::executePixelSampled(float output[4], float x, float y, Pixe
 	}
 }
 
+bool MapUVOperation::read_uv(float x, float y, float &r_u, float &r_v, float &r_alpha)
+{
+	float width = m_inputUVProgram->getWidth();
+	float height = m_inputUVProgram->getHeight();
+	if (x < 0.0f || x >= width || y < 0.0f || y >= height) {
+		r_u = 0.0f;
+		r_v = 0.0f;
+		r_alpha = 0.0f;
+		return false;
+	}
+	else {
+		float col[4];
+		m_inputUVProgram->readSampled(col, x, y, COM_PS_BILINEAR);
+		r_u = col[0] * width;
+		r_v = col[1] * height;
+		r_alpha = col[2];
+		return true;
+	}
+}
+
+void MapUVOperation::pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2], float &r_alpha)
+{
+	float uv[2], alpha; /* temporary variables for derivative estimation */
+	int num;
+
+	read_uv(xy[0], xy[1], r_uv[0], r_uv[1], r_alpha);
+
+	/* Estimate partial derivatives using 1-pixel offsets */
+	const float epsilon[2] = { 1.0f, 1.0f };
+
+	zero_v2(r_deriv[0]);
+	zero_v2(r_deriv[1]);
+
+	num = 0;
+	if (read_uv(xy[0] + epsilon[0], xy[1], uv[0], uv[1], alpha)) {
+		r_deriv[0][0] += uv[0] - r_uv[0];
+		r_deriv[1][0] += uv[1] - r_uv[1];
+		++num;
+	}
+	if (read_uv(xy[0] - epsilon[0], xy[1], uv[0], uv[1], alpha)) {
+		r_deriv[0][0] += r_uv[0] - uv[0];
+		r_deriv[1][0] += r_uv[1] - uv[1];
+		++num;
+	}
+	if (num > 0) {
+		float numinv = 1.0f / (float)num;
+		r_deriv[0][0] *= numinv;
+		r_deriv[1][0] *= numinv;
+	}
+
+	num = 0;
+	if (read_uv(xy[0], xy[1] + epsilon[1], uv[0], uv[1], alpha)) {
+		r_deriv[0][1] += uv[0] - r_uv[0];
+		r_deriv[1][1] += uv[1] - r_uv[1];
+		++num;
+	}
+	if (read_uv(xy[0], xy[1] - epsilon[1], uv[0], uv[1], alpha)) {
+		r_deriv[0][1] += r_uv[0] - uv[0];
+		r_deriv[1][1] += r_uv[1] - uv[1];
+		++num;
+	}
+	if (num > 0) {
+		float numinv = 1.0f / (float)num;
+		r_deriv[0][1] *= numinv;
+		r_deriv[1][1] *= numinv;
+	}
+}
+
 void MapUVOperation::deinitExecution()
 {
 	this->m_inputUVProgram = NULL;
diff --git a/source/blender/compositor/operations/COM_MapUVOperation.h b/source/blender/compositor/operations/COM_MapUVOperation.h
index fe8bfd2a9ac..796ee952607 100644
--- a/source/blender/compositor/operations/COM_MapUVOperation.h
+++ b/source/blender/compositor/operations/COM_MapUVOperation.h
@@ -46,7 +46,9 @@ public:
 	 * the inner loop of this program
 	 */
 	void executePixelSampled(float output[4], float x, float y, PixelSampler sampler);
-	
+
+	void pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2], float &r_alpha);
+
 	/**
 	 * Initialize the execution
 	 */
@@ -58,5 +60,9 @@ public:
 	void deinitExecution();
 
 	void setAlpha(float alpha) { this->m_alpha = alpha; }
+
+private:
+	bool read_uv(float x, float y, float &r_u, float &r_v, float &r_alpha);
 };
+
 #endif
diff --git a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp
index dba3a6f3505..7780023c465 100644
--- a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp
+++ b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp
@@ -36,60 +36,17 @@ extern "C" {
 	#include "BKE_tracking.h"
 }
 
-BLI_INLINE bool isPointInsideQuad(const float x, const float y, const float corners[4][2])
-{
-	float point[2];
-
-	point[0] = x;
-	point[1] = y;
-
-	return isect_point_tri_v2(point, corners[0], corners[1], corners[2]) ||
-	       isect_point_tri_v2(point, corners[0], corners[2], corners[3]);
-}
-
-BLI_INLINE void warpCoord(float x, float y, float matrix[3][3], float uv[2])
+BLI_INLINE void warpCoord(float x, float y, float matrix[3][3], float uv[2], float deriv[2][2])
 {
 	float vec[3] = {x, y, 1.0f};
 	mul_m3_v3(matrix, vec);
-	vec[0] /= vec[2];
-	vec[1] /= vec[2];
+	uv[0] = vec[0] / vec[2];
+	uv[1] = vec[1] / vec[2];
 
-	copy_v2_v2(uv, vec);
-}
-
-BLI_INLINE void resolveUVAndDxDy(const float x, const float y, float matrix[3][3],
-                                 float *u_r, float *v_r, float *dx_r, float *dy_r)
-{
-	float inputUV[2];
-	float uv_a[2], uv_b[2];
-
-	float dx, dy;
-	float uv_l, uv_r;
-	float uv_u, uv_d;
-
-	warpCoord(x, y, matrix, inputUV);
-
-	/* adaptive sampling, red (U) channel */
-	warpCoord(x - 1, y, matrix, uv_a);
-	warpCoord(x + 1, y, matrix, uv_b);
-	uv_l = fabsf(inputUV[0] - uv_a[0]);
-	uv_r = fabsf(inputUV[0] - uv_b[0]);
-
-	dx = 0.5f * (uv_l + uv_r);
-
-	/* adaptive sampling, green (V) channel */
-	warpCoord(x, y - 1, matrix, uv_a);
-	warpCoord(x, y + 1, matrix, uv_b);
-	uv_u = fabsf(inputUV[1] - uv_a[1]);
-	uv_d = fabsf(inputUV[1] - uv_b[1]);
-
-	dy = 0.5f * (uv_u + uv_d);
-
-	*dx_r = dx;
-	*dy_r = dy;
-
-	*u_r = inputUV[0];
-	*v_r = inputUV[1];
+	deriv[0][0] = (matrix[0][0] - matrix[0][2] * uv[0]) / vec[2];
+	deriv[1][0] = (matrix[0][1] - matrix[0][2] * uv[1]) / vec[2];
+	deriv[0][1] = (matrix[1][0] - matrix[1][2] * uv[0]) / vec[2];
+	deriv[1][1] = (matrix[1][1] - matrix[1][2] * uv[1]) / vec[2];
 }
 
 PlaneTrackWarpImageOperation::PlaneTrackWarpImageOperation() : PlaneTrackCommonOperation()
@@ -98,9 +55,6 @@ PlaneTrackWarpImageOperation::PlaneTrackWarpImageOperation() : PlaneTrackCommonO
 	this->addOutputSocket(COM_DT_COLOR);
 	this->m_pixelReader = NULL;
 	this->setComplex(true);
-
-	/* Currently hardcoded to 8 samples. */
-	this->m_osa = 8;
 }
 
 void PlaneTrackWarpImageOperation::initExecution()
@@ -109,8 +63,6 @@ void PlaneTrackWarpImageOperation::initExecution()
 
 	this->m_pixelReader = this->getInputSocketReader(0);
 
-	BLI_jitter_init(this->m_jitter[0], this->m_osa);
-
 	const int width = this->m_pixelReader->getWidth();
 	const int height = this->m_pixelReader->getHeight();
 	float frame_corners[4][2] = {{0.0f, 0.0f},
@@ -127,41 +79,32 @@ void PlaneTrackWarpImageOperation::deinitExecution()
 	this->m_pixelReader = NULL;
 }
 
-void PlaneTrackWarpImageOperation::executePixelSampled(float output[4], float x_, float y_, PixelSampler sampler)
+void PlaneTrackWarpImageOperation::executePixelSampled(float output[4], float x, float y, PixelSampler sampler)
 {
-	float color_accum[4];
+	float xy[2] = {x, y};
+	float uv[2];
+	float deriv[2][2];
 
-	zero_v4(color_accum);
-	for (int sample = 0; sample < this->m_osa; sample++) {
-		float current_x = x_ + this->m_jitter[sample][0],
-		      current_y = y_ + this->m_jitter[sample][1];
-		if (isPointInsideQuad(current_x, current_y, this->m_frameSpaceCorners)) {
-			float current_color[4];
-			float u, v, dx, dy;
+	pixelTransform(xy, uv, deriv);
 
-			resolveUVAndDxDy(current_x, current_y, m_perspectiveMatrix, &u, &v, &dx, &dy);
-
-			/* derivatives are to be in normalized space.. */
-			dx /= this->m_pixelReader->getWidth();
-			dy /= this->m_pixelReader->getHeight();
-
-			this->m_pixelReader->readFiltered(current_color, u, v, dx, dy, COM_PS_NEAREST);
-			add_v4_v4(color_accum, current_color);
-		}
-	}
+	m_pixelReader->readFiltered(output, uv[0], uv[1], deriv[0], deriv[1], COM_PS_BILINEAR);
+}
 
-	mul_v4_v4fl(output, color_accum, 1.0f / this->m_osa);
+void PlaneTrackWarpImageOperation::pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2])
+{
+	warpCoord(xy[0], xy[1], m_perspectiveMatrix, r_uv, r_deriv);
 }
 
 bool PlaneTrackWarpImageOperation::determineDependingAreaOfInterest(rcti *input, ReadBufferOperation *readOperation, rcti *output)
 {
 	float UVs[4][2];
+	float deriv[2][2];
 
 	/* TODO(sergey): figure out proper way to do this. */
-	warpCoord(input->xmin - 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[0]);
-	warpCoord(input->xmax + 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[1]);
-	warpCoord(input->xmax + 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[2]);
-	warpCoord(input->xmin - 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[3]);
+	warpCoord(input->xmin - 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[0], deriv);
+	warpCoord(input->xmax + 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[1], deriv);
+	warpCoord(input->xmax + 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[2], deriv);
+	warpCoord(input->xmin - 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[3], deriv);
 
 	float min[2], max[2];
 	INIT_MINMAX2(min, max);
diff --git a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h
index ed9cc2fe5aa..ceb002c8039 100644
--- a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h
+++ b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h
@@ -36,8 +36,6 @@
 class PlaneTrackWarpImageOperation : public PlaneTrackCommonOperation {
 protected:
 	SocketReader *m_pixelReader;
-	int m_osa;
-	float m_jitter[32][2];
 	float m_perspectiveMatrix[3][3];
 
 public:
@@ -47,6 +45,7 @@ public:
 	void deinitExecution();
 
 	void executePixelSampled(float output[4], float x, float y, PixelSampler sampler);
+	void pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2]);
 
 	bool determineDependingAreaOfInterest(rcti *input, ReadBufferOperation *readOperation, rcti *output);
 };
diff --git a/source/blender/compositor/operations/COM_ReadBufferOperation.cpp b/source/blender/compositor/operations/COM_ReadBufferOperation.cpp
index b93e338ad6c..47d5fc6bcca 100644
--- a/source/blender/compositor/operations/COM_ReadBufferOperation.cpp
+++ b/source/blender/compositor/operations/COM_ReadBufferOperation.cpp
@@ -81,14 +81,16 @@ void ReadBufferOperation::executePixelExtend(float output[4], float x, float y,
 	}
 }
 
-void ReadBufferOperation::executePixelFiltered(float output[4], float x, float y, float dx, float dy, PixelSampler sampler)
+void ReadBufferOperation::executePixelFiltered(float output[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler)
 {
 	if (m_single_value) {
 		/* write buffer has a single value stored at (0,0) */
 		m_buffer->read(output, 0, 0);
 	}
 	else {
-		m_buffer->readEWA(output, x, y, dx, dy, sampler);
+		const float uv[2] = { x, y };
+		const float deriv[2][2] = { {dx[0], dx[1]}, {dy[0], dy[1]} };
+		m_buffer->readEWA(output, uv, deriv, sampler);
 	}
 }
 
diff --git a/source/blender/compositor/operations/COM_ReadBufferOperation.h b/source/blender/compositor/operations/COM_ReadBufferOperation.h
index cdc6e50f6d4..cce519c6eb3 100644
--- a/source/blender/compositor/operations/COM_ReadBufferOperation.h
+++ b/source/blender/compositor/operations/COM_ReadBufferOperation.h
@@ -43,7 +43,7 @@ public:
 	void executePixelSampled(float output[4], float x, float y, PixelSampler sampler);
 	void executePixelExtend(float output[4], float x, float y, PixelSampler sampler,
 	                        MemoryBufferExtend extend_x, MemoryBufferExtend extend_y);
-	void executePixelFiltered(float output[4], float x, float y, float dx, float dy, PixelSampler sampler);
+	void executePixelFiltered(float output[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler);
 	const bool isReadBufferOperation() const { return true; }
 	void setOffset(unsigned int offset) { this->m_offset = offset; }
 	unsigned int getOffset() const { return this->m_offset; }