Archimedean Spiral Algorithm

[fooo]

Hello there,

I'm writing a macro which calculates the toolpath for an archimedean spiral. Some controls do have a G code for that others don't.

My idea is to divide the spiral into circular arcs. For each arc I have to calculate the end point located on the spiral's path. The start point is simply the end point of the former arc. Start and end point are connected by a circular arc.

The problem is to find the radius for each particular circular arc with the smallest possible deviation to the path of the perfect spiral.

Does somebody have a formular for that?

Thanks
Frank

[Geof]

If you are writing a macro why don't you simply calculate the spiral there is no reason to approximate it.

[fooo]

Hmm I don't think this is going to work as I think there could be a significant computational load depending on the resolution of the spiral's curvature you want to achieve.

For example if you want a resolution of lets say 0.001" and the length of the curvature is 1" then you end up with 1000 iterations on a set of calculations.

Thank you anyway for your thoughts.

-Frank

[_Eduardo_]

I use this method for a general approximation with circle arcs.

"Given two points and its slopes, find two tangent circle arcs tangents to each point with the minimun radius difference".

The algebraic steps are a bit hard, the solution written in visual basic is:

Input values:
x1,y1,dx1,dy1: Point 1 + Slope
x2,y2,dx2,dy2: Point 2 + Slope

Return values:
xc,yc: Intersection point for circle arcs
r1: Radius at P1
r2: Radius at P2

'--------------------------------------------------------------------------
Function IncenterRR(x1, y1, dx1, dy1, x2, y2, dx2, dy2, xc, yc, r1, r2)

IncenterRR = 0

p = Sqr(dx1 ^ 2 + dy1 ^ 2)
If p < 0.00000001 Then Exit Function
nx1 = dx1 / p
ny1 = dy1 / p

p = Sqr(dx2 ^ 2 + dy2 ^ 2)
If p < 0.00000001 Then Exit Function
nx2 = dx2 / p
ny2 = dy2 / p

xc = x2 - x1
yc = y2 - y1

p = nx1 * ny2 - ny1 * nx2 ' Note: if p>0 => G03
' if p<0 => G02

If Abs(p) < 0.00000001 Then Exit Function

l1 = (nx1 * yc - ny1 * xc) / p ' Distance P1-Po
If l1 <= 0 Then Exit Function

l2 = (ny2 * xc - nx2 * yc) / p ' Distance P2-Po
If l2 <= 0 Then Exit Function

l3 = Sqr(xc ^ 2 + yc ^ 2) ' Distance P1-P2

p = l1 + l2 + l3
q = l2 / p

xc = x1 + q * (xc + l3 * nx1)
yc = y1 + q * (yc + l3 * ny1) ' return value xc,yc: Incenter

p = p / 2 ' Semiperimeter
l1 = p - l1
l2 = p - l2

q = l1 * l2 * (p - l3) / p ' q = Inradius^2

p = 2 * Sqr(q)
r1 = (q + l1 ^ 2) / p ' return value r1: Radius @ P1
r2 = (q + l2 ^ 2) / p ' return value r2: Radius @ P2

IncenterRR = 1

End Function
'--------------------------------------------------------------------------



Applying to Archimedean spiral.

R = Rini + k*Angle (k=pitch/2PI , Angle in radians)

X = R*cos(Angle)
Y = R*sin(Angle)
dX1 = k*cos(Angle) - Y
dY1 = k*sin(Angle) + X






An example in Excel using 45deg steps (calling the VBA function IncenterRR() )


Pitch = 0.5000
Rini = 4.0000


Ang Radius X Y dX dY Xc Yc R1 R2
----------- ------------------------------- -------------------------------
0 4.0000 4.0000 0.0000 0.0796 4.0000 3.7236 1.5446 4.0097 4.0513
45 4.0625 2.8726 2.8726 -2.8164 2.9289 1.5647 3.7830 4.0722 4.1138
90 4.1250 0.0000 4.1250 -4.1250 0.0796 -1.5925 3.8391 4.1348 4.1763
135 4.1875 -2.9610 2.9610 -3.0173 -2.9047 -3.8985 1.6125 4.1973 4.2388
180 4.2500 -4.2500 0.0000 -0.0796 -4.2500 -3.9546 -1.6403 4.2598 4.3013
225 4.3125 -3.0494 -3.0494 2.9931 -3.1057 -1.6604 -4.0139 4.3223 4.3639
270 4.3750 0.0000 -4.3750 4.3750 -0.0796 1.6881 -4.0701 4.3848 4.4264
315 4.4375 3.1378 -3.1378 3.1941 3.0815 4.1294 -1.7082 4.4473 4.4889
360 4.5000 4.5000 0.0000 0.0796 4.5000 ..................
.........................


And the Gcode looks like:

G00 X4.0000 Y0.0000
G03 X3.7236 Y1.5446 R4.0097
X2.8726 Y2.8726 R4.0513
X1.5647 Y3.7830 R4.0722
X0.0000 Y4.1250 R4.1138
X-1.5925 Y3.8391 R4.1348
X-2.9610 Y2.9610 R4.1763
X-3.8985 Y1.6125 R4.1973
X-4.2500 Y0.0000 R4.2388
X-3.9546 Y-1.6403 R4.2598
X-3.0494 Y-3.0494 R4.3013
X-1.6604 Y-4.0139 R4.3223
X0.0000 Y-4.3750 R4.3639
X1.6881 Y-4.0701 R4.3848
X3.1378 Y-3.1378 R4.4264
X4.1294 Y-1.7082 R4.4473
X4.5000 Y0.0000 R4.4889
..............



