Welcome to
Joachim Breitner's Homepage

As I see it there is only one sensible projection that does what you want. For each pixel (x,y) you determine the distance d to the center (0,0). Then in the transformed image this pixel will be at (x* d / sqrt(x^2+y^2), y * d / sqrt(x^2+y^2)). It is easier in polar coordinates, you transform (r,theta) to (d,theta).

If your function has an inverse this also suggests a possible implementation. Simply walk over the pixels (d,theta) in the output image, and select the color from (f^-1(d),theta).

Thanks for your suggestion, that’s what I’ve been thinking too. But I don’t think it’s the only sensible one, since I can move the pixels with same distance freely on the circle around the center, maybe producing a less distortet picture.

But I’ll implement that and see what happens.

I'm not sure how you could measure or evaluate the distorsion level of your resulting image...

What's for sure is that the straightforward polar projection is the projection which minimizes the distance between a given point in the original picture and its distorted image.
A drawback of this transform is that it is well-defined only if the contour lines are star-convex with respect to the center.

Another way to do things could be to preserve the "density" of points along a given contour line. I'm not sure how to state this formally, but here is a general idea :
take two points a and b on the same contour of distance r to the center and let d be the distance between a and b on the contour. If you call a' and b' the respective images of a and b by the transform, then you could place a' and b' on the circle of radius r, in such a fashion that d' = 2 pi r d / l, where l is the length of the original contour and d' is the distance between a' and b' on the circle.
I'm not sure if I'm clear :/

Yes, that’s an idea I’ve been having as well. It would require that I can measure lengths along the contour line (which ATM is just a set of points), but that should be solvable.

Another problem arising this way would be keeping to adjacent contour lines “similar”, so that close points with different distance arrive at close points in the morphed image.

Anyway, the polar projection yields surprisingly nice results for now, so I’ll make the GUI more pleasant and maybe put it online, before trying out other algorithms.

I miss-clicked and deleted this very good comment by Peter Pöschl:

=====================

I have a different solution:

For each point P_i find the transformed point W_I by adding a displacement D_i:

P_i + D_i = W_i

with the constraint |W_i| = w(P_i), where w is the prescribed walking distance.
Now vector D_i can be written as

D_i = D_ir + D_it

where D_ir is a radial component parallel to W_i and D_it is tangential
(orthogonal to W_i).
We can also think of D_ir and D_it as force vectors, trying to distort the
picture.
D_ir is kind of a stretching force mandated by the walking distance constraint.
D_it tries to rotate the pixel around the origin.

Considering the immediate neighbours j,k,l,m of P_i, there should be no tangential forces in P_i, thus

((D_jt - D_it)+(D_kt - D_it)+(D_lt - Dit)+(D_mt - D_it)) = 0

This is a huge nonlinear equation system to be solved iteratively.

=====================

Very nice idea Peter, I’ll definately try it out once I finshed the GUI to work comfortably with that stuff.

<strong>Trackback:</strong> <a href="http://www.joachim-breitner.de/blog/archives/288-Distance-Map-Morpher.html">Distance Map Morpher</a><br />
In my last blog entry, I was asking for algorithms to morph a picture so that contour lines would become circles. I then started to write some python/gtk/cairo code to actually visualize my thought. It came out relatively pretty, so I made the GUI aroun

You may be able to use level set methods for curvature flow. See for example

http://math.berkeley.edu/~sethian/

and also the work by Guillermo Sapiro.

It's so easy, after I remembered my Analysis lectures :-)

The contour lines in the plot of your walking-distance function should be
interpreted as _equipotential lines_ of a potential function w(x,y).
This function better be well behaved, so I assume exactly one maximum with

w(x_0,y_0) = 0 (the home point)

(an appropriate example would be w(x,y) = - sqrt(x^2 + y^2) for a level meadow)

Perpendicular to the contour lines are _flow lines_, which all should
(for a well posed problem) converge on the home point. The angle with which the
flow line hits the home point uniquely identifies the flow line.

Thus each point in the original plane can be identified through the potential
at the point and the home-hitting angle of the flow line through the point.

How to calculate the angle alpha for a given point (u,v)?

F(x,y) = ( d_x w(x,y), d_y w(x,y) )

