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.

enter image description here

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.

  • 2
    $\begingroup$ I have included your figure into your (well written) question. I have also simplified the tags (no linear algebra here, just analytic geometry). $\endgroup$
    – Jean Marie
    May 3 at 5:53
  • 1
    $\begingroup$ 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$. $\endgroup$
    – dxiv
    May 3 at 6:26
  • 1
    $\begingroup$ @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$. $\endgroup$
    – Blue
    May 3 at 7:29
  • 1
    $\begingroup$ @ Barcus: You want to define a mirror curve or profile that reflects between fixed points $v$ to $a?$ $\endgroup$
    – Narasimham
    May 3 at 18:56
  • $\begingroup$ @Narasimham: The circle of radius $r$ represents the mirror surface. I want to find $p$ for any given $v$ and $a$. $\endgroup$
    – Barcus
    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):

enter image description here

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
      z=z0*conj(conj(z0)*(h*i-z0)); % symmetry using conjugation
      plot([i*h,z0,z0+3*z],LS,'on'); % ray tracing

Part B): (close to your approach)

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

enter image description here

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$:


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:

enter image description here

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


($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
      plot(x,y);hold on;
  • 1
    $\begingroup$ @Blue In fact the quartic equation I find (I am not sure it is the same as yours) is decomposable into 2 quadratics. $\endgroup$
    – Jean Marie
    May 3 at 13:55
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jean Marie
    May 4 at 20:42
  • 1
    $\begingroup$ I have converted my second program with complex numbers into a program with real coordinates, easier to understand. $\endgroup$
    – Jean Marie
    May 4 at 22:03
  • 1
    $\begingroup$ 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). $\endgroup$
    – Jean Marie
    May 5 at 13:10
  • 1
    $\begingroup$ Connected: researchgate.net/publication/… $\endgroup$
    – Jean Marie
    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).

enter image description here

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

  • 1
    $\begingroup$ Mathematica code is available if required. $\endgroup$
    – Narasimham
    May 4 at 19:12
  • $\begingroup$ 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. $\endgroup$
    – Barcus
    May 4 at 19:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.