Find centre of spiral similarity of two lines

Question

I am trying to find a formula for the centre of spiral similarity (see here) of two lines, lineAB and lineCD. My problem comes in the very last step, but I think it might help if I outline the process that got me there.

Start with the following lines and points:

pointA = {a, \[Alpha]};
pointB = {b, \[Beta]};
pointC = {c, \[Gamma]};
pointD = {d, \[Delta]};
lineAB = Line[{pointA, pointB}];
lineAC = Line[{pointA, pointC}];
lineBD = Line[{pointB, pointD}];
lineCD = Line[{pointC, pointD}];

The first step is to find the intersection point pointP1 of lineAC and lineBD. So, I find the functions describing these lines, and solve:

fac[x_] := (\[Gamma] - \[Alpha])/(c - a) x - (a (\[Gamma] - \[Alpha]))/(c - a) + \[Alpha];
fbd[x_] := (\[Delta] - \[Beta])/(d - b) x - (b (\[Delta] - \[Beta]))/(d - b) + \[Beta];
pointP1 = 
  Flatten[{
     x /. Solve[fac[x] == fbd[x], x],
     fac[x] /. x -> x /. Solve[fac[x] == fbd[x], x]
  }];

This works fine.

Next, I find the circumcircles for the triangles {pointP1 ,pointA, pointB} and {pointP1 ,pointC, pointD}:

circlePAB = Circumsphere[{pointP1, pointA, pointB}];
circlePCD = Circumsphere[{pointP1, pointC, pointD}];

This also works, as you can see from evaluating this:

Clear["Global`*"];
Manipulate[
fac[x_] := (\[Gamma] - \[Alpha])/(c - a) x - (
  a (\[Gamma] - \[Alpha]))/(c - a) + \[Alpha];
fbd[x_] := (\[Delta] - \[Beta])/(d - b) x - (
  b (\[Delta] - \[Beta]))/(d - b) + \[Beta];
pointA = {a, \[Alpha]};
pointB = {b, \[Beta]};
pointC = {c, \[Gamma]};
pointD = {d, \[Delta]};
pointP1 = 
 Flatten[{x /. Solve[fac[x] == fbd[x], x], 
fac[x] /. x -> x /. Solve[fac[x] == fbd[x], x]}];
circlePAB = Circumsphere[{pointP1, pointA, pointB}];
circlePCD = Circumsphere[{pointP1, pointC, pointD}];

graphicPointA = {{White, Disk[pointA, 0.5]}, Text["A", pointA]}; 
graphicPointB = {{White, Disk[pointB, 0.5]}, Text["B", pointB]}; 
graphicPointC = {{White, Disk[pointC, 0.5]}, Text["C", pointC]}; 
graphicPointD = {{White, Disk[pointD, 0.5]}, Text["D", pointD]}; 
graphicPointP1 = {{White, Disk[pointP1, 0.5]}, 
  Text["\!\(\*SubscriptBox[\(P\), \(1\)]\)", pointP1]}; 
lineAB = Line[{pointA, pointB}];
lineAC = Line[{pointA, pointC}];
lineBD = Line[{pointB, pointD}];
lineCD = Line[{pointC, pointD}];

Graphics[{{Lighter[Gray], lineAB}, {Lighter[Gray], 
   lineAC}, {Lighter[Gray], lineBD}, {Lighter[Gray], 
   lineCD}, {Darker[Blue], circlePAB}, {Darker[Red], circlePCD}, 
  graphicPointA, graphicPointB, graphicPointC, graphicPointD, 
  graphicPointP1}, ImageSize -> Large]
,
{{a, -7}, -10, 10}, {{\[Alpha], 9}, -10, 10},
{{b, -9}, -10, 10}, {{\[Beta], -7}, -10, 10},
{{c, 9}, -10, 10}, {{\[Gamma], -8}, -10, 10},
{{d, 2}, -10, 10}, {{\[Delta], 8}, -10, 10}
]

A screenshot gives the idea:

The centre of spiral similarity is the second point of intersection (i.e., not pointP1) of circlePAB and circlePCD. So:

RegionIntersection[circlePAB, circlePCD]

But this doesn't evaluate. I just get

