home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch2_6 / conmat.c
Encoding:
C/C++ Source or Header  |  1995-04-04  |  21.3 KB  |  698 lines
  1. /* CONMAT.C - Ellipse tranformation and intersection functions */
  2. /* Written by Kenneth J. Hill, June, 24, 1994  */
  3.  
  4. #include <stdlib.h>
  5. #include <stdio.h>
  6. #include <math.h>
  7. #include <values.h>
  8.  
  9. #ifndef    M_PI
  10. #define M_PI    3.14159265358979323846
  11. #endif
  12.  
  13. #ifndef    M_PI_2
  14. #define M_PI_2    1.57079632679489661923
  15. #endif
  16.  
  17. #ifndef max
  18. #define max(a,b) ((a)>=(b) ? (a) : (b) )
  19. #endif
  20.  
  21. /* Type definitions */
  22.  
  23. /* Point structure */
  24. typedef struct PointTag
  25.    {
  26.    double X,Y;
  27.    } Point;
  28.  
  29. /* Line structure */
  30. typedef struct LineTag
  31.    {
  32.    Point P1,P2;
  33.    } Line;
  34.  
  35. /* Circle structure */
  36. typedef struct CircleTag
  37.    {
  38.    Point Center;
  39.    double Radius;
  40.    } Circle;
  41.  
  42. /* Ellipse structure */
  43. typedef struct EllipseTag
  44.    {
  45.    Point Center;         /* ellipse center     */
  46.    double MaxRad,MinRad;     /* major and minor axis */
  47.    double Phi;             /* major axis rotation  */
  48.    } Ellipse;
  49.  
  50. /* Conic coefficients structure */
  51. typedef struct ConicTag
  52.    {
  53.    double A,B,C,D,E,F;
  54.    } Conic;
  55.  
  56. /* transformation matrix type */
  57. typedef struct TMatTag
  58.    {
  59.    double a,b,c,d;         /* tranformation coefficients */
  60.    double m,n;             /* translation coefficients   */
  61.    } TMat;
  62.  
  63. /* prototypes */
  64.  
  65. /* Functions with names beginning with a O:
  66.  *
  67.  * 1. Ignore ellipse center coordinates.
  68.  *
  69.  * 2. Ignore any translation components in tranformation matrices.
  70.  *
  71.  * 3. Assume that conic coeficients were generated by "O-name" functions.
  72.  *
  73.  * This keeps the computations relatively simple. The ellipse centers
  74.  * can then be transformed separately as points.
  75.  */
  76.  
  77. /* OTransformConic - Transform conic coefficients about the origin */
  78. void OTransformConic(Conic *ConicP,TMat *TMatP);
  79.  
  80. /* OGenEllipseCoefs - Generate conic coefficients of an ellipse */
  81. void OGenEllipseCoefs(Ellipse *EllipseP,Conic *ConicP);
  82.  
  83. /* OGenEllipseGeom - Generates ellipse geometry from conic coefficients */
  84. void OGenEllipseGeom(Conic *ConicP,Ellipse *EllipseP);
  85.  
  86. /* TransformPoint - Transform a point using a tranformation matrix */
  87. void TransformPoint(Point *PointP,TMat *TMatP);
  88.  
  89. /* TransformEllipse - Transform an ellipse using a tranformation matrix */
  90. void TransformEllipse(Ellipse *EllipseP,TMat *TMatP);
  91.  
  92. /* The following functions are defined in cubicand.c, on
  93.    the Graphic Gems I diskette. */
  94. int SolveCubic(double c[4],double s[3]);
  95. int SolveQuadic(double c[3],double dest[2]);
  96.  
  97.  
  98. /* Identity matrix */
  99. static const TMat IdentMat={1.0,0.0,0.0,1.0,0.0,0.0};
  100.  
  101. /* Transformation matrix routines */
  102.  
  103. /* Translate a matrix by m,n */
  104. void TranslateMat(TMat *Mat,double m,double n)
  105.    {
  106.    Mat->m += m;
  107.    Mat->n += n;
  108.    }
  109.  
  110. /* Rotate a matrix by Phi */
  111. void RotateMat(TMat *Mat,double Phi)
  112.    {
  113.    double SinPhi=sin(Phi);
  114.    double CosPhi=cos(Phi);
  115.    TMat temp=*Mat;        /* temporary copy of Mat */
  116.  
  117.    /* These are just the matrix operations written out long hand */
  118.    Mat->a = temp.a*CosPhi - temp.b*SinPhi;
  119.    Mat->b = temp.b*CosPhi + temp.a*SinPhi;
  120.    Mat->c = temp.c*CosPhi - temp.d*SinPhi;
  121.    Mat->d = temp.d*CosPhi + temp.c*SinPhi;
  122.    Mat->m = temp.m*CosPhi - temp.n*SinPhi;
  123.    Mat->n = temp.n*CosPhi + temp.m*SinPhi;
  124.    }
  125.  
  126. /* Scale a matrix by sx, sy */
  127. void ScaleMat(TMat *Mat,double sx,double sy)
  128.    {
  129.    Mat->a *= sx;
  130.    Mat->b *= sy;
  131.    Mat->c *= sx;
  132.    Mat->d *= sy;
  133.    Mat->m *= sx;
  134.    Mat->n *= sy;
  135.    }
  136.  
  137. /* TransformPoint - Transform a point using a tranformation matrix */
  138. void TransformPoint(Point *PointP,TMat *TMatP)
  139.    {
  140.    Point TempPoint;
  141.  
  142.    TempPoint.X = PointP->X*TMatP->a + PointP->Y*TMatP->c + TMatP->m;
  143.    TempPoint.Y = PointP->X*TMatP->b + PointP->Y*TMatP->d + TMatP->n;
  144.    *PointP=TempPoint;
  145.    }
  146.  
  147. /* Conic routines */
  148.  
  149. /* near zero test */
  150. #define EPSILON 1e-9
  151. #define IsZero(x) (x > -EPSILON && x < EPSILON)
  152.  
  153. /* GenEllipseCoefs - Generate conic coefficients of an ellipse */
  154. static void GenEllipseCoefs(Ellipse *elip,Conic *M)
  155.    {
  156.    double sqr_r1,sqr_r2;
  157.    double sint,cost,sin2t,sqr_sint,sqr_cost;
  158.    double cenx,ceny,sqr_cenx,sqr_ceny,invsqr_r1,invsqr_r2;
  159.  
  160.    /* common coeficients */
  161.    sqr_r1 = elip->MaxRad;
  162.    sqr_r2 = elip->MinRad;
  163.    sqr_r1 *= sqr_r1;
  164.    sqr_r2 *= sqr_r2;
  165.    sint = sin(elip->Phi);
  166.    cost = cos(elip->Phi);
  167.    sin2t = 2.0*sint*cost;
  168.    sqr_sint = sint*sint;
  169.    sqr_cost = cost*cost;
  170.    cenx = elip->Center.X;
  171.    sqr_cenx = cenx*cenx;
  172.    ceny = elip->Center.Y;
  173.    sqr_ceny = ceny*ceny;
  174.    invsqr_r1 = 1.0/sqr_r1;
  175.    invsqr_r2 = 1.0/sqr_r2;
  176.  
  177.    /* Compute the coefficients. These formulae are the transformations
  178.       on the unit circle written out long hand */
  179.    M->A = sqr_cost/sqr_r1 + sqr_sint/sqr_r2;
  180.    M->B = (sqr_r2-sqr_r1)*sin2t/(2.0*sqr_r1*sqr_r2);
  181.    M->C = sqr_cost/sqr_r2 + sqr_sint/sqr_r1;
  182.    M->D = -ceny*M->B-cenx*M->A;
  183.    M->E = -cenx*M->B-ceny*M->C;
  184.    M->F = -1.0 + (sqr_cenx + sqr_ceny)*(invsqr_r1 + invsqr_r2)/2.0 +
  185.       (sqr_cost - sqr_sint)*(sqr_cenx - sqr_ceny)*(invsqr_r1 - invsqr_r2)/2.0 +
  186.       cenx*ceny*(invsqr_r1 - invsqr_r2)*sin2t;
  187.    }
  188.  
  189. /* Compute the transformation which turns an ellipse into a circle */
  190. void Elp2Cir(Ellipse *Elp,TMat *CirMat)
  191.    {
  192.    /* Start with identity matrix */
  193.    *CirMat = IdentMat;
  194.    /* Translate to origin */
  195.    TranslateMat(CirMat,-Elp->Center.X,-Elp->Center.Y);
  196.    /* Rotate into standard position */
  197.    RotateMat(CirMat,-Elp->Phi);
  198.    /* Scale into a circle. */
  199.    ScaleMat(CirMat,1.0/Elp->MaxRad,1.0/Elp->MinRad);
  200.    }
  201.  
  202. /* Compute the inverse of the transformation
  203.    which turns an ellipse into a circle */
  204. void InvElp2Cir(Ellipse *Elp,TMat *InvMat)
  205.    {
  206.    /* Start with identity matrix */
  207.    *InvMat = IdentMat;
  208.    /* Scale back into an ellipse. */
  209.    ScaleMat(InvMat,Elp->MaxRad,Elp->MinRad);
  210.    /* Rotate */
  211.    RotateMat(InvMat,Elp->Phi);
  212.    /* Translate from origin */
  213.    TranslateMat(InvMat,Elp->Center.X,Elp->Center.Y);
  214.    }
  215.  
  216. /* OTransformConic - Transform conic coefficients about the origin     */
  217. /* This routine ignores the translation components of *TMatP and
  218.    assumes the conic is "centered" at the origin (i.e. D,E=0, F=-1)     */
  219. /* The computations are just the matrix operations written out long hand */
  220. /* This code assumes that the transformation is not degenerate         */
  221. void OTransformConic(Conic *ConicP,TMat *TMatP)
  222.    {
  223.    double A,B,C,Denom;
  224.  
  225.    /* common denominator for transformed cooefficients */
  226.    Denom = TMatP->a*TMatP->d - TMatP->b*TMatP->c;
  227.    Denom *= Denom;
  228.  
  229.    A = (ConicP->C*TMatP->b*TMatP->b - 2.0*ConicP->B*TMatP->b*TMatP->d +
  230.       ConicP->A*TMatP->d*TMatP->d)/Denom;
  231.  
  232.    B = (-ConicP->C*TMatP->a*TMatP->b + ConicP->B*TMatP->b*TMatP->c +
  233.       ConicP->B*TMatP->a*TMatP->d - ConicP->A*TMatP->c*TMatP->d)/Denom;
  234.  
  235.    C = (ConicP->C*TMatP->a*TMatP->a - 2.0*ConicP->B*TMatP->a*TMatP->c +
  236.       ConicP->A*TMatP->c*TMatP->c)/Denom;
  237.  
  238.    ConicP->A=A;
  239.    ConicP->B=B;
  240.    ConicP->C=C;
  241.    }
  242.  
  243. /* OGenEllipseCoefs - Generate conic coefficients of an ellipse */
  244. /* The ellipse is assumed to be centered at the origin. */
  245. void OGenEllipseCoefs(Ellipse *EllipseP,Conic *ConicP)
  246.    {
  247.    double SinPhi = sin(EllipseP->Phi); /* sine of ellipse rotation   */
  248.    double CosPhi = cos(EllipseP->Phi); /* cosine of ellipse rotation */
  249.    double SqSinPhi = SinPhi*SinPhi;    /* square of sin(phi)         */
  250.    double SqCosPhi = CosPhi*CosPhi;    /* square of cos(phi)         */
  251.    double SqMaxRad = EllipseP->MaxRad*EllipseP->MaxRad;
  252.    double SqMinRad = EllipseP->MinRad*EllipseP->MinRad;
  253.  
  254.    /* compute coefficients for the ellipse in standard position */
  255.    ConicP->A = SqCosPhi/SqMaxRad + SqSinPhi/SqMinRad;
  256.    ConicP->B = (1.0/SqMaxRad - 1.0/SqMinRad)*SinPhi*CosPhi;
  257.    ConicP->C = SqCosPhi/SqMinRad + SqSinPhi/SqMaxRad;
  258.    ConicP->D = ConicP->E = 0.0;
  259.    ConicP->F = -1.0;
  260.    }
  261.  
  262. /* OGenEllipseGeom - Generates ellipse geometry from conic coefficients */
  263. /* This routine assumes the conic coefficients D=E=0, F=-1 */
  264. void OGenEllipseGeom(Conic *ConicP,Ellipse *EllipseP)
  265.    {
  266.    double Numer,Denom,Temp;
  267.    TMat DiagTransform;         /* transform diagonalization */
  268.    Conic ConicT = *ConicP;     /* temporary copy of conic coefficients */
  269.    double SinPhi,CosPhi;
  270.  
  271.    /* compute new ellipse rotation */
  272.    Numer = ConicT.B + ConicT.B;
  273.    Denom = ConicT.A - ConicT.C;
  274.    /* Phi = 1/2 atan(Numer/Denom) = 1/2 (pi/2 - atan(Denom/Numer)
  275.       We use the form that keeps the argument to atan between -1 and 1 */
  276.  
  277.    EllipseP->Phi = 0.5*(fabs(Numer) < fabs(Denom)?
  278.       atan(Numer/Denom):M_PI_2-atan(Denom/Numer));
  279.  
  280.    /* diagonalize the conic */
  281.    SinPhi = sin(EllipseP->Phi);
  282.    CosPhi = cos(EllipseP->Phi);
  283.    DiagTransform.a = CosPhi;     /* rotate by -Phi */
  284.    DiagTransform.b = -SinPhi;
  285.    DiagTransform.c = SinPhi;
  286.    DiagTransform.d = CosPhi;
  287.    DiagTransform.m = DiagTransform.n = 0.0;
  288.    OTransformConic(&ConicT,&DiagTransform);
  289.  
  290.    /* compute new radii from diagonalized coefficients */
  291.    EllipseP->MaxRad = 1.0/sqrt(ConicT.A);
  292.    EllipseP->MinRad = 1.0/sqrt(ConicT.C);
  293.  
  294.    /* be sure MaxRad >= MinRad */
  295.    if (EllipseP->MaxRad < EllipseP->MinRad)
  296.       {
  297.       Temp = EllipseP->MaxRad;     /* exchange the radii */
  298.       EllipseP->MaxRad = EllipseP->MinRad;
  299.       EllipseP->MinRad = Temp;
  300.       EllipseP->Phi += M_PI_2;     /* adjust the rotation */
  301.       }
  302.    }
  303.  
  304. /* TransformEllipse - Transform an ellipse using a tranformation matrix */
  305. void TransformEllipse(Ellipse *EllipseP,TMat *TMatP)
  306.    {
  307.    Conic EllipseCoefs;
  308.  
  309.    /* generate the ellipse coefficients (using Center=origin) */
  310.    OGenEllipseCoefs(EllipseP,&EllipseCoefs);
  311.  
  312.    /* transform the coefficients */
  313.    OTransformConic(&EllipseCoefs,TMatP);
  314.  
  315.    /* turn the transformed coefficients back into geometry */
  316.    OGenEllipseGeom(&EllipseCoefs,EllipseP);
  317.  
  318.    /* translate the center */
  319.    TransformPoint(&EllipseP->Center,TMatP);
  320.    }
  321.  
  322. /* MultMat3 - Multiply two 3x3 matrices */
  323. void MultMat3(double *Mat1,double *Mat2, double *Result)
  324.    {
  325.    int i,j;
  326.  
  327.    for (i = 0;i < 3;i++)
  328.       for (j = 0;j < 3;j++)
  329.      Result[i*3+j] = Mat1[i*3+0]*Mat2[0*3+j] +
  330.         Mat1[i*3+1]*Mat2[1*3+j] +
  331.         Mat1[i*3+2]*Mat2[2*3+j];
  332.    }
  333.  
  334. /* Transform a conic by a tranformation matrix */
  335. void TransformConic(Conic *ConicP,TMat *TMatP)
  336.    {
  337.    double InvMat[3][3],ConMat[3][3],TranInvMat[3][3];
  338.    double Result1[3][3],Result2[3][3];
  339.    double D;
  340.    int i,j;
  341.  
  342.    /* Compute M' = Inv(TMat).M.Transpose(Inv(TMat))
  343.  
  344.    /* compute the transformation using matrix muliplication */
  345.    ConMat[0][0] = ConicP->A;
  346.    ConMat[0][1] = ConMat[1][0] = ConicP->B;
  347.    ConMat[1][1] = ConicP->C;
  348.    ConMat[0][2] = ConMat[2][0] = ConicP->D;
  349.    ConMat[1][2] = ConMat[2][1] = ConicP->E;
  350.    ConMat[2][2] = ConicP->F;
  351.  
  352.    /* inverse transformation */
  353.    D = TMatP->a*TMatP->d - TMatP->b*TMatP->c;
  354.    InvMat[0][0] = TMatP->d/D;
  355.    InvMat[0][1] = -TMatP->b/D;
  356.    InvMat[0][2] = 0.0;
  357.    InvMat[1][0] = -TMatP->c/D;
  358.    InvMat[1][1] = TMatP->a/D;
  359.    InvMat[1][2] = 0.0;
  360.    InvMat[2][0] = (TMatP->c*TMatP->n - TMatP->d*TMatP->m)/D;
  361.    InvMat[2][1] = (TMatP->b*TMatP->m - TMatP->a*TMatP->n)/D;
  362.    InvMat[2][2] = 1.0;
  363.  
  364.    /* compute transpose */
  365.    for (i = 0;i < 3;i++)
  366.       for (j = 0;j < 3;j++)
  367.      TranInvMat[j][i] = InvMat[i][j];
  368.  
  369.    /* multiply the matrices */
  370.    MultMat3((double *)InvMat,(double *)ConMat,(double *)Result1);
  371.    MultMat3((double *)Result1,(double *)TranInvMat,(double *)Result2);
  372.    ConicP->A = Result2[0][0];           /* return to conic form */
  373.    ConicP->B = Result2[0][1];
  374.    ConicP->C = Result2[1][1];
  375.    ConicP->D = Result2[0][2];
  376.    ConicP->E = Result2[1][2];
  377.    ConicP->F = Result2[2][2];
  378.    }
  379.  
  380. /* Compute the intersection of a circle and a line */
  381. /* See Graphic Gems Volume 1, page 5 for a description of this algorithm */
  382. int IntCirLine(Point *IntPts,Circle *Cir,Line *Ln)
  383.    {
  384.    Point G,V;
  385.    double a,b,c,d,t,sqrt_d;
  386.  
  387.    G.X = Ln->P1.X - Cir->Center.X;         /* G = Ln->P1 - Cir->Center */
  388.    G.Y = Ln->P1.Y - Cir->Center.Y;
  389.    V.X = Ln->P2.X - Ln->P1.X;           /* V = Ln->P2 - Ln->P1 */
  390.    V.Y = Ln->P2.Y - Ln->P1.Y;
  391.    a = V.X*V.X + V.Y*V.Y;        /* a = V.V */
  392.    b = V.X*G.X + V.Y*G.Y;b += b;     /* b = 2(V.G) */
  393.    c = (G.X*G.X + G.Y*G.Y) -        /* c = G.G + Circle->Radius^2 */
  394.       Cir->Radius*Cir->Radius;
  395.    d = b*b - 4.0*a*c;            /* discriminant */
  396.  
  397.    if (d <= 0.0)
  398.       return 0;                /* no intersections */
  399.  
  400.    sqrt_d = sqrt(d);
  401.    t = (-b + sqrt_d)/(a + a);           /* t = (-b +/- sqrt(d))/2a */
  402.    IntPts[0].X = Ln->P1.X + t*V.X;      /* Pt = Ln->P1 + t V */
  403.    IntPts[0].Y = Ln->P1.Y + t*V.Y;
  404.    t = (-b - sqrt_d)/(a + a);
  405.    IntPts[1].X = Ln->P1.X + t*V.X;
  406.    IntPts[1].Y = Ln->P1.Y + t*V.Y;
  407.    return 2;
  408.    }
  409.  
  410. /* compute all intersections of two ellipses */
  411. /* E1 and E2 are the two ellipses */
  412. /* IntPts points to an array of twelve points
  413.     (some duplicates may be returned) */
  414. /* The number of intersections found is returned */
  415. /* Both ellipses are assumed to have non-zero radii */
  416. int Int2Elip(Point *IntPts,Ellipse *E1,Ellipse *E2)
  417.    {
  418.    TMat ElpCirMat1,ElpCirMat2,InvMat,TempMat;
  419.    Conic Conic1,Conic2,Conic3,TempConic;
  420.    double Roots[3],qRoots[2];
  421.    static Circle TestCir = {{0.0,0.0},1.0};
  422.    Line TestLine[2];
  423.    Point TestPoint;
  424.    double PolyCoef[4];        /* coefficients of the polynomial */
  425.    double D;            /* discriminant: B^2 - AC */
  426.    double Phi;            /* ellipse rotation */
  427.    double m,n;            /* ellipse translation */
  428.    double Scl;            /* scaling factor */
  429.    int NumRoots,NumLines;
  430.    int CircleInts;        /* intersections between line and circle */
  431.    int IntCount = 0;        /* number of intersections found */
  432.    int i,j,k;
  433.  
  434.    /* compute the transformations which turn E1 and E2 into circles */
  435.    Elp2Cir(E1,&ElpCirMat1);
  436.    Elp2Cir(E2,&ElpCirMat2);
  437.  
  438.    /* compute the inverse transformation of ElpCirMat1 */
  439.    InvElp2Cir(E1,&InvMat);
  440.  
  441.    /* Compute the characteristic matrices */
  442.    GenEllipseCoefs(E1,&Conic1);
  443.    GenEllipseCoefs(E2,&Conic2);
  444.  
  445.    /* Find x such that Det(Conic1 + x Conic2) = 0 */
  446.    PolyCoef[0] = -Conic1.C*Conic1.D*Conic1.D + 2.0*Conic1.B*Conic1.D*Conic1.E -
  447.       Conic1.A*Conic1.E*Conic1.E - Conic1.B*Conic1.B*Conic1.F +
  448.       Conic1.A*Conic1.C*Conic1.F;
  449.    PolyCoef[1] = -(Conic2.C*Conic1.D*Conic1.D) -
  450.       2.0*Conic1.C*Conic1.D*Conic2.D + 2.0*Conic2.B*Conic1.D*Conic1.E +
  451.       2.0*Conic1.B*Conic2.D*Conic1.E - Conic2.A*Conic1.E*Conic1.E +
  452.       2.0*Conic1.B*Conic1.D*Conic2.E - 2.0*Conic1.A*Conic1.E*Conic2.E -
  453.       2.0*Conic1.B*Conic2.B*Conic1.F + Conic2.A*Conic1.C*Conic1.F +
  454.       Conic1.A*Conic2.C*Conic1.F - Conic1.B*Conic1.B*Conic2.F +
  455.       Conic1.A*Conic1.C*Conic2.F;
  456.    PolyCoef[2] = -2.0*Conic2.C*Conic1.D*Conic2.D - Conic1.C*Conic2.D*Conic2.D +
  457.       2.0*Conic2.B*Conic2.D*Conic1.E + 2.0*Conic2.B*Conic1.D*Conic2.E +
  458.       2.0*Conic1.B*Conic2.D*Conic2.E - 2.0*Conic2.A*Conic1.E*Conic2.E -
  459.       Conic1.A*Conic2.E*Conic2.E - Conic2.B*Conic2.B*Conic1.F +
  460.       Conic2.A*Conic2.C*Conic1.F - 2.0*Conic1.B*Conic2.B*Conic2.F +
  461.       Conic2.A*Conic1.C*Conic2.F + Conic1.A*Conic2.C*Conic2.F;
  462.    PolyCoef[3] = -Conic2.C*Conic2.D*Conic2.D + 2.0*Conic2.B*Conic2.D*Conic2.E -
  463.       Conic2.A*Conic2.E*Conic2.E - Conic2.B*Conic2.B*Conic2.F +
  464.       Conic2.A*Conic2.C*Conic2.F;
  465.    NumRoots = SolveCubic(PolyCoef,Roots);
  466.  
  467.    if (NumRoots == 0)
  468.       return 0;
  469.  
  470.    /* we try all the roots, even though it's redundant, so that we
  471.       avoid some pathological situations */
  472.    for (i=0;i<NumRoots;i++)
  473.       {
  474.       NumLines = 0;
  475.  
  476.       /* Conic3 = Conic1 + mu Conic2 */
  477.       Conic3.A = Conic1.A + Roots[i]*Conic2.A;
  478.       Conic3.B = Conic1.B + Roots[i]*Conic2.B;
  479.       Conic3.C = Conic1.C + Roots[i]*Conic2.C;
  480.       Conic3.D = Conic1.D + Roots[i]*Conic2.D;
  481.       Conic3.E = Conic1.E + Roots[i]*Conic2.E;
  482.       Conic3.F = Conic1.F + Roots[i]*Conic2.F;
  483.  
  484.       D = Conic3.B*Conic3.B - Conic3.A*Conic3.C;
  485.       if (IsZero(Conic3.A) && IsZero(Conic3.B) && IsZero(Conic3.C))
  486.      {
  487.      /* Case 1 - Single line */
  488.      NumLines = 1;
  489.      /* compute endpoints of the line, avoiding division by zero */
  490.      if (fabs(Conic3.D) > fabs(Conic3.E))
  491.         {
  492.         TestLine[0].P1.Y = 0.0;
  493.         TestLine[0].P1.X = -Conic3.F/(Conic3.D + Conic3.D);
  494.         TestLine[0].P2.Y = 1.0;
  495.         TestLine[0].P2.X = -(Conic3.E + Conic3.E + Conic3.F)/
  496.            (Conic3.D + Conic3.D);
  497.         }
  498.      else
  499.         {
  500.         TestLine[0].P1.X = 0.0;
  501.         TestLine[0].P1.Y = -Conic3.F/(Conic3.E + Conic3.E);
  502.         TestLine[0].P2.X = 1.0;
  503.         TestLine[0].P2.X = -(Conic3.D + Conic3.D + Conic3.F)/
  504.            (Conic3.E + Conic3.E);
  505.         }
  506.      }
  507.       else
  508.      {
  509.      /* use the espresion for Phi that takes atan of the
  510.         smallest argument */
  511.      Phi = (fabs(Conic3.B + Conic3.B) < fabs(Conic3.A-Conic3.C)?
  512.         atan((Conic3.B + Conic3.B)/(Conic3.A - Conic3.C)):
  513.         M_PI_2 - atan((Conic3.A - Conic3.C)/(Conic3.B + Conic3.B)))/2.0;
  514.      if (IsZero(D))
  515.         {
  516.         /* Case 2 - Parallel lines */
  517.         TempConic = Conic3;
  518.         TempMat = IdentMat;
  519.         RotateMat(&TempMat,-Phi);
  520.         TransformConic(&TempConic,&TempMat);
  521.         if (IsZero(TempConic.C))   /* vertical */
  522.            {
  523.            PolyCoef[0] = TempConic.F;
  524.            PolyCoef[1] = TempConic.D;
  525.            PolyCoef[2] = TempConic.A;
  526.            if ((NumLines=SolveQuadic(PolyCoef,qRoots))!=0)
  527.           {
  528.           TestLine[0].P1.X = qRoots[0];
  529.           TestLine[0].P1.Y = -1.0;
  530.           TestLine[0].P2.X = qRoots[0];
  531.           TestLine[0].P2.Y = 1.0;
  532.           if (NumLines==2)
  533.              {
  534.              TestLine[1].P1.X = qRoots[1];
  535.              TestLine[1].P1.Y = -1.0;
  536.              TestLine[1].P2.X = qRoots[1];
  537.              TestLine[1].P2.Y = 1.0;
  538.              }
  539.           }
  540.            }
  541.         else            /* horizontal */
  542.            {
  543.            PolyCoef[0] = TempConic.F;
  544.            PolyCoef[1] = TempConic.E;
  545.            PolyCoef[2] = TempConic.C;
  546.            if ((NumLines=SolveQuadic(PolyCoef,qRoots))!=0)
  547.           {
  548.           TestLine[0].P1.X = -1.0;
  549.           TestLine[0].P1.Y = qRoots[0];
  550.           TestLine[0].P2.X = 1.0;
  551.           TestLine[0].P2.Y = qRoots[0];
  552.           if (NumLines==2)
  553.              {
  554.              TestLine[1].P1.X = -1.0;
  555.              TestLine[1].P1.Y = qRoots[1];
  556.              TestLine[1].P2.X = 1.0;
  557.              TestLine[1].P2.Y = qRoots[1];
  558.              }
  559.           }
  560.            }
  561.         TempMat = IdentMat;
  562.         RotateMat(&TempMat,Phi);
  563.         TransformPoint(&TestLine[0].P1,&TempMat);
  564.         TransformPoint(&TestLine[0].P2,&TempMat);
  565.         if (NumLines==2)
  566.            {
  567.            TransformPoint(&TestLine[1].P1,&TempMat);
  568.            TransformPoint(&TestLine[1].P2,&TempMat);
  569.            }
  570.         }
  571.      else
  572.         {
  573.         /* Case 3 - Crossing lines */
  574.         NumLines = 2;
  575.  
  576.         /* translate the system so that the intersection of the lines
  577.            is at the origin */
  578.         TempConic = Conic3;
  579.         m = (Conic3.C*Conic3.D - Conic3.B*Conic3.E)/D;
  580.         n = (Conic3.A*Conic3.E - Conic3.B*Conic3.D)/D;
  581.         TempMat = IdentMat;
  582.         TranslateMat(&TempMat,-m,-n);
  583.         RotateMat(&TempMat,-Phi);
  584.         TransformConic(&TempConic,&TempMat);
  585.  
  586.         /* Compute the line endpoints */
  587.         TestLine[0].P1.X = sqrt(fabs(1.0/TempConic.A));
  588.         TestLine[0].P1.Y = sqrt(fabs(1.0/TempConic.C));
  589.         Scl = max(TestLine[0].P1.X,TestLine[0].P1.Y);  /* adjust range */
  590.         TestLine[0].P1.X /= Scl;
  591.         TestLine[0].P1.Y /= Scl;
  592.         TestLine[0].P2.X = - TestLine[0].P1.X;
  593.         TestLine[0].P2.Y = - TestLine[0].P1.Y;
  594.         TestLine[1].P1.X = TestLine[0].P1.X;
  595.         TestLine[1].P1.Y = - TestLine[0].P1.Y;
  596.         TestLine[1].P2.X = - TestLine[1].P1.X;
  597.         TestLine[1].P2.Y = - TestLine[1].P1.Y;
  598.  
  599.         /* translate the lines back */
  600.         TempMat = IdentMat;
  601.         RotateMat(&TempMat,Phi);
  602.         TranslateMat(&TempMat,m,n);
  603.         TransformPoint(&TestLine[0].P1,&TempMat);
  604.         TransformPoint(&TestLine[0].P2,&TempMat);
  605.         TransformPoint(&TestLine[1].P1,&TempMat);
  606.         TransformPoint(&TestLine[1].P2,&TempMat);
  607.         }
  608.      }
  609.  
  610.       /* find the ellipse line intersections */
  611.       for (j = 0;j < NumLines;j++)
  612.      {
  613.      /* transform the line endpts into the circle space of the ellipse */
  614.      TransformPoint(&TestLine[j].P1,&ElpCirMat1);
  615.      TransformPoint(&TestLine[j].P2,&ElpCirMat1);
  616.  
  617.      /* compute the number of intersections of the transformed line
  618.         and test circle */
  619.      CircleInts = IntCirLine(&IntPts[IntCount],&TestCir,&TestLine[j]);
  620.      if (CircleInts>0)
  621.         {
  622.         /* transform the intersection points back into ellipse space */
  623.         for (k = 0;k < CircleInts;k++)
  624.            TransformPoint(&IntPts[IntCount+k],&InvMat);
  625.         /* update the number of intersections found */
  626.         IntCount += CircleInts;
  627.         }
  628.      }
  629.       }
  630.    /* validate the points */
  631.    j = IntCount;
  632.    IntCount = 0;
  633.    for (i = 0;i < j;i++)
  634.       {
  635.       TestPoint = IntPts[i];
  636.       TransformPoint(&TestPoint,&ElpCirMat2);
  637.       if (TestPoint.X < 2.0 && TestPoint.Y < 2.0 &&
  638.      IsZero(1.0 - sqrt(TestPoint.X*TestPoint.X +
  639.      TestPoint.Y*TestPoint.Y)))
  640.      IntPts[IntCount++]=IntPts[i];
  641.       }
  642.    return IntCount;
  643.    }
  644.  
  645. /* Test routines */
  646.  
  647. /* Ellipse with center at (1,2), major radius 2, minor radius 1,
  648.    and rotation 0. */
  649. Ellipse TestEllipse={{2.0,1.0},2.0,1.0,0.0};
  650.  
  651. /* Transform matrix for shear of 45 degrees from vertical */
  652. TMat TestTransform={1.0,0.0,1.0,1.0,0.0,0.0};
  653.  
  654. /* Display an ellipse. This version lists the structure values. */
  655. void DisplayEllipse(Ellipse *EllipseP)
  656.    {
  657.    printf("\tCenter at (%6.3f,%6.3f)\n"
  658.       "\tMajor radius: %6.3f\n"
  659.       "\tMinor radius: %6.3f\n"
  660.       "\tRotation: %6.3f degrees\n\n",
  661.       EllipseP->Center.X,EllipseP->Center.Y,
  662.       EllipseP->MaxRad,EllipseP->MinRad,
  663.       180.0*EllipseP->Phi/M_PI);
  664.    }
  665.  
  666.  
  667. /* test ellipses for intersection */
  668. Ellipse Elp1={{5.0,4.0},1.0,0.5,M_PI/3.0};
  669. Ellipse Elp2={{4.0,3.0},2.0,1.0,0.0};
  670. Ellipse Elp3={{1.0,1.0},2.0,1.0,M_PI_2};
  671. Ellipse Elp4={{1.0,1.0},2.0,0.5,0.0};
  672.  
  673. void main(void)
  674.    {
  675.    Point IntPts[12];
  676.    int IntCount;
  677.    int i;
  678.  
  679.    /* find ellipse intersections */
  680.    printf("Intersections of ellipses 1 & 2:\n\n");
  681.    IntCount=Int2Elip(IntPts,&Elp1,&Elp2);
  682.    for (i = 0;i < IntCount;i++)
  683.       printf("   %f, %f\n",IntPts[i].X,IntPts[i].Y);
  684.    printf("Intersections of ellipses 3 & 4:\n");
  685.    IntCount=Int2Elip(IntPts,&Elp3,&Elp4);
  686.    for (i = 0;i < IntCount;i++)
  687.       printf("   %f, %f\n",IntPts[i].X,IntPts[i].Y);
  688.  
  689.    /* transform ellipses */
  690.    printf("\n\nBefore transformation:\n");
  691.    DisplayEllipse(&TestEllipse);
  692.  
  693.    TransformEllipse(&TestEllipse,&TestTransform);
  694.  
  695.    printf("After transformation:\n");
  696.    DisplayEllipse(&TestEllipse);
  697.    }
  698.