16/PCGPE10/CONIC.CC

   1 \r
   2            ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿\r
   3            ³ A General Conics Sections Scan Line Algorithm ³\r
   4            ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ\r
   5 \r
   6 \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
  11 \r
  12 \r
  13 //\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
  18 \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
  23 //         \r
  24 // Date: 31-Mar-94\r
  25 \r
  26 #include <stdlib.h>\r
  27 #include <stdio.h>\r
  28 #include <math.h>\r
  29 \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
  35 \r
  36 int   debugging = 1;\r
  37 \r
  38 struct ConicPlotter {\r
  39   virtual void plot(int x, int y);\r
  40 };\r
  41 \r
  42 struct DebugPlotter : public ConicPlotter {\r
  43   int xs;\r
  44   int ys;\r
  45   int xe;\r
  46   int ye;\r
  47   int A;\r
  48   int B;\r
  49   int C;\r
  50   int D;\r
  51   int E;\r
  52   int F;      \r
  53 \r
  54   int octant;\r
  55   int d;\r
  56 \r
  57   void plot(int x, int y);\r
  58 };\r
  59 \r
  60 void DebugPlotter::plot(int x, int y)\r
  61 {\r
  62   printf("%3d %3d\n",x,y);\r
  63 \r
  64   if (debugging) {\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
  69   \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
  72     // Calculate F\r
  73     \r
  74     float td = 4*(A*tx*tx + B*tx*ty + C*ty*ty + tD*tx + tE*ty + tF);\r
  75     \r
  76     fprintf(stderr,"O%d ", octant);\r
  77     if (d<0)\r
  78       fprintf(stderr," Inside "); \r
  79     else \r
  80       fprintf(stderr,"Outside "); \r
  81     float err = td - d;\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
  85       abort();\r
  86   }\r
  87   \r
  88 }\r
  89 \r
  90 inline int odd(int n)\r
  91 {\r
  92   return n&1;\r
  93 }\r
  94 \r
  95 inline int abs(int a)\r
  96 {\r
  97   if (a > 0)\r
  98     return a;\r
  99   else\r
 100     return -a;\r
 101 }\r
 102     \r
 103 int getoctant(int gx, int gy)\r
 104 {\r
 105   // Use gradient to identify octant.\r
 106   int upper = abs(gx)>abs(gy);\r
 107   if (gx>=0)                            // Right-pointing\r
 108     if (gy>=0)                          //    Up\r
 109       return 4 - upper;\r
 110     else                                //    Down\r
 111       return 1 + upper;\r
 112   else                                  // Left\r
 113     if (gy>0)                           //    Up\r
 114       return 5 + upper;\r
 115     else                                //    Down\r
 116       return 8 - upper;\r
 117 }\r
 118 \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
 122 {\r
 123   A *= 4;\r
 124   B *= 4;\r
 125   C *= 4;\r
 126   D *= 4;\r
 127   E *= 4;\r
 128   F *= 4;\r
 129   \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
 134   \r
 135   // Work out starting octant\r
 136   int octant = getoctant(D,E);\r
 137   \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
 142 \r
 143   int bsign = BSIGNS[octant];\r
 144   int d,u,v;\r
 145   switch (octant) {\r
 146   case 1:\r
 147     d = A + B/2 + C/4 + D + E/2 + F;\r
 148     u = A + B/2 + D;\r
 149     v = u + E;\r
 150     break;\r
 151   case 2:\r
 152     d = A/4 + B/2 + C + D/2 + E + F;\r
 153     u = B/2 + C + E;\r
 154     v = u + D;\r
 155     break;\r
 156   case 3:\r
 157     d = A/4 - B/2 + C - D/2 + E + F;\r
 158     u = -B/2 + C + E;\r
 159     v = u - D;\r
 160     break;\r
 161   case 4:\r
 162     d = A - B/2 + C/4 - D + E/2 + F;\r
 163     u = A - B/2 - D;\r
 164     v = u + E;\r
 165     break;\r
 166   case 5:\r
 167     d = A + B/2 + C/4 - D - E/2 + F;\r
 168     u = A + B/2 - D;\r
 169     v = u - E;\r
 170     break;\r
 171   case 6:\r
 172     d = A/4 + B/2 + C - D/2 - E + F;\r
 173     u = B/2 + C - E;\r
 174     v = u - D;\r
 175     break;\r
 176   case 7:\r
 177     d = A/4 - B/2 + C + D/2 - E + F;\r
 178     u =  -B/2 + C - E;\r
 179     v = u + D;\r
 180     break;\r
 181   case 8:\r
 182     d = A - B/2 + C/4 + D - E/2 + F;\r
 183     u = A - B/2 + D;\r
 184     v = u - E;\r
 185     break;\r
 186   default:\r
 187     fprintf(stderr,"FUNNY OCTANT\n");\r
 188     abort();\r
 189   }\r
 190   \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
 196 \r
 197   // Work out gradient at endpoint\r
 198   int gxe = xe - xs;\r
 199   int gye = ye - ys;\r
 200   int gx = 2*A*gxe +   B*gye + D;\r
 201   int gy =   B*gxe + 2*C*gye + E;\r
 202   \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
 207   \r
 208   int x = xs;\r
 209   int y = ys;\r
 210   \r
 211   while (octantcount > 0) {\r
 212     if (debugging)\r
 213       fprintf(stderr,"-- %d -------------------------\n", octant); \r
 214     \r
 215     if (odd(octant)) {\r
 216       while (2*v <= k2) {\r
 217         // Plot this point\r
 218         ((DebugPlotter*)plotterdata)->octant = octant;\r
 219         ((DebugPlotter*)plotterdata)->d = d;\r
 220         plotterdata->plot(x,y);\r
 221         \r
 222         // Are we inside or outside?\r
 223         if (d < 0) {                    // Inside\r
 224           x = x + dxS;\r
 225           y = y + dyS;\r
 226           u = u + k1;\r
 227           v = v + k2;\r
 228           d = d + u;\r
 229         }\r
 230         else {                          // outside\r
 231           x = x + dxD;\r
 232           y = y + dyD;\r
 233           u = u + k2;\r
 234           v = v + k3;\r
 235           d = d + v;\r
 236         }\r
 237       }\r
 238     \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
 242       v = v - k2 + k3/2;\r
 243       k1 = k1 - 2*k2 + k3;\r
 244       k2 = k3 - k2;\r
 245       int tmp = dxS; dxS = -dyS; dyS = tmp;\r
 246     }\r
 247     else {                              // Octant is even\r
 248       while (2*u < k2) {\r
 249         // Plot this point\r
 250         ((DebugPlotter*)plotterdata)->octant = octant;\r
 251         ((DebugPlotter*)plotterdata)->d = d;\r
 252         plotterdata->plot(x,y);\r
 253         \r
 254         // Are we inside or outside?\r
 255         if (d > 0) {                    // Outside\r
 256           x = x + dxS;\r
 257           y = y + dyS;\r
 258           u = u + k1;\r
 259           v = v + k2;\r
 260           d = d + u;\r
 261         }\r
 262         else {                          // Inside\r
 263           x = x + dxD;\r
 264           y = y + dyD;\r
 265           u = u + k2;\r
 266           v = v + k3;\r
 267           d = d + v;\r
 268         }\r
 269       }\r
 270       int tmpdk = k1 - k2;\r
 271       d = d + u - v + tmpdk;\r
 272       v = 2*u - v + tmpdk;\r
 273       u = u + tmpdk;\r
 274       k3 = k3 + 4*tmpdk;\r
 275       k2 = k1 + tmpdk;\r
 276       \r
 277       int tmp = dxD; dxD = -dyD; dyD = tmp;\r
 278     }\r
 279     \r
 280     octant = (octant&7)+1;\r
 281     octantcount--;\r
 282   }\r
 283 \r
 284   // Draw final octant until we reach the endpoint\r
 285   if (debugging)\r
 286     fprintf(stderr,"-- %d (final) -----------------\n", octant); \r
 287     \r
 288   if (odd(octant)) {\r
 289     while (2*v <= k2 && x != xe && y != ye) {\r
 290       // Plot this point\r
 291       ((DebugPlotter*)plotterdata)->octant = octant;\r
 292       ((DebugPlotter*)plotterdata)->d = d;\r
 293       plotterdata->plot(x,y);\r
 294       \r
 295       // Are we inside or outside?\r
 296       if (d < 0) {                      // Inside\r
 297         x = x + dxS;\r
 298         y = y + dyS;\r
 299         u = u + k1;\r
 300         v = v + k2;\r
 301         d = d + u;\r
 302       }\r
 303       else {                            // outside\r
 304         x = x + dxD;\r
 305         y = y + dyD;\r
 306         u = u + k2;\r
 307         v = v + k3;\r
 308         d = d + v;\r
 309       }\r
 310     }\r
 311   }\r
 312   else {                        // Octant is even\r
 313     while ((2*u < k2) && (x != xe) && (y != ye)) {\r
 314       // Plot this point\r
 315       ((DebugPlotter*)plotterdata)->octant = octant;\r
 316       ((DebugPlotter*)plotterdata)->d = d;\r
 317       plotterdata->plot(x,y);\r
 318       \r
 319       // Are we inside or outside?\r
 320       if (d > 0) {                      // Outside\r
 321         x = x + dxS;\r
 322         y = y + dyS;\r
 323         u = u + k1;\r
 324         v = v + k2;\r
 325         d = d + u;\r
 326       }\r
 327       else {                            // Inside\r
 328         x = x + dxD;\r
 329         y = y + dyD;\r
 330         u = u + k2;\r
 331         v = v + k3;\r
 332         d = d + v;\r
 333       }\r
 334     }\r
 335   }\r
 336 \r
 337 \r
 338 \r
 339   return 1;\r
 340 }\r
 341 \r
 342 main(int argc, char ** argv)\r
 343 {\r
 344   DebugPlotter db;\r
 345   db.xs = -7;\r
 346   db.ys = -19;\r
 347   db.xe = -8;\r
 348   db.ye = -8;\r
 349   db.A = 1424;\r
 350   db.B = -964;\r
 351   db.C = 276;\r
 352   db.D = 0;\r
 353   db.E = 0;\r
 354   db.F = -40000;\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
 356 }\r