How does 'arc' work?

[luser- -droog]

This may seem like a silly question to some. But how does 'arc'
construct curves to approximate the arc section?

I know you can make a full circle with 4 Bezier curves
whose control points are arranged on a square. (ASCII ART:

* o *
* *
o o
* *
* o * )

So for smaller intervals, do you just trim-up this template?
Constructing as many full quadrants as necessary and smaller
pieces for the ends?

Say you want to do: (0,0) [45..180] r=1.
You want the -x,+y quadrant in full,

* o
*
o

But then, for the partial piece in +x,+y,
hm, I suppose you actually have to compute
the Bezier approximation, huh?

x = r cos t
y = r sin t
45 <= t <= 90

So the end points, of course, are

x = 1 cos 45 = sqrt 2
y = 1 sin 45 = sqrt 2

x = 1 cos 90 = 0
y = 1 sin 90 = 1

and the control points are
...

[bugbear]

luser- -droog wrote:
> This may seem like a silly question to some. But how does 'arc'
> construct curves to approximate the arc section?
>
> I know you can make a full circle with 4 Bezier curves
> whose control points are arranged on a square. (ASCII ART:
>
> * o *
> * *
> o o
> * *
> * o * )
>

You can *approximate* a full circle with 4 Bezier curves.

Depending on the final size of the circle (and thus the pixel
error it engenders) you might want to use more than 4 beziers.

BugBear

[jdaw1]

For a quarter circle of unit radius, PostScript has the internal control points at (0,0.552) and (0.552,0). The number is likely to be an approximation to (√957 - 21)÷18 ≈ 0.55196759, which minimises ∫(r²–1)² wrt angle for a 90° cubic Bézier.

0.552 has the radius too large by a factor of 1 +0.000212 at 17.39° and 90–this, and too small by a factor of 1 –0.000151 at 45°. So the radius is correct to one part in at least ≈4714. At 1200 d.p.i., that’s correct to within a pixel for a diameters ≤ 7.857″ ≈ 19.96 cm.

[luser- -droog]

On Saturday, November 23, 2013 5:06:27 PM UTC-6, jdaw1 wrote:
> For a quarter circle of unit radius, PostScript has the internal control points at (0,0.552) and (0.552,0). The number is likely to be an approximation to (√957 - 21)÷18 ≈ 0.55196759, which minimises ∫(r²–1)² wrt angle for a 90° cubic Bézier.
>
>
>
> 0.552 has the radius too large by a factor of 1 +0.000212 at 17.39° and 90–this, and too small by a factor of 1 –0.000151 at 45°. So the radius is correct to one part in at least ≈4714. At 1200 d.p.i., that’s correct to within a pixel for a diameters ≤ 7.857″ ≈ 19.96 cm.

Ok. Thank you. That's very useful. I suppose for the remaining
segments I can find the points by an application of deCasteljau,
splitting until the desired angle and then discarding half of the
list. Might not be terribly efficient. But this is pass 1:
make it work.

[luser- -droog]

I have found quite a few pdfs online that go into excruciating
detail. But I haven't read them enough times to really
understand it all yet. :)

-
sleepy

[luser- -droog]

Thank you, that is useful.
I could perform the calculations in device coordinates
with a little extra effort. I'm almost more afraid of overkill
than deficient precision. The later can more easily detected
by eye. :)

I must confess, I'm usually afraid of your posts, BugBear!
A few times in the past, they have totally blown my design
apart, and caused lots of extra (ultimately-well-worthwhile)
hard work. Keep it up, brother! :)

[luser- -droog]

Haven't quite solved it yet, but I've got some partial results worth
sharing, I think. Here are the resources I've found.

A gentle introduction to what's going on:
http://whizkidtech.redprince.net/bezier/circle/

Don Lancaster's cubic spline page:
http://www.tinaja.com/cubic01.shtml
Particularly: http://www.tinaja.com/glib/bezarc1.pdf

Arcs to Beziers and vice versa:
http://itc.ktu.lt/itc354/Riskus354.pdf

http://hansmuller-flex.blogspot.fr/2011/04/approximating-circular-arc-with-cubic.html

Implementation in processing:
http://www.flong.com/blog/2009/bezier-approximation-of-a-circular-arc-in-processing/

So, what I've learned so far is, given as input the center point
(x,y), the radius r, and starting and ending angles a1 a2; the end-
points are pretty easy P1 = (x+r*cos(a1),y+r*sin(a1)), P4 =
(x+r*cos(a2),y+r*sin(a2)). But the control points are trickier.

You can get the correct direction easily enough by taking the first
derivative of the function, but the magnitudes are trickier.

I tried skimming through the above resources and using the magic
constant often mentioned, but it's only good for 90 degrees.

So I read further and discovered the way to get the control points
was to select them such that the center point of the curve lies
on the circle. Easier said than done.

I tried using the deCasteljau subdivision algorithm (backwards),
but after pages and pages of algebra, I could only get symmetrical
equations of the two control points, so upon substituting all the
crucial pieces disappear.

So I think the trick is to translate to the center, scale by the
radius, and rotate so the center of the angles is on the x axis.

Then x2 == x3 and y2 == -y3. And that simplifies everything.

I've still got some math to do, but this crucial insight (already
mentioned in the Riskus paper), appears to be quite crucial.

[luser- -droog]

