Programmer 7500

home *** CD-ROM | disk | FTP | other *** search

/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / CNC11TP.ZIP / CONIC.DOC < prev next >

Wrap

Text File | 1989-01-23 | 92.7 KB | 2,509 lines

`Minimum' CONICS Version 1.1 01/19/89 Adam Fritz 133 Main Street Afton, New York 13730 Permission is hereby granted to copy and distribute this document and associated rou- tines for any non-commercial purpose. Any use of this material for commercial advantage without prior written consent of the author is prohibited. INTRODUCTION ============ This writeup presents and discusses a collection of rou- tines for drawing conic sections with points and line segments whose coordinates are generated using rotation of coordinates methods. Timing performance comparisons are made to routines using offset error methods, aka Bresenham. This writeup is distributed with two different versions, viz. CNC11TP.ARC and CNC11TM.ARC where the first contains TURBO Pascal (c) 4.0 routines and the second contains TopSpeed MODULA-2 (c) 1.14 routines. See DISTRIBUTION below. DISCUSSION ========== Polygons Vertex coordinates for a regular polygon with radius r might be computed: 2 Pi x = r cos(i ----) i n 2 Pi y = r sin(i ----) i n Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 An algorithm to use these formulae can be outlined: move(r,0) da = 2*Pi/n for i = 1 to n do begin a = i*da draw(r*cos(a),r*sin(a)) end where i is vertex index and n is number of sides. An algorithm based on these formulae would require n sine and cosine evalua- tions or n-1 if the identity of starting and ending points is noted. Suppose, instead, that the sine and cosine at index i are computed and those at i+1 are desired. Let `a' be angle at i and `da' be step angle = 2Pi/n. One way is to use sum angle formulae where: cos(a + da) = cos(a) * cos(da) - sin(a) * sin(da) sin(a + da) = sin(a) * cos(da) + cos(a) * sin(da) and sine and cosine of `da' need be computed just once at first vertex as do sine and cosine of initial vertex angle. Each subsequent vertex requires just four multiplies and two addi- tions. If vertex angle functions are initialized to include an offset rotation angle `ta' then an algorithm to draw a polygon can be outlined: da = 2*Pi/n cosda = cos(da) sinda = sin(da) cosa = cos(ta) sina = sin(ta) move(r*cosa,r*sina) for i = 1 to n do begin t = cosa * cosda - sina * sinda sina = sina * cosda + cosa * sinda cosa = t draw(r*cosa,r*sina) end Routine POLYRA.PAS is a TURBO Pascal 4.0 floating point routine to apply this algorithm and includes the effects of aspect ratio, a specified center, and inverted vertical coordi- nate. Procedure PolyRA in POLYGON.MOD is a TopSpeed Modula-2 1.14 routine with the same features. In general, routines with `*RA*.x' names use this `rotation angle' approach. Computations can be done in fixed point arithmetic rather than floating point. Sines and cosines vary from -1.0 to 1.0 so that if a 16 bit integer is used to represent this range it is convenient to scale sine or cosine at B14. That is, the imagined binary point is to the right of bit 14 when the bits are labeled right to left. For example, 1@B12 = b0001.0000 0000 0000, where the binary point is imagined, 1@B2 = $0004, etc. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 If values are coded fixed point then it is more convenient to make them longint rather than integer for formation of products and dividends. The algorithm above can be rewritten so that the loop is done in fixed point: da = 2*Pi/n icosda = round(cos(da)*16384) isinda = round(sin(da)*16384) icosa = round(cos(ta)*16384) isina = round(sin(ta)*16384) ix = longhi(r*icosa shl 2) iy = longhi(r*isina shl 2) move(ix,iy) for i = 1 to n do begin it = longhi((icosa * icosda - isina * isinda) shl 2) isina = longhi((isina * icosda + icosa * isinda) shl 2) icosa = it ix = longhi(r*icosa shl 2) iy = longhi(r*isina shl 2) draw(ix,iy) end where angle functions are scaled B14 and radius r and display coordinates ix and iy are scaled B0. LongHi is a function that swaps the high order 16 bits of a 32 bit operand to low order 16 bits and does sign extend - see file CONIC which includes several utility values and routines. Similar routines for Modula-2 are found in MathFun.MOD. Obviously, r greater than 32767 will produce overflow. Routine POLYRAF is a fixed point routine using this method that includes aspect correction, offset (discussed below), and coordinate reflection. It also computes rotation on radiused circle scaled B6 rather than unit circle scaled B14. A generic driver is found in DEMOPOLY that enables selection of algo- rithm, and specification of parameters for radius, number of sides, and rotation angle. A driver to estimate average time is presented in BNCHPOLY where timing is compared to TURBO's BGI Circle routine. Timing results with and without a numeric co- processor are tabulated in the PERFORMANCE section below. It is possible to generate sines and cosines, among other things, using second order linear difference equations. For example, the form for such a generator is: z = c * z - z n n-1 n-2 where, for circular functions, `c' is computed 2 * cos(da) and z and z are initialized to any two successive function n-1 n-2 values separated by da. Using this relationship for a polygon generator can be outlined: Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 da = 2*Pi/n cosda2 = 2*cos(da) cos1 = cos(ta) cos2 = cos(ta-da) sin1 = sin(ta) sin2 = sin(ta-da) move(r*cos1,r*sin1) for i = 1 to n do begin cos0 = cosda2 * cos1 - cos2 sin0 = cosda2 * sin1 - sin2 draw(r*cos0,r*sin0) cos2 = cos1 cos1 = cos0 sin2 = sin1 sin1 = sin0 end where a generator is assigned to each sine and cosine which requires just two multiplies and two additions for each vertex. Routine POLYDA is a floating point and POLYDAF is a fixed point routine to apply this `difference' angle algorithm with DEMOPOLY and BNCHPOLY provisions similar to those for the rotation angle routines. Similar routines for circles, ellipses, and hyperbolae use a `*DA*.x' syntax. Incremental Angles If the above routines are used to generate polygons on a graphics plane with finite resolution there is some number of sides n such that the polygon with n+1 sides is not readily distinguished. For a Hercules (tm) graphics adapter with 720x348 resolution this happens for n between 25 and 30. Further, for this value there is no marked distinction between a polygon and a circle generated with an offset error routine, e.g., Bresenham's circle algorithm. This suggests the idea of using polygons to approximate circles. If a regular polygon with n sides is inscribed in a circle with radius r then the maximum distance between circumference and chord is: da dr = r [1 - cos(--)] 2 where da = 2Pi/n is the step angle. If dr is required to be no greater than 0.5, i.e. half a pixel, then: da 1 cos(--) >= 1 - -- 2 2r Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 and the step angle can be estimated, via series expansions: 1 da <= 2 sqrt( - ) r An integer number of sides is convenient, but not essential, but requires: 2 Pi n >= Round( ---- ) da with da readjusted: 2 Pi da = ---- n These relations claim that a polygon with n sides can fit a circle with radius r to within half a pixel if da is computed as indicated. The following table shows r, n, and da for several values. Routine DELCANG generates these values. DELTA ANGLES ================================================== | exact | approximate | adjusted -------------------------------------------------- r | da n | da n | da | (deg) | (deg) | (deg) -------------------------------------------------- 1 120.00000 3 114.59156 3 120.00000 2 82.81924 4 81.02847 4 90.00000 4 57.91005 6 57.29578 6 60.00000 8 40.72827 9 40.51423 9 40.00000 16 28.72302 13 28.64789 13 27.69231 32 20.28358 18 20.25712 18 20.00000 64 14.33329 25 14.32395 25 14.40000 128 10.13186 36 10.12856 36 10.00000 256 7.16314 50 7.16197 50 7.20000 512 5.06469 71 5.06428 71 5.07042 1024 3.58113 101 3.58099 101 3.56436 ================================================== Ignoring symmetry, the number of chords grows with the square root of radius. The number of pixels using Bresenham circle grows linearly with radius. This suggests chord genera- tion is more efficient than Bresenham. This is misleading because the chords are presumably drawn with Bresenham line to produce the same effective number of pixels. Therefore, the comparison is between Bresenham line and Bresenham circle. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 The code for Bresenham line, for both x and y increasing and dx greater than dy looks like: ie := (dx+1) DIV 2 ; LOOP Plot(x1,y1,c) ; IF x1 = x2 THEN EXIT END ; INC(x1) ; INC(ie,dy) ; IF ie > dx THEN DEC(ie,dx) ; INC(y1) ; END ; END ; and the code for Bresenham circle, using octant symmetry and ignoring both aspect correction and symmetry looks like: WHILE ix < iy DO IF ie < 0 THEN ie := ie + 2*iy - 1 ; DEC(iy) ; END ; ie := ie - 2*ix - 1 ; INC(ix) ; PLOT(xc+ix,yc+iy,c) ; END ; If details about symmetry and setup are ignored, then the differences between these loops are the error offset terms. The offset terms can be computed by summing constant terms to form the linear offsets. Routine CircMB shows how this loop struc- ture is only slightly complicated by incremental offset and aspect correction. In effect the aspect corrected routine for circle about doubles the work per pixel over that for line. It is this difference that makes it useful to consider stroke over pixel algorithms. This difference grows with resolution. Rou- tines named `*MB*.x' use offset error methods where MB denotes modified Bresenham - the first published examples of this class of algorithm. See Offset Error Methods somewhere below. Rou- tines named `*MBR*.x' include both rotation and aspect correction. Incremental Angle Generating circles using polygons with two multiplies and two adds is interesting. Is there anything faster? Suppose that the sines and cosines of a step angle were very simple binary numbers, that is, had only a few bits set or not set. Such values permit multiplication to be done with a few shifts and adds. If a sine or cosine had just one bit set then a shift corresponding to the bit position would be a multiply. Similar- Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 ly, floating point multiply might be done with exponent augmenta- tion. Let angles be computed 360/n and their sines and cosines be computed and scaled Bxx. Those angles listed below have minimum distinct set bit counts. Routine INCCANG generates these numbers. There are no such angles scaled B16 for radii less than 2000 pixels using rounding to generate the scaled values. There are some using truncation or a larger set bit count threshold. INCREMENTAL ANGLE B14 ======================================================= n ang r sin cos number of bits (deg) sin cos total ------------------------------------------------------- 50 7.20000 253 $0805 $3F7F 3 -1 4 134 2.68657 1819 $0300 $3FEE 2 -2 4 ======================================================= INCREMENTAL ANGLE B15 ======================================================= n ang r sin cos number of bits (deg) sin cos total ------------------------------------------------------- 100 3.60000 1013 $080A $7FBF 3 -1 4 ======================================================= For example, if the parameters for case n = 50 at B14 are used, then sine and cosine scaled B14 can be generated with the rotation angle algorithm by: it = icosa * icosda - isina * isinda isina = isina * icosda + icosa * isinda icosa = it or, using hex, it = icosa * $3F7F - isina * $0805 isina = isina * $3F7F + icosa * $0805 icosa = it or, using binary, it = icosa * b0011111101111111 - isina * b000010000000 101 isina = isina * b0011111101111111 + icosa * b0000100000000101 icosa = it or, using complementarity, i.e. $3F7F = $4000 - $0080, and sub- stituting ix = icosa, iy = isina, with operands assumed to be longint, Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 it = ix shl 14 - ix shl 7 - iy shl 11 - iy shl 2 - iy iy = iy shl 14 - iy shl 7 + ix shl 11 + ix shl 2 + ix ix = it and, factoring, it = (((ix shl 3 - iy) shl 4 - ix) shl 5 - iy) shl 2 - iy iy = (((iy shl 3 + ix) shl 4 - iy) shl 5 + ix) shl 2 + ix ix = it where the operands are now scaled B28. A further shift left two and capture of the high order 16 bits of a 32 bit longint operand with sign extension properly rescales - see RoundScale16 and LongHi functions in CONIC or MathFun. That is, B28 + B2 - B16 = B14. If instead, operands are scaled B30 then the bit operations can be to the right with `div' rather than `shr' which does not preserve sign. Routine CIRCRAIR implements the rotation angle algorithm with right shifting for the same parameters used above. Routine CIRCRAI2 illustrates application of the parameters for n = 100 with scaling at B15. Routines with names `*I*.x' denote usage of incremental angle values. Incremental Cosine The limited selection of angles above suggests considera- tion of difference equation generators where step angles with cosines having few distinct bits are sought. Because only the cosine, rather than sine and cosine together, is considered, it can be expected there are more step angles and, therefore, more good fitting polygons available. If angles are again computed 360/n and their cosines are scaled Bxx then those angles listed below have minimum set bit counts. Routine INCCOS genrates these tableau. INCREMENTAL COSINE B14 ===================================================================== n ang sin cos 2*cos atan r nbit |nbit| (deg) (deg) 2cos 2cos --------------------------------------------------------------------- 6 60.00000 $376D $2000 $4000 60.00000 4 1 1 100 3.60000 $03FF $3FE0 $7FBF 3.58157 1024 -1 1 139 2.58993 $02EA $3FEF $7FDF 2.61030 1927 -1 1 140 2.57143 $02D4 $3FF0 $7FDF 2.53235 2048 -1 1 ===================================================================== Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 INCREMENTAL COSINE B15 ===================================================================== n ang sin cos 2*cos atan r nbit |nbit| (deg) (deg) 2cos 2cos --------------------------------------------------------------------- 6 60.00000 $6EDA $4000 $8000 60.00000 4 1 1 100 3.60000 $080F $7FBF $FF7F 3.60945 1008 -1 1 141 2.55319 $05BE $7FDF $FFBF 2.57162 1986 -1 1 197 1.82741 $041F $7FEF $FFDF 1.84568 3855 -1 1 ===================================================================== INCREMENTAL COSINE B16 ======================================================================= n ang sin cos 2*cos atan r nbit |nbit| (deg) (deg) 2cos 2cos ----------------------------------------------------------------------- 6 60.00000 $DDB4 $8000 $00010000 60.00000 4 1 1 71 5.07042 $169B $FF00 $0001FDFF 5.06593 512 -1 1 199 1.80905 $081F $FFDF $0001FFBF 1.81833 3972 -1 1 ======================================================================= For example, take the case for n = 100 at B15 and compute sine and cosine scaled, say, B30. Then, ix0 = i2cos * ix1 - ix2 iy0 = i2cos * iy1 - iy2 or, with hex, ix0 = $0FF7F * ix1 - ix2 iy0 = $0FF7F * iy1 - iy2 or, in binary, ix0 = b0000 1111 1111 0111 1111 * ix1 - ix2 iy0 = b0000 1111 1111 0111 1111 * iy1 - iy2 or using complementarity, i.e., $0FF7F = $10000 - $00080, ix0 = (ix - ix div 512) shl 1 - ix2 iy0 = (iy - iy div 512) shl 1 - iy2 Therefore, for this case, where the operands remain scaled at B30, sine and cosine are generated with two shifts and two adds. If TURBO or TopSpeed had a `sar' it would be used rather than `div'. There are many nearby angles with duplicate bit patterns and the tabulated data has been edited. There are many useful angles with two distinct bits, etc. Routine CIRCDAI2 illus- trates application for this example with difference angle scaled Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 B15 and minimum distinct bit count ideas for n = 100 which corre- sponds to a circle with radius about 1000. `Minimum' Angle The selection of angles above was limited to 360/n and not many angles were found. These angles correspond to regular polygons so that integral, or nearly integral, numbers of sides occur. Suppose, instead, that the search is over the sine and cosine bit patterns themselves with only the requirement that sine squared plus cosine squared be about one. The number of sides may not be integral but this proves not to be essential. Suppose one asks whether there are such angles whose sine and cosine taken together contain few distinct bits. Such angles permit coordinate rotations with a correspondingly small number of shift and add operations. For each scaling of the circular functions, e.g. B14, B15, etc., there are angles whose distinct set bit count is minimum. Generally, these angles change with scaling. Routine MINCANG searches for such angles where sine and cosine are scaled Bxx and produces the table below for 3 or fewer distinct bits. Note that each pair has complementary angles. The value for 89.55238 degrees contains just one distinct bit. The values for 7.18076 and and 82.81936 scaled at B15 are also interesting because they each contain only two distinct bits. The first is close to 7.2 degrees or exactly 1/50 of a circle. Routine CIRCRAM uses this angle to draw circles. Its utility is limited because fifty sides is about the right number for a circle with radius about 255. Smaller circles are using too many sides and larger circles, not enough. See DELTA ANGLE table above. MINIMUM ANGLE B14 ================================================== ang r sin cos number of bits (deg) sin cos total -------------------------------------------------- 89.10474 2 $3FFE $0100 -1 1 2* 82.80538 2 $3F7F $0804 -1 2 3 7.19485 254 $0804 $3F7F 2 -1 3 7.18781 254 $0802 $3F7F 2 -1 3 7.18428 254 $0801 $3F7F 2 -1 3 7.18076 255 $0800 $3F7F 1 -1 2* 3.63939 991 $0410 $3FDF 2 -1 3 3.61135 1007 $0408 $3FDF 2 -1 3 1.90275 3627 $0220 $3FF7 2 -1 3 ================================================== Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 MINIMUM ANGLE B15 ================================================== ang r sin cos number of bits (deg) sin cos total -------------------------------------------------- 89.55238 2 $7FFF $0100 0 1 1* 89.10474 2 $7FFC $0200 -2 1 3 88.65710 2 $7FF7 $0300 -1 2 3 82.81936 2 $7EFF $1000 -1 1 2* 14.47751 63 $2000 $7BEF 1 -2 3 7.18428 254 $1002 $7EFF 2 -1 3 7.18252 255 $1001 $7EFF 2 -1 3 7.18076 255 $1000 $7EFF 1 -1 2* 3.61135 1007 $0810 $7FBF 2 -1 3 3.59734 1015 $0808 $7FBF 2 -1 3 1.84677 3850 $0420 $7FEF 2 -1 3 ================================================== MINIMUM ANGLE B16 ================================================== ang r sin cos number of bits (deg) sin cos total -------------------------------------------------- 89.55238 2 $FFFE $0200 -1 1 2 82.81936 2 $FDFE $2000 -2 1 3 75.52263 2 $F7DF $4000 -2 1 3 14.47751 63 $4000 $F7DF 1 -2 3 7.18076 255 $2000 $FDFE 1 -2 3 3.59734 1015 $1010 $FF7F 2 -1 3 3.59033 1019 $1008 $FF7F 2 -1 3 1.81878 3970 $0820 $FFDF 2 -1 3 ================================================== Routine CIRCRAM uses parameters at B15 for angle 7.18076 and ELIPRAM uses parameters at B14 for the same angle to draw ellipses. Routines with `*M*.x' names use this minimum angle idea. Note that one way of thinking about ellipses is in terms of two circles, viz., major and minor axes circles. An ellipse can be generated by taking the x coordinate from the major circle and the y coordinate from the minor circle. That is; 2 2 cos (A) + sin (A) = 1 Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 or 2 2 2 2 a cos (A) b sin (A) --------- + --------- = 1 2 2 a b or, the usual equation of an ellipse, 2 2 x y -- + -- = 1 2 2 a b where a and b are major and minor semi-axes and x = a cos(A) and y = b sin(A). Routine CIRCRAS uses 89.10474 at B15 to draw points on a circular constellation. The angular step is not commensurate with 360 so that the pattern of points precesses slightly with each revolution. This routine is included just to show what some of the values in these tables mean. The algorithm for this value is unstable, i.e., if more and more points are drawn they no longer lie on a circle but spiral in or out depending upon effective computational modulus. Any of the `minimum' angle values - incommensurate with 360 - can be applied to generate constellations this way. There are similar elliptic constellations. For the elliptic routines there are drivers DEMOELIP and BNCHELIP to show and time the ellipse drawing routines distributed with this release. See DISTRIBUTION section below. `Minimum' Cosine Routine MINCOS tabulates rounded number of chord seg- ments, angle, and cosine and sine for angles whose cosine bit pattern contains only one distinct bit. MINIMUM COSINE B14 ========================================================== n ang r cos sin 2*cos nbit |nbit| (deg) 2cos 2cos ---------------------------------------------------------- 4 89.10471 2 $3FFE $0100 $0200 1 1 4 88.20921 2 $3FF8 $0200 $0400 1 1 4 86.41668 2 $3FE0 $0400 $0800 1 1 4 82.81924 2 $3F7F $0800 $1000 1 1 5 75.52249 2 $3DF8 $1000 $2000 1 1 6 60.00202 4 $376D $2000 $3FFF -1 1 Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 6 60.00000 4 $376D $2000 $4000 1 1 9 41.41227 8 $2A56 $3000 $5FFF -1 1 12 28.95864 16 $1EFD $3800 $6FFF -1 1 18 20.36916 32 $1647 $3C00 $77FF -1 1 25 14.36856 64 $0FE2 $3E00 $7BFF -1 1 35 10.15172 127 $0B48 $3F00 $7DFF -1 1 50 7.18065 255 $0800 $3F80 $7EFF -1 1 71 5.08569 508 $05AC $3FC0 $7F7F -1 1 100 3.60945 1008 $0407 $3FE0 $7FBF -1 1 140 2.57162 1986 $02DF $3FF0 $7FDF -1 1 195 1.84568 3855 $0210 $3FF8 $7FEF -1 1 ========================================================== MINIMUM COSINE B15 ========================================================== n ang r cos sin 2*cos nbit |nbit| (deg) 2cos 2cos ---------------------------------------------------------- 4 89.55238 2 $7FFF $0100 $0200 1 1 4 89.10471 2 $7FFC $0200 $0400 1 1 4 88.20921 2 $7FF0 $0400 $0800 1 1 4 86.41668 2 $7FC0 $0800 $1000 1 1 4 82.81924 2 $7EFF $1000 $2000 1 1 5 75.52249 2 $7BEF $2000 $4000 1 1 6 60.00101 4 $6EDA $4000 $7FFF -1 1 6 60.00000 4 $6EDA $4000 $8000 1 1 9 41.41095 8 $54AB $6000 $BFFF -1 1 12 28.95683 16 $3DF9 $7000 $DFFF -1 1 18 20.36665 32 $2C8C $7800 $EFFF -1 1 25 14.36504 64 $1FC2 $7C00 $F7FF -1 1 35 10.14676 128 $168D $7E00 $FBFF -1 1 50 7.17365 255 $0FFC $7F00 $FDFF -1 1 71 5.07582 510 $0B53 $7F80 $FEFF -1 1 100 3.59554 1016 $0807 $7FC0 $FF7F -1 1 141 2.55206 2016 $05B3 $7FE0 $FFBF -1 1 ========================================================== MINIMUM COSINE B16 ============================================================== n ang r cos sin 2*cos nbit |nbit| (deg) 2cos 2cos -------------------------------------------------------------- 4 89.55238 2 $FFFE $0200 $00000400 1 1 4 89.10471 2 $FFF8 $0400 $00000800 1 1 4 88.20921 2 $FFE0 $0800 $00001000 1 1 4 86.41668 2 $FF80 $1000 $00002000 1 1 4 82.81924 2 $FDFE $2000 $00004000 1 1 5 75.52249 2 $F7DF $4000 $00008000 1 1 6 60.00051 4 $DDB4 $8000 $0000FFFF -1 1 6 60.00000 4 $DDB4 $8000 $00010000 1 1 9 41.41028 8 $A955 $C000 $00017FFF -1 1 Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 12 28.95593 16 $7BF0 $E000 $0001BFFF -1 1 18 20.36539 32 $5917 $F000 $0001DFFF -1 1 25 14.36327 64 $3F81 $F800 $0001EFFF -1 1 35 10.14428 128 $2D17 $FC00 $0001F7FF -1 1 50 7.17015 255 $1FF4 $FE00 $0001FBFF -1 1 71 5.07088 511 $16A1 $FF00 $0001FDFF -1 1 100 3.58856 1020 $1006 $FF80 $0001FEFF -1 1 142 2.54222 2032 $0B5B $FFC0 $0001FF7F -1 1 ============================================================== Routine CIRCDAM2 implements the algorithm for scaling at B15 and two times cosine equal to $FEFF with 71 segments corresponding to a circle with radius about 510. Routine CIRCDAS uses the values for 75.52249 to execute a difference equation approach to constellation generation where computation of sine or cosine requires one shift and one add. This selection is arbitrary and most of the entries in the table above will work as well or better as long as the effective step angle is incommensurate to 360. All of those that this author has tried, for both rotation and difference angles, are unstable. Delta Hyperbolic Angles The ideas outlined above, excepting incremental angles, all apply to hyperbolae. A hyperbola cannot be said to have a period in the same sense that circles and ellipses can, so the idea of dividing infinity by n to get a step interval is not immediately useful. Given; 2 2 cosh (A) - sinh (A) = 1 then, 2 2 x y -- - -- = 1 . 2 2 a b Hyperbolic coordinates are unbounded. If display coordi- nates are arbitrarily restriced to 16 bit 2's complement numbers then there is a relation between hyperbolic function argument step and number of steps to drive a coordinate asymptote to display limit, i.e., r cosh x < 32768. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 Geometrically, the step angle is estimated from the rela- tion: dx r + 0.5 cosh(--) = ------- 2 r which is approximated by, 1 dx = 2 sqrt( - ) r much like the relations for circular functions. The sum angle formula for cosh gives: cosh (x+dx) = cosh(x) * cosh(dx) + sinh(x) * sinh(dx) which can be iterated to produce the exact number of rotations required to reach a display limit. An approximation is derived by observing that: cosh (x+dx) = cosh(x) (cosh(dx) + dx) for asymptotically large r and x. If r = 1 then cosh(dx) = 2 so that: r cosh(x) < 16384 is the iteration limit because the next step would produce over- flow. Therefore, the number of rotations can be estimated: 32768 ln( ----- ) r ndxe = ----------------- ln(cosh(dx) + dx) Routine DELHANG computes the exact and estimated num- ber of steps for radii from 1 to 1024 tabulated below. DELTA HYPERBOLIC ANGLES ============================ r dx ndx ndxe ---------------------------- 1 2.00000 6 5 2 1.41421 7 7 4 1.00000 10 9 8 0.70711 12 12 16 0.50000 16 15 32 0.35355 20 19 64 0.25000 25 25 128 0.17678 32 31 Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 256 0.12500 39 38 512 0.08839 48 47 1024 0.06250 56 55 ---------------------------- The number of strokes computed above is worst case and necessary only to a display coordinate space scaled to 16 bits. For example, a Hercules display adapter with resolution 720 x 348 requires, typically, 10-12 strokes per asymptote for a complete, i.e., two branches or four asymptotes, display centered hyperbola. A more realistic stroke count can be estimated by requir- ing that the hyperbolic coordinates be bounded by offset display window limits. That is, if the conic is centered at (0,0) rather than at display center then the maximum coordinate is display width. Therefore, if coordinate limit is computed: r cosh(x) <= abs(xc) + xd <= 32767 where xd is display width, then stroke count is reduced. The table below shows this relation, computed by DELHANGD, for a Hercules graphics adapter with xc = 360: DELTA HYPERBOLIC ANGLES DISPLAY COORDINATES ============================ r dx ndx ndxe ---------------------------- 1 2.00000 4 4 2 1.41421 5 5 4 1.00000 7 6 8 0.70711 8 8 16 0.50000 10 10 32 0.35355 12 12 64 0.25000 15 14 128 0.17678 16 16 256 0.12500 17 17 512 0.08839 16 16 1024 0.06250 6 11 ---------------------------- These estimates are for an unrotated hyperbola. The esti- mate for number of steps used in all hyperbolic routines is based on the two dimensional extension of this idea where an effective radial coordinate limit is computed from horizontal and vertical coordinate limits. There is no need to adjust dx for rounded ndx because closure is not an issue. This estimate is still as much as two times too high depending upon actual center and eccentricity. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 `Minimum' Hyperbolic Angles Routine MINHANG searches the bit patterns for hyperbolic sines and cosines scaled Bxx and lists those arguments for which the sum of distinct set bits in sinh and cosh is less than 5. Also listed is the radius for a corresponding equilateral hyper- bola. Routine HYPRRAM demonstrates application of the parameters for the case r equals 257. MINIMUM HYPERBOLIC ANGLE B14 ================================================== arg r sinh cosh number of bits sinh cosh total -------------------------------------------------- 0.69315 8 $3000 $5000 2 2 4 0.69310 8 $2FFF $4FFF -1 -2 3* 0.69305 8 $2FFE $4FFF -2 -2 4 0.24936 64 $1020 $4200 2 2 4 0.12516 255 $0808 $4080 2 2 4 0.12492 256 $0804 $4080 2 2 4 0.12480 257 $0802 $4080 2 2 4 0.12474 257 $0801 $4080 2 2 4 0.12468 257 $0800 $4080 1 2 3* 0.06295 1010 $0408 $4020 2 2 4 0.06270 1017 $0404 $4020 2 2 4 0.06258 1021 $0402 $4020 2 2 4 0.06252 1023 $0401 $4020 2 2 4 0.06246 1025 $0400 $4020 1 2 3* 0.03173 3972 $0208 $4008 2 2 4 ================================================== MINIMUM HYPERBOLIC ANGLE B15 ================================================== arg r sinh cosh number of bits sinh cosh total -------------------------------------------------- 0.69315 8 $6000 $A000 2 2 4 0.69312 8 $5FFF $9FFF -1 -2 3* 0.69310 8 $5FFE $9FFF -2 -2 4 0.24936 64 $2040 $8400 2 2 4 0.12492 256 $1008 $8100 2 2 4 0.12480 257 $1004 $8100 2 2 4 0.06270 1017 $0808 $8040 2 2 4 0.06258 1021 $0804 $8040 2 2 4 0.06252 1023 $0802 $8040 2 2 4 0.06249 1024 $0801 $8040 2 2 4 0.06246 1025 $0800 $8040 1 2 3* 0.03173 3972 $0410 $8010 2 2 4 ================================================== Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 MINIMUM HYPERBOLIC ANGLE B16 ====================================================== arg r sinh cosh number of bits sinh cosh total ------------------------------------------------------ 0.69315 8 $C000 $00014000 2 2 4 0.69313 8 $BFFF $00013FFF -1 -2 3* 0.69312 8 $BFFE $00013FFF -2 -2 4 0.24936 64 $4080 $00010800 2 2 4 0.12492 256 $2010 $00010200 2 2 4 0.06258 1021 $1008 $00010080 2 2 4 0.06252 1023 $1004 $00010080 2 2 4 0.06249 1024 $1002 $00010080 2 2 4 0.06247 1025 $1001 $00010080 2 2 4 0.06246 1025 $1000 $00010080 1 2 3* 0.03149 4034 $0810 $00010020 2 2 4 ====================================================== *minima `Minimum' Hyperbolic Cosines The difference equations for sinh or cosh generation are just like those for the circular functions with a coefficient c = 2 * cosh(da). Routine MINCOSH searches the interval (1.0 < cosh <= 1.25) and lists the arguments whose hyperbolic cosines scaled Bxx contain fewer than three distinct bits. Routine HYPRDAM implements the B16, r = 512 parameters. MINIMUM HYPERBOLIC COSINE B14 ====================================================== arg r sinh cosh 2*cosh nbit |nbit| 2cosh 2cosh ------------------------------------------------------ 0.69311 8 $2FFF $5000 $9FFF -2 2 0.49493 16 $20FC $4800 $9000 2 2 0.35174 32 $16FA $4400 $8800 2 2 0.24935 64 $1020 $4200 $8400 2 2 0.17655 128 $0B5C $4100 $8200 2 2 0.12492 256 $0804 $4080 $8100 2 2 0.08836 512 $05AA $4040 $8080 2 2 0.06249 1024 $0400 $4020 $8040 2 2 0.04419 2048 $02D4 $4010 $8020 2 2 ====================================================== Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 MINIMUM HYPERBOLIC COSINE B15 ========================================================== arg r sinh cosh 2*cosh nbit |nbit| 2cosh 2cosh ---------------------------------------------------------- 0.69313 8 $5FFF $A000 $00013FFF -2 2 0.49493 16 $41F8 $9000 $00012000 2 2 0.35174 32 $2DF5 $8800 $00011000 2 2 0.24935 64 $2040 $8400 $00010800 2 2 0.17655 128 $16B7 $8200 $00010400 2 2 0.12492 256 $1008 $8100 $00010200 2 2 0.08836 512 $0B53 $8080 $00010100 2 2 0.06249 1024 $0801 $8040 $00010080 2 2 0.04419 2048 $05A9 $8020 $00010040 2 2 ========================================================== MINIMUM HYPERBOLIC COSINE B16 ============================================================== arg r sinh cosh 2*cosh nbit |nbit| 2cosh 2cosh -------------------------------------------------------------- 0.69314 8 $BFFF $00014000 $00027FFF -2 2 0.49493 16 $83F0 $00012000 $00024000 2 2 0.35174 32 $5BEA $00011000 $00022000 2 2 0.24935 64 $4080 $00010800 $00021000 2 2 0.17655 128 $2D6E $00010400 $00020800 2 2 0.12492 256 $2010 $00010200 $00020400 2 2 0.08836 512 $16A6 $00010100 $00020200 2 2 0.06249 1024 $1002 $00010080 $00020100 2 2 0.04419 2048 $0B51 $00010040 $00020080 2 2 ============================================================== Parabola Routines to generate parabolas are included just to cover all conics. They do not demonstrate incremental or minimum angle ideas. A parabola can be defined, among many ways, as a curve equidistant between a point and a line. That distance is usually denoted by p - the parameter. The point is called the focus and the line is called the directrix. The point on the curve midway between the focus and directrix is the vertex. If the focus is at (-p/2,0) and the vertex is at (0,0) then the equation of the parabola is y*y = -2*x*p. The figure is open to the left rather than open upwards and this orientation is defined, here, to be at zero rotation angle. The curve can be generated by stepping the polar angle from 0 at the vertex to 180 at the asymptotic limit and drawing both asymptotes by reflection. The discussion regarding hyperbolic Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 limits also applies here and limits the polar radius to abs(xf)+xd where (xf,yf) is the focus and xd is display width in pixels. The step angle can be computed from the local radius of curvature by evaluating da = 2.0*sqrt(1.0/R) where R is radius of curvature estimated by: 2 r 3/2 R = p ( 1 + -- ) 2 p where r is polar radius. Routine PARASA generates parabolas with specified focal coordinates, parameters, and rotation angles. Routine DEMOPARA is a driver. Routines ParaRA and ParaRAF present floating and fixed, respectively, rotation angle algorithms and routines ParaMB and ParaMBR present offset error algortihms. Offset Error Methods Let a circle be defined: 2 2 2 x + y = r and pick a point, say, (ix,iy). If that point is picked to satisfy the equation then the error at that point can be taken to be zero. If a positive step of one unit in x is taken then the equation becomes: 2 2 2 (ix+1) + iy = r + ie, or 2 2 2 (ix + 2 ix + 1) + iy = r + ie, where ie is the offset error for a unit step in x: +x: ie = 2 ix + 1. Similarly, for a negative step in x or steps in y: -x: ie = -2 ix + 1, +y: ie = 2 iy + 1, -y: ie = -2 iy + 1. The idea of the offset error method is to pick a sequence of x and y steps to attempt to reduce the cumulative offset error. For example, take the starting point for a circle at (0,r) and step positive in x and negative in y until x = y. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 This segment is one octant of a circle. Copying this segment with suitable reflections produces a complete circle. Obviously, the starting point can be anywhere on the circle with self consistent arrangements for stepping and copying to produce a complete figure. The algorithm can be outlined: ix = 0 iy = r ie = 0 while ix < iy do begin ie = ie + 2*ix + 1 ix = ix + 1 if ie > 0 then begin ie = ie - 2*iy + 1 iy = iy - 1 end pixel(ix,iy) ... end where the copying is denoted by the ellipsis. With suitable consideration of quadrant and octant points this is Bresenham's circle algorithm. In a certain sense, the offset error is the error in area of the figure. If another property like radius or curvature is used then another kind of offset error algorithm results. In two dimensions area offset error reduction produces esthetically satisfactory results. Offset error methods are not restricted to conics or to two dimensions. There are many modifications to this algorithm. For example, if a step in y is taken only if the absolute error is reduced then the algorithm appears: ix = 0 iy = r ie = 0 while ix < iy do begin ie = ie + 2*ix + 1 ix = ix + 1 iey = ie - 2*iy + 1 if abs(ie) > abs(iey) then begin ie = iey iy = iy - 1 end pixel(ix,iy) ... end This algorithm is trying to minimize the offset error. If aspect ratio is considered then the equation becomes: 2 2 2 2 x + y c = r Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 where c is a vertical scale factor. This aspect ratio parameter is computed from the effective display dimensions and resolutions. A standard broadcast video monitor has an aspect of 4 units horizontally and 3 units vertically. For a Hercules graphics adaptor the resolutions are 720 by 348. Ideally, such an adaptor driving such a monitor would produce vertical pixel distances of 3/347 and horizontal distances of 4/719. The parameter c would be: c = 3/347 / 4/719 = 1.55... For computer displays an effective aspect of 3 to 2 is more common with: c = 2/347 / 3/719 = 1.38..., The wrong value makes circles into ellipses. The particu- lar monitor on which this writeup is prepared has an actual c = 1.43. Let w denote width in pixels and h denote height in pixels. Then: c = 2/(h-1) / 3/(w-1), and ignore the -1 terms: c = 2/h / 3/w. The circle equation may be written - by clearing the fractions: 2 2 2 (3h ix) + (2w iy) = (3h r) and the relevant unit offset errors are: 2 +x: ie = (2 ix + 1) 9h , 2 -y: ie = -(2 iy - 1) 4w , and the algorithm becomes: ix = 0 iy = r px2 = sqr(3h) py2 = sqr(2w) ie = 0 while ix*px2 < iy*py2 do begin ie = ie + (2*ix + 1) px2 ix = ix + 1 iey = ie - (2*iy - 1) py2 if abs(ie) > abs(iey) then begin ie = iey Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 iy = iy - 1 end pixel(ix,iy) ... end while iy > 0 do begin ie = ie - (2*iy - 1) dey2 iy = iy - 1 iex = ie + (2*ix + 1) dex2 if abs(ie) > abs(iex) then begin ie = iex ix = ix + 1 end pixel(ix,iy) ... end where a couple things have changed. First, octant symmetry is lost and a quadrant is drawn in two segments with stepping in x dominating the first loop and stepping in y dominating the second. The changeover in stepping in x to stepping in y is at the point where the absolute slope is about one. The ellipsis denote quadrant copying operations. This algorithm is called modified Bresenham. Second, the loop control condition for the second segment, viz. iy > 0, is correct but misleading. Area is invariant - not coordinate. The loop control should be iy*py2, i.e. the area weighted coordinate. Note the loop control condition for the first segment. If the figure is drawn with the condition ix < iy then something looking more like sixteen sided polygons appears. Let the equation for an ellipse be written: 2 2 u + v -- -- = 1 2 2 a b and suppose a coordinate rotation with equations: u = x cos(A) + y c sin(A) v = -x sin(A) + y c cos(A) where A is rotation angle an c is the vertical scale factor. Substituting, expanding, and clearing fractions produces the equation: 2 2 2 2 2 2 x (3h) [b cos (A) + a sin (A)] + 2 2 2 x y (3h) (2w) (b - a ) sin(A) cos(A) + Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 2 2 2 2 2 2 2 2 2 y (2w) [b sin (A) + a cos (A)] = a b (3h) Offsetting x by one produces an error offset: 2 2 2 2 2 +x: ie = (2x+1) (3h) [b cos (A) + a sin (A)] + 2 2 2 y (3h) (2w) (b - a ) sin(A) cos(A) with similar expressions for negative offset in x and offsets in y. These expressions are only slightly more complicated than those for a circle. Note that the coefficients are constants: +x: ie = (2x+1) pxx + 2 y pxy -x: ie = -(2x-1) pxx - 2 y pxy +y: ie = (2y+1) pyy + 2 x pxy -y: ie = -(2y-1) pyy - 2 x pxy If the point on a ellipse where plus y is an extremum is taken as a starting point then an ellipse may be produced by four segments and semi-elliptic reflection. The control outline of this algorithm is: initialize ix = ix0 iy = iy0 ie = 0 step clockwise while ... do begin ix = ix + 1 if abs(ie) > abs(iey) then begin iy = iy - 1 end end while ... do begin iy = iy - 1 if abs(ie) > abs(iex) then begin ix = ix + 1 end end reinitialize ix = ix0 iy = iy0 ie = 0 step counter-clockwise while ... do begin ix = ix - 1 if abs(ie) > abs(iey) then begin iy = iy - 1 end end Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 while ... do begin iy = iy - 1 if abs(ie) > abs(iex) then begin ix = ix - 1 end end where the operations to compute initial coordinates, offset errors and increments, turning points, and plot are omitted. This algorithm is modified Bresenham with rotation and routines named `*MBR*.x' apply it, e.g. ElipMBR in Ellipse.MOD. Obvious- ly, there is a structurally identical algorithm which starts at an extremum in x with the roles of ix and iy interchanged. Similar logic is applied to produce parabolas, cf. PARAMBR.PAS. Rotated hyperbolae are generated by using a vertex as starting point rather than an extremum point which may be at infinity. This requires up to three rather that two segments per asymptote, cf. HYPRMBR.PAS or procedure HyprMBR in Hyprbola.MOD. Note that there is nothing that requires offset error to be computed fixed point. Also, for TURBO, a floating point version of ELIPMBR is somewhat faster - if a coprocessor is present; --------------------------------------------------- TURBO Pascal TopSpeed Modula-2 Routine Mode Timing (sec) Timing (sec) w/ w/o w/ w/o --------------------------------------------------- ELIPMBR Fix 0.334 0.304 0.150 0.199 ELIPMBR Float 0.249 0.547 0.262 2.352 --------------------------------------------------- It would diminish the scaling problem for small parameter figures to compute it floating point. For example, an ellipse with parameters a = 10 and b = 2 is unsatisfactory at most rotation angles. Roughly speaking, this problem arises because the projection of small parameter radii upon the lower resolution vertical axis is less than one vertical pixel unit which is typically about one and one half horizontal pixels. Parameter `s' in parabolic and hyperbolic `*MBR.x' routines tries to extend the useful parameter range. Radius Scaling Fixed point arithmetic requires that operands be scaled. Scaling depends on word length and desired or required operand range. For example, if radii 0-511 are desired then scaling at B6 is possible for coordinate variables encoded using 16 bit integers or B22 using 32 bit integers. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 Radius Offset The inscribed polygon is required to have a chord-tangent distance not greater than 0.5. But the ideal and chord circle must both be plotted on a discrete grid. If the ideal discrete circle is taken to be that produced by Bresenham circle then it is possible to consider offsetting the chord circle radius to get closer agreement. Routine DFCTCIRC displays a stroke circle, clears pixels on the corresponding Bresenham circle, and counts the remaining pixels. This procedure is performed for radius offsets from 0.25 to 0.75 in steps of about 0.031. For circles with radii 100 and 200 using routine CIRCRA the results are: DEFECT COUNT VRS. OFFSET ========================================= r = 100 r = 200 offset offset number of number of fixed float pixels pixels ----------------------------------------- 16 0.25 124 472 18 0.28 124 460 20 0.31 104 406 22 0.34 104* 406 24 0.38 104* 406 26 0.41 104 406 28 0.44 148 406 30 0.47 148 348 32 0.50 148 268 34 0.53 178 254* 36 0.56 178 312 38 0.59 178 312 40 0.63 200 274 42 0.66 200 368 44 0.69 212 378 46 0.72 212 424 48 0.75 212 424 ========================================= *minima This table suggests an offset of 0.3-0.5 for floating point generators. In the algorithms using floating point to generate coordinates an offset of 0.5 will generally be used, e.g. CIRCRA. If rotation angle arithmetic is done fixed point then there is interaction between scaling, rounding, and offset, and resulting closure, overflow, and defect count. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 Step variables may be truncated or rounded and display variables may be truncated or rounded; ROUNDING STRATEGIES ============================ case step display variables variables ---------------------------- 1 round round 2 round truncate 3 truncate round 4 truncate truncate ============================ Truncation of display variables, i.e. cases 2 and 4, produces large defect counts for either float or fixed arithmetic. For case 1, with both step and display rounding, the relation between offset scaled B6 and defect counts for radii 100 and 200 are tabulated below. Routine DFCTCIRC using CIRCRAF4 with step and display variable rounding enabled produces these counts. DEFECT COUNT VRS. OFFSET ========================================= r = 100 r = 200 offset offset number of number of fixed float pixels pixels ----------------------------------------- 16 0.25 124 460 18 0.28 124 460 20 0.31 104 406 22 0.34 104* 406 24 0.38 104* 406 26 0.41 104 406 28 0.44 148 406 30 0.47 148 348 32 0.50 148 268* 34 0.53 148 268* 36 0.56 148 270 38 0.59 148 270 40 0.63 170 326 42 0.66 170 326 44 0.69 182 382 46 0.72 182 382 48 0.75 182 382 ========================================= *minima Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 where the minimum for r = 100 is at 22-24 or about 0.34-0.38 scaled B6, and for r = 200 it is at 32-34 or about 0.50-0.53. This means that the radius offsets are essentially the same for floating point and fixed point arithmetic - but that there is a radius dependence. Radius Constraint Using fixed point arithmetic there is a problem with closure and overflow. That is, a complete circle produced by rotation of coordinates may not close, or worse, arithmetic overflow due to roundoff accuulation occurs. Failure to close produces leaks on flood fill while overflow both leaks and pro- duces spurious graphics lines. For each algorithm and each implementation it is necessary to test for closure and overflow. Routine LEAKCRAF illustrates this test specifically instrumented for a fixed point rotation angle algorithm. The test is run for the same range of offsets as that used in DFCTCIRC and radii from 500 through 511 are tested because it is only at the limits of the scaling range that this problem arises. LEAK AND OVERFLOW =============================== offset leaks overflow FIX FLOAT ------------------------------- 16 0.25 0 0 18 0.28 0 1 20 0.31 0 1 22 0.34 0 1 24 0.38 0 1 26 0.41 0 1 28 0.44 0 1 30 0.47 1 1 32 0.50 0 0 34 0.53 0 0 36 0.56 1 0 38 0.59 1 1 40 0.63 1 1 42 0.66 1 1 44 0.69 1 1 46 0.72 0 1 48 0.75 0 1 =============================== The earlier defect testing suggested fixed point offset of 32 at B6 or 0.5 floating point. There are no leaks or overflow Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 for this selection for this algorithm and implementation. The tables below show a radius parameter constraint when the above testing produced leak or overflow. For example, routine ELIPDAI which uses the difference angle algorithm for increment angle at B15 and n = 50 is restricted to major and minor radii less than 508. The leak problem can be overcome by forcing closure on starting point - with or without symmetry - but the overflow problem can be best overcome by restricting radii. The rounding strategies noted above interact with closure and overflow. For each routine there is a possibility of round- ing or truncating coordinate or display variables. Display vari- able truncation aggravates the bumpy appearance of stroked fig- ures at low resolution or small radii. Truncation rather than rounding of coordinate variables produces slightly asymmetric figures but can solve the overflow problem. This is an interest- ing problem pending availability of shift arithmetic right imple- mentations because the rounding would be asymmetric. The discussion above about restricting radii to something less than 512 rather than 512 is for a worst case where the algorithm generates chord coordinates for a full circle. If semi-circular, quadrant, or octant symmetry is exploited, then overflow does not occur. Further, there are several cases where closure fails with- out overflow for specified radius range or where overflow occurs on final stroke. For those cases a final closing stroke can be drawn to close the figure rather than computing a final stroke. USAGE ===== TURBO drivers have an include header with all but one of the routines commented out. For example DEMOCIRC begins: ... {-$I circmb.pas } ... {$I circra.pas } ... {-$I circras.pas } ... where CIRCRA is selected. Compile and execute the driver with the appropriate routine include enabled. Each routine for a particular kind of figure has the same procedure name. For example, circle routines are declared with the name StrokeCircle independent of whether they are stroke. TopSpeed drivers are organized with all the procedures in Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 one module and each procedure is distinctly named. `Make' the driver and execute it with a command line argument specifying the desired routine. Remember to force initialization of the appropriate graphics driver - all routines default to InitHerc. For example: bnchcirc circmb would perform and report timing estimates for routine CircMB. Note that TopSpeed requires that decimal numbers be entered with decimal points and fractions. That is, 30 # 30.0. PARAMETERS ========== POLYGON ROUTINES ===================================================== Routine Mode Parameters Constraint n sin cos 2*cos ----------------------------------------------------- POLYRA FLT POLYRAF FIX r<512 POLYDA FLT POLYDAF FIX r<511 ----------------------------------------------------- CIRCLE ROUTINES ===================================================== Routine Mode Parameters Constraint n sin cos 2*cos ----------------------------------------------------- CIRCDA2 FLT CIRCDAF FIX r<510 CIRCDAI FIX B16 71 $169B $FF00 $1FDFF r<511 CIRCDAI2 FIX B15 100 $080F $7FBF $FF7F r<511 CIRCDAM2 FIX B15 71 $0B53 $7F80 $FEFF r<511 CIRCDAS FIX B15 $7BEF $2000 $4000 r<512 CIRCRA FLT CIRCRAF4 FIX r<512 CIRCRAI2 FIX B15 100 $080A $7FBF r<511 CIRCRAIR FIX B14 50 $0805 $3F7F r<512 CIRCRAM FIX B15 50 $1000 $7EFF r<510 CIRCRAS FIX B15 $7EFF $1000 r<504 ----------------------------------------------------- Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 ELLIPSE ROUTINES ===================================================== Routine Mode Parameters Constraint n sin cos 2*cos ----------------------------------------------------- ELIPDA2 FLT ELIPDAF FIX r<510 ELIPDAI FIX B15 100 $080F $7F8F $FF7F r<508 ELIPDAM FIX B16 50 $1FF4 $FE00 $1FBFF r<512 ELIPRA FLT ELIPRA2 FLT ELIPRAF FIX r<511 ELIPRAF2 FIX r<511 ELIPRAI2 FIX B14 50 $0805 $3F7F r<510 ELIPRAM FIX B14 50 $0800 $3F7F r<511 ----------------------------------------------------- PARABOLA ROUTINES ===================================================== Routine Mode Parameters Constraint n sin cos 2*cos ----------------------------------------------------- PARARA FLT PARARAF FIX r<511 ----------------------------------------------------- HYPERBOLA ROUTINES ===================================================== Routine Mode Parameters Constraint n sinh cosh 2*cosh ----------------------------------------------------- HYPRDA FLT HYPRDAF FIX HYPRDAM FIX B16 $16A6 $10100 $20200 HYPRDAMG FIX B16 $2010 $10200 $20400 (SAR) " " " " " HYPRRA FLT HYPRRAF FIX HYPRRAM FIX B14 $0800 $4080 ----------------------------------------------------- PERFORMANCE =========== All timings for both TURBO Pascal 4.0 (c) and TopSpeed Modula-2 1.14 (c) are on an 8MHz 8086/8087 driving a Hercules graphics adapter with the usual 720x348 resolution. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 Times are average time per figure for a more-or-less random parameter selection. For example, TopSpeed's GRAPH Ellipse routine exercised by BnchElip produced a set of trials; 0.150, 0.153, 0.144, 0.137, 0.145, 0.138, 0.147, 0.152, 0.143, 0.138, 0.147, and 0.139. This sequence has a mean of 0.144 so that plus or minus ten percent of the timing numbers is not significant. Comparison numbers are shown for TURBO's BGI and TopSpeed's GRAPH circle and ellipse routines. BGI circles correct aspect - GRAPH circles don't. The differences in timing between Versions 1.0 and 1.1 arises mostly because the centers are randomized for 1.1 while they are screen centered for 1.0. Remember to disable compiler check code generation, i.e. no stack checking, I/O errors, etc., - makes about ten percent difference. POLYGON ROUTINES =================================================== TURBO Pascal TopSpeed Modula-2 Routine Mode Timing (sec) Timing (sec) w/ w/o w/ w/o --------------------------------------------------- POLYRA FLT 0.076 0.177 0.119 0.220 POLYRAF FIX 0.089 0.118 0.115 0.176 POLYDA FLT 0.074 0.158 0.118 0.258 POLYDAF FIX 0.082 0.135 0.120 0.176 --------------------------------------------------- CIRCLE ROUTINES =================================================== TURBO Pascal TopSpeed Modula-2 Routine Mode Timing (sec) Timing (sec) w/ w/o w/ w/o --------------------------------------------------- Circle 0.381 0.359 0.073 0.053 CIRCMB 0.249 0.280 0.098 0.080 CIRCDA2 FLT 0.086 0.153 0.086 0.210 CIRCDAF FIX 0.119 0.134 0.110 0.127 CIRCDAI FIX 0.256 0.253 0.154 0.126 CIRCDAI2 FIX 0.234 0.214 0.143 0.118 CIRCDAM2 FIX 0.177 0.163 0.121 0.097 CIRCDAS FIX 0.433 0.394 0.074 0.069 CIRCRA FLT 0.099 0.280 0.097 0.418 CIRCRAF4 FIX 0.086 0.108 0.083 0.108 CIRCRAI2 FIX 0.164 0.171 0.152 0.120 CIRCRAIR FIX 0.236 0.224 0.112 0.101 CIRCRAM FIX 0.263 0.270 0.134 0.116 CIRCRAS FIX 0.671 0.650 0.111 0.097 --------------------------------------------------- Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 ELLIPSE ROUTINES =================================================== TURBO Pascal TopSpeed Modula-2 Routine Mode Timing (sec) Timing (sec) w/ w/o w/ w/o --------------------------------------------------- Ellipse 0.316 0.334 0.150 0.101 ELIPMB 0.218 0.293 0.146 0.126 ELIPMBR 0.281 0.304 0.163 0.199 ELIPDA2 FLT 0.110 0.265 0.143 0.459 ELIPDAF FIX 0.160 0.175 0.187 0.215 ELIPDAI FIX 0.505 0.466 0.310 0.252 ELIPDAM FIX 0.277 0.253 0.223 0.186 ELIPRA FLT 0.140 0.511 0.182 0.845 ELIPRA2 FLT 0.105 0.282 0.157 0.527 ELIPRAF FIX 0.175 0.186 0.207 0.225 ELIPRAF2 FIX 0.120 0.154 0.182 0.203 ELIPRAI2 FIX 0.135 0.150 0.191 0.165 ELIPRAM FIX 0.308 0.302 0.210 0.190 --------------------------------------------------- PARABOLA ROUTINES =================================================== TURBO Pascal TopSpeed Modula-2 Routine Mode Timing (sec) Timing (sec) w/ w/o w/ w/o --------------------------------------------------- PARAMB 0.337 0.377 0.163 0.146 PARAMBR 0.410 0.642 0.272 0.299 PARARA FLT 0.136 0.530 0.216 1.052 PARARAF FIX 0.229 0.259 0.253 0.274 --------------------------------------------------- HYPERBOLA ROUTINES =================================================== TURBO Pascal TopSpeed Modula-2 Routine Mode Timing (sec) Timing (sec) w/ w/o w/ w/o --------------------------------------------------- HYPRMB 0.332 0.400 0.172 0.158 HYPRMBR 0.376 0.466 0.259 0.272 HYPRDA FLT 0.087 0.276 0.197 0.504 HYPRDAF FIX 0.097 0.185 0.226 0.248 HYPRDAM FIX 0.156 0.227 0.248 0.250 HYPRDAMG FIX 0.148 0.198 0.220 0.229 (SAR) 0.109 0.151 0.207 HYPRRA FLT 0.087 0.338 0.218 0.585 HYPRRAF FIX 0.107 0.202 0.227 0.247 HYPRRAM FIX 0.209 0.272 0.267 0.228 --------------------------------------------------- Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 DISTRIBUTION ============ Below are listed the files and functions with some com- ments and performance parameters that are distributed with this release. FILENAME DESCRIPTION -------- ----------- TURBO TopSpeed Pascal Modula-2 CNC11TP CNC11TM CONIC.DOC CONIC.DOC This write-up. DELCANG n/a Delta circular angle. INCCANG n/a Incremental angle. INCCOS n/a Incremental cosine. MINCANG n/a Minimum angle. MINCOS n/a Minimum cosine. DELHANG n/a Delta hyperbolic angle. DELHANGD n/a Delta hyperbolic angle - display. MINHANG n/a Minimum hyperbolic angle. MINCOSH n/a Minimum hyperbolic cosine. CONIC MathFun Utility routines. n/a GraphFun Utility routine. n/a Timing Utility routine. Polygon Routines DEMOPOLY DemoPoly Demonstrate polygon routines. BNCHPOLY BnchPoly Timing. POLYRA PolyRA POLYRAF PolyRAF POLYDA PolyDA POLYDAF PolyDAF Circle Routines DEMOCIRC DemoCirc Demonstrate circle routines. BNCHCIRC BnchCirc Timing. DFCTCIRC DfctCirc Geometric comparison. LEAKCRAF LeakCRAF Closure and overflow. CIRCMB CircMB CIRCDA2 CircDA2 CIRCDAF CircDAF CIRCDAI CircDAI CIRCDAI2 CircDAI2 Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 CIRCDAM2 CircDAM2 CIRCDAS CircDAS CIRCRA CircRA CIRCRAF4 CircRAF4 CIRCRAI2 CircRAI2 CIRCRAIR CircRAIR CIRCRAM CircRAM CIRCRAS CircRAS Ellipse Routines DEMOELIP DemoElip Demonstrate ellipse routines. BNCHELIP BnchElip Timing. ELIPMB ElipMB ELIPMBR ElipMBR ELIPDA2 ElipDA2 ELIPDAF ElipDAF ELIPDAI ElipDAI ELIPDAM ElipDAM ELIPRA ElipRA ELIPRA2 ElipRA2 ELIPRAF ElipRAF ELIPRAF2 ElipRAF2 ELIPRAI2 ElipRAI2 ELIPRAM ElipRAM Parabola Routines DEMOPARA DemoPara Demonstrate parabola routines. BNCHPARA BnchPara Timing. PARAMB ParaMB PARAMBR ParaMBR PARARA ParaRA PARARAF ParaRAF Hyperbola Routines DEMOHYPR DemoHypr Demonstrate hyperbola routines. BNCHHYPR BnchHypr Timing. HYPRMB HyprMB HYPRMBR HyprMBR HYPRDA HyprDA HYPRDAF HyprDAF HYPRDAM HyprDAM HYPRDAMG HyprDAMG (SAR) n/a HYPRRA HyprRA HYPRRAF HyprRAF HYPRRAM HyprRAM NOTES Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 ===== 1. For most centers only one hyperbolic branch or asymptote, at most, is required. Likewise for parabolae. 2. TURBO Pascal and Topspeed Modula-2 lack a shift arithmetic right capability. Many of the routines use div to get arith- metic right shifts. This is typically one third again as slow. For example, in the HYPERBOLA table the entry (SAR) is a special test of routine HYPRDAMG where the rotation angle is restriced to first quadrant and shr is substituted for div to show sar timing effect. Note that div and sar would not be equivalent and each application would have to be tested. For example, -7 div 2 = -3, while -7 sar 1 = -4, etc. Likewise, neither TURBO nor TopSpeed has a facility to direct- ly manipulate a floating point exponent. An IncExp and DecExp function would enable rapid power of two operations. 3. It is surprising that TURBO's BGI and TopSpeed's GRAPH circle and ellipse routines are so slow - especially when it is noted that all the routines presented here are higher order language. 4. Note that the timing comparisons are for a random selection of parameters. If the tests are restriced to the parameter ranges for which the respective routines are ideal, then the minimum angle routines are 5-10 times faster. For example, if radii for CIRCDAM2 are set at 510 then TURBO Circle requires 0.903 sec while the fixed point `minimum' difference angle routine requires about the same 0.192 sec per figure. 5. A high performance implementation of minimum angle ideas would justify a polyalgorithm where the parameter set is matched to specified radius for several radius ranges. This algorithm would also economically compensate the radius dependent offset noted in defect testing for smaller radii. A high performance implementation of rotation angle or difference angle ideas would allow overlap of numeric coprocessor operation with graphics operations. That is, the coordinates for the next stroke could be generated while the pixels on the present stroke are toggled. 6. Another kind of fixed point routine is possible where the scaling shift necessary to normalize radius is computed and used to scale for coordinate computation and descale for display variable computation. For the anticipated range of display resolutions and radii the B6 scaling for circles and ellipses and B0 scaling for hyperbolae is probably adequate. 7. Two natural logarithms is a lot of arithmetic to get a useful stroke count for hyperbolae. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 8. The defect counts indicate disagreement between TURBO's Circle and stroked circles of about one pixel in seven or eight when using TURBO's GetAspectRatio parameters. 9. Overhead in calls to Line, for example, moderates use of higher than quadrant symmetry. That is, routine performance is not greatly improved by using octant symmetry over quadrant symmetry because stroke computation is reduced but argument fetch, etc. is not. 10. Neither TURBO nor TopSpeed appear to be using the numeric coprocessor for integer and long integer operations. Note the anomalous TopSpeed results where routines run faster without numeric coprocessor, i.e. Circle and Ellipse. TopSpeed does not appear to have an option to disable coprocessor usage and has a double or nothing floating operand restriction. 11. If increment and minimum angle tableaus for scalings above about B20 are required then floating point arithmetic in INCxx and MINxx routines should be done double rather than single. 12. These same incremental and minimum angle ideas are applicable to fast Fourier and fast Hartley transforms as well as to function generation, particularly, variable precision function generation. 13. Using extremum point geometry for parabolas means that parameters greater than about 400 don't draw anything because the extremum point is outside the coordinate limits as given. The routines can be reformulated to use vertex point geometry, cf. HyprMBR, which requires an additional loop per asymptote. Coordinate limits may be imposed with a Cartesian norm, e.g., abs(ix) + abs(ix) < irx, which relieves the problem slightly with some penalty in running time. 14. TopSpeed's GRAPH restriction of Circle, etc., centers to cardinal values is unnecessarily restrictive. Further, the restriction of Line coordinates to cardinal values is unaccept-able, particularly because the compiler accepts integer arguments. What happens is that GRAPH's Line routine treats the negative values as large cardinal values and tests the whole range for display. Such arguments make the display appear to be inoperable for 10's of seconds. Use the Line routine from GraphFun.MOD. This routine could be further improved for closed figure application by noting that both end points do not have to be displayed. Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730 COMMENTS ======== 1. It would be useful to have a simple graphics standard for display scaling and aspect. This might not require an inter- national conference. Just a set of conventions so that a procedure call, e.g. Circle(200,200,100), produced a round circle with the same display relative size and positioning independent of display adapter. 2. Neither TURBO's nor TopSpeed's graphics libraries can be recommended for speed. It would help to combine line and pixel routines and drop color as an argument for each pixel. 3. As a rule, TURBO 4.0 does better than TopSpeed 1.14 with floating point while TopSpeed does better than TURBO with fixed point. OBSERVATIONS ============ 1. To this author, the displays produced when parabolae with random parameters are displayed vaguely but often suggest a network or lattice with nodes and correlated asymptotes. Ran- dom ellipses and hyperbolae do dot produce the same effect. I wonder if Lowell was observing something similar. NOTICES ======= CP/M (c) Digital Research Corporation MS-DOS (c) Microsoft Corporation TURBO Pascal (c) Borland International TopSpeed Modula-2 (c) Jensen Partners, International Hercules (tm) Hercules, Inc. iAPX8086/8087 (tm) Intel Corporation Copyright (C) 1989 Adam Fritz, 133 Main St., Afton, N.Y. 13730