home *** CD-ROM | disk | FTP | other *** search
/ HomeWare 14 / HOMEWARE14.bin / prog / pcgpe10.arj / CONIC.CC < prev    next >
C/C++ Source or Header  |  1994-05-10  |  9KB  |  357 lines
  1.  
  2.            ┌───────────────────────────────────────────────┐
  3.            │ A General Conics Sections Scan Line Algorithm │
  4.            └───────────────────────────────────────────────┘
  5.  
  6.  
  7. The following code is the complete algorithm for the general conic
  8. drawer as mentioned in Foley/VanDam. It is included here with the
  9. permission of Andrew W. Fitzgibbon, who derived the remaining code
  10. sections not included in the book.
  11.  
  12.  
  13. //
  14. // CONIC  2D Bresenham-like conic drawer.
  15. //       CONIC(Sx,Sy, Ex,Ey, A,B,C,D,E,F) draws the conic specified
  16. //       by A x^2 + B x y + C y^2 + D x + E y + F = 0, between the
  17. //       start point (Sx, Sy) and endpoint (Ex,Ey).
  18.  
  19. // Author: Andrew W. Fitzgibbon (andrewfg@ed.ac.uk),
  20. //         Machine Vision Unit,
  21. //         Dept. of Artificial Intelligence,
  22. //         Edinburgh University,
  23. //         
  24. // Date: 31-Mar-94
  25.  
  26. #include <stdlib.h>
  27. #include <stdio.h>
  28. #include <math.h>
  29.  
  30. static int DIAGx[] = {999, 1,  1, -1, -1, -1, -1,  1,  1};
  31. static int DIAGy[] = {999, 1,  1,  1,  1, -1, -1, -1, -1};
  32. static int SIDEx[] = {999, 1,  0,  0, -1, -1,  0,  0,  1};
  33. static int SIDEy[] = {999, 0,  1,  1,  0,  0, -1, -1,  0};
  34. static int BSIGNS[] = {99, 1,  1, -1, -1,  1,  1, -1, -1};
  35.  
  36. int   debugging = 1;
  37.  
  38. struct ConicPlotter {
  39.   virtual void plot(int x, int y);
  40. };
  41.  
  42. struct DebugPlotter : public ConicPlotter {
  43.   int xs;
  44.   int ys;
  45.   int xe;
  46.   int ye;
  47.   int A;
  48.   int B;
  49.   int C;
  50.   int D;
  51.   int E;
  52.   int F;      
  53.  
  54.   int octant;
  55.   int d;
  56.  
  57.   void plot(int x, int y);
  58. };
  59.  
  60. void DebugPlotter::plot(int x, int y)
  61. {
  62.   printf("%3d %3d\n",x,y);
  63.  
  64.   if (debugging) {
  65.     // Translate start point to origin...
  66.     float tF = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
  67.     float tD = D + 2 * A * xs + B * ys;
  68.     float tE = E + B * xs + 2 * C * ys;
  69.   
  70.     float tx = x - xs + ((float)DIAGx[octant] + SIDEx[octant])/2;
  71.     float ty = y - ys + ((float)DIAGy[octant] + SIDEy[octant])/2;
  72.     // Calculate F
  73.     
  74.     float td = 4*(A*tx*tx + B*tx*ty + C*ty*ty + tD*tx + tE*ty + tF);
  75.     
  76.     fprintf(stderr,"O%d ", octant);
  77.     if (d<0)
  78.       fprintf(stderr," Inside "); 
  79.     else 
  80.       fprintf(stderr,"Outside "); 
  81.     float err = td - d;
  82.     fprintf(stderr,"Real(%5.1f,%5.1f) = %8.2f Recurred = %8.2f err = %g\n", 
  83.             tx, ty, td/4, d/4.0f, err);
  84.     if (fabs(err) > 1e-14)
  85.       abort();
  86.   }
  87.   
  88. }
  89.  
  90. inline int odd(int n)
  91. {
  92.   return n&1;
  93. }
  94.  
  95. inline int abs(int a)
  96. {
  97.   if (a > 0)
  98.     return a;
  99.   else
  100.     return -a;
  101. }
  102.     
  103. int getoctant(int gx, int gy)
  104. {
  105.   // Use gradient to identify octant.
  106.   int upper = abs(gx)>abs(gy);
  107.   if (gx>=0)                            // Right-pointing
  108.     if (gy>=0)                          //    Up
  109.       return 4 - upper;
  110.     else                                //    Down
  111.       return 1 + upper;
  112.   else                                  // Left
  113.     if (gy>0)                           //    Up
  114.       return 5 + upper;
  115.     else                                //    Down
  116.       return 8 - upper;
  117. }
  118.  
  119. int conic(int xs, int ys, int xe, int ye,
  120.           int A, int B, int C, int D, int E, int F,
  121.           ConicPlotter * plotterdata)
  122. {
  123.   A *= 4;
  124.   B *= 4;
  125.   C *= 4;
  126.   D *= 4;
  127.   E *= 4;
  128.   F *= 4;
  129.   
  130.   // Translate start point to origin...
  131.   F = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
  132.   D = D + 2 * A * xs + B * ys;
  133.   E = E + B * xs + 2 * C * ys;
  134.   
  135.   // Work out starting octant
  136.   int octant = getoctant(D,E);
  137.   
  138.   int dxS = SIDEx[octant]; 
  139.   int dyS = SIDEy[octant]; 
  140.   int dxD = DIAGx[octant];
  141.   int dyD = DIAGy[octant];
  142.  
  143.   int bsign = BSIGNS[octant];
  144.   int d,u,v;
  145.   switch (octant) {
  146.   case 1:
  147.     d = A + B/2 + C/4 + D + E/2 + F;
  148.     u = A + B/2 + D;
  149.     v = u + E;
  150.     break;
  151.   case 2:
  152.     d = A/4 + B/2 + C + D/2 + E + F;
  153.     u = B/2 + C + E;
  154.     v = u + D;
  155.     break;
  156.   case 3:
  157.     d = A/4 - B/2 + C - D/2 + E + F;
  158.     u = -B/2 + C + E;
  159.     v = u - D;
  160.     break;
  161.   case 4:
  162.     d = A - B/2 + C/4 - D + E/2 + F;
  163.     u = A - B/2 - D;
  164.     v = u + E;
  165.     break;
  166.   case 5:
  167.     d = A + B/2 + C/4 - D - E/2 + F;
  168.     u = A + B/2 - D;
  169.     v = u - E;
  170.     break;
  171.   case 6:
  172.     d = A/4 + B/2 + C - D/2 - E + F;
  173.     u = B/2 + C - E;
  174.     v = u - D;
  175.     break;
  176.   case 7:
  177.     d = A/4 - B/2 + C + D/2 - E + F;
  178.     u =  -B/2 + C - E;
  179.     v = u + D;
  180.     break;
  181.   case 8:
  182.     d = A - B/2 + C/4 + D - E/2 + F;
  183.     u = A - B/2 + D;
  184.     v = u - E;
  185.     break;
  186.   default:
  187.     fprintf(stderr,"FUNNY OCTANT\n");
  188.     abort();
  189.   }
  190.   
  191.   int k1sign = dyS*dyD;
  192.   int k1 = 2 * (A + k1sign * (C - A));
  193.   int Bsign = dxD*dyD;
  194.   int k2 = k1 + Bsign * B;
  195.   int k3 = 2 * (A + C + Bsign * B);
  196.  
  197.   // Work out gradient at endpoint
  198.   int gxe = xe - xs;
  199.   int gye = ye - ys;
  200.   int gx = 2*A*gxe +   B*gye + D;
  201.   int gy =   B*gxe + 2*C*gye + E;
  202.   
  203.   int octantcount = getoctant(gx,gy) - octant;
  204.   if (octantcount <= 0)
  205.     octantcount = octantcount + 8;
  206.   fprintf(stderr,"octantcount = %d\n", octantcount);
  207.   
  208.   int x = xs;
  209.   int y = ys;
  210.   
  211.   while (octantcount > 0) {
  212.     if (debugging)
  213.       fprintf(stderr,"-- %d -------------------------\n", octant); 
  214.     
  215.     if (odd(octant)) {
  216.       while (2*v <= k2) {
  217.         // Plot this point
  218.         ((DebugPlotter*)plotterdata)->octant = octant;
  219.         ((DebugPlotter*)plotterdata)->d = d;
  220.         plotterdata->plot(x,y);
  221.         
  222.         // Are we inside or outside?
  223.         if (d < 0) {                    // Inside
  224.           x = x + dxS;
  225.           y = y + dyS;
  226.           u = u + k1;
  227.           v = v + k2;
  228.           d = d + u;
  229.         }
  230.         else {                          // outside
  231.           x = x + dxD;
  232.           y = y + dyD;
  233.           u = u + k2;
  234.           v = v + k3;
  235.           d = d + v;
  236.         }
  237.       }
  238.     
  239.       d = d - u + v/2 - k2/2 + 3*k3/8; 
  240.       // error (^) in Foley and van Dam p 959, "2nd ed, revised 5th printing"
  241.       u = -u + v - k2/2 + k3/2;
  242.       v = v - k2 + k3/2;
  243.       k1 = k1 - 2*k2 + k3;
  244.       k2 = k3 - k2;
  245.       int tmp = dxS; dxS = -dyS; dyS = tmp;
  246.     }
  247.     else {                              // Octant is even
  248.       while (2*u < k2) {
  249.         // Plot this point
  250.         ((DebugPlotter*)plotterdata)->octant = octant;
  251.         ((DebugPlotter*)plotterdata)->d = d;
  252.         plotterdata->plot(x,y);
  253.         
  254.         // Are we inside or outside?
  255.         if (d > 0) {                    // Outside
  256.           x = x + dxS;
  257.           y = y + dyS;
  258.           u = u + k1;
  259.           v = v + k2;
  260.           d = d + u;
  261.         }
  262.         else {                          // Inside
  263.           x = x + dxD;
  264.           y = y + dyD;
  265.           u = u + k2;
  266.           v = v + k3;
  267.           d = d + v;
  268.         }
  269.       }
  270.       int tmpdk = k1 - k2;
  271.       d = d + u - v + tmpdk;
  272.       v = 2*u - v + tmpdk;
  273.       u = u + tmpdk;
  274.       k3 = k3 + 4*tmpdk;
  275.       k2 = k1 + tmpdk;
  276.       
  277.       int tmp = dxD; dxD = -dyD; dyD = tmp;
  278.     }
  279.     
  280.     octant = (octant&7)+1;
  281.     octantcount--;
  282.   }
  283.  
  284.   // Draw final octant until we reach the endpoint
  285.   if (debugging)
  286.     fprintf(stderr,"-- %d (final) -----------------\n", octant); 
  287.     
  288.   if (odd(octant)) {
  289.     while (2*v <= k2 && x != xe && y != ye) {
  290.       // Plot this point
  291.       ((DebugPlotter*)plotterdata)->octant = octant;
  292.       ((DebugPlotter*)plotterdata)->d = d;
  293.       plotterdata->plot(x,y);
  294.       
  295.       // Are we inside or outside?
  296.       if (d < 0) {                      // Inside
  297.         x = x + dxS;
  298.         y = y + dyS;
  299.         u = u + k1;
  300.         v = v + k2;
  301.         d = d + u;
  302.       }
  303.       else {                            // outside
  304.         x = x + dxD;
  305.         y = y + dyD;
  306.         u = u + k2;
  307.         v = v + k3;
  308.         d = d + v;
  309.       }
  310.     }
  311.   }
  312.   else {                        // Octant is even
  313.     while ((2*u < k2) && (x != xe) && (y != ye)) {
  314.       // Plot this point
  315.       ((DebugPlotter*)plotterdata)->octant = octant;
  316.       ((DebugPlotter*)plotterdata)->d = d;
  317.       plotterdata->plot(x,y);
  318.       
  319.       // Are we inside or outside?
  320.       if (d > 0) {                      // Outside
  321.         x = x + dxS;
  322.         y = y + dyS;
  323.         u = u + k1;
  324.         v = v + k2;
  325.         d = d + u;
  326.       }
  327.       else {                            // Inside
  328.         x = x + dxD;
  329.         y = y + dyD;
  330.         u = u + k2;
  331.         v = v + k3;
  332.         d = d + v;
  333.       }
  334.     }
  335.   }
  336.  
  337.  
  338.  
  339.   return 1;
  340. }
  341.  
  342. main(int argc, char ** argv)
  343. {
  344.   DebugPlotter db;
  345.   db.xs = -7;
  346.   db.ys = -19;
  347.   db.xe = -8;
  348.   db.ye = -8;
  349.   db.A = 1424;
  350.   db.B = -964;
  351.   db.C = 276;
  352.   db.D = 0;
  353.   db.E = 0;
  354.   db.F = -40000;
  355.   conic(db.xs,db.ys,db.xe,db.ye,db.A,db.B,db.C,db.D,db.E,db.F, &db);
  356. }
  357.