On Wednesday, December 11, 2013 1:02:08 AM UTC-6, luser- -droog wrote:
> On Sunday, December 1, 2013 5:58:18 AM UTC-6, luser- -droog wrote:
>
> > On Sunday, December 1, 2013 5:48:50 AM UTC-6, luser- -droog wrote:
>
> >
> > > On Saturday, November 23, 2013 5:06:27 PM UTC-6, jdaw1 wrote:
> > > > For a quarter circle of unit radius, PostScript has the internal control points at (0,0.552) and (0.552,0). The number is likely to be an approximation to (√957 - 21)÷18 ≈ 0.55196759, which minimises ∫(r²–1)² wrt angle for a 90° cubic Bézier.
>
> > > >
> > > > 0.552 has the radius too large by a factor of 1 +0.000212 at 17.39° and 90–this, and too small by a factor of 1 –0.000151 at 45°. So the radius is correct to one part in at least ≈4714. At 1200 d.p.i., that’s correct to within a pixel for a diameters ≤ 7.857″ ≈ 19.96 cm.
>
> > >
> > > Ok. Thank you. That's very useful. I suppose for the remaining
> > > segments I can find the points by an application of deCasteljau,
> > > splitting until the desired angle and then discarding half of the
> > > list. Might not be terribly efficient. But this is pass 1:
> > > make it work.
> >
> >
> >
> > I have found quite a few pdfs online that go into excruciating
> > detail. But I haven't read them enough times to really
> > understand it all yet. :)
> >
>
> Haven't quite solved it yet, but I've got some partial results worth
> sharing, I think. Here are the resources I've found.
>
> A gentle introduction to what's going on:
> http://whizkidtech.redprince.net/bezier/circle/
>
> Don Lancaster's cubic spline page:
> http://www.tinaja.com/cubic01.shtml
> Particularly: http://www.tinaja.com/glib/bezarc1.pdf

Yes! This is the one. Last page, the second to last equations.
x1 = (4 - cos(a)) / 3
y1 = ((1 - cos(a))(cos(a) - 3)) / (3*sin(a))

So, repeating the revelation above:

x1 = (4 - cos(a)) / 3
y1 = ((1 - cos(a))(cos(a) - 3)) / (3*sin(a))

and of course

x2 = x1
y2 = -y1
x0 = cos(a)
y0 = cos(a)
x3 = x0
y3 = -y0

Woohoo!

So here's a portable postscript implementation for angles
from 0 to 90.

/arcbez { % draw single bezier % x y r angle1 angle2
5 dict begin
/mat matrix currentmatrix def
5 3 roll translate % r angle1 angle2
3 2 roll dup scale % angle1 angle2
2 copy exch sub /da exch def % da=a2-a1
add 2 div rotate
/da_2 da 2 div def
/sin_a da_2 sin def
/cos_a da_2 cos def
4 cos_a sub 3 div % x1
1 cos_a sub cos_a 3 sub mul
3 sin_a mul div % x1 y1
neg
1 index % x1 y1 x2(==x1)
1 index neg % x1 y1 x2 y2(==-y1)
cos_a sin_a neg % x1 y1 x2 y2 x3 y3
cos_a sin_a % ... x0 y0
4 { 8 2 roll transform } repeat
mat setmatrix
4 { 8 2 roll itransform } repeat
{ currentpoint pop pop } stopped { moveto }{ lineto } ifelse
curveto
end
} def

And for full circles:

% x y r angle1 angle2 arc -
% append a circular arc to the current path
/arc {
2 { dup 360 gt { dup 360 div truncate 360 mul sub } if exch } repeat
2 { dup 0 lt { 360 add } if exch } repeat
2 copy exch sub dup 90 gt { % recurse
.5 mul % x y r a1 a2 da/2
6 copy sub arc
3 2 roll add exch arc
}{ pop
% draw single bezier
arcbez
} ifelse
} def

[luser- -droog]

Dammit. y0 = sin(a) of course.

[luser- -droog]

> > y0 = [sin](a)

It can be done with half the transforms by using 3 matrices.

> > And for full circles:
> >
> > % x y r angle1 angle2 arc -
> > % append a circular arc to the current path
> > /arc {
> > 2 { dup 360 gt { dup 360 div truncate 360 mul sub } if exch } repeat
> >

There's a bug here. Trimming 360 means 360 become 0 and a full
circle won't draw anything. So bump the 360s up to 720, or just
first one.

> > 2 { dup 0 lt { 360 add } if exch } repeat
> > 2 copy exch sub dup 90 gt { % recurse
> > .5 mul % x y r a1 a2 da/2
> > 6 copy sub arc
> > 3 2 roll add exch arc

These are drawn in the wrong order. It results in
criss-crosses through the center when drawing a circle.
S.b.:
6 copy

3 2 roll add exch arc

sub arc

[tlvp]

On Wed, 11 Dec 2013 00:20:11 -0800 (PST), luser- -droog wrote:

>> x0 = cos(a)
>>
>> y0 = cos(a)
>
> Dammit. y0 = sin(a) of course.

Ah, the perils of "cut, paste, and edit", when we forget the "edit" part
:-) . Cheers, -- tlvp
--
Avant de repondre, jeter la poubelle, SVP.

[luser- -droog]

On Friday, February 28, 2014 2:47:58 AM UTC-6, tlvp wrote:
> On Wed, 11 Dec 2013 00:20:11 -0800 (PST), luser- -droog wrote:
>
> >> x0 = cos(a)
> >> y0 = cos(a)
> >
> > Dammit. y0 = sin(a) of course.
>
> Ah, the perils of "cut, paste, and edit", when we forget the "edit" part
>
> :-) . Cheers, -- tlvp
>

Very generous, but I'm pretty sure I typed that. I can get a little
"lyxdexic" in the wee hours.

--
call me the space cowboy.

[jdaw1]

> So I read further and discovered the way to get the control points
> was to select them such that the center point of the curve lies
> on the circle. Easier said than done.

If you do that for the 90° case, you’d get a result of (√2–1)×4/3 ≈ 0.55228475, also close the PostScript value of 0.552.

For the 90° case I have a really good spreadsheet with all the cases you might want. My contact details are at http://www.jdawiseman.com/author.html — send me an email and I’ll send .xls.

[jdaw1]

You have forced me to do some simple work on this, and the answer is actually quite easy.

http://www.jdawiseman.com/2014/Bezier_part_circle_20140502.nb