BooleanRegion[#1 && #2 &, {Sphere[{-((-b^2 [Alpha] - b c [Alpha] + b d [Alpha]... (etc, etc, etc)

How do I take the final step to find the second point of intersection?

Answer 1

Assume we have 2 straight lines between 2 points: l1 and l2 and want to know the spiral similarity between them. Further assume that line l2 was created from line l1 by first rotating around an unknown center and then stretching relative to the origin.

As an example we choose 2 lines:

SeedRandom[224]
l1 = RandomReal[{-1, 1}, {2, 2}];
l2 = RandomReal[{-1, 1}, {2, 2}];

The stretch factor must be the relation between the length of l1 and l2. And the rotated l1, call it l10, must be the scale factor times l2:

scal = Norm[Subtract @@ l1]/Norm[Subtract @@ l2];
l10 = scal  l2;

The rotation turns the first point of l1 into the first point of l10 and similar with the second points. Therefore the center must be on the perpendicular bisector of l1[[1]] and l10[[1]] and also on the perpendicular bisector of l1[[2]] and l10[[2]]. Therefore, the center is at the cut of both bisectors. To get this point, we define a function for the bisectors "perbi":

perpbi[p1_, p2_, lam_] = (p1 + p2)/2 + lam {{0, -1}, {1, 0}}.(p1 - p2);

p1 and p2 are the two points and lam is some variable, that defines the points on the line. With this the center is:

center = perpbi[l1[[1]], l10[[1]], lam1] /. 
  Solve[perpbi[l1[[1]], l10[[1]], lam1] == 
     perpbi[l1[[2]], l10[[2]], lam2], {lam1, lam2}][[1]]

Finally we can make a nice plot to illustrate all this:

Graphics[{Line[{l1, l2}], Green, Line[l10], Red, PointSize[0.03], 
  Point[center], Lighter[Blue], 
  Line[{{l10[[1]], {0, 0}}, {l10[[2]], {0, 0}}}], Darker[Yellow], 
  Circle[center, Norm[l1[[1]] - center]], 
  Circle[center, Norm[l1[[2]] - center]], Black, 
  Text["l1", {0.2, -0.2}], Text["l2", {-0.8, -0.23}], 
  Text["l10", {-1, -0.6}]}, Axes -> True]

Addendum

If we assume that the rotation and dilation have the same center and we first rotate around this center and then dilate relative to this center, we may proceed like the following:

Again we specify 3 points: p1,..p4 and 2 lines: l1={p1,p2}, l2={p3,p4}. The dilation is again given by the ratio of the lengths of l2 to l1:

Clear["Global`*"]
SeedRandom[224];
{p1, p2, p3, p4} = RandomReal[{-1, 1}, {4, 2}];
{l1, l2} = {{p1, p2}, {p3, p4}};
dilation = Norm[p3 - p4]/Norm[p1 - p2];

We define symbolically the center and the rotation matrix as a function of the angle:

rot[phi_] := {{Cos[phi], -Sin[phi]}, {Sin[phi], Cos[phi]}};
center = {cx, cy};

We know, that l2 is constructed by first rotate l1 and then dilate, all relative to the center. This gives the equations:

eq = l2 == 
   dilation (rot[phi].# & /@ (l1 - {center, center})) + {center, center};

Because of the circular functions "Solve" is not able to solve this, we need the more powerful Reduce. To be able to use the output of Reduce as rules, we use "ToRule". With this we obtain the center point: pc, and the rotation angle: angle:

{pc, angle} = {{cx, cy}, ArcTan[Cos[phi], Sin[phi]]} /. 
  ToRules[Reduce[eq, {cx, cy}]]

For the graphics we calculate the rotation matrix and the rotated points of l1: p10, p20:

rotmat = RotationMatrix[ArcTan[Cos[angle], Sin[angle]]];
{p10, p20} = rotmat.# & /@ ({p1, p2} - {pc, pc}) + {pc, pc};

With this we may draw the graphics:

Graphics[{
  Line[{{p1, p2}, {p3, p4}}], Green, Line[t = ({p10, p20})], Blue, 
  Line[{{t[[1]], pc}, {t[[2]], pc}}], Darker[Yellow], 
  Circle[pc, Norm[p1 - pc]], Circle[pc, Norm[p2 - pc]], Red, 
  PointSize[0.03], Point[pc], Black, Text["l1", {0.2, -0.2}], 
  Text["l2", {-0.8, -0.23}], Text["l10", {-1, -0.5}]
  }, Axes -> True]

Better

But a much more elegant method is to use complex numbers as described e.g. in: Wiki:

Define the following complex numbers:

im = {1, I};
{a, b, c, d} = im.# & /@ {p1, p2, p3, p4};

Then the center is given by:

pcenter = ReIm[(a d - b c)/(a + d - b - c)]