/****************************************************************************** * Cagd_Cnc.c - Curve representation of arcs and circles * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Jun. 90. * ******************************************************************************/ #include #include #include #include "irit_sm.h" #include "extra_fn.h" #include "cagd_loc.h" #define CONIC_MAX_VAL -10 #define CONIC_MIN_VAL 10 #define CONIC_DELTA_VAL (CONIC_MAX_VAL - CONIC_MIN_VAL) #define CAGD_QUADRIC_EXTENT (CONIC_MAX_VAL) #define QUADRIC_INVERT(x) ((x) == 0 ? IRIT_INFNTY : 1.0 / sqrt(x)) /***************************************************************************** * DESCRIPTION: M * Construct rational quadratic curve out of the 6 coefficients of the M * conic section: A x^2 + B xy + C y^2 + D x + E y + F = 0. M * Based on: M * Bezier Curves and Surface Patches on Quadrics, by Josef Hoschek, M * Mathematical methods in Computer aided Geometric Design II, Tom Lyche and M * Larry L. Schumaker (eds.), pp 331-342, 1992. M * * * PARAMETERS: M * A, B, C, D, E, F: The six coefficients of the conic curve. M * ZLevel: Sets the Z level of this XY parallel conic curve. M * * * RETURN VALUE: M * CagdCrvStruct *: A quadratic curve representing the conic. M * * * SEE ALSO: M * BspCrvCreateCircle, BspCrvCreateUnitCircle, CagdCrvCreateArc, M * BzrCrvCreateArc, CagdCreateConicCurve2, CagdCreateQuadricSrf M * * * KEYWORDS: M * CagdCreateConicCurve M *****************************************************************************/ CagdCrvStruct *CagdCreateConicCurve(CagdRType A, CagdRType B, CagdRType C, CagdRType D, CagdRType E, CagdRType F, CagdRType ZLevel) { VectorType Trans; CagdRType **Points, RotAngle = APX_EQ(B, 0.0) ? 0.0 : atan2(B, A - C) * 0.5, A1 = ((A + C) + B * sin(2 * RotAngle) + (A - C) * cos(2 * RotAngle)) * 0.5, B1 = B * cos(2 * RotAngle) - (A - C) * sin(2 * RotAngle), C1 = ((A + C) - B * sin(2 * RotAngle) - (A - C) * cos(2 * RotAngle)) * 0.5, D1 = D * cos(RotAngle) + E * sin(RotAngle), E1 = -D * sin(RotAngle) + E * cos(RotAngle), F1 = F; MatrixType Mat; CagdCrvStruct *PwrCrv, *TCrv, *Crv = NULL; if (!APX_EQ(B1, 0.0) || (APX_EQ(A1, 0.0) && APX_EQ(C1, 0.0))) { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } if (APX_EQ(A1, 0.0) || APX_EQ(C1, 0.0)) { /* It is a parabola. */ return CagdCreateConicCurve2(A, B, C, D, E, F, ZLevel, NULL, NULL); } /* We now have conic: A1 x^2 + C1 y^2 + D1 x + E1 y + F1 = 0. */ /* Compute the translation factors by converting to, */ /* A1(x + D1/(2A1))^2 + C1(y + E1/(2C1))^2 + */ /* F1 - D1^2/4A1 - E1^2/4C1 = 0. */ Trans[0] = -D1 / (2.0 * A1); Trans[1] = -E1 / (2.0 * C1); Trans[2] = ZLevel; F1 -= (C1 * SQR(D1) + A1 * SQR(E1)) / (4 * A1 * C1); if (A1 < 0.0) { /* Make sure A1 is positive. */ A1 = -A1; C1 = -C1; F1 = -F1; } if (A1 * C1 > 0.0) { /* Conic is an ellipse A1 x^2 + C1 y^2 = -F1, A1 > 0, C1 > 0. */ if (F1 <= 0) { CagdRType R[2]; Crv = BspCrvCreateUnitCircle(); R[0] = sqrt(-F1 / A1); R[1] = sqrt(-F1 / C1); /* Scale the circular shape into an ellipse. */ MatGenMatScale(R[0], R[1], 1.0, Mat); TCrv = CagdCrvMatTransform(Crv, Mat); CagdCrvFree(Crv); Crv = TCrv; } else { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } } else { CagdBType SwapXY; if (F1 > 0.0) { SWAP(CagdRType, A1, C1); A1 = -A1; C1 = -C1; F1 = -F1; SwapXY = TRUE; } else SwapXY = FALSE; C1 = -C1; /* Conic is now A1 x^2 - C1 y^2 = -F1, A1 > 0, C1 > 0, -F1 > 0. */ /* The above hyperbola can be written as a parametrization of, */ /* x = (1/w) * (a^2 + b^2) / sqrt(A1), */ /* y = (1/w) * 2ab / sqrt(C1), */ /* and w equals w = (a^2 - b^2) / sqrt(-F1). */ /* The only (large) remaining question is how to parametrize a, b. */ /* For no better selection, we select a = t, b = 1 - t, or */ /* x = (1/w) * (1 - 2 * t + 2 * t^2) / sqrt(A1), */ /* y = (1/w) * (2 * t - 2 * t^2) / sqrt(C1), */ /* and w equals w = (-1 + 2 * t) / sqrt(-F1). */ A1 = 1.0 / sqrt(A1); C1 = 1.0 / sqrt(C1); F1 = F1 == 0.0 ? IRIT_INFNTY : 1.0 / sqrt(-F1); PwrCrv = CagdCrvNew(CAGD_CPOWER_TYPE, CAGD_PT_P2_TYPE, 3); Points = PwrCrv -> Points; Points[0][0] = -F1; /* The w term. */ Points[0][1] = 2.0 * F1; Points[0][2] = 0.0; Points[1][0] = A1; /* The x term. */ Points[1][1] = -2.0 * A1; Points[1][2] = 2.0 * A1; Points[2][0] = 0.0; /* The y term. */ Points[2][1] = 2.0 * C1; Points[2][2] = -2.0 * C1; if (SwapXY) { SWAP(CagdRType, Points[1][0], Points[2][0]); SWAP(CagdRType, Points[1][1], Points[2][1]); SWAP(CagdRType, Points[1][2], Points[2][2]); } Crv = CnvrtPower2BezierCrv(PwrCrv); CagdCrvFree(PwrCrv); } CagdCrvTransform(Crv, Trans, 1.0); MatGenMatRotZ1(RotAngle, Mat); TCrv = CagdCrvMatTransform(Crv, Mat); CagdCrvFree(Crv); Crv = TCrv; CAGD_SET_GEOM_TYPE(Crv, CAGD_GEOM_CONIC_SEC); return Crv; } /***************************************************************************** * DESCRIPTION: M * Construct rational quadratic curve out of the 6 coefficients of the M * conic section: A x^2 + B xy + C y^2 + D x + E y + F = 0. M * * * PARAMETERS: M * A, B, C, D, E, F: The six coefficients of the conic curve. M * ZLevel: Sets the Z level of this XY parallel conic curve. M * PStartXY, PEndXY: Domain of conic section - starting/end points, in M * the XY plane. If NULL, the most complete conic M * possible is created. M * * * RETURN VALUE: M * CagdCrvStruct *: A quadratic curve representing the conic. M * * * SEE ALSO: M * BspCrvCreateCircle, BspCrvCreateUnitCircle, CagdCrvCreateArc, M * BzrCrvCreateArc, CagdCreateConicCurve, CagdCreateQuadricSrf M * * * KEYWORDS: M * CagdCreateConicCurve2 M *****************************************************************************/ CagdCrvStruct *CagdCreateConicCurve2(CagdRType A, CagdRType B, CagdRType C, CagdRType D, CagdRType E, CagdRType F, CagdRType ZLevel, CagdRType *PStartXY, CagdRType *PEndXY) { VectorType Trans; CagdRType Desc, x, x1, x2, dx, y1, y2, dy, RotAngle = APX_EQ(B, 0.0) ? 0.0 : atan2(B, A - C) * 0.5, A1 = ((A + C) + B * sin(2 * RotAngle) + (A - C) * cos(2 * RotAngle)) * 0.5, B1 = B * cos(2 * RotAngle) - (A - C) * sin(2 * RotAngle), C1 = ((A + C) - B * sin(2 * RotAngle) - (A - C) * cos(2 * RotAngle)) * 0.5, D1 = D * cos(RotAngle) + E * sin(RotAngle), E1 = -D * sin(RotAngle) + E * cos(RotAngle), F1 = F; MatrixType Mat; CagdCrvStruct *TCrv, *Crv = NULL; if (!APX_EQ(B1, 0.0) || (APX_EQ(A1, 0.0) && APX_EQ(C1, 0.0))) { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } /* Computes the translation factors. */ Trans[0] = Trans[1] = 0.0; Trans[2] = ZLevel; /* We now have conic: A1 x^2 + C1 y^2 + D1 x + E1 y + F1 = 0. */ Desc = A1 * C1; if (APX_EQ_EPS(Desc, 0.0, IRIT_UEPS)) { CagdCrvStruct *PwrCrv = CagdCrvNew(CAGD_CPOWER_TYPE, CAGD_PT_E3_TYPE, 3); CagdRType **Points= PwrCrv -> Points; /* Its a parabola. */ PwrCrv -> Order = PwrCrv -> Length = 3; ZAP_MEM(PwrCrv -> Points[1], 3 * sizeof(CagdRType)); ZAP_MEM(PwrCrv -> Points[2], 3 * sizeof(CagdRType)); ZAP_MEM(PwrCrv -> Points[3], 3 * sizeof(CagdRType)); if (!APX_EQ(A1, 0.0) && !APX_EQ(E1, 0.0)) { /* Update the translation factors. */ Trans[0] = -D1 / (2.0 * A1); /* Update the control points. */ if (PStartXY != NULL && PEndXY != NULL) { x1 = PStartXY[0] - Trans[0]; x2 = PEndXY[0] - Trans[0]; y1 = PStartXY[1] - Trans[1]; y2 = PEndXY[1] - Trans[1]; dx = x2 - x1; dy = y2 - y1; } else { x1 = CONIC_MIN_VAL - Trans[0]; x2 = CONIC_MAX_VAL - Trans[0]; y1 = CONIC_MIN_VAL - Trans[1]; y2 = CONIC_MAX_VAL - Trans[1]; dx = CONIC_DELTA_VAL; dy = CONIC_DELTA_VAL; } x = x1 * cos(RotAngle) - y1 * sin(RotAngle); y1 = x1 * sin(RotAngle) + y1 * cos(RotAngle); x1 = x; Points[1][0] = x1; Points[1][1] = dx; Points[2][0] = SQR(x1) * (- A1 / E1); Points[2][1] = 2 * dx * x1 * (- A1 / E1); Points[2][2] = SQR(dx) * (- A1 / E1); } else if (!APX_EQ(C1, 0.0) && !APX_EQ(D1, 0.0)) { /* Update the translation factors. */ Trans[1] = -E1 / (2.0 * C1); /* Update the control points. */ if (PStartXY != NULL && PEndXY != NULL) { x1 = PStartXY[0] - Trans[0]; x2 = PEndXY[0] - Trans[0]; y1 = PStartXY[1] - Trans[1]; y2 = PEndXY[1] - Trans[1]; dx = x2 - x1; dy = y2 - y1; } else { x1 = CONIC_MIN_VAL - Trans[0]; x2 = CONIC_MAX_VAL - Trans[0]; y1 = CONIC_MIN_VAL - Trans[1]; y2 = CONIC_MAX_VAL - Trans[1]; dx = CONIC_DELTA_VAL; dy = CONIC_DELTA_VAL; } x = x1 * cos(RotAngle) - y1 * sin(RotAngle); y1 = x1 * sin(RotAngle) + y1 * cos(RotAngle); x1 = x; Points[1][0] = SQR(y1) * (- C1 / D1); Points[1][1] = 2 * dy * y1 * (- C1 / D1); Points[1][2] = SQR(dy) * (- C1 / D1); Points[2][0] = y1; Points[2][1] = dy; } else { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } Crv = CnvrtPower2BezierCrv(PwrCrv); CagdCrvFree(PwrCrv); } else { /* Update the translation factors. */ Trans[0] = -D1 / (2.0 * A1); Trans[1] = -E1 / (2.0 * C1); if (PStartXY != NULL && PEndXY != NULL) { x1 = PStartXY[0] - Trans[0]; x2 = PEndXY[0] - Trans[0]; y1 = PStartXY[1] - Trans[1]; y2 = PEndXY[1] - Trans[1]; x = x1 * cos(RotAngle) - y1 * sin(RotAngle); y1 = x1 * sin(RotAngle) + y1 * cos(RotAngle); x1 = x; dx = x2 - x1; dy = y2 - y1; } F1 -= (C1 * SQR(D1) + A1 * SQR(E1)) / (4 * A1 * C1); /* We now have conic in canonic form of: A1 x^2 + C1 y^2 + F1 = 0. */ if (Desc > 0) { CagdRType R[2]; /* Its an ellipse. */ if (F1 / A1 < 0 && F1 / C1 < 0) { STATIC_DATA CagdPtStruct POrigin = { NULL, NULL, { 0.0, 0.0, 0.0 } }; R[0] = sqrt(-F1 / A1); R[1] = sqrt(-F1 / C1); if (PStartXY != NULL && PEndXY != NULL) { CagdRType StartAngle, EndAngle; /* Create an ellipse from PStartXY to PEndXY! */ StartAngle = atan2(y1 / R[1], x1 / R[0]); if (StartAngle < 0.0) StartAngle = RAD2DEG(StartAngle) + 360.0; else StartAngle = RAD2DEG(StartAngle); EndAngle = atan2(y2 / R[1], x2 / R[0]); if (EndAngle < 0.0) EndAngle = RAD2DEG(EndAngle) + 360.0; else EndAngle = RAD2DEG(EndAngle); Crv = CagdCrvCreateArc(&POrigin, 1.0, StartAngle, EndAngle); } else { /* Create a full ellipse! */ Crv = BspCrvCreateUnitCircle(); } } else { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } /* Scale the circular shape into an ellipse. */ MatGenMatScale(R[0], R[1], 1.0, Mat); TCrv = CagdCrvMatTransform(Crv, Mat); CagdCrvFree(Crv); Crv = TCrv; } else { /* t < 0 */ /* Its an hyperbola. */ CAGD_FATAL_ERROR(CAGD_ERR_HYPERBOLA_NO_SUPPORT); return NULL; } } CagdCrvTransform(Crv, Trans, 1.0); MatGenMatRotZ1(RotAngle, Mat); TCrv = CagdCrvMatTransform(Crv, Mat); CagdCrvFree(Crv); Crv = TCrv; CAGD_SET_GEOM_TYPE(Crv, CAGD_GEOM_CONIC_SEC); return Crv; } /***************************************************************************** * DESCRIPTION: M * Constructs an ellipse in the XY plane through the given 3 points of M * minimal area. The A,B,C,D,E,F coefficients of the bounding ellipse as in M * A x^2 + B xy + C y^2 + D x + E y + F = 0. M * are returned. M * Algorithm: V * 1. Compute center of mass, C := (Pt1 + Pt2 + Pt3) / 3 V * V * 3 V * 2. Computer a 2x2 matrix N = 1/3 Sum (Pti - C) (Pti - C)^T V * i=1 V * V * 3. M = N^{-1} V * V * 4. The ellipse E: (P - C)^T M (P - C) - Z = 0, Z constant, P = (x, y). V * M * See also: "Exact Primitives for Smallest Enclosing Ellipses", M * by Bernd Gartner and Sven Schonherr, Proceedings of the 13th annual M * symposium on Computational geometry, 1997. M * * * PARAMETERS: M * Pt1, Pt2, Pt3: The 3 input points. Assumed non-linear. M * A, B, C, D, E, F: Coefficients of the computed bounding ellipse. M * * * RETURN VALUE: M * int: TRUE if succesful, FALSE otherwise. M * * * SEE ALSO: M * CagdEllipseOffset, CagdCreateConicCurve, CagdCreateConicCurve2 M * * * KEYWORDS: M * CagdEllipse3Points, ellipse M *****************************************************************************/ int CagdEllipse3Points(CagdPType Pt1, CagdPType Pt2, CagdPType Pt3, CagdRType *A, CagdRType *B, CagdRType *C, CagdRType *D, CagdRType *E, CagdRType *F) { CagdPType Center, DPt1, DPt2, DPt3; CagdRType Det, Nrml, N[2][2], M[2][2]; /* Compute the center of mass the three vectors from it to points. */ PT_COPY(Center, Pt1); PT_ADD(Center, Center, Pt2); PT_ADD(Center, Center, Pt3); PT_SCALE(Center, 1.0 / 3.0); PT_SUB(DPt1, Pt1, Center); PT_SUB(DPt2, Pt2, Center); PT_SUB(DPt3, Pt3, Center); N[0][0] = (SQR(DPt1[0]) + SQR(DPt2[0]) + SQR(DPt3[0])) / 3; N[1][0] = N[0][1] = (DPt1[0] * DPt1[1] + DPt2[0] * DPt2[1] + DPt3[0] * DPt3[1]) / 3; N[1][1] = (SQR(DPt1[1]) + SQR(DPt2[1]) + SQR(DPt3[1])) / 3; Det = N[0][0] * N[1][1] - N[0][1] * N[1][0]; if (APX_EQ_EPS(Det, 0.0, IRIT_UEPS)) return FALSE; M[0][0] = N[1][1] / Det; M[1][1] = N[0][0] / Det; M[0][1] = M[1][0] = -N[0][1] / Det; *A = M[0][0]; *B = M[0][1] + M[1][0]; *C = M[1][1]; Nrml = 1.0 / MAX(MAX(FABS(*A), FABS(*B)), FABS(*C)); if (APX_EQ_EPS(Nrml, 0.0, IRIT_UEPS)) return FALSE; *A *= Nrml; *B *= Nrml; *C *= Nrml; *D = (-2 * M[0][0] * Center[0] - Center[1] * (M[1][0] + M[0][1])) * Nrml; *E = (-2 * M[1][1] * Center[1] - Center[0] * (M[1][0] + M[0][1])) * Nrml; *F = -(*A * SQR(Pt1[0]) + *B * Pt1[0] * Pt1[1] + *C * SQR(Pt1[1]) + *D * Pt1[0] + *E * Pt1[1]); return TRUE; } /***************************************************************************** * DESCRIPTION: M * Update the implicit form of the given ellipse with some offset Offset. M * * * PARAMETERS: M * A, B, C, D, E, F: The six coefficients of the ellipse. M * Offset: Offset amount. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * CagdEllipse3Points, CagdCreateConicCurve, CagdCreateConicCurve2 M * * * KEYWORDS: M * CagdEllipseOffset, ellipse M *****************************************************************************/ int CagdEllipseOffset(CagdRType *A, CagdRType *B, CagdRType *C, CagdRType *D, CagdRType *E, CagdRType *F, CagdRType Offset) { VectorType Trans; CagdRType Desc, RotAngle = APX_EQ(*B, 0.0) ? 0.0 : atan2(*B, *A - *C) * 0.5, A1 = ((*A + *C) + *B * sin(2 * RotAngle) + (*A - *C) * cos(2 * RotAngle)) * 0.5, B1 = *B * cos(2 * RotAngle) - (*A - *C) * sin(2 * RotAngle), C1 = ((*A + *C) - *B * sin(2 * RotAngle) - (*A - *C) * cos(2 * RotAngle)) * 0.5, D1 = *D * cos(RotAngle) + *E * sin(RotAngle), E1 = -*D * sin(RotAngle) + *E * cos(RotAngle), F1 = *F; if (!APX_EQ(B1, 0.0) || (APX_EQ(A1, 0.0) && APX_EQ(C1, 0.0))) { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return FALSE; } /* We now have conic: A1 x^2 + C1 y^2 + D1 x + E1 y + F1 = 0. */ if ((Desc = A1 * C1) <= 0.0) { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return FALSE; } /* Update the translation factors. */ Trans[0] = -D1 / (2.0 * A1); Trans[1] = -E1 / (2.0 * C1); F1 -= (C1 * SQR(D1) + A1 * SQR(E1)) / (4 * A1 * C1); /* We now have conic in canonic form of: A1 x^2 + C1 y^2 + F1 = 0. */ /* Update the major and minor axes with the desired offset. */ A1 = -F1 / SQR(sqrt(-F1 / A1) + Offset); C1 = -F1 / SQR(sqrt(-F1 / C1) + Offset); D1 = -2.0 * A1 * Trans[0]; E1 = -2.0 * C1 * Trans[1]; F1 += (C1 * SQR(D1) + A1 * SQR(E1)) / (4 * A1 * C1); /* And rotate back. */ RotAngle = -RotAngle; *A = ((A1 + C1) + B1 * sin(2 * RotAngle) + (A1 - C1) * cos(2 * RotAngle)) * 0.5; *B = B1 * cos(2 * RotAngle) - (A1 - C1) * sin(2 * RotAngle); *C = ((A1 + C1) - B1 * sin(2 * RotAngle) - (A1 - C1) * cos(2 * RotAngle)) * 0.5; *D = D1 * cos(RotAngle) + E1 * sin(RotAngle); *E = -D1 * sin(RotAngle) + E1 * cos(RotAngle); *F = F1; return TRUE; } /***************************************************************************** * DESCRIPTION: M * Transform given conic form A x^2 + B xy + C y^2 + D x + E y + F = 0, M * using Mat, in the XY plane. Algorithm: M * M * 1. Convert the implicit conic to a matrix form as: V * [ A B/2 0 D/2] [x] V * V * [x, y, 0, 1] [B/2 C 0 E/2] [y] = P^T M P = 0 V * V * [ 0 0 1 0 ] [0] V * V * [D/2 E/2 0 F ] [1] V * V * 2. Compute N = Mat^{-1} the inverse of the desired transformation. V * V * 3. Compute K = N M N^T and decompose K back to A-F coefficients. V * * * PARAMETERS: M * A, B, C, D, E, F: The six coefficients of the conic. Updated in place. M * Mat: Transformation matrix in the XY plane. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * CagdQuadricMatTransform, CagdSrfTransfrom, CagdCrvTransform M * * * KEYWORDS: M * CagdConicMatTransform, implicit M *****************************************************************************/ int CagdConicMatTransform(CagdRType *A, CagdRType *B, CagdRType *C, CagdRType *D, CagdRType *E, CagdRType *F, CagdMType Mat) { MatrixType M, K, InvMat, TransposeMat; /* Construct the M matrix. */ MatGenUnitMat(M); M[0][0] = *A; M[0][1] = M[1][0] = *B * 0.5; M[1][1] = *C; M[0][3] = M[3][0] = *D * 0.5; M[1][3] = M[3][1] = *E * 0.5; M[3][3] = *F; MatInverseMatrix(Mat, InvMat); /* Compute N M N^T. */ MatMultTwo4by4(K, InvMat, M); MatTranspMatrix(InvMat, TransposeMat); MatMultTwo4by4(K, K, TransposeMat); *A = K[0][0]; *B = K[0][1] * 2.0; *C = K[1][1]; *D = K[0][3] * 2.0; *E = K[1][3] * 2.0; *F = K[3][3]; return TRUE; } /***************************************************************************** * DESCRIPTION: M * Transform given quadric form M * A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G x + H y + I z + J = 0, M * using Mat. Algorithm: M * M * 1. Convert the implicit quadric to a matrix form as: V * [ A D/2 E/2 G/2] [x] V * V * [x, y, z, 1] [D/2 B F/2 H/2] [y] = P^T M P = 0 V * V * [E/2 F/2 C I/2] [z] V * V * [G/2 H/2 I/2 J ] [1] V * V * 2. Compute N = Mat^{-1} the inverse of the desired transformation. V * V * 3. Compute K = N M N^T and decompose K back to A-J coefficients. V * * * PARAMETERS: M * A, B, C, D, E, F, G, H, I, J: The ten coefficients of the quadric. M 8 Updated in place. M * Mat: Transformation matrix. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * CagdConicMatTransform, CagdSrfTransfrom, CagdCrvTransform M * * * KEYWORDS: M * CagdQuadricMatTransform, implicit M *****************************************************************************/ int CagdQuadricMatTransform(CagdRType *A, CagdRType *B, CagdRType *C, CagdRType *D, CagdRType *E, CagdRType *F, CagdRType *G, CagdRType *H, CagdRType *I, CagdRType *J, CagdMType Mat) { MatrixType M, K, InvMat, TransposeMat; /* Construct the M matrix. */ M[0][0] = *A; M[1][1] = *B; M[2][2] = *C; M[0][1] = M[1][0] = *D * 0.5; M[0][2] = M[2][0] = *E * 0.5; M[1][2] = M[2][1] = *F * 0.5; M[0][3] = M[3][0] = *G * 0.5; M[1][3] = M[3][1] = *H * 0.5; M[2][3] = M[3][2] = *I * 0.5; M[3][3] = *J; MatInverseMatrix(Mat, InvMat); /* Compute N M N^T. */ MatMultTwo4by4(K, InvMat, M); MatTranspMatrix(InvMat, TransposeMat); MatMultTwo4by4(K, K, TransposeMat); *A = K[0][0]; *B = K[1][1]; *C = K[2][2]; *D = K[0][1] * 2.0; *E = K[0][2] * 2.0; *F = K[1][2] * 2.0; *G = K[0][3] * 2.0; *H = K[1][3] * 2.0; *I = K[2][3] * 2.0; *J = K[3][3]; return TRUE; } /***************************************************************************** * DESCRIPTION: M * Construct rational quadric surface above (and below) the given conic M * in the XY plane of Z height. The conic is given in the A-F coefficients M * as A x^2 + B xy + C y^2 + D x + E y + F = 0 and the quadric is returned M * in the A-J coefficients as: M * A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G x + H y + I z + J = 0. M * * * PARAMETERS: M * A, B, C, D, E, F, G, H, I, J: Input: in A-F the conic and in J the M * Z height. Output: the new 10 coefficients of the quadric. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * CagdCreateQuadricSrf, CagdEllipse3Points, CagdEllipseOffset M * * * KEYWORDS: M * CagdConic2Quadric M *****************************************************************************/ int CagdConic2Quadric(CagdRType *A, CagdRType *B, CagdRType *C, CagdRType *D, CagdRType *E, CagdRType *F, CagdRType *G, CagdRType *H, CagdRType *I, CagdRType *J) { CagdRType RotAngle = APX_EQ(*B, 0.0) ? 0.0 : atan2(*B, *A - *C) * 0.5, A1 = ((*A + *C) + *B * sin(2 * RotAngle) + (*A - *C) * cos(2 * RotAngle)) * 0.5, B1 = *B * cos(2 * RotAngle) - (*A - *C) * sin(2 * RotAngle), C1 = ((*A + *C) - *B * sin(2 * RotAngle) - (*A - *C) * cos(2 * RotAngle)) * 0.5, D1 = *D * cos(RotAngle) + *E * sin(RotAngle), E1 = -*D * sin(RotAngle) + *E * cos(RotAngle), F1 = *F, ZHeight = *J; if (!APX_EQ(B1, 0.0) || (APX_EQ(A1, 0.0) && APX_EQ(C1, 0.0))) { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return FALSE; } F1 -= (C1 * SQR(D1) + A1 * SQR(E1)) / (4 * A1 * C1); /* Update the modified coefficients. */ *J = *F; *I = 0; *H = *E; *G = *D; *E = *F = 0; *D = *B; *B = *C; *C = -F1 / SQR(ZHeight); return TRUE; } /***************************************************************************** * DESCRIPTION: M * Construct rational quadric surface out of the 9 coefficients of: M * A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + G x + H y + I z + J = 0. M * Based on: M * Bezier Curves and Surface Patches on Quadrics, by Josef Hoschek, M * Mathematical methods in Computer aided Geometric Design II, Tom Lyche and M * Larry L. Schumaker (eds.), pp 331-342, 1992. M * * * PARAMETERS: M * A, B, C, D, E, F, G, H, I, J: The ten coefficients of the quadric. M * * * RETURN VALUE: M * CagdSrfStruct *: A quadric surface representing the given form. M * * * SEE ALSO: M * BspCrvCreateCircle, BspCrvCreateUnitCircle, CagdCrvCreateArc, M * BzrCrvCreateArc, CagdCreateConicCurve, CagdCreateConicCurve2, M * CagdConic2Quadric M * * * KEYWORDS: M * CagdCreateQuadricSrf M *****************************************************************************/ CagdSrfStruct *CagdCreateQuadricSrf(CagdRType A, CagdRType B, CagdRType C, CagdRType D, CagdRType E, CagdRType F, CagdRType G, CagdRType H, CagdRType I, CagdRType J) { CagdBType SwitchXY = FALSE, SwitchXZ = FALSE, SwitchYZ = FALSE; int i; CagdRType **Points, N[9], b[3]; VectorType Trans; MatrixType Mat, Mat2, M, U, d, V; CagdSrfStruct *PwrSrf, *TSrf, *BzrSrf = NULL; /* Computes the translation factors. Substituting (x-X0) for x (and */ /* same for y and z), one gets the following constraints for the linear */ /* terms to vanish (G, H, I in the original general quadric): */ /* [2A D E] [X0] = [G] */ /* [ D 2B F] [Y0] = [H] */ /* [ E F 2C] [Z0] = [I] */ /* The solution for X0, Y0, Z0 zeros these linear terms and yields new */ /* constant J term: J + A X0^2 + B Y0^2 + C Z0^2 + */ /* D X0 Y0 + E X0 Z0 + F Y0 Z0 - G X0 - H Y0 - I Z0. */ N[0] = 2 * A; N[3] = D; N[6] = E; N[1] = D; N[4] = 2 * B; N[7] = F; N[2] = E; N[5] = F; N[8] = 2 * C; b[0] = G; b[1] = H; b[2] = I; IritQRUnderdetermined(N, NULL, NULL, 3, 3); IritQRUnderdetermined(NULL, Trans, b, 3, 3); /* Find the translational factors. */ J += A * SQR(Trans[0]) + B * SQR(Trans[1]) + C * SQR(Trans[2]) + D * Trans[0] * Trans[1] + E * Trans[0] * Trans[2] + F * Trans[1] * Trans[2] - G * Trans[0] - H * Trans[1] - I * Trans[2]; G = H = I = 0.0; /* Compute the rotation factors. Now we have a quadrics of the form: */ /* A x^2 + B y^2 + C z^2 + D xy + E xz + F yz + J = 0. */ /* Writing the quadric form as: */ /* [x y z 1] [A D/2 E/2 G/2] [x] = [x y z 1] M [x] = 0 */ /* [D/2 B F/2 H/2] [y] [y] */ /* [E/2 F/2 C I/2] [z] [z] */ /* [G/2 H/2 I/2 J ] [1] [1] */ /* we seek a rotation of M that make it diagonal. Toward this end we */ /* employ SVD to decompose M into U N V, N diagonal and U and V are two */ /* rotational matrices that are also transpose (==inverse) of each other.*/ M[0][0] = A; M[1][1] = B; M[2][2] = C; M[3][3] = J; M[0][1] = M[1][0] = D * 0.5; M[0][2] = M[2][0] = E * 0.5; M[1][2] = M[2][1] = F * 0.5; M[0][3] = M[3][0] = G * 0.5; M[1][3] = M[3][1] = H * 0.5; M[2][3] = M[3][2] = I * 0.5; JacobiMatrixDiag4x4(M, U, d, V); A = d[0][0]; B = d[1][1]; C = d[2][2]; J = d[3][3]; /* We now have quadric: A x^2 + B y^2 + C z^2 + J = 0. */ if (J == 0.0) { if (A < 0.0) { A = -A; B = -B; C = -C; } } else { if (J > 0.0) { /* Make sure J is not positive. */ A = -A; B = -B; C = -C; J = -J; } if (A * J >= 0.0) { if (B * J >= 0.0) { if (C * J >= 0.0) { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } else { SwitchXZ = TRUE; SWAP(CagdRType, A, C); } } else { SwitchXY = TRUE; SWAP(CagdRType, A, B); } } } /* Now J is not positive and A is not negative. */ if (B >= 0.0 && C >= 0.0) { /* Quadrics is ellipsoid A x^2 + B y^2 + C z^2 + J w^2 = 0, */ /* A >= 0, B >= 0, C >= 0, J <= 0. */ STATIC_DATA CagdVType Center = { 0.0, 0.0, 0.0 }; CagdPType R; BzrSrf = CagdPrimSphereSrf(Center, 1.0, TRUE); R[0] = sqrt(-J / A); R[1] = sqrt(-J / B); R[2] = sqrt(-J / C); /* Scale the circular shape into an ellipse. */ MatGenMatScale(R[0], R[1], R[2], Mat); TSrf = CagdSrfMatTransform(BzrSrf, Mat); CagdSrfFree(BzrSrf); BzrSrf = TSrf; } else if (B <= 0.0 && C <= 0.0) { /* Quadrics is now A x^2 + B y^2 + C z^2 + J w^2 = 0, */ /* A >= 0, B <= 0, C <= 0, J <= 0. */ /* The above hyperboloid of two sheets can be written as, */ /* x = (1/w) * (a^2 + b^2 + c^2) / sqrt(A), */ /* y = (1/w) * 2ab / sqrt(B), */ /* z = (1/w) * 2ac / sqrt(C), */ /* and w equals w = (a^2 - b^2 - c^2) / -sqrt(-J). */ /* The only (large) remaining question is how to parametrize a, b, c.*/ /* For no better selection, we select a = u, b = 1 - v, c = v, or */ /* x = (1/w) * (u^2 + 1 - 2 * v + 2 * v^2) / sqrt(A), */ /* y = (1/w) * (2 * u - 2 * uv) / sqrt(B), */ /* z = (1/w) * (2 * u * v) / sqrt(C), */ /* and w equals w = (u^2 - 1 + 2 * v - 2 * v^2) / sqrt(-J) */ A = QUADRIC_INVERT(A); B = QUADRIC_INVERT(-B); C = QUADRIC_INVERT(-C); J = -QUADRIC_INVERT(-J); PwrSrf = CagdSrfNew(CAGD_SPOWER_TYPE, CAGD_PT_P3_TYPE, 3, 3); Points = PwrSrf -> Points; for (i = 0; i <= 3; i++) ZAP_MEM(Points[i], 9 * sizeof(CagdRType)); /* Coefficients of the biquadratic surface are orders as follows */ /* 0) 1 1) u 2) u^2 */ /* 3) v 4) u * v 5) u^2 * v */ /* 6) v^2 7) u * v^2 8) u^2 * v^2 */ Points[0][0] = -J; /* The w term. */ Points[0][2] = J; Points[0][3] = 2.0 * J; Points[0][6] = -2.0 * J; Points[1][0] = A; /* The x term. */ Points[1][2] = A; Points[1][3] = -2.0 * A; Points[1][6] = 2.0 * A; Points[2][1] = 2.0 * B; /* The y term. */ Points[2][4] = -2.0 * B; Points[3][4] = 2.0 * C; /* The z term. */ BzrSrf = CnvrtPower2BezierSrf(PwrSrf); CagdSrfFree(PwrSrf); } else if ((B >= 0.0 && C <= 0.0) || (B <= 0.0 && C >= 0.0)) { if (B <= 0.0 && C >= 0.0) { /* Exchange the role of the B and C coefficients. */ SWAP(CagdRType, B, C); SwitchYZ = TRUE; } /* Quadrics is now A x^2 + B y^2 + C z^2 + J w^2 = 0, */ /* A >= 0, B >= 0, C <= 0, J <= 0. */ /* The above hyperboloid of one sheet can be written as, */ /* x = (1/w) * (a^2 - b^2 + c^2) / sqrt(A), */ /* y = (1/w) * 2ab / sqrt(B), */ /* z = (1/w) * 2ac / -sqrt(-C), */ /* and w equals w = (a^2 + b^2 - c^2) / -sqrt(-J). */ /* The only (large) remaining question is how to parametrize a, b, c.*/ /* For no better selection, we select a = u, b = 1 - v, c = v, or */ /* x = (1/w) * (u^2 - 1 + 2 * v) / sqrt(A), */ /* y = (1/w) * (2 * u - 2 * uv) / sqrt(B), */ /* z = (1/w) * (2 * u * v) / -sqrt(C), */ /* and w equals w = (1 + u^2 - 2 * v) / -sqrt(-J). */ A = QUADRIC_INVERT(A); B = QUADRIC_INVERT(B); C = -QUADRIC_INVERT(-C); J = -QUADRIC_INVERT(-J); PwrSrf = CagdSrfNew(CAGD_SPOWER_TYPE, CAGD_PT_P3_TYPE, 3, 3); Points = PwrSrf -> Points; for (i = 0; i <= 3; i++) ZAP_MEM(Points[i], 9 * sizeof(CagdRType)); /* Coefficients of the biquadratic surface are orders as follows */ /* 0) 1 1) u 2) u^2 */ /* 3) v 4) u * v 5) u^2 * v */ /* 6) v^2 7) u * v^2 8) u^2 * v^2 */ Points[0][0] = J; /* The w term. */ Points[0][2] = J; Points[0][3] = -2.0 * J; Points[1][0] = -A; /* The x term. */ Points[1][2] = A; Points[1][3] = 2.0 * A; Points[2][1] = 2.0 * B; /* The y term. */ Points[2][4] = -2.0 * B; Points[3][4] = 2.0 * C; /* The z term. */ BzrSrf = CnvrtPower2BezierSrf(PwrSrf); CagdSrfFree(PwrSrf); } else { CAGD_FATAL_ERROR(CAGD_ERR_INVALID_CONIC_COEF); return NULL; } if (SwitchYZ) { MatGenMatRotX1(M_PI * 0.5, Mat); MatGenMatScale(1, 1, -1, Mat2); MatMultTwo4by4(Mat, Mat2, Mat); TSrf = CagdSrfMatTransform(BzrSrf, Mat); CagdSrfFree(BzrSrf); BzrSrf = TSrf; } if (SwitchXZ) { MatGenMatRotY1(-M_PI * 0.5, Mat); MatGenMatScale(1, 1, -1, Mat2); MatMultTwo4by4(Mat, Mat2, Mat); TSrf = CagdSrfMatTransform(BzrSrf, Mat); CagdSrfFree(BzrSrf); BzrSrf = TSrf; } if (SwitchXY) { MatGenMatRotZ1(M_PI * 0.5, Mat); MatGenMatScale(1, -1, 1, Mat2); MatMultTwo4by4(Mat, Mat2, Mat); TSrf = CagdSrfMatTransform(BzrSrf, Mat); CagdSrfFree(BzrSrf); BzrSrf = TSrf; } if (BzrSrf != NULL) { TSrf = CagdSrfMatTransform(BzrSrf, V); CagdSrfFree(BzrSrf); BzrSrf = TSrf; PT_SCALE(Trans, -1.0); CagdSrfTransform(BzrSrf, Trans, 1.0); CAGD_SET_GEOM_TYPE(BzrSrf, CAGD_GEOM_CONIC_SEC); } return BzrSrf; }