http://www.jdawiseman.com/2014/Bezier_part_circle_20140502.pdf

Justify that the inner control points are a distance away from the ends of (4/3)×Tan[Angle÷4].

[jdaw1]

Slightly improved, showing the size of the maximum radius error (varying, for small angles, as the sixth power of that angle).
http://www.jdawiseman.com/2014/Bezier_part_circle_20140508.nb
http://www.jdawiseman.com/2014/Bezier_part_circle_20140508.pdf

[jdaw1]

All the above reasoning was about 90° angles. What about smaller angles?

For the Bézier curve, obviously the two outer control points are at the y=sin() x=cos() ends of the segment of the circle. Obviously the two inner control points are at an angle that is tangent to that segment of the circle. This leaves one missing parameter: how far away are the inner control points?

For 90° of the unit circle, 0.552, as discussed above. I didn’t realise until recently that 90° is a special case: extrapolating the distances for smaller angles to 90° would give 0.552285, the 0.552 crossover happening near 89.9582°.

For code to generate numbers, and subsequent analysis, see
http://www.jdawiseman.com/2016/20160401_bezier_arcs.ps
http://www.jdawiseman.com/2016/20160401_bezier_arcs.xlsx

Comment and improvements welcomed.

Do any readers of this forum know the actual angle→distance formula?

[luser- -droog]

You did see this message of mine, right?
https://groups.google.com/d/msg/comp.lang.postscript/YiuM6lj5ngY/4Lo9AW0M2xgJ

Aside from the cos/sin bungle, what do you need that these formulas
don't get you?

x0 = cos(a)
y0 = sin(a)

x1 = (4 - cos(a)) / 3
y1 = ((1 - cos(a))(cos(a) - 3)) / (3*sin(a))

x2 = x1
y2 = -y1

[Scott Hemphill]

jdaw1 <jdawi...@gmail.com> writes:

> All the above reasoning was about 90° angles. What about smaller angles?

What above reasoning?

For real PostScript? Feed this to your printer. Modify it as desired.

========================================================================
%!PS

% What is the distance of an inner Bezier curve control point to the outer
% point as a function of angle in an "arc" operator

/dist
{
/angle exch def
gsave
newpath
0 0 1000 0 angle arc
{pop pop} {pop pop} {pop pop pop pop exch pop} {} pathforall
grestore
1000 div
}
bind def

/tmpstr 20 string def

% I use paper size "letter"

/in { 72 mul } bind def
/Times-Roman findfont 10 scalefont setfont
/lineskip 12 def
/nl { 1 in y moveto /y y lineskip sub def } bind def
/y 10 in def
nl

10 10 90 { dup tmpstr cvs show ( ) show dist tmpstr cvs show nl } for

showpage
========================================================================

Scott
--
Scott Hemphill hemp...@alumni.caltech.edu
"This isn't flying. This is falling, with style." -- Buzz Lightyear

[jdaw1]

On Sunday, 3 April 2016 05:35:24 UTC+1, luser- -droog wrote:

> Aside from the cos/sin bungle, what do you need that these formulas
> don't get you?
>
> x0 = cos(a)
> y0 = sin(a)
> x1 = (4 - cos(a)) / 3
> y1 = ((1 - cos(a))(cos(a) - 3)) / (3*sin(a))
> x2 = x1
> y2 = -y1
> x3 = x0
> y3 = -y0

Are they correct? I am expecting, with an angle of π/2 = 90°, that Sqrt[(x0-x1)^2+(y0-y1)^2] should be about 0.552. I think your formula returns 2×Sqrt[13]/3 ≈ 2.4037. Even with a = π/4 = 45°, it returns a distance of 1.09565.

Or have I made an error?

[luser- -droog]

Very cool. There may indeed be something wrong with my implementation.

gs gives:
10 0.0582147
20 0.11665
30 0.175536
40 0.235099
50 0.295591
60 0.357264
70 0.420398
80 0.485293
90 0.552282

whereas xpost gives:
10 0.115506
20 0.231321
30 0.346987
40 0.461987
50 0.575741
60 0.687515
70 0.796415
80 0.901347
90 1.000998

[Scott Hemphill]

I changed the above program to make an arc with radius 1024 instead of
1000 so that presumably the processor can do the division with no
additional roundoff error.

My HP LaserJet Pro M452dw gives:

10 0.0582143
20 0.116652
30 0.175537
40 0.235103
50 0.295593
60 0.357266
70 0.420398
80 0.485294
90 0.55228

There is something funky going on with small radii, which is why I used
a larger radius to begin with. If I change the radius from 1024 to 1, I
get this:

10 0.0582886
20 0.116699
30 0.175537
40 0.235107
50 0.295593
60 0.357239
70 0.42041
80 0.485291
90 0.552307

[luser- -droog]

On Tuesday, April 5, 2016 at 11:48:46 PM UTC-5, luser- -droog wrote:
> On Sunday, April 3, 2016 at 9:40:55 AM UTC-5, Scott Hemphill wrote:
>
> > For real PostScript? Feed this to your printer. Modify it as desired.
> >
> > ========================================================================
> > %!PS
> >
> > % What is the distance of an inner Bezier curve control point to the outer
> > % point as a function of angle in an "arc" operator
> >

<snip>

> Very cool. There may indeed be something wrong with my implementation.
>
> gs gives:
> 10 0.0582147
> 20 0.11665
> 30 0.175536
> 40 0.235099
> 50 0.295591
> 60 0.357264
> 70 0.420398
> 80 0.485293
> 90 0.552282
>
> whereas xpost gives:
> 10 0.115506
> 20 0.231321
> 30 0.346987
> 40 0.461987
> 50 0.575741
> 60 0.687515
> 70 0.796415
> 80 0.901347
> 90 1.000998

I've double-checked my code and my source material:
http://www.tinaja.com/glib/bezarc1.pdf

I don't see where the problem is (nor indeed what to
do about if it is located). My output looks like arcs
to me.

