Interpolating a rectangle [Resolved]

I have a rectangle defined by it's bottom left coordinates, and a width and height.

I have a destination rectangle, and i want to interpolate between the 2.

I've done this in a linear fasion by doing linear interpolation on the 4 values (x,y,w,h).

Unfortunately this is no good for what i want as the rectangle represents a 2D viewport.

If the destination rectangle is much smaller than the initial size (which is often is) then to start with the width/height changes by a small percentage, but by the end it's changing by a much larger percentage.
As a result the animation looks like it starts of slowly zooming in, gradually speeding up.

What i want is a smooth zooming in. I guess that would meen for each time step i want the rectangle to change by a constant percentage.

I've got no idea of the maths involved though, so any help would be much appreciated.

If you want the have the constant percentage of change just use the desired percentage on each step and do this until until the desired value is reached. That way you haver to specify the percentage of change.
For example you want a 10&#37; change, calculate the next step by NewValue=LatValue*.9. Continue until the NewValue >= DesiredValue (of course the check must be <= if the DesiredValue is smaller then the originalValue).

Thanks for the reply,
Ideally i want a function that can give the value as a function of start value, end value, and time (between 0 and 1), rather than calculating it as an iterative process.

With your method i'm not sure how to make it reach the end value after a given number of steps.

Sorry, your setup is somewhat unclear to me. You have a viewport like this

Code:

---

_____________________________
 |                          |
 |        _______          |
 |        |    |          |
y|        h|    |          |
 |        |_____|          |
 |          w              |
 |___________________________|
              x

---

where the larger rectangle given by x and y is the initial rectangle and the smaller rectangle given by w and h is the final rectangle. You're transforming the large rectangle into the small rectangle using some sort of animation. Is the display showing a 2D or 3D world? ("Viewport" is context sensitive which is why I ask.) What exactly is doing the "zooming"? Is it the world inside the viewport, or is it the rectangle?


Opus' method may well be what you want, but I can't tell from the information given. However,

Quote:

---

Ideally i want a function that can give the value as a function of start value, end value, and time (between 0 and 1)

---

can be done. Take a start value I, a final value F, and a ratio R where each step multiplies the current value by a factor of R. For instance, if I=1000, F=512, and R=0.8, the steps give widths 1000, 1000*0.8=800, 800*0.8=640, 640*0.8=512.

The ith step is computed as
C(i) = C(i-1)*R = C(i-2)*R*R = ... = C(i-k)*R^k
C(0) = I
=> C(i) = C(i-i)*R^i = C(0)*R^i = I*R^i

Suppose you want N total steps. That is, C(N) = F. Then we can find the ratio required:

C(N) = F = I*R^N
=> F/I = R^N
=> (F/I)^(1/N) = R

In the example above with I=1000, F=512, suppose we wanted 3 steps (4 including the initial value, but 3 multiplications). The ratio we want is then

(512/1000)^(1/3) = 0.8

which gave 3 computation steps above. You can easily make this description into a function of time instead of step number. Assigning time 1 to step N and time 0 to step 0, you can convert from time t to step number n by multiplying t by N. Your function of time is then

Code:

---

f(t) = C(n) = C(t*N) = I*R^(t*N)

---

Since R = (F/I)^(1/N), this simplifies to

Code:

---

f(t) = I*R^(t*N) = I*((F/I)^(1/N))^(t*N) = I*(F/I)^(tN/N) = I*(F/I)^t

---

You can see that f(0) = I*(F/I)^0 = I*1 = I, f(1) = I*(F/I)^1 = I*F/I = F, so f indeed interpolates between I and F as t goes from 0 to 1.

Thanks for steping in, jemidiah, I havn't had a coffee by the time I was replying.

@SLH

Quote:

---

Ideally i want a function that can give the value as a function of start value, end value, and time (between 0 and 1)

---

Your time-value will be anything between 0 and 1? In that case you need to define useable timesteps in order to use the function mentioned by jemidiah.

Thanks for that huge reply jemidiah. My intial explanation wasn't quite clear, although i think i can use your explanation/formula to achieve what i'm after.

I'm rendering a fractal animation, zooming in on the fractal. So my 'viewport' is really a rectangle on the complex plane, with x being the bottom-left corner position along the real axis and y the bottom-left corner position along the imaginary. Similarly w is the width of this rectangle along the real axis, and y the height of the rectangle, along the imaginary axis.

