Correcting for a curved film plane?

[Mark Sirota]

Is it possible to use ImageMagick to correct the distortion resulting from a cylindrical film plane?

Specifically, imagine a camera like this one -- it's a pinhole camera with the pinhole at the geometric center of a cylinder, with the film along the surface of the cylinder. I'd like to scan the film, then correct the distortion to a rectilinear projection (as though the film were flat).

I don't think any of ImageMagick's current projection corrections can do this, right?

[anthony]

No a distortion like that is not yet available. But with the right mathematics you can DIY it.

Once you have the distortion you can create a distortion map and process lots of images very very fast.

Not many examples of distortion map has been put online, but some are.
http://www.imagemagick.org/Usage/mapping/

[Mark Sirota]

Thanks for the reply. I'll look into distortion maps.

[fmw42]

test your mathematics with -fx first. it will be slow, but you can perfect it. then move on to the distortion map.

Note my cylinderize script below does just the opposite. It wraps a planar image onto a cylinder. But I have not given your situation much thought.

[anthony]

Mathematically you have rays going through a pinhole to intersect a cylinder. Now you want to map those rays to intersect a flat plane.

However distortion mappings map from destination image (the flat plane) to the source image (a cylinder).
So you need to work out a formula to map the coordinate i,j on a flat plane to a x,y coodinate on the cylinder.

NOTE i,j and x,y are in pixel coordinates. that is a integer value is the center of a pixel, as such you should first add 0.5 to i,j do the transformation then subtract 0.5 from x,y before looking up the pixel value in the source image.

the simplist example of such a distortion is a direct 1 to 1 mapping. You can see this in the FX 'verbose' output example in
http://www.imagemagick.org/Usage/distor ... rt_verbose

That is an direct FX mapping.

Distortion mapping (using the compose method "Distort" or the FX equivelent) is the same except each pixel coordinate is pre-calculated and stored in two gray scale maps - X coordinate in a red channel, and Y coordinate in the green channel.

I really should make some time to rewrite the examples of using distortion and displacement mapping.

[fmw42]

Here is what I believe is the formulae. Taking the first image from your reference on the camera.



infile="P90achop.jpg"
rr=`convert $infile -ping -format "%[fx:w*57/90]" info:`
hh=`convert $infile -ping -format "%h" info:`
w2=`convert $infile -ping -format "%[fx:w/2]" info:`
ww=`convert xc: -format "%[fx:2*$rr*tan(0.5*90/57)]" info:`
ww2=`convert xc: -format "%[fx:$ww/2]" info:`
convert -size ${ww}x${hh} xc: $infile -monitor \
-fx "xx=(i-$ww2); v.p{$rr*atan(xx/$rr)+$w2,j}" +monitor \
P90achop_corrected.jpg

and this will be faster (10 sec vs 47 sec), since the horizontal formula is one-dimensional; that is, it is the same for every line.

infile="P90achop.jpg"
rr=`convert $infile -ping -format "%[fx:w*57/90]" info:`
hh=`convert $infile -ping -format "%h" info:`
w2=`convert $infile -ping -format "%[fx:w/2]" info:`
ww=`convert xc: -format "%[fx:2*$rr*tan(0.5*90/57)]" info:`
ww2=`convert xc: -format "%[fx:$ww/2]" info:`
max=`convert xc: -format "%[fx:$rr*atan($ww2/$rr)+$w2]" info:`
convert -size ${ww}x1 xc: -monitor \
-fx "xx=(i-$ww2); ($rr*atan(xx/$rr)+$w2)/$max" +monitor \
-scale ${ww}x${hh}! tmp.png
time convert tmp.png $infile -monitor -fx "v.p{$max*u,j}" +monitor \
P90achop_corrected.jpg





However, FX is slow and the next step would be to convert this to using a center relative displacement map (i.e. convert ... -compose distort ... )

[Mark Sirota]

Fantastic, thank you! It looks like you've corrected for the X axis, but not the Y.

I've worked out the math for the X axis somewhat differently than you did, at least I think it's different. I need to work out your formula. Assuming the image center is at 0,0 and x1,y1 is the location of an original pixel, x2 = cotan(arccos(x1/r)), where r is the radius of the film plane in pixels. That doesn't take into account the half-pixel error that Anthony mentioned either. I hadn't figured out how to make that formula work in ImageMagick, though. It's possible that your formula and mine are equivalent.

The Y axis needs to be stretched as well, and as with the X axis the amount varies with the horizontal distance from center. But this time it's a simple multiplicative factor rather than trigonometric. I have an answer but haven't yet verified it, so I'm not yet confident enough to post it here.

In any case, I greatly appreciate the interest you're all showing in this little challenge, and I'm learning a lot more about ImageMagick in the process. Thanks!

[fmw42]

The Y axis needs to be stretched as well, and as with the X axis the amount varies with the horizontal distance from center. But this time it's a simple multiplicative factor rather than trigonometric. I have an answer but haven't yet verified it, so I'm not yet confident enough to post it here.

Mark,

