Groups	Advanced Groups Search    Preferences    Groups Help

In article <13656@ztivax.UUCP>, rencken@ztivax.UUCP (Dr. Wolfgang Rencken) writes:
|> Does anyone out there know where I can obtain a library function or 
|> C source code of a routine, which is capable of drawing a rotated ellipse
|> under X11. The Xlib routine XDrawArc can only render ellipses which
|> are parallel to the x-axis!!
|> Any help will be greatly appreciated. Thanks in advance.

Here is a nice routine written by James Tough.  It will be incorporated into
xfig 2.1.5 and is therefore freely copyable as long as James' name and the
xfig copyright notice is kept intact.

First are his original routines unmodified for X.  There are three routines
with their descriptions.  The lines generated produce a filled ellipse.

Next, is the one I used in xfig, modified for X.


-- 
Brian V. Smith    (bvsmith@lbl.gov)
Lawrence Berkeley Laboratory
I don't speak for LBL; they don't pay me enough for that.


James Tough's original routines:

/*
 *  An Ellipse Generator.
 *  Written by James Tough   7th May 92
 *
 *  The following routines displays a filled ellipse on the screen from the
 *    semi-minor axis 'a', semi-major axis 'b' and angle of rotation
 *    'phi'.
 *
 *  It works along these principles .....
 *
 *        The standard ellipse equation is
 *
 *             x*x     y*y
 *             ---  +  ---
 *             a*a     b*b
 *
 *
 *        Rotation of a point (x,y) is well known through the use of
 *
 *            x' = x*COS(phi) - y*SIN(phi)
 *            y' = y*COS(phi) + y*COS(phi)
 *
 *        Taking these to together, this gives the equation for a rotated
 *      ellipse centered around the origin.
 *
 *           [x*COS(phi) - y*SIN(phi)]^2   [x*SIN(phi) + y*COS(phi)]^2
 *           --------------------------- + ---------------------------
 *                      a*a                           b*b
 *
 *        NOTE -  some of the above equation can be precomputed, eg,
 *
 *              i = COS(phi)/a        and        j = SIN(phi)/b
 *       
 *        NOTE -  y is constant for each line so,
 *
 *              m = -yk*SIN(phi)/a    and     n = yk*COS(phi)/b
 *              where yk stands for the kth line ( y subscript k)
 *
 *        Where yk=y, we get a quadratic,
 *
 *              (i*x + m)^2 + (j*x + n)^2 = 1
 *
 *        Thus for any particular line, y, there is two corresponding x
 *      values. These are the roots of the quadratic. To get the roots, 
 *      the above equation can be rearranged using the standard method, 
 *
 *          -(i*m + j*n) +- sqrt[i^2 + j^2 - (i*n -j*m)^2]
 *      x = ----------------------------------------------
 *                           i^2 + j^2 
 *
 *        NOTE -  again much of this equation can be precomputed.
 *     
 *           c1 = i^2 + j^2 
 *           c2 = [COS(phi)*SIN(phi)*(a-b)]
 *           c3 = [b*b*(COS(phi)^2) + a*a*(SIN(phi)^2)]
 *           c4 = a*b/c3
 * 
 *      x = c2*y +- c4*sqrt(c3 - y*y),    where +- must be evaluated once
 *                                      for plus, and once for minus. 
 *
 *        We also need to know how large the ellipse is. This condition
 *      arises when the sqrt of the above equation evaluates to zero.
 *      Thus the height of the ellipse is give by 
 * 
 *              sqrt[ b*b*(COS(phi)^2) + a*a*(SIN(phi)^2) ]
 *
 *       which just happens to be equal to sqrt(c3).
 *
 *         It is now possible to create a routine that will scan convert 
 *       the ellipse on the screen.
 *
 *        NOTE -  c2 is the gradient of the new ellipse axis.
 *                c4 is the new semi-minor axis, 'a'.
 *           sqr(c3) is the new semi-major axis, 'b'.
 *
 *         These values could be used in a 4WS or 8WS ellipse generator
 *       that does not work on rotation, to give the feel of a rotated
 *       ellipse. These ellipses are not very accurate and give visable
 *       bumps along the edge of the ellipse. However, these routines
 *       are very quick, and give a good approximation to a rotated ellipse.
 *
 *       NOTES on the code given.
 * 
 *           All the routines take there parameters as ( x, y, a, b, phi ),
 *           where x,y are the center of the ellipse ( relative to the 
 *           origin ), a and b are the vertical and horizontal axis, and
 *           phi is the angle of rotation in RADIANS.
 *
 *           The 'moveto(x,y)' command moves the screen cursor to the 
 *               (x,y) point.
 *           The 'lineto(x,y)' command draws a line from the cursor to
 *               the point (x,y).
 *
 *
 *   Examples,
 *
 *       NOTE, all angles must be given in radians.
 *
 *       quickellipse(0,0,100,20,PI/4)         an ellipse at the origin,
 *                                             major axis 100, minor 20
 *                                             at angle 45 degrees.
 *
 *       ellipse4ws(100,100,100,20,PI/8)       an ellipse at (100,100)
 *                                             major axis 20, minor 100
 *                                             at angle 22.5 degrees.
  *                                 
 *
 */