Given this rectangle and my desired image size, i can workout what complex number each pixel in the resulting image relates to. I then calculate the pixel colour using a function of this complex number.


I think i can apply the I*(F/I)^t formula to my setup, thanks a lot for the working out; very good breakdown, i just about understand it!
I'll let you know how i get on and see if i can post a little gif or something of the results (if successfull!).

Ok, i've sort of got it working.

jemidiah's formula is perfect for the width and height.

I've changed it so that the rectangle is defined by a point in the center of the rectangle its the width and height.

The width and height i've got interpolating as per jemidiah's formula.

Interpolating the x/y values in a linear way, which is causing the problem.

The change in x/y values get larger in respect to the width and height changes as t gets closer to 1.

As a result the image starts panning a lot, to the extent that sometimes the change in x or y is more than the width/height!

My feeling is that i've somehow got to do the same sort of interpolation for X and Y, however it's not quite as straight forward as the start / finish x and y values can be +ve or -ve, meaning the formula I * (F/I)^t falls over if I and F don't have the same sign.

My thought is that i have get them both into the +ve by adding the lower number to both, then calculate the new value using the formula, then subtract what was added originally.

Any ideas?

The interpolation I derived in my previous post doesn't seem to make sense while panning. For instance, the interpolation has a very different shape if you just move the origin. However, moving the origin shouldn't change pan speed at all.

I think I understand your setup now. You have one target image size, but you have two rectangular areas on the complex plane. You want to move smoothly (visually) from the first rectangle to the second. The width and height change between the two, and the center of the rectangles change as well, resulting in both a pan and a zoom. Ignoring the pan, you're happy using the exponential interpolation to get the zoom working. You want to get the pan working as well.

I think instead of absolute position you should be thinking in terms of velocity. If you do a simple linear interpolation, your initial velocity will be magnified as you zoom in. Say you pan straight down. When your viewing height is half of its original height, the velocity of the pan will appear twice as fast. Without going into detail, you can generalize this to say if the initial width and height are wI and hI, and the current width and height are w(t) and h(t), the velocity of a linear interpolation will be magnified by a factor of M(t) = Sqrt((wI/w(t))^2 + (hI/h(t))^2). You can counteract this by scaling the linear interpolation's velocity down by that same factor.

A linear interpolation can be thought of as moving from one point to another at a constant velocity, V. To keep the apparent speed constant, modify this interpolation by making it time dependent, moving instead at velocity V/M(t) with M(t) as above. At this point, you theoretically have enough information to code a smooth pan. However, you'd have to precompute the constant V, which is fixed by the condition that you start at your initial position at frame 0 and end at your final position at frame N. However, you'd prefer a non-iterative interpolation function of time, which might well work. Let's see if a nice form exists.


Try interpolating the initial rectangle center x coordinate xI to xF. We want to find the interpolated x coordinate at time t, x(t). At the same time, the initial rectangle width wI is interpolating to wF according to w(t) = wI*(wF/wI)^t, and the initial rectangle height hI is interpolating to hF according to h(t) = hI*(hF/hI)^t. Under these conditions the scale factor M(t) can be found:

wI/w(t) = wI/(wI*(wF/wI)^t) = (wI/wF)^t = (Sqrt(W))^t
hI/h(t) = (hI/hF)^t = (Sqrt(H))^t
M(t) = Sqrt(W^t + H^t)

with convenient substitutions H = (hI/hF)^2 and W = (wI/wF)^2.

The velocity is then V/M(t) = V/Sqrt(W^t + H^t). Velocity is the rate of change of position, which is also called the derivative. Going from velocity to position, position is the antiderivative of velocity, also called the integral. Again without going into details, we then have

v(t) = V/M(t)
=> x(t) = xI + Integral(v(T), T, 0, t)
=> x(t) = xI + Integral(V/M(t), T, 0, t)
=> x(t) = xI + V*Integral(1/M(T), T, 0, t)

where T is a dummy variable. The velocity scale constant, V, is given by the condition that

x(1) = xF = xI + V*Integral(1/M(T), T, 0, 1)
=> xF-xI = V*Integral(1/M(T), T, 0, 1)
=> V = (xF-xI)/Integral(1/M(T), T, 0, 1)

Evaluating this integral analytically is unfortunately not possible (at least, Mathematica can't do it, and it's hideously good at integrals). At this point you have a few options.

