Cylindrical Mirror Anamorphosis – Development of a System of Equations

I have a personal project to develop an application in Python to generate a Cylindrical Mirror Anamorphosis from a bitmap. I’ve developed working code that involves iteration but is slow as a result. I’m trying to develop some equations to replace my iteration routines.

With reference to Figure 1, I’m trying to find the coordinates of $$p$$ given $$v$$, $$a$$, $$o$$ and $$r$$.

For my project $$x_o$$, $$y_o$$ and $$x_v$$ equal zero.

The gradients of $$A$$, $$B$$ and $$C$$ are therefore: $$gradA = \frac{y_p-y_a}{x_p-x_a} \; \; \; \; \; gradB = \frac{y_p}{x_p} \; \; \; \; \; gradC = \frac{y_p-y_v}{x_p}$$

Given the angle of incidence and reflection Ø are equal then: $$\frac{gradB-gradA}{1+gradB.gradA}=\frac{gradC-gradB}{1+gradC.gradB}$$

Gives: $$\frac{\frac{y_p}{x_p} - \frac{y_p-y_a}{x_p-x_a}}{1+\frac{y_p}{x_p}.\frac{y_p-y_a}{x_p-x_a}}=\frac{\frac{y_p-y_v}{x_p}-\frac{y_p}{x_p}}{1+\frac{y_p-y_v}{x_p}.\frac{y_p}{x_p}}$$

Which condenses down to:

$$x_p^3(y_a+y_v)-x_p^2(y_px_a+y_vx_a)-x_p2y_vy_py_a+y_p^2(x_py_a+y_vx_a+y_vx_p)-y_p^3x_a=0$$

Given $$p$$ is on a circle of radius $$r$$:

$$y_p=\sqrt{r^2-x_p^2}=\sqrt{(r+x_p)(r-x_p)} = \sqrt{r+x_p}\sqrt{r-x_p}$$

I now want to substitute $$y_p$$ into the above equation and develop an expression for $$x_p$$ and this is where I could do with some guidance from this community. How, if at all possible, do I free $$x_p$$ from the square root term?

I’ll apologies now if my question appears naive.

• I have included your figure into your (well written) question. I have also simplified the tags (no linear algebra here, just analytic geometry). May 3 at 5:53
• You have two bivariate polynomial equations of (total) degree $3$ and $2$ respectively. Eliminating one variable such as $y_p$ between the two equations can be done using resultants, but is computationally expensive, and will produce an equation of degree $6$ in $x_p$.
– dxiv
May 3 at 6:26
• @Barcus: Your first relation has the form $A_3 y_p^3+A_2 y_p^2+A_1 y_p+A_0=0$; your second, $y_p^2=B$. Typically, non-linear elimination is done with resultants or Grobner bases, but that's not needed here; rather, "all you have to do" is finagle the first equation to have only even powers of $y_p$. You can do that by re-writing it as $$y_p(A_3y_p^2+A_1)=-(A_2 y_p^2+A_0)\quad\to\quad y_p(A_3B+A_1)=-(A_2 B+A_0)$$ then squaring: $$y_p^2(A_3B+A_1)^2=(A_2B+A_0)^2\quad\to\quad B(A_3B+A_1)^2=(A_2B+A_0)^2$$ Now expand. In this case, cancellation leaves you with a "mere" quartic in $x_p$.
– Blue
May 3 at 7:29
• @ Barcus: You want to define a mirror curve or profile that reflects between fixed points $v$ to $a?$ May 3 at 18:56
• @Narasimham: The circle of radius $r$ represents the mirror surface. I want to find $p$ for any given $v$ and $a$. May 4 at 6:35

Two very different answers: A) and B)

Part A) As you will digitalize the scene (your reference to a bitmap) consider preprocessing once for all (for a fixed point $$v$$) the "plane" in this way (a kind of ray tracing):

Fig. 1.

not stopping there, but building out of this information a rectangular array (corresponding to pixels) where each entry is the value corresponding to the angular value "transported" by (one of the) rays crossing it (assuming a certain density and/or using interpolation, with adequate rounding of course).

Do you see what I mean ?

Here is the little Matlab program that I wrote for the generation of figure 1. It uses complex numbers geometry (I am used to it for issues involving angles) with $$r=1$$ WLOG.

   clear all;close all;hold on;
axis equal;axis([-3,3,0,3]);
h=3; % position of point v
plot(exp(i*(0:0.01:pi)),'r'); % half circle
b=asin(1/h);% limit angle
for a=b:pi/100:pi-b
z0=exp(i*a);
z=z0*conj(conj(z0)*(h*i-z0)); % symmetry using conjugation
plot([i*h,z0,z0+3*z],LS,'on'); % ray tracing
end;


Part B): (close to your approach)

Take a look at the following figure (We assume that, up to a change of scale, $$r=1$$.)

Fig. 2.

Indeed, we can as well WLOG, assume that, up to rotation, line $$(L)=BE$$ (where $$B,E$$ are resp. the begin and end point) is horizontal with equation $$y=h$$.

I now assume you know what an harmonic division and a crossratio are.

In this case, you must know that the interior line bisector and the exterior line bissector issued from the solution point $$M$$ (this exterior bissector being the tangent to the circle in $$M$$) define on line $$(L)$$ an harmonic division $$(B,E;X_1,X_2)=-1$$, i.e., with notations of Fig. 2, $$B(a,h),E(b,h),X_1(x_1,h),X_2(x_2,h)$$, a "cross-ratio" of abscissas equal to $$-1$$:

$$\frac{\frac{x_1-a}{x_1-b}}{\frac{x_2-a}{x_2-b}}=-1\tag{(a)}$$

