transforming coordinates from one distorted coordinate system to another

the problem is best explained with an example:

http://dl.dropbox.com/u/1013446/distortedcoordinatespace.exe

drag and drop the little red square inside the small square on the right. it corresponds to the red square in the big quadrilateral on the left. you can also drag 开发者_如何学编程the 4 corners of the big quadrilateral on the left to see how it occupies a distorted version of the space within the square.

given the absolute coordinates for the 4 points of a square and the coordinates of an arbitrary point within the square, it's simple matter to remap the point's coordinates to an arbitrary quadrilateral.

what I want is to be able to start off with an arbitrary quadrilateral, and be able to do the same thing, transforming the quadrilateral to any other 4 sided shape, but maintaining the relative distorted position of the point,

so given the 4 absolute coordinates of each of 2 irregular quadrilaterals, A and B, how can I convert the coordinates of point C given it's absolute coordinates?

also helpful, would be any terminology that I'm missing here for what these transformations would be called, because I'd like to look into them more

ok, I'm attempting to implement btilly's solution, and here's what I have so far:

   #include<complex>
#define cf complex<float>

cf i=sqrt(complex<float>(-1));

cf GetZ(float x,float y)
{
    return cf(x)+(cf(y)*i);
}

cf GetPathIntegral(cf p1,cf p2,cf q1,cf q2, int n)
{
    cf sum;
    for (int index=0;index<=n;index++)
    {
        cf s=cf(float(index)/float(n));
        cf weight;
        if (index==0||index==n)
            weight=1;
        else if(index%2)
            weight=4;
        else weight =2;


        sum+=(((cf(1)-s)*q1)+(s*q2))*(p2-p1)*weight;
    }
    return sum/cf((3.0*(n-1.0)));
}

before I move on from here, I want to make sure I'm right so far...

also, this paragraph confused me a bit:

OK, so we can do path integrals. What is the value of that? Well suppose we take a random point z0 = x + iy somewhere in our region. Suppose that f(z) is defined on the path. Then the Cauchy Integral Formula says that the integral around our region (which is the sum of 4 piecewise integrals that we know how to do) of f(z)/(2 * π * i * (z - z0)) is a really nice function that is going to match our original function on the boundary.

what does the function do exactly?

(My first pass was a complicated derivation of a natural seeming formula for this. But then I realized that there is a far, far better solution. Which I would have remembered earlier if I had used Complex Analysis in the last 20 years.)

The right way to do this is to apply the Cauchy Integral Formula. With this you can map any polygon to any other polygon. If the polygons don't self-intersect, it will send the boundary to the boundary and the interior to the interior. The mapping will also have the excellent property that it is conformal, meaning that angles are preserved. By that I mean that if a pair of curves intersect in your region, then they will be mapped to a pair of curves that intersect at the same angle. (Many of Escher's drawings are based on conformal mappings.)

Enough hype. How do you do it? I'll explain it, assuming that you know absolutely nothing about complex analysis. I'll use some Calculus terms, but you should be able to follow my directions even if you don't know any Calculus at all. Since I am assuming so little, the explanation has to be kind of long. I'm sorry for that.

Every point (x, y) in the real plane can be viewed as a complex number z = x + iy. We can add and multiply complex numbers using the usual rules of algebra and the fact that i * i = -1. Furthermore note that 1 = (x + iy) * (x - iy)/(x2 + y2) so we can divide if we let 1/z = (x - iy)/(x2 + y2). We therefore have all of the usual rules of arithmetic.

But we can do better than that. We can do Calculus. In particular we can do path integrals around curves. An integral of a function along a curve is a kind of weighted average of that function over the points in that curve. You can read up on how to do it in general. But here is how to do it in this case.

Suppose that the starting region has corners P1, P2, P3, P4. The path around the region is defined by the four line segments (P1, P2), (P2, P3), (P3, P4), (P4, P1). I'll talk about how to handle the first line segment. The others are similar.