#include <math.h>

void ellipse(float, float, float, float, float);
void quickellipse(float, float, float, float, float);
void ellipse4ws(float, float, float, float, float);


/*
 *  Ellipse,  uses floating point arithmetic to impliment the discussed
 *    method from above. Note, that the arithmetic has not been optimised
 *    in anyway.
 *
 *  Written by James Tough
 *  7th May 1992
 * 
 */


void ellipse(float xcen, float ycen, float a, float b, float phi)  {

 double c1, c2, c3, c4, c5, cphi, sphi, cphisqr, sphisqr;
 double xleft, xright, d, asqr, bsqr;
 int ymax, i=0;

 cphi = cos(phi);
 sphi = sin(phi);
 cphisqr = cphi*cphi;
 sphisqr = sphi*sphi;
 asqr = a*a;
 bsqr = b*b;
 
 c1 = (cphisqr/asqr)+(sphisqr/bsqr);
 c2 = ((cphi*sphi/bsqr)-(cphi*sphi/asqr))/c1;
 c3 = (bsqr*cphisqr) + (asqr*sphisqr);
 ymax = sqrt(c3);
 c4 = a*b/c3;
 c5 = 0;

 if (ymax % 2) {
  d = sqrt(c3);
  moveto(xcen+d,ycen+i);
  lineto(xcen+d,ycen+i);
  c5 = c2;
  i=1;
 }
 
 for (; i<ymax; i++) {
  d = c4*sqrt(c3-(i*i));
  xleft = c5-d;
  xright = c5+d;                        
  moveto(xcen+xleft,ycen+i);
  lineto(xcen+xright,ycen+i);
  moveto(xcen-xleft,ycen-i);
  lineto(xcen-xright,ycen-i);
  c5+=c2;
 }

 return;
}


/*
 *  QuickEllipse, uses the same method as Ellipse, but uses incremental
 *    methods to reduce the amount of work that has to be done inside
 *    the main loop. The speed increase is very noticeable.
 *
 *  Written by James Tough
 *  7th May 1992
 * 
 */


void quickellipse(float xcen, float ycen, float a, float b, float phi)  {

 double c1, c2, c3, c4, c5, c6, v1, cphi, sphi, cphisqr, sphisqr;
 double xleft, xright, d, asqr, bsqr;
 int ymax, y=0;

 cphi = cos(phi);
 sphi = sin(phi);
 cphisqr = cphi*cphi;
 sphisqr = sphi*sphi;
 asqr = a*a;
 bsqr = b*b;
 
 c1 = (cphisqr/asqr)+(sphisqr/bsqr);
 c2 = ((cphi*sphi/asqr)-(cphi*sphi/bsqr))/c1;
 c3 = (bsqr*cphisqr) + (asqr*sphisqr);
 ymax = sqrt(c3);
 c4 = a*b/c3;
 c5 = 0;
 v1 = c4*c4;
 c6 = 2*v1;
 c3 = c3*v1-v1;

 if (ymax % 2) {
  d = sqrt(c3);
  moveto(xcen+d,ycen);
  lineto(xcen+d,ycen);
  c5 = c2;
  y=1;
 }
 
 while (c3>=0) {
  d = sqrt(c3);
  xleft = c5-d;
  xright = c5+d;                        
  moveto(xcen+xleft,ycen+y);
  lineto(xcen+xright,ycen+y);
  moveto(xcen-xleft,ycen-y);
  lineto(xcen-xright,ycen-y);
  c5+=c2;
  v1+=c6;
  c3-=v1;
  y=y+1;
 }

 return;
}



/*
 *  Ellipse4WS,  uses a four way symmetry algorithm ( not discussed ) to
 *    approximate the rotated ellipse. The method is to draw the non
 *    rotated ellipse as usual, but to add a gradient ('c2') to an
 *    offset every iteration. The drawn ellipse then appears rotated.
 *    The speed of this algorithm is achieved because all of the ellipse

 *    routine can be handled in integer arithmetic with only additions 
 *    and subtractions. It is therefore a good method to draw ellipses
 *    in assembler so long as accuracy is not a high priority.
 *
 *  NOTE - There is a lot that can be done to speed up the running of this
 *       routine. The implimentation is meant to be instructive.
 *
 *  Written by James Tough
 *  7th May 1992
 * 
 */