[luser- -droog]

I'm not sure! Referring to this file of yours:
http://www.jdawiseman.com/2014/Bezier_part_circle_20140508.pdf
I'm not sure I understand it. Does the "Knots" stuff accurately
model the Bezier interpolation?

Empirically, my code seems to draw things that look like circular
arcs to me.

And my math-guru spent a few days with the equations and
came up with the same results. The derivation only worries
about the middle point of the curve, where the naïve recursive
deCasteljau subdivision will later chop the curve in twain and
select the point as the join. We then just assume by symmetry
that further subdivision points will be on the same arc.

I also don't see in your notes any easy formula that I can grab
and apply to the problem of calculating the control points given
theta.

On the other hand, based on the outputs from Scott's program,
it appears my lengths are fairly consistent in being close
to double the "desired" measure. So it's probably worthwhile
to try scaling my tangent lengths and compare outputs. Maybe my
circles are not as pretty as they might be...

[jdaw1]

On Wednesday, 6 April 2016 05:48:46 UTC+1, luser- -droog wrote:
> Very cool. There may indeed be something wrong with my implementation.
>
> gs gives:
> 10 0.0582147
> 20 0.11665
> 30 0.175536
> 40 0.235099
> 50 0.295591
> 60 0.357264
> 70 0.420398
> 80 0.485293
> 90 0.552282
>
> whereas xpost gives:
> 10 0.115506
> 20 0.231321
> 30 0.346987
> 40 0.461987
> 50 0.575741
> 60 0.687515
> 70 0.796415
> 80 0.901347
> 90 1.000998

If xpost really had the 90° inner control points more than a radius away from the outer control points, the Bézier curve would look more square than circular. You say they look circular. So that isn't the problem. Is this a reporting error: are your curves good, but the numbers in these tables wrong?

[luser- -droog]

On Monday, April 18, 2016 at 11:40:24 PM UTC-5, luser- -droog wrote:
> On Sunday, April 3, 2016 at 3:54:55 PM UTC-5, jdaw1 wrote:
> > On Sunday, 3 April 2016 05:35:24 UTC+1, luser- -droog wrote:
> >
> > > Aside from the cos/sin bungle, what do you need that these formulas
> > > don't get you?
> > >
> > > x0 = cos(a)
> > > y0 = sin(a)
> > > x1 = (4 - cos(a)) / 3
> > > y1 = ((1 - cos(a))(cos(a) - 3)) / (3*sin(a))
> > > x2 = x1
> > > y2 = -y1
> > > x3 = x0
> > > y3 = -y0
> >
> > Are they correct? I am expecting, with an angle of π/2 = 90°, that Sqrt[(x0-x1)^2+(y0-y1)^2] should be about 0.552. I think your formula returns 2×Sqrt[13]/3 ≈ 2.4037. Even with a = π/4 = 45°, it returns a distance of 1.09565.
> >
> > Or have I made an error?
>
> I'm not sure! Referring to this file of yours:
> http://www.jdawiseman.com/2014/Bezier_part_circle_20140508.pdf
> I'm not sure I understand it. Does the "Knots" stuff accurately
> model the Bezier interpolation?
>
> Empirically, my code seems to draw things that look like circular
> arcs to me.
>
> And my math-guru spent a few days with the equations and
> came up with the same results.

This is not correct. The result for y1 is different.
The new result is:

y1 = (1 - x1*cos(a)) / sin(a)

or fully expanded:

y1 = ((3 - cos(a))(1 - cos(a)))/ sin(a)

> The derivation only worries
> about the middle point of the curve, where the naïve recursive
> deCasteljau subdivision will later chop the curve in twain and
> select the point as the join. We then just assume by symmetry
> that further subdivision points will be on the same arc.
>
> I also don't see in your notes any easy formula that I can grab
> and apply to the problem of calculating the control points given
> theta.
>
> On the other hand, based on the outputs from Scott's program,
> it appears my lengths are fairly consistent in being close
> to double the "desired" measure. So it's probably worthwhile
> to try scaling my tangent lengths and compare outputs. Maybe my
> circles are not as pretty as they might be...

The output does appear better. The edges seem smoother.
Have not yet checked Scott's measurements.

[jdaw1]

I think you want the half-way point to have exactly the correct radius. Please confirm that is so.

If so, that which I call ‘z’, the distance of the inner control points from the outer, of course along the tangent to the circle, satisfies z == Tan[a/4] * 4/3.

With a=90°, Tan[a/4] = Tan[22½°] = √2-1, so z = (√2-1)*4/3 ≈ 0.55228474983. Adobe’s value of 0.552 won’t have exactly the correct radius at 45°. Indeed, at 45° Adobe’s radius is Sqrt[499849/500000] ≈ 0.999848988597778.

There are several sensible grounds for choosing z, all of which give values close to 0.552. Correct radius at half-way? Correct fill area (z = 2-Sqrt[(22-5π)/3])? Max and min radii symmetrical around ½ (z = root of a order-12 polynomial: 2097152, -23592960, 117129216, -342814720, 668905728, -921504768, 883891456, -515232000, 96817248, 63288000, -28977264, 466560, 18225, so z ≈ 0.551915)? Minimal absolute dRadius/dAngle (z ≈ 0.552031)? They’re all close to 0.552.

But what is Adobe’s formula for non-90° angles?

[jdaw1]

Droog: I think that I have found your error.
http://www.jdawiseman.com/2016/20160420_Arc_Radius_Droog.nb
http://www.jdawiseman.com/2016/20160420_Arc_Radius_Droog.pdf

[luser- -droog]