is a vector valued function, the gradient flow, with d_x/d_y denoting partial
derivatives.
(this implies, that you have to approximate your distance function with
something continuously differentiable, e.g. parabolic splines.
Beware of the home point, the function might not be integrable there, as in
the example: F= - 1 / w * (x, y), w looks like a cone!).

Integration of the ordinary differential equation

(x', y') = F(x,y)

with starting point (u,v) with any ODE solver(e.g. Runge Kutta) yields the
flow line. The integration has to be stopped, when it is sufficiently near the
home point. The polar angle of the last calculated point is our alpha.

Nice, that’s a good starting point for me.

Unfortunately, the function is not that well-behaved (sometimes a “bubble” moves out of a equipotential line when going farther away, and appears shortly after). Maybe these regions can be detected well and then ignored. The have always have a center with F'=0, maybe that helps.

For calculating F I can approximate it by looking at surrounding pixels and just take apropriate differences, right?

> ... bubble ...
This I think depends, if the 'tip' of the buble represents a local maximum or a minimum.
If it is a maximum, you have a real problem, because it will attract all flow lines in the neighbourhood, thus acting as a false home point. You would have to determine the 'basin of attraction' of this maximum, as there will be no path from any point in this area to your home point.
If it is a minimum it should be no problem because then it is a repulsive fix point, so if you aren't starting exactly on the minimum you will drift away.

I was wondering from the beginning how you are constructing your distance function. There are infinitely many paths to go from A to B and each one can have a different length. So which one do you choose?

> ... calculating F
(Hesitant) yes (I have to admit that I have next to no experience in numerics). You would have to ensure, that your numeric differencials are continuous, otherwise your ODE stepper will have problems.
That's why I mentioned splines. With them you can approximate the w-function with piecewise polynomials of order 2 which also have continuus derivatives.

Should minima be problematic and maxima ok? Since the potential is a distance funtion, the home point ist the global minimum and attracts all lines. Anyways, I think these bubbles will then not be problems as they are clearly local maxima and repulse the flow lines. (I’m assuming I start the from a point in the image and find my way to the home point).

The distance funktion works like this: For the edges of the given graph, I calculate the minimum distance from the home point along the graph. Then for all other points, I consider all edges of the facet they are on, calculate the orthogonal footpoint on these edges and take the one that’s the fastest.

In other words: I walk along streets as long as I can, and then go orthogonal to a street to the points that I want to reach.

I didn't say it too explicitly:

My potential function is the negative of your distance function.

Thus a maximum of my potential is an attractive fixpoint, while a minimum is a repulsive fixpoint. You can see this if you evaluate the F of my example potential near the origin: The result vector points towards the origin.

If you change the sign of the potential, attractive fixpoints become repulsive, and vice versa.

Ok, then it fits.

> The polar angle of the last calculated point is our alpha.

This can be simplified:
Let F(eps) be the last point before you hit the home point.

Then

- w(x,y) / | _F_(eps) | * _F_(eps)

is the point where you put your transformed pixel

All this requires floating point, right (otherwise I will only have a few directions in which are points drawn at all). Sounds slow :-)

BTW, what do I do next, assuming I do it numerically? I then have the map T which assigns every point in the old image it’s new position – what do I do best with those points in the targed that are not hit by T? (Another reason to interpolate and then work with “real” functions I guess.)

> what do I do best with those points in the targed that are not hit by T?

Depends, how you want to color target areas:

Solution 1:

Use as starting coordinates of the integrations the _corners_ of each pixel.
As the transformation preserves neighbourhood you end up with four points forming a distorted quadrangle. Paint this with the color of the original pixel.

Solution 2:

(Not sure, if this always works. I guess, this would impose some invertibility conditions on the potential function)
o In the target system, calculate the length of the vector B into the blank area pixel and store as b.
o multiply B by a tiny factor, so that the resulting vector B' has a length of several times eps, the smallest representable number not 0.
o use B' as starting point for an ODE solver for gradient field -F (i.e. start in the vicinity of the home point, which is now repulsive, because we turned the potential upside down and integrate away until we hit a region, where the potential is equal b.
o color target pixel with color of pixel where we stopped integrating.

Nice, solution 2 is what I had in mind as well :-)

I just hope that rounding errors don’t make too much problems (i.e. that large parts of the resulting images were painted with information from the same pixel). But I think I’ll try that next.