(1) You could just use the iterative form of this solution--i.e. divide your time up into frames, compute the velocity scale factor at each frame in turn, and move the point between frames according to that velocity. You'd have to compute the velocity scale factor V by first running this iteration with a guessed value for V until you reach your target final position xF. You'd see how long it took and scale V accordingly so that, for instance, instead of 20 frames you'd double the speed at each point so the interpolation only takes 10 frames.

(2) You could use a numeric integration routine to evaluate V and then x(t) using the expressions above. There have to be a zillion libraries that do this. These functions aren't very bad, either, and your own home-brew integration routine might well work fine. I can give more details if you want to pursue this, though an internet search for "numeric integration" would also be very fruitful.

(3) You could use an approximate form of the function 1/M(t), which could be relatively easily integrated. Using a 2nd order Maclaurin expansion gives the approximation (using Mathematica so the algebra is correct)

1/M(t) ~= 1/Sqrt[2] - ( ((Log[H] + Log[W])) t)/(4 Sqrt[2]) - (((Log[H]^2 - 6 Log[H] Log[W] + Log[W]^2)) t^2)/(32 Sqrt[2])

I won't give you the details of where this comes from unless you're interested. Anyway, using this quadratic approximation, the definite integral of M is simply

I(t) ~= Integral(1/M(T), T, 0, t) ~= t/Sqrt[2] - ((Log[H] + Log[W]) t^2)/(8 Sqrt[2]) - ((Log[H]^2 - 6 Log[H] Log[W] + Log[W]^2) t^3)/(96 Sqrt[2])

V and x(t) are then simply

V=(xF-xI)/I(1)
x(t) = xI + V*I(t)


There's no guarantee this approximation is any good for this particular function. It looked like the terms were decreasing pretty quickly, or I wouldn't have suggested it at all.


The upshot is, using the exponential width and height interpolation, you can indeed get a smooth pan, but it's not particularly easy. This post contains the exact method you were asking about, but it's a bit complicated, especially if you're not used to calculus. Sorry :(. I've bolded important equations or substitutions to try and make it a little easier to digest.

Wow, thanks for the very detailed explanation!

I've read it through a few times, and am still not 100&#37;, but it's starting to make sense.

I've got to refactor my code a bit so that instead of interpolating the x,y,w,h values seperately it interpolates the rectangle as a whole (so the w/h values are available when interpolating x/y) so i'll get that done first, then have a go at implimenting the calculations.

I think numerical integration (option 2) is probably the way you'll go. In a brief search, ALGLIB offers code to do so in its autogk module, which is available as VB6 code. If you download the whole thing, testautogk gives some examples of how to use it, and I'm sure there's more elsewhere. It's not very idiomatic code for VB6 since it was ported from C.

There are also many other sources of such a routine. You could write your own easily enough:

Code:

---

Function Integrate(f as function, a as number, b as number)
    size = (b-a)/10000
    For i = a to b Step size
        Integrate = Integrate + f(i) * size

---

is pseudocode for a very rudimentary, but probably accurate enough, numerical integration routine. The above pseudocode is some strange mix of VB6, C, and Python. Sorry :P.

Right, i think i've implimented option 2, however i'm still getting the same issue so obviously i've done something slightly wrong.

It looks like up until about half way through the animation the panning looks ok, but then it very rapidly seems to speed up.

Here is my c# version of the integration routine:

Code:

---

private double Integrate(Func<double, double> f, double Start, double End)
{
    double Size = (End - Start) / 10000.0;

    double Result = 0.0;

    for (double i = Start; i < End; i += Size)
    {
        Result += f(i) * Size;
    }

    return Result;
}

---

With an small wrapper function that calls it passing in the M(x) function.
iR is the initial rectangle and fR is the final rectangle.

Code:

---

private double I(RectangleD iR, RectangleD fR, double t)
{
    double W = Math.Pow(iR.Width / fR.Width, 2);
    double H = Math.Pow(iR.Height / fR.Height, 2);

    return Integrate(
        (T) => 1 / Math.Sqrt(Math.Pow(W, T) + Math.Pow(H, T)),
        0,
        t);
}

---

Finally, here's my code that calculates the x,y coordinates as well as the width and height (variables tX, tY, tW, tH respectively):

Code:

---

double double tX, tY, tW, tH;

//Calculate Width
double i = iR.Width;
double f = fR.Width;
tW = i * Math.Pow(f / i, t);

