/*
 * Refinement plugin
 */

/****************************************************************************
 *                                                                          *
 *  Module  :   bbtp.c                                                      *
 *                                                                          *
 *  Purpose :   Tools for Bernstein Bezier Triangular Patch                 *
 *                                                                          *
 ****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "rpplugin.h"
#include "time.h"              /* pick up struct tm; */
#include <rpdbgerr.h>
#include <rwcore.h>
#include <rpworld.h>
#include "rprefine.h"

/* Check for SSE. */

#if (defined(__ICL) && defined(RWSIMD))
#if (400 <= __ICL)

#include <rtintel.h>

#endif /* (400 <= __ICL) */
#endif /* (defined(__ICL) && defined(RWSIMD)) */

#include "bbtptool.h"

#if (!defined(DOXYGEN_SHOULD_SKIP_THIS))
static const char   rcsid[] __RWUNUSED__ =
    "@@(#)$Id: bbtptool.c,v 1.18 2001/01/26 11:05:19 johns Exp $";
#endif /* (!defined(DOXYGEN_SHOULD_SKIP_THIS)) */

/****************************************************************************
 Local types
 */

/****************************************************************************
 Public variables
 */

/****************************************************************************
 External variables
 */

/****************************************************************************
 Local variables
 */

/****************************************************************************
 Local functions for the Bernstein-Bezier Triangular Patch toolkit
 */

/*****************************************************************************
 PATCH code
 */

/* generate ordinates from triangle vertices and vertex normals */
RwReal
_rtbbtpPatchEvaluate(BBTPOrdinates * ord, RwReal a, RwReal b, RwReal c)
{
    RwReal              res;

    RWFUNCTION(RWSTRING("_rtbbtpPatchEvaluate"));

    res = (ord->b021 * b * b * c) + (ord->b012 * b * c * c);
    res += (ord->b120 * a * b * b) + (ord->b102 * a * c * c);
    res += (ord->b210 * a * a * b) + (ord->b201 * a * a * c);
    res += ord->b111 * a * b * c * 2.0f;
    res *= 3.0f;

    RWRETURN(res);
}

RwReal
_rtbbtpComputeOrdinates(RwV3d * vpV0, RwV3d * vpV1, RwV3d * vpN0,
                        RwV3d * __RWUNUSED__ normalDir)
{
    RwV3d               edge1Dir;
    RwReal              alpha, dp2;

    RWFUNCTION(RWSTRING("_rtbbtpComputeOrdinates"));

    edge1Dir.x = vpV1->x - vpV0->x;
    edge1Dir.y = vpV1->y - vpV0->y;
    edge1Dir.z = vpV1->z - vpV0->z;
    /* RwV3dSub(&edge1Dir, vpV1, vpV0); */

    /* compute scaling factor */
    dp2 = edge1Dir.x * vpN0->x +
        edge1Dir.y * vpN0->y + edge1Dir.z * vpN0->z;

    /* dp2=RwV3dDotProduct(&edge1Dir, vpN0); */
    /* alpha =- (dp2) / (3.0f * dp1); */

    alpha = -dp2 * (RwReal) (0.333333);

    RWRETURN(alpha);
}

RwReal
_rt_rtbbtpComputeOrdinatesSlow(RwV3d * vpV0, RwV3d * vpV1, RwV3d * vpN0,
                               RwV3d * normalDir)
{
    RwV3d               edge1Dir;
    RwReal              alpha, dp1, dp2;

    RWFUNCTION(RWSTRING("_rt_rtbbtpComputeOrdinatesSlow"));

    edge1Dir.x = vpV1->x - vpV0->x;
    edge1Dir.y = vpV1->y - vpV0->y;
    edge1Dir.z = vpV1->z - vpV0->z;

    dp1 = normalDir->x * vpN0->x +
        normalDir->y * vpN0->y + normalDir->z * vpN0->z;

    /* compute scaling factor */
    dp2 = edge1Dir.x * vpN0->x +
        edge1Dir.y * vpN0->y + edge1Dir.z * vpN0->z;

    alpha = -(dp2) / (3.0f * dp1);

    RWRETURN(alpha);
}

/*****************************************************************************
 SURFACE code
 */