Ok, I tried this idea (solution 2), but the result was not promising. It turned out that the local maxima of the height function do indeed attract the integration path (which follows the gradient showing up the path), so most paths diverge fast from the edges of the graph and get stuck in the middle of the facets.

I like the idea of following a path, but the gradient of the height function is not the right thing it seems. It should be some vector field that encourages paths that have a good walking time to euclidian distance ration, so that the paths follow the edges for a while (without all going the same way, of course). Not sure how thats easily possible.

The algorithm is in the code (although only in a rough approxmation that does not avoid rounding errors very hard), you can try it out if you get the program to run.

Or maybe the transformation should have some notion of density, so that it avoids mapping a lot of points to the same area and insteads spreads them out more even on the circles of the output image.

I’ll now improve the radial layout by calculating it backwards, e.g. finding the source pixel for each pixel of the target image, which allows for aliasing and should improve image quality a lot.

I've thought a lot over the weekend and found several nasty problems.

> ... that the local maxima ...

I believe, that local maxima should not happen and are an artifact of the
interpolation scheme.


> I like the idea of following a path ...
It was a good idea, but it does not work with your problem.

Why? I haven't looked at your code, but I suspect that you modell streets as
a graph, where nodes are connected by straight lines.
What would be my expectation for the flow lines be? It's per your definition
the shortest path I would take from any given point to the home point.
Beeing on-road allows me to travel at a high speed, while off-road I have to
crawl at a snail's pace. There being a potential, there is a force field
(the gradient field) which basically tells me where the next street is.
Thus I go on the shortest path perpendicular to the next street, take a
90 degree turn onto the street and zip along at high speed on the shortest
combination of streets to the home point.
The result: As soon as I turn onto a street, I loose any information about
the point where I started.
==> you're busted, because your whole picture is transformed
into some lines going out from the home point.


Can you do better?
Yes, use streets with finite width. But now interpolation of the potential
becomes much, much harder. Using an analogy with electric potentials it
turns out, that you have to solve a
http://en.wikipedia.org/wiki/Poisson's_equation.

And the result? Starting from an off-road point I would go roughly perpendicular
to the next street, making a sharp turn at the rim of the road and zoom along,
while more and more flow lines are entering at the rim and pushing me toward
the middle of the street. No passing and crossing of flow lines, if you please,
so everything stays strictly in order when we arrive at the home point.
In reality, I predict that roundoff errors will have a field day and drown
any meaningfull result.


I had another insight:
Consider an equilateral triangle ABC with the home point h dividing _AB_ in
two equal parts and m the middle of the triangle. This is the simplest case
of streets forming a closed loop.

How should, qualitatively, the potential look like?
Mostly, the equipotential lines should be parallel to the sides of the triangle.

But along _Am_ and _Bm_ the potential must be discontinuous, because south
of them you can go directly onto _AB_ while otherwise you have to take the long
way via A or B. As a consequence, neighbouring points across this rift lines
will be placed far away from each other.

There is no discontinuity along _mC_, but this could be a special case,
valid in an equilateral triangle only.
Anyway, I think that this should be treated as another rift.
If the triangle is part of a bigger graph, this rift will propagate into
neighbouring cells until it reaches the edge of the image.
Every other loop in the graph will have to be cut open, too, and introduce
another rift.

The result will be a simply connected graph, or multi-tree.
Maybe you could let this tree 'relax' in the sense, that whenever n edges
meet at a node, set the angle between adjacent edges to 360°/n.
Then stretch the edges, such that the distance of each node to the home point
equals the minimum walking distance. Now apply this stretch+rotate transform
plus the original polar projection to the off-road pixels.

Very interesting thoughts.

In my implementation, I don’t have the problem with the triangle, as it takes the fastest paths via any of the three edges, not just via the closest.

I do have rifts along edges though, when you can reach them faster directly than via the graph, but are not allowed to cross them.

I have someone clean up the radial layout (now does antialiasing and actually looks quite nice), so for now I’m not pursuing the idea. But that should not have to stop from having more ideas, and you are welcome to put them here. Maybe I’ll come back later to the idea.

s/someone clean/somewhat cleaned/

I do not employ people for my fun hacks :-)

> - w(x,y) / | F(eps) | * F(eps)

wrong sign, must be '+ w(x,y)' because w is negative.