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.
 A: 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;


A: 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.
