The weighting function I used was a cubic approx. to a bell-shape. it turns out that the function f(x) = x*x*(3-2*x) maps the interval [0,1] onto itself and has a derivative of zero at each endpoint. The weight I use for points is f(distance from center of image).
My projection scheme didn't use the standard 8x8 matrix, instead i mapped both quadrilaterals onto a unit square using 4x4 matrices, and then took the inverse of one of them. A translate, followed by a rotate, followed by a skew puts the quadrilateral in a position where its easy to solve for how to project it onto the unit square. Multiplying these projections toghether yields the final 4x4 transform onto the unit square. This is probably not any faster then the 8x8 matrix method, but its a different approach, and I found it intuitive.