void ellipse4ws(float xcen, float ycen, float a, float b, float phi)  {

        double  c, g, o=0, cphi, sphi, cphisqr, sphisqr;
        double  xleft, xright, asqr, bsqr;
 long x, y=0, t1, t2, t4, t5, t6, t7=0, d, ymax;
        int     i=0, loop;

        cphi = cos(phi);
        sphi = sin(phi);
        cphisqr = cphi*cphi;
        asqr = a*a;
        bsqr = b*b;
       
        c = (bsqr*cphisqr) + (asqr*sphisqr);
        g = (cphi*sphi*(asqr-bsqr))/c;
 ymax = sqrt(c);
 x = (a*b)/ymax;

 t1 = x*x;
 t2 = 2*t1;
 t4 = 2*c;
 t5 = x*c;
 t6 = 2*t5;
 d = ((t1+c)/4)-t5;
 loop = x+ymax;

 for (i=0; i<loop; i++) 
  if (d<0) {
   xleft = o-x;
   xright = o+x;                        
   moveto(xcen+xleft,ycen+y);
   lineto(xcen+xright,ycen+y);
   moveto(xcen-xleft,ycen-y);
   lineto(xcen-xright,ycen-y);
   o=o+g;
   t7=t7+t2;
   d=d+t7;
   y++;
  } else {
   x=x-1; 
   d=d-t6;
   t6=t6-t4; 
  }
 
 return;
}

***************************************************
As adapted for xfig:

/*
 *  An Ellipse Generator.
 *  Written by James Tough   7th May 92
 *
 *  The following routines displays a filled ellipse on the screen from the
 *    semi-minor axis 'a', semi-major axis 'b' and angle of rotation
 *    'phi'.
 *
 *  It works along these principles .....
 *
 *        The standard ellipse equation is
 *
 *             x*x     y*y
 *             ---  +  ---
 *             a*a     b*b
 *
 *
 *        Rotation of a point (x,y) is well known through the use of
 *
 *            x' = x*COS(phi) - y*SIN(phi)
 *            y' = y*COS(phi) + y*COS(phi)
 *
 *        Taking these to together, this gives the equation for a rotated
 *      ellipse centered around the origin.
 *
 *           [x*COS(phi) - y*SIN(phi)]^2   [x*SIN(phi) + y*COS(phi)]^2
 *           --------------------------- + ---------------------------
 *                      a*a                           b*b
 *
 *        NOTE -  some of the above equation can be precomputed, eg,
 *
 *              i = COS(phi)/a        and        j = SIN(phi)/b
 *       
 *        NOTE -  y is constant for each line so,
 *
 *              m = -yk*SIN(phi)/a    and     n = yk*COS(phi)/b
 *              where yk stands for the kth line ( y subscript k)
 *
 *        Where yk=y, we get a quadratic,
 *
 *              (i*x + m)^2 + (j*x + n)^2 = 1
 *
 *        Thus for any particular line, y, there is two corresponding x
 *      values. These are the roots of the quadratic. To get the roots, 
 *      the above equation can be rearranged using the standard method, 
 *
 *          -(i*m + j*n) +- sqrt[i^2 + j^2 - (i*n -j*m)^2]
 *      x = ----------------------------------------------
 *                           i^2 + j^2 
 *
 *        NOTE -  again much of this equation can be precomputed.
 *     
 *           c1 = i^2 + j^2 
 *           c2 = [COS(phi)*SIN(phi)*(a-b)]
 *           c3 = [b*b*(COS(phi)^2) + a*a*(SIN(phi)^2)]
 *           c4 = a*b/c3
 * 
 *      x = c2*y +- c4*sqrt(c3 - y*y),    where +- must be evaluated once
 *                                      for plus, and once for minus. 
 *
 *        We also need to know how large the ellipse is. This condition
 *      arises when the sqrt of the above equation evaluates to zero.
 *      Thus the height of the ellipse is give by 
 * 
 *              sqrt[ b*b*(COS(phi)^2) + a*a*(SIN(phi)^2) ]
 *
 *       which just happens to be equal to sqrt(c3).
 *
 *         It is now possible to create a routine that will scan convert 
 *       the ellipse on the screen.
 *
 *        NOTE -  c2 is the gradient of the new ellipse axis.
 *                c4 is the new semi-minor axis, 'a'.
 *           sqr(c3) is the new semi-major axis, 'b'.
 *
 *         These values could be used in a 4WS or 8WS ellipse generator
 *       that does not work on rotation, to give the feel of a rotated
 *       ellipse. These ellipses are not very accurate and give visable
 *       bumps along the edge of the ellipse. However, these routines
 *       are very quick, and give a good approximation to a rotated ellipse.
 *
 *       NOTES on the code given.
 * 
 *           All the routines take there parameters as ( x, y, a, b, phi ),
 *           where x,y are the center of the ellipse ( relative to the 
 *           origin ), a and b are the vertical and horizontal axis, and
 *           phi is the angle of rotation in RADIANS.
 *
 *           The 'moveto(x,y)' command moves the screen cursor to the 
 *               (x,y) point.
 *           The 'lineto(x,y)' command draws a line from the cursor to
 *               the point (x,y).
 *
 *
 *   Examples,
 *
 *       NOTE, all angles must be given in radians.
 *
 *       angle_ellipse(0,0,100,20,PI/4)         an ellipse at the origin,
 *                                             major axis 100, minor 20
 *                                             at angle 45 degrees.
 *
 *
 */


