/******************************************************************************
* EvalCurv.c - Evaluate curvature properties of curves and surfaces.	      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, Nov. 96.					      *
******************************************************************************/

#include <string.h>
#include "symb_loc.h"
#include "geomat3d.h"
#include "miscattr.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Preprocess a given surface so we can evaluate curvature properties from  M
* it efficiently, at every point.  See SymbEvalCurvature for actual          M
* curvature at surface point evaluations				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Surface to preprocess.                                       M
*   Init:	TRUE for initializing, FALSE for clearing out.               M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbEvalCurvature                                                        M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbEvalSrfCurvPrep, curvature                                           M
*****************************************************************************/
void SymbEvalSrfCurvPrep(CagdSrfStruct *Srf, CagdBType Init)
{
    CagdSrfStruct **DSrfs;

    if (Init) {
        DSrfs = (CagdSrfStruct **) IritMalloc(6 * sizeof(CagdSrfStruct *));
	DSrfs[0] = SymbSrfNormalSrf(Srf);
	DSrfs[1] = CagdSrfDerive(Srf, CAGD_CONST_U_DIR);
	DSrfs[2] = CagdSrfDerive(Srf, CAGD_CONST_V_DIR);
	DSrfs[3] = CagdSrfDerive(DSrfs[1], CAGD_CONST_U_DIR);
	DSrfs[4] = CagdSrfDerive(DSrfs[2], CAGD_CONST_V_DIR);
	DSrfs[5] = CagdSrfDerive(DSrfs[2], CAGD_CONST_U_DIR);

	AttrSetPtrAttrib(&Srf -> Attr, "_EvalCurv", DSrfs);
    }
    else {
        int i;

      	if ((DSrfs = (CagdSrfStruct **) AttrGetPtrAttrib(Srf -> Attr,
							 "_EvalCurv")) == NULL)
	    return;

	/* Release the three preprocessed derivative surfaces. */
	for (i = 0; i < 6; i++)
	    CagdSrfFree(DSrfs[i]);
	IritFree(DSrfs);
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Evaluate a given surface's curvature properties.  Returns the principal  M
* curvatures and directions for the given surface location.		     M
*   This function must be invoked after SymbEvalSrfCurvPrep was called to    M
* initialize the proper data structures, for fast curvature at a point       M
* evaluation.  SymbEvalSrfCurvPrep should be called at the end to release    M
* these data structures.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Surface to evaluate its curvature properties.                M
*   U, V:	Location of evaluation.					     M
*   DirInUV:	If TRUE principal directions are given in UV, otherwise in   M
*		Euclidean 3-space.					     M
*   K1, K2:	Principal curvatures.					     M
*   D1, D2:	Principal directions.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:	TRUE if succesful, FALSE otherwise.                          M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbEvalSrfCurvPrep                                                      M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbEvalSrfCurvature, curvature                                          M
*****************************************************************************/
int SymbEvalSrfCurvature(CagdSrfStruct *Srf,
			 CagdRType U,
			 CagdRType V,
			 CagdBType DirInUV,
			 CagdRType *K1,
			 CagdRType *K2,
			 CagdVType D1,
			 CagdVType D2)
{
    CagdRType *R, E, F, G, L, M, N, A, B, C, Disc, v1, v2;
    CagdVType Nrml, DSrfU, DSrfV, DSrfUU, DSrfVV, DSrfUV, Vec;

    CagdSrfStruct **DSrfs;

    if ((DSrfs = (CagdSrfStruct **) AttrGetPtrAttrib(Srf -> Attr, "_EvalCurv"))
								   == NULL)
        return FALSE;

    R = CagdSrfEval(DSrfs[0], U, V);
    CagdCoerceToE3(Nrml, &R, -1, DSrfs[0] -> PType);
    PT_NORMALIZE(Nrml);

    R = CagdSrfEval(DSrfs[1], U, V);
    CagdCoerceToE3(DSrfU, &R, -1, DSrfs[1] -> PType);

    R = CagdSrfEval(DSrfs[2], U, V);
    CagdCoerceToE3(DSrfV, &R, -1, DSrfs[2] -> PType);

    R = CagdSrfEval(DSrfs[3], U, V);
    CagdCoerceToE3(DSrfUU, &R, -1, DSrfs[3] -> PType);

    R = CagdSrfEval(DSrfs[4], U, V);
    CagdCoerceToE3(DSrfVV, &R, -1, DSrfs[4] -> PType);

    R = CagdSrfEval(DSrfs[5], U, V);
    CagdCoerceToE3(DSrfUV, &R, -1, DSrfs[5] -> PType);

    /* Compute the coefficients of the first and second fundamental forms.   */
    E = DOT_PROD(DSrfU, DSrfU);
    F = DOT_PROD(DSrfU, DSrfV);
    G = DOT_PROD(DSrfV, DSrfV);

    L = DOT_PROD(DSrfUU, Nrml);
    M = DOT_PROD(DSrfUV, Nrml);
    N = DOT_PROD(DSrfVV, Nrml);

    /* Solve the quadratic equation for the principal curvature:             */
    /* (EG - F^2) k^2 - (GL + EN - 2FM) k + (LN - M^2) = 0.		     */
    A = E * G - F * F;
    B = -(G * L + E * N - 2 * F * M);
    C = L * N - M * M;

    Disc = sqrt(B * B - 4 * A * C);
    *K1 = (-B - Disc) / (2 * A);
    *K2 = (-B + Disc) / (2 * A);

    /* Compute the principal directions - K1. */
    A = L - *K1 * E;
    B = M - *K1 * F;
    C = N - *K1 * G;
    if (fabs(A) > fabs(C)) {
        v1 = B;
        v2 = -A;
    }
    else {
        v1 = C;
	v2 = -B;
    }
    if (DirInUV) {
	D1[0] = v1;
	D1[1] = v2;
	D1[2] = 0.0;
    }
    else {
	PT_COPY(D1, DSrfU);
	PT_SCALE(D1, v1);
	PT_COPY(Vec, DSrfV);
	PT_SCALE(Vec, v2);
	PT_ADD(D1, D1, Vec);
	PT_NORMALIZE(D1);
    }

    /* Compute the principal directions - K2. */
    A = L - *K2 * E;
    B = M - *K2 * F;
    C = N - *K2 * G;
    if (fabs(A) > fabs(C)) {
        v1 = B;
        v2 = -A;
    }
    else {
        v1 = C;
	v2 = -B;
    }
    if (DirInUV) {
	D2[0] = v1;
	D2[1] = v2;
	D2[2] = 0.0;
    }
    else {
	PT_COPY(D2, DSrfU);
	PT_SCALE(D2, v1);
	PT_COPY(Vec, DSrfV);
	PT_SCALE(Vec, v2);
	PT_ADD(D2, D2, Vec);
	PT_NORMALIZE(D2);
    }

    return TRUE;
}