//Calculate Height
i = iR.Height;
f = fR.Height;
tH = i * Math.Pow(f / i, t);

//Calculate X
double V = (fR.X - iR.X) / I(iR, fR, 1);
tX = iR.X + V * I(iR, fR, t);

//Calculate Y
V = (fR.Y - iR.Y) / I(iR, fR, 1);
tY = iR.Y + V * I(iR, fR, t);

---

I'm hoping that i've just misunderstood the mathematical notation in your post and so got something wrong with the translation into code. I can follow the workings in your post (although the ideas etc. are certainly beyond what i could come up with!) so i'm a bit unsure where i've gone wrong.

I haven't been able to find mistakes in either your code or my previous reasoning. Perhaps posting a sample of the current output would help.

I implemented my method in a throw-away VB6 project. The attached video zooms and pans along the Mandelbrot Set. The width and height interpolation is the exponential one derived in my first post. The center position interpolation is the integral one derived in my second post. It seems to be nice and smooth to me.

The viewport's initial center is (-1, 0) with width 1, height 1. The final center is (-2, 0) with width 0.003, height 0.003. The video has been heavily resized and so is very pixelated, but there's enough detail to see what's going on. It's encoded with H.264.

I've just ran my animation routine using your parameters and i get an identical result.

I'm trying it for a julia fractal with formula and starting coordinates as below:

Code:

---

<Fractal name="Julia" maxiterations="512" resolution="1" colourscheme="NewScheme">
  <Formulas>
    <Init formula="Z = Pixel;C = {-0.8,0.156}" />
    <Iteration formula="Z = Z ^ 2 + C" />
    <Bailout formula="Z.MagSq() >= 4" />
    <Weighting formula="" />
  </Formulas>
  <Dimensions>
    <Origin>
      <ComplexNumber r="0" i="0" />
    </Origin>
    <Size>
      <ComplexNumber r="3.0" i="3.0" />
    </Size>
  </Dimensions>
</Fractal>

---

The init function runs one per pixel, the iteration formula runs once per iteration per pixel and the bailout formula is tested once per interation.

I'm trying to interpolate over 500 frames to:

Code:

---

<Rectangle x="-1.0254925351942"
y="0.477417634828581"
w="7.73634685202751E-08"
h="7.73634685202751E-08" />

---

I'm wondering if it's as simple as precision problems, although i can zoom in further than that manually with no problem.

By the way, is that a fractal in your avatar? If so, what's the formula?