/*
 *  QuickEllipse, uses the same method as Ellipse, but uses incremental
 *    methods to reduce the amount of work that has to be done inside
 *    the main loop. The speed increase is very noticeable.
 *
 *  Written by James Tough
 *  7th May 1992
 * 
 */

static int x[MAXNUMPTS/4][4],y[MAXNUMPTS/4][4];
static int nump[4];
static int totpts,i,j;
static int order[4]={0,1,3,2};

angle_ellipse(center_x, center_y, radius_x, radius_y, angle,
       op, thickness, style, style_val, fill_style, color)
    int      center_x, center_y;
    int      radius_x, radius_y;
    float     angle;
    int      op,thickness,style,fill_style,color;
    float     style_val;
{
 float xcen, ycen, a, b; 

 double c1, c2, c3, c4, c5, c6, v1, cphi, sphi, cphisqr, sphisqr;
 double xleft, xright, d, asqr, bsqr;
 int ymax, yy=0;
 int k,m,n,dir;
 XPoint *pnts,*ppnts;

 if (radius_x == 0 || radius_y == 0)
  return;

 xcen = center_x;
 ycen = center_y;
 a = radius_x;
 b = radius_y;

 cphi = cos((double)angle);
 sphi = sin((double)angle);
 cphisqr = cphi*cphi;
 sphisqr = sphi*sphi;
 asqr = a*a;
 bsqr = b*b;
 
 c1 = (cphisqr/asqr)+(sphisqr/bsqr);
 c2 = ((cphi*sphi/asqr)-(cphi*sphi/bsqr))/c1;
 c3 = (bsqr*cphisqr) + (asqr*sphisqr);
 ymax = sqrt(c3);
 c4 = a*b/c3;
 c5 = 0;
 v1 = c4*c4;
 c6 = 2*v1;
 c3 = c3*v1-v1;
 totpts = 0;
 for (i=0; i<=3; i++)
  nump[i]=0;
 i=0; j=0;
 /* odd first points */
 if (ymax % 2) {
  d = sqrt(c3);
  newpoint(xcen-d,ycen);
  newpoint(xcen+d,ycen);
  c5 = c2;
  yy=1;

 }
 while (c3>=0) {
  d = sqrt(c3);
  xleft = c5-d;
  xright = c5+d;                        
  newpoint(xcen+xleft,ycen+yy);
  newpoint(xcen+xright,ycen+yy);
  newpoint(xcen-xright,ycen-yy);
  newpoint(xcen-xleft,ycen-yy);
  c5+=c2;
  v1+=c6;
  c3-=v1;
  yy=yy+1;
 }
 dir=0; n=0;
 totpts++; /* add another point to join with first */
 pnts = (XPoint *) malloc(totpts * sizeof(XPoint));
 ppnts = pnts;
 /* now go down the 1st column, up the 2nd, down the 4th 
    and up the 3rd to get the points in the correct order */
 for (k=0; k<=3; k++)
     {
     if (dir==0)
  for (m=0; m<nump[k]; m++)
      {
      ppnts->x=x[m][order[k]];
      ppnts->y=y[m][order[k]];
      ppnts++;
      n++;
      }
     else
  for (m=nump[k]-1; m>=0; m--)
      {
      ppnts->x=x[m][order[k]];
      ppnts->y=y[m][order[k]];
      ppnts++;
      n++;
      }
     dir = 1-dir;
     } /* next k */
 /* add another point to join with first */
 ppnts->x = pnts->x;
 ppnts->y = pnts->y;
 n++;
 /* use either DrawLines or FillPolygon */
*** XDrawLines(display, window, gc, pnts, n, CoordModeOrigin);
*** XFillPolygon(display, window, gc, points, npoints,
***                         Convex, CoordModeOrigin);
 free(pnts);
 return;

}

/* store the points across (row-wise in) the matrix */

newpoint(xp,yp)
    float    xp,yp;
{
    x[i][j]=round(xp);
    y[i][j]=round(yp);
    nump[j]++;
    totpts++;
    if (++j > 3)
 {
 j=0; 
 i++;
 }
}