void
_rtbbtpSurfaceEvaluate(RwV3d * res, BBTPControlPoints * cps,
                       RwReal a, RwReal b, RwReal c)
{
    RwReal              a2 = a * a, a3 = a2 * a;
    RwReal              b2 = b * b, b3 = b2 * b;
    RwReal              c2 = c * c, c3 = c2 * c;

    RWFUNCTION(RWSTRING("_rtbbtpSurfaceEvaluate"));

    res->x = (cps->b021.x * b2 * c) + (cps->b012.x * b * c2);
    res->x += (cps->b120.x * a * b2) + (cps->b102.x * a * c2);
    res->x += (cps->b210.x * a2 * b) + (cps->b201.x * a2 * c);
    res->x += cps->b111.x * a * b * c * 2.0f;
    res->x *= 3.0f;
    res->x += cps->b300.x * a3 + cps->b030.x * b3 + cps->b003.x * c3;

    res->y = (cps->b021.y * b2 * c) + (cps->b012.y * b * c2);
    res->y += (cps->b120.y * a * b2) + (cps->b102.y * a * c2);
    res->y += (cps->b210.y * a2 * b) + (cps->b201.y * a2 * c);
    res->y += cps->b111.y * a * b * c * 2.0f;
    res->y *= 3.0f;
    res->y += cps->b300.y * a3 + cps->b030.y * b3 + cps->b003.y * c3;

    res->z = (cps->b021.z * b2 * c) + (cps->b012.z * b * c2);
    res->z += (cps->b120.z * a * b2) + (cps->b102.z * a * c2);
    res->z += (cps->b210.z * a2 * b) + (cps->b201.z * a2 * c);
    res->z += cps->b111.z * a * b * c * 2.0f;
    res->z *= 3.0f;
    res->z += cps->b300.z * a3 + cps->b030.z * b3 + cps->b003.z * c3;

    RWRETURNVOID();
}

/*****************************************************************************
 Intersect the tangent plane at v0 with the line from third of the way along
 the edge in the direction of the interpolated normal
 */
RwV3d
_rt_rtbbtpComputeControlPointsOld(RwV3d * v0, RwV3d * n0, RwV3d * v1,
                                  RwV3d * n1)
{
    RwV3d               N, tmp, V, E;
    RwReal              dp1, dp2;

    RWFUNCTION(RWSTRING("_rt_rtbbtpComputeControlPointsOld"));

    /* compute a unit length average normal */
    RwV3dScale(&N, n0, (RwReal) 0.666667);
    RwV3dScale(&tmp, n1, (RwReal) 0.333333);
    RwV3dAdd(&N, &N, &tmp);

    RwV3dNormalize(&N, &N);

    /* edge dir */
    RwV3dSub(&E, v1, v0);

    dp1 = RwV3dDotProduct(&E, n0);
    dp2 = RwV3dDotProduct(&N, n0) * -0.3f;
    RwV3dScale(&N, &N, dp1 / dp2);
    RwV3dScale(&N, &N, 0.95f);

    /* third way (TM, New Labour) along the edge */
    RwV3dScale(&V, v0, (RwReal) 0.666667);
    RwV3dScale(&tmp, v1, (RwReal) 0.333333);
    RwV3dAdd(&V, &V, &tmp);

    RwV3dAdd(&V, &V, &N);

    RWRETURN(V);
}

/*****************************************************************************
 Given an edge rotate the third-length edge vector about the vertex such that
 it is mutually perpendicular to n0 and (n0 cross edge).
 */
RwV3d
_rtbbtpComputeControlPoints(RwV3d * v0, RwV3d * n0, RwV3d * v1,
                            RwV3d * __RWUNUSED__ n1)
{
    RwV3d               A, E, B, V;
    RwReal              length;

    RWFUNCTION(RWSTRING("_rtbbtpComputeControlPoints"));

    RwV3dSub(&E, v1, v0);
    length = RwV3dLength(&E);
    RwV3dScale(&E, &E, 1.0f / length);
    RwV3dCrossProduct(&A, n0, &E);
    /*RwV3dNormalize(&A, &A); */
    RwV3dCrossProduct(&B, &A, n0);
    RwV3dNormalize(&B, &B);
    RwV3dScale(&B, &B, (RwReal) 0.333333 * length);
    RwV3dAdd(&V, v0, &B);

    RWRETURN(V);
}