Thanks for this! It looks like everything I might need to
dig further. But, I'm kinda burnt-out on xpost and taking
a few days on other projects. I've spent the last few weeks
debugging memory problems raised by a change to the memory
tables. The good news is it's almost 50% faster and suddenly
collecting credible quantities of garbage (for a long time,
the gc results have been suspiciously low). The bad news is
there's still (at least) one bug that's still proving
difficult to diagnose. If the gc runs at certain unfortuitous
moments, it hits a bunch of errors. As a result of this week's
efforts, this is now a regular postscript VMerror and not a
crash. But it's still bad, but can be avoided by compiling
with -DXPOST_NO_GC or by tweaking the
XPOST_GARBAGE_COLLECTION_THRESHOLD so it doesn't happen in
just that wrong spot. Even weirder, running on different
OSes produce slightly different errors. It's always
"can't find table for ent XXXXXX" with some out-of-bounds
number. But mys2, Cygwin, and Linux all report a different
ent number.

I discovered that part of the problem in the numbers
revealed by Scott's program is that xpost is erroneously
producing 2 beziers for the 90-degree arc, for some reason.
Re-reading the code (quoted below), I still don't see why it's
doing that.

One apparent difference is that you're taking a tangent to
find the z value, whereas the idea I've been following
was to solve that equation for y1. I'm not sure I see where
the calculation for z comes from, but perhaps that will be
more clear when my brain is less mushy.

And I could always peek at gs or even ralpage for further
reference (and stealing).

For reference, my current implementation (with un-updated
postscript original in comments). This version does not (yet)
incorporate the suggested calculations in the above pdf.
https://github.com/luser-dr00g/xpost/blob/19ace0b18a946ddd95a9e489a2b3759dd14f1221/src/lib/xpost_op_path.c#L511

/*
% packs the center-point, radius and center-angle in a matrix
% then performs the simpler task of calculating a bezier
% for the arc that is symmetrical about the x-axis
% formula derived from http://www.tinaja.com/glib/bezarc1.pdf
/arcbez { % draw single bezier % x y r angle1 angle2 . x1 y1 x2 y2 x3 y3 x0 y0
DICT
%5 dict
begin
%/mat matrix def
5 3 roll mat translate pop % r angle1 angle2
3 2 roll dup mat1 scale mat mat concatmatrix pop % angle1 angle2

2 copy exch sub /da exch def % da=a2-a1

add 2 div mat1 rotate mat mat concatmatrix pop

/da_2 da 2 div def
/sin_a da_2 sin def
/cos_a da_2 cos def
4 cos_a sub 3 div % x1
1 cos_a sub cos_a 3 sub mul
3 sin_a mul div % x1 y1
neg
1 index % x1 y1 x2(==x1)
1 index neg % x1 y1 x2 y2(==-y1)
cos_a sin_a neg % x1 y1 x2 y2 x3 y3
cos_a sin_a % ... x0 y0

4 { 8 2 roll mat transform } repeat
%pstack()=
end
}
dup 0 10 dict
dup /mat matrix put
dup /mat1 matrix put
put
bind
def
*/


static
void _transform(Xpost_Matrix mat, real x, real y, real *xres, real *yres)
{
*xres = mat.xx * x + mat.xy * y + mat.xz;
*yres = mat.yx * x + mat.yy * y + mat.yz;
}

static
Xpost_Object _arc_start_proc;

static
int _arcbez(Xpost_Context *ctx,
Xpost_Object x, Xpost_Object y, Xpost_Object r,
Xpost_Object angle1, Xpost_Object angle2)
{
Xpost_Matrix mat1, mat2, mat3;
real da_2, sin_a, cos_a;
real x0, y0, x1, y1, x2, y2, x3, y3;

xpost_matrix_scale(&mat1, r.real_.val, r.real_.val);
xpost_matrix_translate(&mat2, x.real_.val, y.real_.val);
xpost_matrix_mult(&mat2, &mat1, &mat3);
xpost_matrix_rotate(&mat2, (real)(((angle1.real_.val + angle2.real_.val) / 2.0) * RAD_PER_DEG));
xpost_matrix_mult(&mat3, &mat2, &mat1);

da_2 = (real)(((angle2.real_.val - angle1.real_.val) / 2.0) * RAD_PER_DEG);
sin_a = (real)sin(da_2);
cos_a = (real)cos(da_2);
x0 = cos_a;
y0 = sin_a;
x1 = (real)((4 - cos_a) / 3.0);
//y1 = - (((1 - cos_a) * (cos_a - 3)) / (3 * sin_a));
y1 = (1 - x1*cos_a) / sin_a;

x2 = x1;
y2 = -y1;

x3 = cos_a;
y3 = -sin_a;
_transform(mat1, x0, y0, &x0, &y0);
_transform(mat1, x1, y1, &x1, &y1);
_transform(mat1, x2, y2, &x2, &y2);
_transform(mat1, x3, y3, &x3, &y3);
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(x1));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(y1));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(x2));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(y2));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(x3));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(y3));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(x0));
xpost_stack_push(ctx->lo, ctx->os, xpost_real_cons(y0));
return 0;
}

static
int _arc(Xpost_Context *ctx,
Xpost_Object x, Xpost_Object y, Xpost_Object r,
Xpost_Object angle1, Xpost_Object angle2)
{
double a1 = angle1.real_.val;
double a2 = angle2.real_.val;
while (a2 < a1)
{
double t;
t = a2 + 360;
a2 = t;
}
if ((a2 - a1) > 90)
{
_arc(ctx, x, y, r, xpost_real_cons(a1), xpost_real_cons((real)(a2 - ((a2 - a1)/2.0))));
_arc(ctx, x, y, r, xpost_real_cons((real)(a1 + ((a2 - a1)/2.0))), xpost_real_cons(a2));
}
else
{
//Xpost_Object path = _cpath(ctx);
//int pathlen = xpost_dict_length_memory(xpost_context_select_memory(ctx, path), path);
_arcbez(ctx, x, y, r, xpost_real_cons(a1), xpost_real_cons(a2));
xpost_stack_push(ctx->lo, ctx->es, xpost_operator_cons_opcode(_curveto_opcode));
xpost_stack_push(ctx->lo, ctx->es, _arc_start_proc);
/*
if (pathlen)
xpost_stack_push(ctx->lo, ctx->es, xpost_operator_cons_opcode(_lineto_opcode));
else
xpost_stack_push(ctx->lo, ctx->es, xpost_operator_cons_opcode(_moveto_opcode));
*/
}
return 0;
}