The path integral of f(z) over (P1, P2) is the integral from 0 to 1 of f((1-s)P1 + sP2)(P2 - P1). To evaluate that integral, the easiest thing to do is numerical integration using Simpson's Rule. To do this pick an odd number n and for the values s = 0, 1/n, 2/n, ..., (n-1)/n, 1 assign them weights in the pattern 1, 4, 2, 4, 2, ..., 2, 4, 1. (The end points are 1, everything else alternates between 4 and 2.) Now for each point calculate f((1-s)P1 + sP2)(P2 - P1), multiply by the weight, and add them all together. Then divide by the magic value 3 * (n-1). The result is approximately your integral. (As n grows, the error in this approximation is O(1/n4). In your case if you take n = 21 then the approximation should wind up good enough to map pixels to the right pixel, except for some pixels near the boundary. Make it a little larger, and the problematic area will get smaller. Right at the edge you'll want some multiple of the number of pixels on a side to make the error small .)

OK, so we can do path integrals. What is the value of that? Well suppose we take a random point z0 = x + iy somewhere in our region. Suppose that f(z) is defined on the path. Then the Cauchy Integral Formula says that the integral around our region (which is the sum of 4 piecewise integrals that we know how to do) of f(z)/(2 * π * i * (z - z0)) is a really nice function that is going to match our original function on the boundary. I won't get into all of the "really nice" things about it, but what I was saying above about conformal is part of it.

Now what function f do we use? Well suppose that our region is being mapped to a region with corners Q1, Q2, Q3, Q4. We want the first path piece to map to the second path piece. So we want f((1-s)P1 + sP2) to be (1-s)Q1 + sQ2. That tells us how to calculate f at all of the points we need to do our integral.

Now, you ask, how do you reverse it? That's simple. Just reverse the role of the two polygons, and calculate the reverse transformation! Which brings a really good unit test. You should define a couple of weird regions, pick a point in the middle, and verify that if you map from the first to the second and back again that you wind up close to where you started. If you pass that test, then you probably have made no mistakes.

And finally what about my general polygon claim that I made? Well we defined our path as four pieces we traversed linearly. A higher degree polygon just has more pieces to its path, but otherwise the calculation is done in exactly the same way.

found the solution. I have to say, it's much more complicated than I had expected:

assuming a square or quadrilateral has the four corners:

AB
CD

you need an interpolation factor: xt for the x axis, and yt for the y axis, so that

if you define a linear interpolation formula:

lerp(j,k,t)
{
     return (t*(k-j))+j;
}

a point p within the ABCD quad is defined as:

p.x=lerp(lerp(a.x,b.x,xt),lerp(c.x,d.x,xt),yt)

and

p.y=lerp(lerp(a.y,c.y,yt),lerp(b.y,d.y,yt),xt)

then the values you need to define are xt and yt

xt= ((2* c.x* a.y) - (d.x* a.y) - (2 *a.x c.y) + (b.x c.y) - (c.x* b.y) + (a.x* d.y) - (a.y* p.x) + (c.y* p.x )+ (b.y p.x) - (d.y p.x) + (a.x p.y) - (b.x p.y) - (c.x* p.y) + (d.x* p.y) - Sqrt(-4* ((c.x* a.y) - (d.x* a.y) - (a.x* c.y) + (b.x* c.y) - (c.x* b.y) + (d.x* b.y) + (a.x d.y) - (b.x d.y))* ((c.x* a.y) - (a.x* c.y) - (a.y* p.x) + (c.y* p.x) + (a.x* p.y) - (c.x* p.y)) + ((-2 *c.x a.y) + (d.x a.y) + (2 *a.x c.y) - (b.x c.y) + (c.x* b.y) - (a.x* d.y) + (a.y* p.x) - (c.y* p.x) - (b.y* p.x) + (d.y* p.x) - (a.x* p.y) + (b.x* p.y) + (c.x* p.y) - (d.x p.y))^2))/(2 ((c.x* a.y) - (d.x* a.y) - (a.x* c.y) + (b.x* c.y) - (c.x* b.y) + (d.x *b.y) + (a.x *d.y) - ( b.x *d.y)))

and once you have that

yt=(p.x-lerp(a.x,b.x,('xt')))/(lerp(c.x,d.x,('xt'))-lerp(a.x,b.x,('xt')))