void
_rtbbtpGenerateOrdinates(BBTPOrdinates * ord, RwV3d * verts,
                         RwV3d * norms)
{
    RwV3d              *vpV0, *vpV1, *vpV2;
    RwV3d              *vpN0, *vpN1, *vpN2;
    RwV3d               tmp;

    RWFUNCTION(RWSTRING("_rtbbtpGenerateOrdinates"));

    vpV0 = verts + 0;
    vpV1 = verts + 1;
    vpV2 = verts + 2;
    vpN0 = norms + 0;
    vpN1 = norms + 1;
    vpN2 = norms + 2;

    RwV3dAdd(&tmp, vpN0, vpN1);
    RwV3dNormalize(&tmp, &tmp);
    /* RwV3dScale(&tmp, 0.5f); */

    ord->b210 = _rtbbtpComputeOrdinates(vpV0, vpV1, vpN0, &tmp);
    ord->b120 = _rtbbtpComputeOrdinates(vpV1, vpV0, vpN1, &tmp);

    RwV3dAdd(&tmp, vpN0, vpN2);
    RwV3dNormalize(&tmp, &tmp);
    /* RwV3dScale(&tmp, 0.5f); */

    ord->b201 = _rtbbtpComputeOrdinates(vpV0, vpV2, vpN0, &tmp);
    ord->b102 = _rtbbtpComputeOrdinates(vpV2, vpV0, vpN2, &tmp);

    RwV3dAdd(&tmp, vpN1, vpN2);
    RwV3dNormalize(&tmp, &tmp);
    /* RwV3dScale(&tmp, 0.5f); */

    ord->b021 = _rtbbtpComputeOrdinates(vpV1, vpV2, vpN1, &tmp);
    ord->b012 = _rtbbtpComputeOrdinates(vpV2, vpV1, vpN2, &tmp);

    ord->b111 = 0.25f * (ord->b201 + ord->b102 + ord->b021 +
                         ord->b012 + ord->b210 + ord->b120);

    RWRETURNVOID();
}

void
_rtbbtpGenerateControlPoints(BBTPControlPoints * cps, RwV3d * v,
                             RwV3d * n)
{
    RwV3d               tmp;

    RWFUNCTION(RWSTRING("_rtbbtpGenerateControlPoints"));

    /* vertex 0 */
    cps->b201 = _rtbbtpComputeControlPoints(v + 0, n + 0, v + 2, n + 2);
    cps->b210 = _rtbbtpComputeControlPoints(v + 0, n + 0, v + 1, n + 1);

    /* vertex 1 */
    cps->b021 = _rtbbtpComputeControlPoints(v + 1, n + 1, v + 2, n + 2);
    cps->b120 = _rtbbtpComputeControlPoints(v + 1, n + 1, v + 0, n + 0);

    /* vertex 2 */
    cps->b102 = _rtbbtpComputeControlPoints(v + 2, n + 2, v + 0, n + 0);
    cps->b012 = _rtbbtpComputeControlPoints(v + 2, n + 2, v + 1, n + 1);

    /* bums */
    cps->b300 = v[0];
    cps->b030 = v[1];
    cps->b003 = v[2];

    RwV3dAdd(&cps->b111, &cps->b210, &cps->b201);
    RwV3dAdd(&cps->b111, &cps->b111, &cps->b120);
    RwV3dAdd(&cps->b111, &cps->b111, &cps->b102);
    RwV3dAdd(&cps->b111, &cps->b111, &cps->b012);
    RwV3dAdd(&cps->b111, &cps->b111, &cps->b021);
    RwV3dScale(&cps->b111, &cps->b111, 0.25f);

    RwV3dAdd(&tmp, &cps->b300, &cps->b030);
    RwV3dAdd(&tmp, &tmp, &cps->b003);
    RwV3dScale(&tmp, &tmp, (-1.0f / 6.0f));

    RwV3dAdd(&cps->b111, &cps->b111, &tmp);

    RWRETURNVOID();
}
