| It sounds as if you have "truncation" errors rather than round-off errors.
x = a + r*cos(T)
y = b + r*sin(T)
With these equations, you can get minor errors from the sin & cos routines,
the multiplication, and the addition. The relative error of 1.5x10^-7 seems
reasonable with single-precision arithmetic. For an angle of 45 degrees, you
ARE dealing with numbers that aren't exactly representable by the hardware
(namely, some square roots of 2).
The fact that GKS uses 2 different pixels for the points (0.4999997,0.5) and
(0.5,0.5) shouldn't be puzzling -- there's always SOME boundary value that
separates two pixels. It seems that GKS is truncating (rather than rounding)
its inputs to three significant digits -- this is reasonable.
Realizing this, the solution is to increase a & b (in the above equations) by
some small amount, so that "near" values (such as 0.4999997) will get truncated
to "nice" values (0.5). Assuming that a & b are fairly small (|a| < 10), just
add 10^-5 to each.
- Gilbert
|
| C Here is a short program that shows how to draw a near circle
C using only multiplications and additions, no square roots or
C trigonometric functions.
C
C Reference:
C ---------
C Beeler, Gosper, and Schroeppel. HAKMEM, Artificial Intelligence
C Memo 239, Massachusetts Institute of Technology,
C A.I. Laboratory. Cambridge: 1972. Page 73, item 149.
C
C The algorithm is due to Marvin Minsky.
C The following program is for a VT125 or VT240.
C Change initial value of X or EPS for different size or
C resolution circles.
C
PROGRAM CIRCLE
CALL CLEAR
EPS=0.01
X=50.0
Y=0.0
DO I=1,10000
X=X-EPS*Y
Y=Y+EPS*X ! Yes, that's not a typo. Use "new x" here not "old x".
CALL PLOT(X,Y)
END DO
TYPE 100, 27
100 FORMAT(1X,A1,'\')
CALL EXIT
END
SUBROUTINE CLEAR
TYPE 100, 27,27,27
100 FORMAT(1X,A1,'[H',A1,'[2J',A1,'P1ps(M0(L0)(AL0))S(I0)S(E)S[0,0]')
END
SUBROUTINE PLOT(X,Y)
IX=X+200
IY=Y+200
TYPE 100, IX,IY
100 FORMAT(' P[',I3,',',I3,']V[]')
END
|
| Rounding to produce pleasing results will depend on how GKS converts
from floating point coordinates to fixed point pixel-addresses.
If GKS adds an offset to (X,Y) so both are positive and then converts to
an integer (via a MODF instruction or similar) then adding a quantity
of about 1/2 a pixel to X and Y will be correct. However, if GKS
first truncates to an integer and then adds an integer offset to convert
to its pixel-address space this won't look good because there will be
a gap near (X,Y) = (0,0). In this case you would probably have to
simulate the rounding somehow to get the result you want. You might try
plotting a few shallow arcs or lines to see what GKS really does if its
documentation is not explicit on this.
Incidentally, there are much better ways of plotting circles in an
integer coordinate space than the ellipse approximation Stan referrs to,
that don't require multiply instructions in the loop. They're based on
always choosing X and Y steps (in a series of kings-moves) to minimize
the difference between (x(k)**2 + y(k)**2) and r**2. Gosper mentions the
technique in HAKMEM. Also, J. Bresenham has summarized the technique
in 'A Linear Algorithm for Incremental Digital Display of Circular Arcs',
CACM Feb-77. Too bad REGIS and GDS under P/OS don't use this technique -
it produces much better looking circular arcs than the 32-sided polygon
garbage they use.
But for your problem of generating regular polygons, a parameteric
equation is correct.
- Jim
|
| The Gosper reference that Jim Roth refers to is:
Item 177 (Drawing Curves Incrementally) in
Beeler, Gosper, and Schroeppel. HAKMEM. Artificial Intelligence Memo 239.
Massachusetts Institute of Technology, A.I. Laboratory.
Cambridge: 1972, page 82.
You should also study Turtle Geometry, and the LOGO computer language
in particular, for techniques for drawing curves incrementally.
|
| The following program compares circles plotted via the REGIS circle plotting
command vs an incremental algorithm. The increased accuracy of the incremental
method can be seen. Naturally, with direct acces to bitmap hardware the
incremental method would be no slower than drawing the polygon.
- Jim
PROGRAM CIRCPLOT
C
C Program to compare circles plotted with a polygon approximation
C vs circles plotted with an incremental algorithm.
C
IMPLICIT INTEGER*4(A-Z)
CALL CLEAR
C
DO 100 RAD=100,700,40
CALL CIRCLE(RAD,-90)
CALL CIRCLE(RAD+3,-90)
CALL CIRCLE(RAD+6,-90)
CALL CIRCLE(RAD+9,-90)
CALL ICIRCL(RAD+20,-90)
CALL ICIRCL(RAD+23,-90)
CALL ICIRCL(RAD+26,-90)
100 CALL ICIRCL(RAD+29,-90)
TYPE 101,27
101 FORMAT (1X,A1,'\')
STOP
END
SUBROUTINE ICIRCL(RAD,ARC)
C
C This subroutine can plot ellipses with proper choice of the constants A and B
C Alternatively, A and B can ba chosen to compensate for a non-isotropic
C display by drawing an ellipse which will be displayed as a circle.
C
C Note that any quadratic function (hyperbolas, etc) can be handled with this
C technique with proper starting values and error function calculation.
C
C Cubics would require keeping second order finite differences as well.
C
IMPLICIT INTEGER*4 (A-Z)
A = 1 ! Use a 1:1 aspec ratio
B = 1
A2 = A*A
B2 = B*B
X = RAD
Y = 0
ERR=0
DDEX = 2*A2 ! Init finite differences
DDEY = 2*B2
DEX=-A2+2*A2*X
DEY=-B2
TYPE 101, X
101 FORMAT(' P['I3,',0]')
C
100 CONTINUE
IF (X.LT.1) RETURN
EX = ERR+DEX ! Calc errors for various step directions
EY = ERR+DEY
EXY = EX+DEY
IF (IABS(EX).LT.IABS(EY)) GOTO 300 ! Choose one which minimizes
IF (IABS(EXY).LT.IABS(EY)) GOTO 400 ! resulting error
200 Y = Y+1
ERR = EY
DEY = DEY-DDEY
TYPE 201
201 FORMAT(' P[,+1]V[]')
GOTO 100
300 IF (IABS(EXY).LT.IABS(EX)) GOTO 400
X = X-1
ERR=EX
DEX = DEX-DDEX
TYPE 301
301 FORMAT(' P[-1,]V[]')
GOTO 100
400 X = X-1
Y = Y+1
ERR = EXY
DEX = DEX-DDEX
DEY = DEY-DDEY
TYPE 401
401 FORMAT(' P[-1,+1]V[]')
GOTO 100
END
SUBROUTINE CIRCLE(RAD,ARC)
IMPLICIT INTEGER*4 (A-Z)
TYPE 100, ARC,RAD,0
100 FORMAT (' P[0,0]C(A-90)[',I3,',]')
RETURN
END
SUBROUTINE CLEAR
IMPLICIT INTEGER*4 (A-Z)
ESC=27
TYPE 100, ESC,ESC,ESC
100 FORMAT(1X,A1,'[H',A1,'[2J',A1,'P1ps(M0(L0)(AL0))S(I0)S(E)S[0,0]')
RETURN
END
SUBROUTINE PLOT(X,Y)
IMPLICIT INTEGER*4 (A-Z)
TYPE 100, X,Y
100 FORMAT(' P[',I3,',',I3,']V[]')
RETURN
END
|