As the equations of line OM and the tangent in $$M$$ are resp.

$$y=x \tan \theta \ \ \text{and} \ \ x \cos \theta+ y \sin \theta=1$$

it is easy to deduce that the abscissas of $$X_1$$ and $$X_2$$ are resp.:

$$x_1=\frac{h}{\tan \theta} \ \ \text{and} \ \ x_2=\frac{1-h \sin \theta}{\cos \theta} \tag{(b)}$$

Plugging expressions (b) into (a) gives equation

$$\frac{h\cos \theta - a\sin \theta}{h\cos \theta - b\sin \theta} = - \ \frac{1-h\sin \theta - a\cos \theta}{1-h\sin \theta - b\cos \theta}\tag{(c)}$$

equivalent to:

$$2h\cos \theta - (a+b)\sin \theta +h(a+b)(\sin^2 \theta-\cos^2 \theta) + 2(ab-h^2)\cos \theta\sin \theta=0 \tag{(d)}$$

yielding a fourth degree polynomial equation in $$t=\tan \frac12 \theta$$:

$$2h(1-t^4)-2(a+b)(t+t^3)+h(a+b)(-t^4+6t^2-1)+4(ab-h^2)(t-t^3)=0 \tag{(e)}$$

when using Weirstrass substitution formulas:

$$\cos \theta = \dfrac{1-t^2}{1+t^2}, \ \ \ \ \sin \theta = \dfrac{2t}{1+t^2}.$$

Appendix: A different approach I had proposed at first is to consider auxiliary curves depicted on fig.3:

Fig. 3.

Indeed, the different ellipses featured on the figure are loci (plural of locus) of points $$M$$ with a given constant optical path length of the form

$$BM+ME=constant,$$

($$B$$ and $$E$$ being their common foci... plural of focus :). We are looking (general principle in optics) to a smallest length optical path with an $$M$$ belonging to the unit circle. Only one of them (see remark 1 below) is candidate: "the" ellipse tangent to the unit circle (in point $$I$$ which is the solution we are looking for).

How is this accessible to computation ? Plainly by the computation of roots of polynomials that are more difficult to explicitate than with the previous solution.

Remark 1: In fact, two ellipses are tangent to the unit circle, the one we have considered and another one, much larger, externally tangent to the unit circle, with no physical meaning.

Remark 2: These concentric ellipses are well understood if one knows the generic parametric equations:

$$x=\cosh(u)\cos(t), \ \ y=\sinh(u)\sin(t)$$

of ellipses $$(E_u)$$ having their foci $$F$$ and $$F'$$ in $$(-1,0)$$ and $$(1,0)$$ resp. Matlab program given below simply maps foci $$F$$ and $$F'$$ onto foci $$B$$ and $$E$$ by a convenient similitude operation (rotation, enlargment, translation):

Matlab program for figure 3:

   close all;axis equal;grid on
t=0:0.01:2*pi;plot(exp(i*t),'r'); % unit circle
xb=0;yb=3; % start point
xe=2;ye=2; % end point
mx=(xb+xe)/2;my=(yb+ye)/2; % midpoint
dx=(xe-xb)/2;dy=(ye-yb)/2; % direction of focal axis
for u=0.1:0.1:1.2
c=cosh(u);s=sinh(u);
co=cos(t);si=sin(t);
x=mx+dy*s*co+dx*c*si;
y=my-dx*s*co+dy*c*si;
plot(x,y);hold on;
end;

• @Blue In fact the quartic equation I find (I am not sure it is the same as yours) is decomposable into 2 quadratics. May 3 at 13:55
• 1) Letter $M$ is for the generic point of a given ellipse, which is the locus of points whose sum of distances to the two foci is (a) constant ("gardener's construction", see animation here). 2) b and e are the complex numbers associated to point B and E, bu the important geometrical thing is that m is the midpoint and d gives the direction of line segment $[EB]$ (analogous to $\frac12\vec{EB}$). If I have time tomorrow, i will try to give a complete solution with real numbers. I must also rectify what I have said about the two quadratics. May 4 at 20:42
• I have converted my second program with complex numbers into a program with real coordinates, easier to understand. May 4 at 22:03
• I have now a simpler explanation using cross ratio (instead of ellipses] giving a rather simple fourth degree equation ((e)) (unfortunately not reducible to 2 quadratics as I thought at first). May 5 at 13:10
• Connected: researchgate.net/publication/… May 5 at 15:05

Sixth degree equation is cumbersome to put it on to pixel work. At the end you want only the P coordinates.

A numerical procedure is suggested here.

Taken input data with unit circle and the origin:

$$\{x_v=0, y_v=1.4, x_a=1.2,y_a=1.2 \}$$

For equal $$\theta$$ inclination dot product divided by modulus products gives a result ( details omitted I am sure you can work it out, you showed your work. Please check for any typos.)

$$\frac{1-y_v \sin \theta}{\sqrt{1-2 \sin \theta\; y_v + y_v^2}}=\frac{1-x_A \cos \theta +y_A\sin \theta \;y_v}{\sqrt{1-2 (x_A \cos \theta + y_A \sin \theta) + x_A^2+y_A^2}}$$

which formula can be used for iteration. A radial line through O gives P. Out of six solutions I get two imaginary, two extraneous, two real points, out of which one point is the required internal bisector and the other an external bisector (not shown).

verifies $$\theta$$ equality reasonably well for assumed data by Geogebra construction.

• Mathematica code is available if required. May 4 at 19:12
• I appreciate your contribution. I already have an iterative solution to my project but I will look into your proposal as it might be more efficient than mine. May 4 at 19:52