2 ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿
\r
3 ³ A General Conics Sections Scan Line Algorithm ³
\r
4 ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ
\r
7 The following code is the complete algorithm for the general conic
\r
8 drawer as mentioned in Foley/VanDam. It is included here with the
\r
9 permission of Andrew W. Fitzgibbon, who derived the remaining code
\r
10 sections not included in the book.
\r
14 // CONIC 2D Bresenham-like conic drawer.
\r
15 // CONIC(Sx,Sy, Ex,Ey, A,B,C,D,E,F) draws the conic specified
\r
16 // by A x^2 + B x y + C y^2 + D x + E y + F = 0, between the
\r
17 // start point (Sx, Sy) and endpoint (Ex,Ey).
\r
19 // Author: Andrew W. Fitzgibbon (andrewfg@ed.ac.uk),
\r
20 // Machine Vision Unit,
\r
21 // Dept. of Artificial Intelligence,
\r
22 // Edinburgh University,
\r
30 static int DIAGx[] = {999, 1, 1, -1, -1, -1, -1, 1, 1};
\r
31 static int DIAGy[] = {999, 1, 1, 1, 1, -1, -1, -1, -1};
\r
32 static int SIDEx[] = {999, 1, 0, 0, -1, -1, 0, 0, 1};
\r
33 static int SIDEy[] = {999, 0, 1, 1, 0, 0, -1, -1, 0};
\r
34 static int BSIGNS[] = {99, 1, 1, -1, -1, 1, 1, -1, -1};
\r
38 struct ConicPlotter {
\r
39 virtual void plot(int x, int y);
\r
42 struct DebugPlotter : public ConicPlotter {
\r
57 void plot(int x, int y);
\r
60 void DebugPlotter::plot(int x, int y)
\r
62 printf("%3d %3d\n",x,y);
\r
65 // Translate start point to origin...
\r
66 float tF = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
\r
67 float tD = D + 2 * A * xs + B * ys;
\r
68 float tE = E + B * xs + 2 * C * ys;
\r
70 float tx = x - xs + ((float)DIAGx[octant] + SIDEx[octant])/2;
\r
71 float ty = y - ys + ((float)DIAGy[octant] + SIDEy[octant])/2;
\r
74 float td = 4*(A*tx*tx + B*tx*ty + C*ty*ty + tD*tx + tE*ty + tF);
\r
76 fprintf(stderr,"O%d ", octant);
\r
78 fprintf(stderr," Inside ");
\r
80 fprintf(stderr,"Outside ");
\r
82 fprintf(stderr,"Real(%5.1f,%5.1f) = %8.2f Recurred = %8.2f err = %g\n",
\r
83 tx, ty, td/4, d/4.0f, err);
\r
84 if (fabs(err) > 1e-14)
\r
90 inline int odd(int n)
\r
95 inline int abs(int a)
\r
103 int getoctant(int gx, int gy)
\r
105 // Use gradient to identify octant.
\r
106 int upper = abs(gx)>abs(gy);
\r
107 if (gx>=0) // Right-pointing
\r
119 int conic(int xs, int ys, int xe, int ye,
\r
120 int A, int B, int C, int D, int E, int F,
\r
121 ConicPlotter * plotterdata)
\r
130 // Translate start point to origin...
\r
131 F = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
\r
132 D = D + 2 * A * xs + B * ys;
\r
133 E = E + B * xs + 2 * C * ys;
\r
135 // Work out starting octant
\r
136 int octant = getoctant(D,E);
\r
138 int dxS = SIDEx[octant];
\r
139 int dyS = SIDEy[octant];
\r
140 int dxD = DIAGx[octant];
\r
141 int dyD = DIAGy[octant];
\r
143 int bsign = BSIGNS[octant];
\r
147 d = A + B/2 + C/4 + D + E/2 + F;
\r
152 d = A/4 + B/2 + C + D/2 + E + F;
\r
157 d = A/4 - B/2 + C - D/2 + E + F;
\r
162 d = A - B/2 + C/4 - D + E/2 + F;
\r
167 d = A + B/2 + C/4 - D - E/2 + F;
\r
172 d = A/4 + B/2 + C - D/2 - E + F;
\r
177 d = A/4 - B/2 + C + D/2 - E + F;
\r
182 d = A - B/2 + C/4 + D - E/2 + F;
\r
187 fprintf(stderr,"FUNNY OCTANT\n");
\r
191 int k1sign = dyS*dyD;
\r
192 int k1 = 2 * (A + k1sign * (C - A));
\r
193 int Bsign = dxD*dyD;
\r
194 int k2 = k1 + Bsign * B;
\r
195 int k3 = 2 * (A + C + Bsign * B);
\r
197 // Work out gradient at endpoint
\r
200 int gx = 2*A*gxe + B*gye + D;
\r
201 int gy = B*gxe + 2*C*gye + E;
\r
203 int octantcount = getoctant(gx,gy) - octant;
\r
204 if (octantcount <= 0)
\r
205 octantcount = octantcount + 8;
\r
206 fprintf(stderr,"octantcount = %d\n", octantcount);
\r
211 while (octantcount > 0) {
\r
213 fprintf(stderr,"-- %d -------------------------\n", octant);
\r
216 while (2*v <= k2) {
\r
218 ((DebugPlotter*)plotterdata)->octant = octant;
\r
219 ((DebugPlotter*)plotterdata)->d = d;
\r
220 plotterdata->plot(x,y);
\r
222 // Are we inside or outside?
\r
223 if (d < 0) { // Inside
\r
239 d = d - u + v/2 - k2/2 + 3*k3/8;
\r
240 // error (^) in Foley and van Dam p 959, "2nd ed, revised 5th printing"
\r
241 u = -u + v - k2/2 + k3/2;
\r
243 k1 = k1 - 2*k2 + k3;
\r
245 int tmp = dxS; dxS = -dyS; dyS = tmp;
\r
247 else { // Octant is even
\r
250 ((DebugPlotter*)plotterdata)->octant = octant;
\r
251 ((DebugPlotter*)plotterdata)->d = d;
\r
252 plotterdata->plot(x,y);
\r
254 // Are we inside or outside?
\r
255 if (d > 0) { // Outside
\r
270 int tmpdk = k1 - k2;
\r
271 d = d + u - v + tmpdk;
\r
272 v = 2*u - v + tmpdk;
\r
277 int tmp = dxD; dxD = -dyD; dyD = tmp;
\r
280 octant = (octant&7)+1;
\r
284 // Draw final octant until we reach the endpoint
\r
286 fprintf(stderr,"-- %d (final) -----------------\n", octant);
\r
289 while (2*v <= k2 && x != xe && y != ye) {
\r
291 ((DebugPlotter*)plotterdata)->octant = octant;
\r
292 ((DebugPlotter*)plotterdata)->d = d;
\r
293 plotterdata->plot(x,y);
\r
295 // Are we inside or outside?
\r
296 if (d < 0) { // Inside
\r
312 else { // Octant is even
\r
313 while ((2*u < k2) && (x != xe) && (y != ye)) {
\r
315 ((DebugPlotter*)plotterdata)->octant = octant;
\r
316 ((DebugPlotter*)plotterdata)->d = d;
\r
317 plotterdata->plot(x,y);
\r
319 // Are we inside or outside?
\r
320 if (d > 0) { // Outside
\r
342 main(int argc, char ** argv)
\r
355 conic(db.xs,db.ys,db.xe,db.ye,db.A,db.B,db.C,db.D,db.E,db.F, &db);
\r