static
int _arcn(Xpost_Context *ctx,
Xpost_Object x, Xpost_Object y, Xpost_Object r,
Xpost_Object angle1, Xpost_Object angle2)
{
real a1 = angle1.real_.val;
real a2 = angle2.real_.val;
while (a2 > a1)
{
double t;
t = a2 - 360;
a2 = t;
}
if ((a1 - a2) > 90)
{
_arcn(ctx, x, y, r, xpost_real_cons(a1), xpost_real_cons(a2 + (real)((a1 - a2)/2.0)));
_arcn(ctx, x, y, r, xpost_real_cons(a1 - (real)((a1 - a2)/2.0)), xpost_real_cons(a2));
}
else
{
//Xpost_Object path = _cpath(ctx);
//int pathlen = xpost_dict_length_memory(xpost_context_select_memory(ctx, path), path);
_arcbez(ctx, x, y, r, xpost_real_cons(a1), xpost_real_cons(a2));
xpost_stack_push(ctx->lo, ctx->es, xpost_operator_cons_opcode(_curveto_opcode));
xpost_stack_push(ctx->lo, ctx->es, _arc_start_proc);
/*
if (pathlen)
xpost_stack_push(ctx->lo, ctx->es, xpost_operator_cons_opcode(_lineto_opcode));
else
xpost_stack_push(ctx->lo, ctx->es, xpost_operator_cons_opcode(_moveto_opcode));
*/
}
return 0;
}

[luser- -droog]

I tried coming at it from a different angle. I re-read chapter 7
http://www.math.ubc.ca/~cass/graphics/manual/pdf/ch7.pdf
and tried to write the arc->Bezier as a general tool for
parametric curves. But when it finally mostly worked, the
tangent vectors are again the wrong lengths. So I visually
tweaked the constant so the first few iterations seem to match up.

$ cat bez.ps
%!

% a r g s [/a /r /g /s] n . -
/argsbegin {
dict begin
dup length 1 sub -1 0 {
2 copy get
4 3 roll def
pop
} for
pop
} def

/2= {
2 copy exch =only (,)print =
gsave
2 copy newpath
currentlinewidth 3 mul 0 360 arc fill
grestore
} def

% t0 t1 f g f' g' N . - (creates path)
/piecewisebezier {
{t0 t1 f g f' g' N} 16 argsbegin
/h t1 t0 sub N div def
(h=)print h =
/h3 h 1.8 div def % h 3 div makes the vectors too short
(h3=)print h3 =
/t t0 def
/updatexy {
(t=)print t =
/x t f def
/y t g def
/x' t f' h3 mul def
/y' t g' h3 mul def
%(x,y=)print x =only (,)print y =
} def
updatexy

x y
2=
moveto % (f(t0), g(t0))

N {
x x' add
y y' add % (f(t)+f'(t)*dt/3, g(t)+g'(t)*dt/3)
2=

/t t h add def
updatexy
x x' sub
y y' sub % (f(t+dt)-f'(t+dt)*dt/3, g(t+dt)-g'(t+dt)*dt/3)
2=

x y % (f(t+dt), g(t+dt))
2=

curveto
%lineto pop pop pop pop
} repeat
end
} def

300 400 translate
currentlinewidth
300 dup dup scale
div setlinewidth

0 1 {90 mul cos}{90 mul sin}{90 mul sin neg}{90 mul cos} 1 piecewisebezier stroke
2 1 3 {
0 1 {90 mul cos}{90 mul sin}{90 mul sin neg}{90 mul cos} 7 6 roll piecewisebezier stroke
} for
showpage



coordinates:
h=1.0
h3=0.555556
t=0
1.0,0.0
1.0,0.555556
t=1.0
0.555556,1.0
0.0,1.0
h=0.5
h3=0.277778
t=0
1.0,0.0
1.0,0.277778
t=0.5
0.903525,0.510688
0.707107,0.707107
0.510688,0.903525
t=1.0
0.277778,1.0
0.0,1.0
h=0.333333
h3=0.185185
t=0
1.0,0.0
1.0,0.185185
t=0.333333
0.958618,0.339625
0.866025,0.5
0.773433,0.660375
t=0.666667
0.660375,0.773433
0.5,0.866025
0.339625,0.958618
t=1.0
0.185185,1.0
0.0,1.0

[luser- -droog]

On Friday, April 22, 2016 at 1:32:22 PM UTC-5, luser- -droog wrote:
> I tried coming at it from a different angle. I re-read chapter 7
> http://www.math.ubc.ca/~cass/graphics/manual/pdf/ch7.pdf
> and tried to write the arc->Bezier as a general tool for
> parametric curves. But when it finally mostly worked, the
> tangent vectors are again the wrong lengths. So I visually
> tweaked the constant so the first few iterations seem to match up.
>

Improved. Uses .552 for 90 degrees and proportionally decreasing
lengths for shorter angles and choosing the number of segments
by the base-90 logarithm of the angle span.



josh@LAPTOP-ILO10OOF ~

$ cat bez.ps
%!

% a r g s [/a /r /g /s] n . -
/argsbegin {
dict begin
dup length 1 sub -1 0 {
2 copy get
4 3 roll def
pop
} for
pop
} def

/2= {

%2 copy exch =only ( )print =

gsave
2 copy newpath
currentlinewidth 3 mul 0 360 arc fill
grestore
}

%pop {} %uncomment to suppress points as little circles

def