Regards.

Eduardo.

[lerman]

Eduardo,

So -- I think:

You compute the intersection point of the lines perpendicular to the tangents at the two points of interest. Then you compute the distance from that intersection to each of the points to get the two radii.

But now we want to approximate the curve between the two points with a circular arc. Ideally, we would like to minimize the maximum error. Since (in principle), we know nothing about the curve between those points, we don't have much to work with.

There are somethings we can do, though. We can require that the arc goes through both of those points. Since there are an infinity of arcs that do that, it is not a problem. The simple thing to do is to use a straight line, but that misses the point of the exercise.

It is clear that at each point, we must start at the previous point. We must also end at the next point. If we don't do that, we risk having unbounded error accumulation.

So, one approach might be to pick a radius between the two that were calculated. Find the center corresponding to it and draw the arc between the two points at that radius. One could pick the larger radius, the smaller radius, the average radius, the geometric mean radius, or something else.

Assuming that the curve is reasonably well behaved (whatever that means), I would think that the error would be less than the difference between the two radii. Hmmm. Perhaps reasonably well behaved means that the curvature is monotonic between the two points.

Tell me more. For instance, do you have a quick calculation for the center, given the two points and the radius. (Give me a break, please, it's been almost fifty years since I studied plane geometry.)

Thanks,

Ken

[fooo]

Hi Kenneth,

I haven't evaluated Eduardo's function yet but your idea corresponds with mine:

I calculate two points on the path of the spiral. Then I take the radius of each point (i.e. the distance to the spiral's center) and calculate the average radius wich interconnects the two points. This results into an arc with a displaced midpoint to the spiral's center. This approach seems to work as long as you don't start the spiral in the center (i.e. r = 0).

I'm looking for other possible approaches. I know there's a formular to calculate the radius of curvature. That might be another way to approximate.

I also found this on the net: http://wotan.liu.edu/docis/dbl/wscgws/2003___PCAOSA.htm

Here's a MuPad function to calculate the midpoint. No warrenty of course ;-) Maybe somebody knows an easier way?



/* Calculates the midpoint of a circle through two points (as polar coordinats) on the circumference and the radius. */

GetMidpoint := proc(r1, a1, r2, a2, r)
local p1, p2, b, l, a, i, k;
begin
p1 := PolarToCartesian(r1, DegreeToRadian(a1));
p2 := PolarToCartesian(r2, DegreeToRadian(a2));

/* Bisector of chord. */
b := sqrt((p1[1] - p2[1])^2 + (p1[2] - p2[2])^2) / 2;

/* I.e. radius minus segment's height. */
l := sqrt(r^2 - b^2);

/* Angle. */
a := arctan((p1[1] - p2[1]) / (p2[2] - p1[2]));

/* Do we need to tweak angle? */
if p2[2] < p1[2] then
a := a + PI;
end_if;

/* What direction counterclockwise or clockwise? */
if a2 > a1 then
/* X distance p1 to center */
i := -(sin(a) * b + cos(a) * l);

/* Y distance p1 to center */
k := cos(a) * b - sin(a) * l;
else
/* X distance p1 to center */
i := -(sin(a) * b - cos(a) * l);

/* Y distance p1 to center */
k := cos(a) * b + sin(a) * l;
end_if;

/* Verify results. */
assert(is(i^2 + k^2 ~= r^2));
assert(is((p2[1] - (p1[1] + i))^2 + (p2[2] - (p1[2] + k))^2 ~= r^2));

/* Return p1, p2 and midpoint as cartesian coordinates. */
(p1, p2, p1[1] + i, p1[2] + k);
end_proc;

[_Eduardo_]

The image attached show the method, the clue is evaluate the incenter of the triangle P1-P2-P3.
The amazing property of the incenter is what the circle arcs has minimun radius difference, that mean what the curvature transition is the smoothest.
Of course, this method is only 'good' when you have the curve equation.


>But now we want to approximate the curve between the two points with a
>circular arc. Ideally, we would like to minimize the maximum error. Since (in
>principle), we know nothing about the curve between those points, we don't
>have much to work with.

Having the equation of the curve you may check the deviation.

For any point of the approximation [x,y] ==> R' = sqrt(x^2+y^2) and Angle' = atan(y/x) + n*PI
Note: n = 0 for 1st quadrant, 1 for 2nd & 3rd, 2 for 4th & 5th ,,etc


With a spiral equation R = Ro + k*Angle

the value dR = R' - Ro - k*Angle' is the radial error of the aproximation

Unless the steps are greater than 45-60 deg the error is far less than 0.0001



>Tell me more. For instance, do you have a quick calculation for the center,
>given the two points and the radius. (Give me a break, please, it's been
>almost fifty years since I studied plane geometry.)


Given Pini = [x1,y1] , Pend = [x2,y2] , R and sign = 1 if G03 , -1 if G02

Then
dx = (x2-x1)/2
dy = (y2-y1)/2
p = sign*R/sqrt(dx^2+dy^2)
h = sign*sqrt(p^2-1)

I = dx - h*dy
J = dy + h*dx

I,J distance from the center to Pini

Then the circle center
Ox = x1 + I
Oy = y1 + J


Also
xm = Ox + p*dy
ym = Oy - p*dx

xm,ym Arc middle point




Regards.
Eduardo.