I did not see any mention of any y stretch on the reference page you mention. Do you have any further documentation? Why should it vary with distance from the center? Why is it not just a uniform stretch? This does not make sense if the film "plane" is a quarter cylinder (fov=90 degrees) and the sides are straight. I can understand a uniform stretch, thought I did not see any mention, but I could have missed it.

Anyway, if you draw a tangent to a quarter-circle ( looking down along y on the quarter-cylinder) and given R as the radius of the quarter-cirle, then any point on the quarter-circle (your input image) will have its x coordinate given by

x=R*theta

where theta is measured from the radius perpendicular to the tangent (ie. divided symmetrically in half, i.e tangent point at image center)

If you draw any ray from the center of the quarter-circle through the quarter-circle and find where it intercepts the tangent (the output image), then its X coordinate is given by

X=R*tan(theta)

For -fx you need the inverse transform, from output back to input, so combining

theta=atan(X/R)

and

theta=x/R

we get

x=R*tan(X/R)

That is the basis of how I worked it out.

Fred

[Mark Sirota]

The Y axis must stretch as well, and it must stretch further as you move away from center. Here's a quick diagram (I did it in Sketchup rather than ImageMagick, which is shameful, but it's past my bedtime and I'm faster in Sketchup):


As you can see, where the ray trace (red) passes through the curved film surface (green), it is shorter than where it intersects the flat plane behind it (blue). At the center, they'd be the same heights since the curved film surface and flat plane intersect there. The further you move from center the greater the distance between film and flat plane, and so the greater the "stretch" in the Y axis. It's a similar triangles problem.

Off to bed, I'll digest your trigonometry tomorrow!

[fmw42]

Let me make sure I understand. Are you saying that the pinhole is at the bottom of the cylinder or in the middle of it (vertically)?

The projection I used converted to a perspective projection in both axes. The y axis (if the pinhole is in the center) is symmetric perspective onto vertical line (cylinder) about its center (by definition since it is a pinhole camera). The x axis was a perspective projection onto a cylinder, but my equations convert that to perspective. So the result should be equivalent to a 90 perspective projection horizontally and whatever fov it is vertically.

Are you wanting to now convert the perspective to an orthographic projection (parallel rays)?

In my opinion, the perspective image is what one should be working towards.

Now if the pinhole is at the bottom of the cylinder, then one has an asymmetric perspective and I have not taken that into account. But it would seem very odd that they would do that and not have the pinhole in the vertical center of the film plane cylinder.

If I am misunderstanding something, please feel free to correct me.

Fred

P.S. I may be missing something that I have to think about further. That is for this camera the vertical perspective is the same for every column, but in a real 2D perspective image, the rays go higher as one moves further horizontally from the center. I have not accounted for this. If corrected the result will be a top and bottom edge that will look somewhat like an hourglass on its side.

[Mark Sirota]

Yes, that last bit (your P.S.) is exactly right. The top and bottom of the corrected image would be concave. (Presumably, this would then be cropped back to a rectangle, as indicated in my diagram.)

[fmw42]

Right. Sorry, I guess I just misunderstood your comments and your drawing. I still have to think about this further as I am not sure these two aspects are separable. Perhaps to first order you could do that, but I am not yet convinced that you don't need a 2D transformation rather than a horizontal followed by a vertical transformation. I have no time estimate, but if you get any further ideas or equations, feel free to run them by me. You can email me directly if you want at fmw at alink dot net.

[whugemann]

Hello everybody,

the transformation we need over here is illustrated by:


The angle in x-direction is just phi = x/r, and x' = r * tan(phi) , which gives
x' = r tan(x/r).

For the y-direction we apply the theorem of intersecting lines (hope I got that right in English) and reach at
y'/y = sqrt(x'^2 + r^2)/r = sqrt(1 + (x'/r)^2) = sqrt (1 + tan(phi)^2) = 1 / cos (phi) = 1 / cos (x/r)

Thus
y' = y / cos (x/r)

But this is just a quick guess; there should be lots of pages in cylindrical projection on the Web.

[fmw42]

For -fx, you need the inverse transformation. But Wolfgang's solution looks correct to me as I had been thinking about a similar approach over the weekend, though was planning on working on the equations this morning.

so inverting the x equation gives:

x=r*atan(x'/r) as I specified earlier

and inverting his y equation:

y' = y / cos (x/r)

y = y' * cos(x/r)

y = y' * cos(atan(x'/r))

which I think is correct. I will post the result of fx using this shortly.

Thanks Wolfgang for jump-starting me again on this. You diagram and equations were very helpful.

[fmw42]

Here is the results of using fx with the formulae above on a half size image keeping the same output size. I will work on expanding the output size next.




rr=`convert P90achop_sm.jpg -ping -format "%[fx:w*57/90]" info:`
w2=`convert P90achop_sm.jpg -ping -format "%[fx:w/2]" info:`
h2=`convert P90achop_sm.jpg -ping -format "%[fx:h/2]" info:`
echo "rr=$rr; w2=$w2; h2=$h2"
convert P90achop_sm.jpg -monitor \
-fx "xx=(i-$w2); aa=atan(xx/$rr); yy=(j-$h2); u.p{$rr*aa+$w2,yy*cos(aa)+$h2}" +monitor \
P90achop_sm_corrected4.jpg