% t0 t1 f g f' g' N . - (creates path)
/piecewisebezier {
{t0 t1 f g f' g' N} 16 argsbegin
/h t1 t0 sub N div def

%(h=)print h =
%/h3 h 3 div def % makes the vectors too short
%/h3 h 1.8 div def % == .5556 mul
%/h3 h .552282 mul def
/h3 h .00613 mul def % .00613 == .552 / 90
%(h3=)print h3 =
/t t0 def
/updatexy {
%(t=)print t =

/draw {
gsave
{moveto}
{lineto}
{6 4 roll lineto 4 2 roll lineto lineto}
{closepath} pathforall
stroke
grestore
stroke
}
%pop {stroke}
def

%180 3.14159 div =

305 400 translate

{
currentlinewidth
300 dup dup scale
div setlinewidth

/A 0 def
/B 90 def

{
%A B {cos}{sin}{sin neg}{cos} 1 piecewisebezier draw
1 1 1 {
%2 mul
A B {cos}{sin}{sin neg}{cos} 7 6 roll piecewisebezier draw
} for
} pop

{
10 10 90 {
A exch {cos}{sin}{sin neg}{cos} 1 piecewisebezier draw
} for
} pop
} pop

/myarc {
{x y r a1 a2} 10 argsbegin
a1 a2
{cos r mul}{sin r mul}{sin neg r mul}{cos r mul}
a2 a1 sub log 90 log div ceiling cvi
piecewisebezier
end
} def

0 0 300 0 90 myarc draw

10 10 90 {
0 0 2 index 2 mul 0 5 4 roll myarc draw
} for

10 10 90 {
0 0 2 index 2 mul 90 5 4 roll 180 add myarc draw
} for

showpage

josh@LAPTOP-ILO10OOF ~
$

[luser- -droog]

On Saturday, April 23, 2016 at 10:49:47 PM UTC-5, luser- -droog wrote:
> On Friday, April 22, 2016 at 1:32:22 PM UTC-5, luser- -droog wrote:
> > I tried coming at it from a different angle. I re-read chapter 7
> > http://www.math.ubc.ca/~cass/graphics/manual/pdf/ch7.pdf
> > and tried to write the arc->Bezier as a general tool for
> > parametric curves. But when it finally mostly worked, the
> > tangent vectors are again the wrong lengths. So I visually
> > tweaked the constant so the first few iterations seem to match up.
> >
>
> Improved. Uses .552 for 90 degrees and proportionally decreasing
> lengths for shorter angles and choosing the number of segments
> by the base-90 logarithm of the angle span.
>

Re-write/cleanup. Drop-in replacement for arc.

josh@LAPTOP-ILO10OOF ~
$ cat arc.ps
%!

/.args {

dup length 1 sub -1 0 {
2 copy get 4 3 roll def pop
} for
pop
} def

/.argsbegin {
dict begin
.args
} def

/arc {
{x y r a1 a2} 20 .argsbegin
matrix currentmatrix x y translate

a1 a2 {cos r mul}{sin r mul}{sin neg r mul}{cos r mul}
a2 a1 sub log 90 log div ceiling cvi

{t0 t1 f g f' g' N} .args

/h t1 t0 sub N div def

/hdelta h .00613 mul def
/updatexy {

/x t f def
/y t g def

/x' t f' hdelta mul def
/y' t g' hdelta mul def
} def

/t t0 def updatexy
x y moveto

N {
x x' add y y' add

/t t h add def updatexy
x x' sub y y' sub

x y curveto
} repeat

setmatrix
end
} def

{
/X 306 def
/Y 400 def

X Y 300 0 90 arc stroke

10 10 90 {
X Y 2 index 2 mul 0 5 4 roll arc stroke
} for

10 10 90 {
X Y 2 index 2 mul 90 5 4 roll 180 add arc stroke
} for

showpage
} pop

(xpost/arcdist.ps) run



And Scott's arcdist program yields:

10 0.0612999
20 0.122599
30 0.183900
40 0.245200
50 0.306500
60 0.367999
70 0.429099
80 0.490400
90 0.551700

[jdaw1]

On Sunday, 24 April 2016 04:49:47 UTC+1, luser- -droog wrote:
> Improved. Uses .552 for 90 degrees and proportionally decreasing
> lengths for shorter angles and choosing the number of segments
> by the base-90 logarithm of the angle span.

I’m sure you don’t mean “base-90 logarithm”. IIRC, Adobe Distiller uses 90° pieces, and if necessary a smaller last piece (so dividing by 90°).

Adobe Distiller doesn’t quite use proportionality either, though that isn’t bad. As far as I can tell, Adobe might be doing something like
0.005825526 a - 4.49157E-07 a^2 + 4.32805E-08 a^3, with an exception at 90° of exactly 0.552. (Adobe employees: please reveal the exact formula.)

Linear differs from Adobe by the most at ≈52.4°, differing by ≈ 0.0111.

[luser- -droog]

I peeked at ghostscript's implementation:
http://git.ghostscript.com/?p=ghostpdl.git;a=blob;f=base/gspath1.c;h=25644c2bae912fad7e9e9206a973b245fe6d67d0;hb=HEAD

and here's a version that more closely copies that behavior. It partitions the
angle span into orthogonal 90-degree segments, possibly with partial initial
and final <90-degree segments.

$ cat arc.ps
%!

/debug false def % set true to print points to stdout

/.args {

dup length 1 sub -1 0 {
2 copy get 4 3 roll def pop
} for
pop
} def

/.argsbegin {
dict begin
.args
} def

% x y . x-y*floor(x/y)
/fmod {
2 copy div floor mul sub
} def
%45.0 90.0 fmod = % 45.0
%927.0 90 fmod = % 27.0

/point {
gsave

newpath
currentlinewidth 3 mul 0 360 arc fill
grestore

} bind def

/debugpoint {
debug { 2 copy exch =only ( )print = } if
} def

/next_arc_curve {
/t0 a1 def
/t1 anext def
debug {
(next_arc_curve)=
t0 t1 exch =only (..)=only =
} if
/h t1 t0 sub def
/hdelta h .00613647 mul def
/updatexy {
/x t cos r mul def
/y t sin r mul def
/x' t sin neg r mul hdelta mul def
/y' t cos r mul hdelta mul def

} def
/t t0 def updatexy
x y

debugpoint
%moveto

{ currentpoint pop pop } stopped { moveto }{ lineto } ifelse

x x' add y y' add

debugpoint

/t t h add def updatexy
x x' sub y y' sub

debugpoint
x y
debugpoint
curveto
/a1 anext def
} def

/next_arc_quadrant {
/t0 a1 def
/t1 anext def
debug {
(next_arc_quadrant)=
t0 t1 exch =only (..)=only =
} if
/h 90 def
/hdelta .552282 def
/updatexy {
/x t cos r mul def
/y t sin r mul def
/x' t sin neg r mul hdelta mul def
/y' t cos r mul hdelta mul def

} def
/t t0 def updatexy
x y

debugpoint
%moveto

{ currentpoint pop pop } stopped { moveto }{ lineto } ifelse

x x' add y y' add

debugpoint

/t t h add def updatexy
x x' sub y y' sub

debugpoint
x y
debugpoint
curveto
/a1 anext def

} def

/arc {
{x y r a1 a2} 20 .argsbegin
matrix currentmatrix x y translate

debug {
(a1=)print a1 =
(a2=)print a2 =
} if

r 0 lt {
/r r neg def
/a1 a1 180 add def
/a2 a2 180 add def
} if

{ a2 a1 lt not { exit } if % while (a2<a1)
(adjust by 360)=
/a2 a2 360 add def
} loop

{
% first part up to 90 deg
a1 90 fmod 0 ne {
/anext a1 90 div floor 90 mul def
anext a2 lt not {exit} if % goto last part
debug { (first part)= } if
next_arc_curve
} if

% multiples of 90 deg
{
/anext a1 90 add def
anext a2 le not {exit} if % while (anext<a2)
next_arc_quadrant
} loop

exit } loop

% last part, if any
a1 a2 lt {
/anext a2 def
debug { (last part)= } if
next_arc_curve
} if

setmatrix
end
} def

/draw {
false { gsave
{point}
{point}
{point point point}
{} pathforall
grestore } if
true { gsave

{moveto}
{lineto}
{6 4 roll lineto 4 2 roll lineto lineto}
{closepath} pathforall
stroke

grestore } if
stroke
}
%pop {stroke} %uncomment to suppress lines through control points

def

{
/X 306 def
/Y 400 def

X Y 300 0 90 arc draw

10 10 90 { X Y 2 index 2 mul 0 5 4 roll arc draw } for

10 10 90 { X Y 2 index 2 mul 90 5 4 roll 180 add arc draw } for

showpage
} exec%pop

%(xpost/arcdist.ps) run

[luser- -droog]

On Saturday, April 30, 2016 at 11:47:32 PM UTC-5, luser- -droog wrote:
> On Sunday, April 24, 2016 at 6:01:56 AM UTC-5, jdaw1 wrote:
> > On Sunday, 24 April 2016 04:49:47 UTC+1, luser- -droog wrote:
> > > Improved. Uses .552 for 90 degrees and proportionally decreasing
> > > lengths for shorter angles and choosing the number of segments
> > > by the base-90 logarithm of the angle span.
> >
> > I’m sure you don’t mean “base-90 logarithm”. IIRC, Adobe Distiller uses 90° pieces, and if necessary a smaller last piece (so dividing by 90°).
> >
> > Adobe Distiller doesn’t quite use proportionality either, though that isn’t bad. As far as I can tell, Adobe might be doing something like
> > 0.005825526 a - 4.49157E-07 a^2 + 4.32805E-08 a^3, with an exception at 90° of exactly 0.552. (Adobe employees: please reveal the exact formula.)
> >
> > Linear differs from Adobe by the most at ≈52.4°, differing by ≈ 0.0111.
>
> I peeked at ghostscript's implementation:
> http://git.ghostscript.com/?p=ghostpdl.git;a=blob;f=base/gspath1.c;h=25644c2bae912fad7e9e9206a973b245fe6d67d0;hb=HEAD
>
> and here's a version that more closely copies that behavior. It partitions the
> angle span into orthogonal 90-degree segments, possibly with partial initial
> and final <90-degree segments.
>

This update suppresses the degenerate lineto between each segment. And it implements the above cubic formula for the
tangent vector lengths. Scott Hemphill's measurement program
yields (gs):

10 0.0582551
20 0.116678
30 0.17553
40 0.235071
50 0.295563
60 0.357264
70 0.420432
80 0.485327
90 0.552282

(xpost yields the same numbers to an unnecessary length
of decimal places).

josh@LAPTOP-ILO10OOF ~
$ cat arc.ps
%!

/debug true def % set true to print points to stdout

%/hdelta h .00613647 mul def
/hdelta
.005825526 h mul
4.49157E-07 h h mul mul sub
4.32805E-08 h h h mul mul mul add

def
/updatexy {
/x t cos r mul def
/y t sin r mul def
/x' t sin neg r mul hdelta mul def
/y' t cos r mul hdelta mul def
} def
/t t0 def updatexy

initialseg {

x y
debugpoint
%moveto
{ currentpoint pop pop } stopped { moveto }{ lineto } ifelse

/initialseg false def
} if

initialseg {

x y
debugpoint
%moveto
{ currentpoint pop pop } stopped { moveto }{ lineto } ifelse

/initialseg false def
} if

/initialseg true def

} pop%exec%pop

(xpost/arcdist.ps) run

[jdaw1]

A function has been written for those wishing to render a mathematical in few Bézier cubics. Details in
https://groups.google.com/forum/#!topic/comp.lang.postscript/3RIq0Jnwrbo