# Find centre of spiral similarity of two lines

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.

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?

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]


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)]

• Thanks @Daniel Huber. I won't be able to look at this properly until tomorrow morning. In the meantime, can I ask (a) if there is a particular reason why you chose a different method, and (b) if you have any insight into why RegionIntersection didn't work for me? May 6, 2021 at 19:46
• Hi, check circlePAB, circlePCD , they are symbolic, I can not see where you evaluate them. Why did I choose another method? Well, I do not understand why the intersection should be the center. May 7, 2021 at 7:45
• Hi @Daniel Huber. I may be mis-using the word "evaluate". What I mean is that I would like the expression RegionIntersection[circlePAB, circlePCD] to produce two sets of coordinates corresponding to the required intersections of the circles. If that means that I need to use some different form for circlePAB and circlePCD, then I guess that was the answer I was seeking. I'll mark your answer as correct once I've gone through it in detail - but I would really like to know what I was doing wrong; that's how I learn! May 7, 2021 at 10:04
• Sorry, @Daniel Huber - a second query, and rather more important! Your stretching transformation appears to be with respect to the origin rather than that centre of rotation. My understanding is that both rotation and stretching have the same origin. Or have I misunderstood your diagram? May 7, 2021 at 10:21
• You are right I assumed that the center of dilation is the origin. I have to think about the case where the center of dilation and rotation is the same. May 7, 2021 at 11:09