There's certainly a possibility of precision problems with the numeric integration routine I provided. It adds up 10,000 small numbers to get its result, which would chop off several significant figures due to rounding errors. The real culprit is probably that for t near zero, I(t) is, say, 1.0, while for t near 1, I(t) starts to change only very small amounts, say by 0.00000001 per frame. This small amount is obtained by adding even smaller numbers to the initial I(t) in the integral summation, which would be another 5 orders of magnitude below the value added, say at 0.0000000000001. At this point you're approaching the limit of double-precision arithmetic (though I don't think you should have exceeded it, at least not before orders of magnitude of rounding errors accumulate). There might be an effect I'm missing which pushes things over the edge, so-to-speak. Perhaps the exponential or power functions in .NET don't provide ~16 decimal digits of precision, so M is slightly incorrect; I dunno.

A few things come to mind for you to try.

(1) You could see if changing the magic number 10,000 in my routine to something smaller, say 10, 100, or 1000, has any noticeable effect. Without testing it, I'd guess that 10 would be pretty horrifically bad.

(2) You could use a different numeric integration routine.

(3) You could do all of the calculations in the velocity domain instead of the position domain. That is, on each frame, you could calculate v(t) and just add v(t)*dt to the current position vector, where dt is the time between frames. This approximation basically says that the pan velocity stays constant for the time it takes to move between two frames before recalculating for the next frame.


Since #3 is the most promising, it's the one I tested. At each step the x and y coordinates need to be updated in a new way by projecting the velocity onto the vector pointing from the initial to final center coordinates. I won't get into details, but (using the notation from my earlier post), you need the following convenient constants:

Code:

---

dx = (xF - xI)
dy = (yF - yI)
D = Sqrt(dx^2 + dy^2)
dt = 1 / (NumFrames + 1) // length of time between frames
px = dx / D              // Projection of unit velocity vector in x direction
py = dy / D              // Projection in y direction

---

You also calculate V using

Code:

---

V = D / I(1)

---

using the same integral as before. However, it turns out V will be calculated *too* exactly using a large number of iterations, so we actually change the magic number in my integral routine from 10,000 to NumFrames - 1 (I haven't thought through if the -1 is strictly necessary, but it made my test work out and was natural using another notation I'm translating from).

At each frame you need to increment the positions using the following:

Code:

---

dpos = V / M(t) * dt
x = x + dpos * px
y = y + dpos * py

---

Using this method and your previous post's constants, I generated the attached video (again, heavily pixelated, H.264 encoded). I only rendered every 5th frame (i.e. 100 frames are encoded in the video) both for size and speed, but the parameters were updated every frame (i.e. 501 times; I happened to render frame numbers 0 to 500). The video looks quite smooth to me. I'm sure it would be beautiful at a high resolution, maybe with the missing frames added in.


I'm curious, how long does it take for your program to compute the animation in your previous post, and at what resolution?


It certainly is a fractal in my avatar. I don't remember the formula, or even if I generated it or found it online. Fractal generation was a hobby of mine years ago; maybe I should update it with another one in my collection :).

Thanks for another very detailed post!

I've change my code around a bit to use a the method described in your post (incrementing at each step).

Unfortunately i'm not quite getting the results as expected.
I'm trying to render the mandelbrot animation again, with X going from -1 to -2.

When i run it X actually only goes to -1.98646453189044, although the animation does look very smooth!

I altered my integrate function, replacing 10000 with EndFrame - StartFrame - 1

Below is my code to calculate tX and tY:

Code:

---

double dx = fR.X - iR.X;
double dy = fR.Y - iR.Y;

double D = Math.Sqrt(dx * dx + dy * dy);

double dt = 1.0 / (EndFrame - StartFrame + 1);

double px = dx / D;
double py = dy / D;

double V = D / I(iR, fR, 1.0);

double dpos = V / M(iR, fR, t) * dt;

tX = (PreviousValue as RectangleD).X + dpos * px;
tY = (PreviousValue as RectangleD).Y + dpos * py;

---

I'm hoping i've done something wrong, rather than it being a precision thing again!

As for how long it takes to render, i've just done the Julia animation with the same parameters as earlier, with an image size of 117 x 117 and it took 10mins, altough i've not done any optomisation.

Your NumFrames calculation looks to be slightly off. If your start frame is 1 and your end frame is 10, there are 10-1 + 1 = 10 frames, not 10-1 frames. That is, NumFrames = EndFrame - StartFrame + 1. You'd have to change the magic constant in the integrate routine and your calculation of dt.

When I rendered the movie attached to my previous post, I specifically checked to make sure the final position was correct, and it was. Since it's a (much) more complicated case that's (much) more sensitive to approximation errors, and it worked, I'm pretty confident in the general method.

Thanks for all your patience!

Hmm, i've tried altering it as you say but there must be something else i'm doing wrong.
Can i take a look at your code to see if there's anything i'm doing differently?

My code is attached. Hopefully it helps debug the problem.

Woohoo! Got it sorted.
The issue was that i was incrementing the current frame number before applying the formula rather than afterwards.

One thing that i don't understand though, is why in the numerical integrate routine the magic number had to be changed to the number of frames.

As promised attached is a little animation (compressed to hell though!).

Even with the incredibly low bitrate parts of that look pretty nice. And it looks like the zoom and pan are very smooth!

To be honest, I actually hacked the numerical integration routine to do the required calculation when computing V while staying in the velocity domain. It just turned out that my approximation was exactly the routine required with the magic number tweak.

I'll run through a cut-down example. Say the magnification level M(t) happens to be 1, 2, and then 4 for a three frame animation. The rate of motion at each frame is V/1, V/2, and V/4, for some constant V, so that the motion appears smooth visually. Say we need to move a total distance of D=5 units during the pan. Since distance = rate * time, and our rate is V/1, V/2, or V/4, for time steps dt between frames, we can calculate the distance panned between each frame as V/n*dt, n=1, 2, or 4. The total distance moved is then V/1*dt + V/2*dt + V/4*dt = D = 5. Thus,

V = D/(1/1*dt + 1/2*dt + 1/4*dt).

In the integral call, Start and End will always be 0 and 1, so Size is basically 1/(NumFrames+1) = dt. [It looks like maybe you'll get slightly more accurate results if you use the magic constant NumFrames+1 instead of NumFrames-1, but meh, I knew it was NumFrames possibly with an off-by-1 error, so I just tested it to see what gave basically correct results, and NumFrames-1 did :).]

The integral call computes the sum of 1/M(t)*Size, with M(t) computed at each frame, with each term multiplied effectively by dt = Size. This happens to be exactly the denominator necessary to compute V in this case. It's really just a convenient coincidence.

Yeah, i've been playing around with colouring. Currently you can place any number of colours on a scale 0..1 and it interpolates between them. If you click the fractal it calculates the iterations/maxitterations, pops up a colourchooser and adds that colour in. It's great for interactively altering a fractal's colour.

Next thing on the list is adding animateable parameters to the fractal's formulas, as such:

Code:

---

<Fractal name="MandelbrotAlt" maxiterations="512" resolution="1" colourscheme="NewScheme">
  <FormulaParameters>
    <FormulaParameter name="SinMultiplier" defaultvalue="{0,0}" />
  </FormulaParameters>
  <Formulas>
    <Init formula="Z = 0;C = Pixel" />
    <Iteration formula="Z = Z ^ 2 + C + Z.Sin() * SinMultiplier" />
    <Bailout formula="Z.MagSq() >= 4" />
    <Weighting formula="" />
  </Formulas>
  <Dimensions>
    <Origin>
      <ComplexNumber r="-0.5" i="0" />
    </Origin>
    <Size>
      <ComplexNumber r="3.0" i="3.0" />
    </Size>
  </Dimensions>
</Fractal>

---

Code:

---

<Animation name="MandelbrotAltAnim" fractal="MandelbrotAlt" frames="500">
  <KeyFrames>
    <KeyFrame objectname="SinMultiplier" property="Value" frame="500" value="{1,0}" >
    </KeyFrame>
  </KeyFrames>
</Animation>

---

If i could figure out some good formula's to apply paramters to i think i should be able to get some pretty neat animations.

I noticed my translation of my code several posts back was slightly off. It should have been dt = 1 / (NumFrames - 1) instead of 1 / (NumFrames + 1). The integration routine looks to have been translated correctly in the first place. Then again, if you're getting pan/zoom results you're happy with maybe you've tweaked it to overcome this error.

One thing you could try is to morph one fractal you like into another between keyframes. That is, you could, say, linearly interpolate between using z |-> f1(z, c) and z |-> f2(z, c) by using a combination of the two, say z |-> t*f1(z, c) + (1-t)*f2(z, c) where t ranges from 0 to 1. You might need to scale and translate one of the fractals [which could be done by letting Pixel be replaced by (Pixel - NewOrigin) * (WidthScale, HeightScale) in your Init formulas; if Pixel is only used to initialize c, you'd use... z |-> t*f1(z, c) + (1-t)*f2(z, (c - NewOrigin) * (WidthScale, HeightScale))] so that the interesting parts match up during a smooth zoom/pan.

There's also just perturbing an existing fractal like you've done by adding a time-varying sine term. I'm sure some of those are quite interesting, but I think guess and check is about all that's possible to see what works and what doesn't. There are many other "tweaks", of course: multiplicative terms, additive terms, exponent terms, base terms (i.e. SinMultiplier^z); each term could contain standard arithmetic (+, -, *, /, \, ^, fractional part), standard special functions (logarithms, trig functions--standard, hyperbolic, inverse, hyperbolic inverse), z and/or c, ....

Of course, calculations in all of the above needs to be complex. For instance, defining cos for complex inputs can be done with cos z = [e^(i z) + e^(-i z)] / 2 where e^(i z) = cos(Re[z]) + i sin(Im[z]).

Going from one fractal formula to another over time sounds like it could produce some good animations; i'll definately have to give that a go!

Cheers for all your help, what started as a pretty simple question has really given me a good example of following logical/mathematical reasoning through (and i've got some pretty good fractal animations too!).

Thanks for all your help!

This a very interesting thread to me as I've been scratching my head over something similar for some time, and all I've gotten so far is dandruff under my fingernails!!

I have two maps of the same geographical area done by two different government agencies. Both maps show the same typical section lines at approximately one mile intervals. A section is of course defined as one square mile. Anomolies exist between the two maps based primarily on the differing surveying and projection methods used. Consequently, when measuring a given one square mile section on one map, you'll find that the section as projected on the map will probably not have real 90 degree corners, will not have four equal length sides, and the sides may not really be parallel. Somewhere inside of this section I have a single point that I'm interested in. I can digitally "project" this section onto a true square cartesian system and find the absolute values of the four corners and the single point. For example, assume that the lower left corner is at x = 0, y = 0; the upper left at 2, 100; the upper right at 99, 101; the lower right at 102, -2, and the single point at 25, 25. I can also project the same section from the second map onto the cartesian system and locate it's four corners, again assuming that the lower left corner is at 0, 0. It is highly unlikely that any of the other corners will match between the two maps.

What I want to do is take the x and y coordinates of the single point on map 1, and through interpolation, find it's corresponding location on map 2. This would be a simple task if both sections at least had all 90 degree corners, but with so many variables involved, it gets real complicated real fast. I've come close a couple of times, but the method used is always so esoteric that it's been nearly impossible to reproduce programmatically.

Just so everyone is talking and responding in the same language, assume the following variable names:

MAP 1
Lower left corner: Xp11, Yp11
Upper left corner: Xp12, Yp12
Upper right corner: Xp13, Yp13
Lower right corner: Xp14, Yp14
Single point: X1, Y1

Lower left corner: Xp21, Yp21
Upper left corner: Xp22, Yp22
Upper right corner: Xp23, Yp23
Lower right corner: Xp24, Yp24
Single point: X2, Y2

Totally confused now?? See what you can do with that!!

Regards,

Chief Brad

As I understand it, you have two 4-sided figures. You want to map each corner of each figure to its corresponding corner in the other figure. Further, for points inside the two figures, you want some sort of reasonable interpolation, where reasonable is in the eye of the beholder.

If your two shapes were parallelograms, this would be a simple linear algebra problem. However, allowing arbitrary angles and side lengths makes it more complicated.

My solution first maps one shape to the unit square. It then maps those points on the unit square to points on the second shape. We hope to make a(n invertible) function taking in points on the shape and spitting out points on the unit square. Say a point on the shape is given by (x,y), and a point on the unit square is given by (u,v). We want to find a way to convert between these two pairs. Say g does this translation: g(x,y) = (u,v). Say it's inverse is f, so f(u,v) = (x,y).

Take the points of the shape to be, going clockwise from lower left, p1 = (Xp1, Yp1), p2 = (Xp2, Yp2), p3 = (Xp3, Yp3), p4 = (Xp4, Yp4). For instance, g(Xp1,Yp1) = (0,0), g(Xp2,Yp2) = (0,1), g(Xp3,Yp3) = (1,1), g(Xp4,Yp4) = (1,0); f(0,0) = (Xp1,Yp1), f(0,1) = (Xp2,Yp2), f(1,1) = (Xp3,Yp3), f(1,0) = (Xp4,Yp4). A reasonable way to find the interior points, to me, is an interpolation linear in u or v when v or u is held constant. Unfortunately a fully linear interpolation isn't possible in general because the shape isn't always a parallelogram.

Take the left and right sides of the shape. As you sweep across the shape starting at the left side going to the right side, imagine dragging a line which transforms from being parallel to the left side at the start to being parallel to the right side at the end. This is something like what you might do if you wanted to erase a box on a chalkboard with one stroke.

I'll give this transformation algebraically:

Define U1 = p2-p1 (left edge), U2 = p3-p2 (top edge), U3 = p3-p4 (right edge), U4 = p4-p1 (bottom edge). Transform U1 into U3 as u moves from 0 to 1 using a simple linear interpolation, U1*(1-u)+U3*u. Halfway across, u=1/2, and our chalkboard eraser is midway between the left and right edge's orientations. The bottom of the eraser is also halfway across the bottom edge at 1/2*U4 = u*U4. If you wanted to get the point halfway up along the eraser (say, at v=1/2) you'd go from the bottom edge u*U4 and add half the chalkboard eraser, 1/2*(U1*(1-u)+U3*u). In general, the 1/2 would just be v. In general, then, continuing to use vector math,

(x,y) = f(u,v) = u*U4 + v*(U1*(1-u)+U3*u) = u*U4 + uv*(U3 - U1) + v*U1

Note that if the left and right edges are parallel and of equal length, U3 = U1, so U3-U1 = 0, and the transformation is linear as expected. Also, taking v constant, this is a linear function of u; taking u constant, this is a linear function of v.

This can be written as a system of equations of the form

x = u*c1 + uv*c2 + v*c3
y = u*d1 + uv*d2 + v*d3

We can translate from (u,v) to (x,y) using these equations. However, to accomplish your interpolation, we need to translate from (x,y) on the first shape into (u,v) on the unit square, and from (u,v) into (x,y) on the second shape. That is, we need to be able to go from (x,y) to (u,v), so we'll need to invert these equations. This could be done with a page or two of algebra, most likely, but I just plugged it into Mathematica's Solve function, which gave the following two solutions for g(x,y) = (u,v):

u = -(1/(2 c2 d1 - 2 c1 d2))(c3 d1 - c1 d3 + d2 x - c2 y + Sqrt[4 (c3 d2 - c2 d3) (d1 x - c1 y) + (c3 d1 - c1 d3 - d2 x + c2 y)^2])
v = -(1/(2 c3 d2 - 2 c2 d3))(c3 d1 - c1 d3 - d2 x + c2 y + Sqrt[4 (c3 d2 - c2 d3) (d1 x - c1 y) + (c3 d1 - c1 d3 - d2 x + c2 y)^2])

or

u = (1/(2 c2 d1 - 2 c1 d2))(-c3 d1 + c1 d3 - d2 x + c2 y + Sqrt[4 (c3 d2 - c2 d3) (d1 x - c1 y) + (c3 d1 - c1 d3 - d2 x + c2 y)^2])
v = (1/(2 c3 d2 - 2 c2 d3))(-c3 d1 + c1 d3 + d2 x - c2 y + Sqrt[4 (c3 d2 - c2 d3) (d1 x - c1 y) + (c3 d1 - c1 d3 - d2 x + c2 y)^2])

(Only one of these solutions will fit the conditions that 0<=u<=1 and 0<=v<=1. This is basically solving a quadratic equation and tossing out a solution which doesn't fit certain conditions. The "correct" solution is probably always the top set or always the bottom set, but I'd have to do a large chunk of algebra to figure it out, and a computer would be happy just executing an "if" statement instead. Note that the two solutions differ only in the sign in front of the Sqrt term.)

So, we can now go from (x,y) to (u,v) with this interpolation as well; how nifty. At this point, you have everything necessary. More explicitly to convert from a point (X1, Y1) on shape 1 to a point (X2, Y2) on shape 2...

From shape 1 to unit square:
1. Calculate U1, U2, U3, and U4 using shape 1's coordinates (Xp1i, Yp1i), i=1...4.
2. Calculate c1, d1, c2, d2, c3, and d3 using the components of U4, (U3-U1), and U1, respectively.
3. Calculate u and v using x=X1, y=Y1 and the horrific formulas above. That is, calculate (u,v) = g(X1,Y1). Pick the solution which lands (u,v) on the unit square [0,1]x[0,1].

From unit square to shape 2:
1. Calculate U1, U2, U3, and U4 using shape 2's coordinates (Xp2i, Yp2i), i=1...4.
2. Calculate c1, d1, c2, d2, c3, and d3 using the components of U4, (U3-U1), and U1, respectively.
3. Calculate x=X2 and y=Y2 using the less horrific formulas above and u and v. That is, calculate (X2,Y2) = f(u,v).

Even more succinctly, (X2,Y2) = f2(g1(X1,Y1)), where g1 is g on shape 1 and f2 is f on shape 2. The reverse, of course, is simply (X1,Y1) = f1(g2(X2,Y2)).


If I were you, I'd double check my algebra (except the nonlinear equation solving step which Mathematica performed) to make sure no terms have been dropped; hunting around for an algebra mistake while implementing this would probably be awful. There's also the distinct possibility that Mathematica ignored some edge cases (i.e. if all the variables are 0, the denominator (c2d1 - c1d2) is zero). These probably can't physically occur, but without significantly more fiddling with the algebra I can't tell. It's easier to deal with these cases as they come up. I would also plug in (Xpi, Ypi) to g and make sure you get (0, 0), (1, 0), etc., as well as the reverse with (0, 0), (1, 0), etc. and f.

The above appears to be entirely correct. I had a vision of distorted zebras I wanted to realize and so briefly implemented it in the attached VB6 project. The attached image shows the project in action. The shape coordinates are hard-coded, though they were chosen so the shapes wouldn't be parallelograms (the most difficult case).

It was generated by first distorting a source picture of zebras by converting from (u,v) to (x,y). The distorted picture was then used as a source shape, undistorted (which is quite pretty to watch), and distorted again in a different way to produce the lower right image.