/******************************************************************************
 *
 * Project:  Microstation DGN Access Library
 * Purpose:  Code to stroke Arcs/Ellipses into polylines.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2001, Avenza Systems Inc, http://www.avenza.com/
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "dgnlibp.h"
#include <cmath>

constexpr double DEG_TO_RAD = M_PI / 180.0;

/************************************************************************/
/*                         ComputePointOnArc()                          */
/************************************************************************/

static void ComputePointOnArc2D(double dfPrimary, double dfSecondary,
                                double dfAxisRotation, double dfAngle,
                                double *pdfX, double *pdfY)

{
    // dfAxisRotation and dfAngle are supposed to be in Radians
    const double dfCosRotation = cos(dfAxisRotation);
    const double dfSinRotation = sin(dfAxisRotation);
    const double dfEllipseX = dfPrimary * cos(dfAngle);
    const double dfEllipseY = dfSecondary * sin(dfAngle);

    *pdfX = dfEllipseX * dfCosRotation - dfEllipseY * dfSinRotation;
    *pdfY = dfEllipseX * dfSinRotation + dfEllipseY * dfCosRotation;
}

/************************************************************************/
/*                            DGNStrokeArc()                            */
/************************************************************************/

/**
 * Generate a polyline approximation of an arc.
 *
 * Produce a series of equidistant (actually equi-angle) points along
 * an arc.  Currently this only works for 2D arcs (and ellipses).
 *
 * @param hFile the DGN file to which the arc belongs (currently not used).
 * @param psArc the arc to be approximated.
 * @param nPoints the number of points to use to approximate the arc.
 * @param pasPoints the array of points into which to put the results.
 * There must be room for at least nPoints points.
 *
 * @return TRUE on success or FALSE on failure.
 */

int DGNStrokeArc(CPL_UNUSED DGNHandle hFile, DGNElemArc *psArc, int nPoints,
                 DGNPoint *pasPoints)
{
    if (nPoints < 2)
        return FALSE;

    if (psArc->primary_axis == 0.0 || psArc->secondary_axis == 0.0)
    {
        CPLError(CE_Warning, CPLE_AppDefined,
                 "Zero primary or secondary axis in DGNStrokeArc().");
        return FALSE;
    }

    const double dfAngleStep = psArc->sweepang / (nPoints - 1);
    for (int i = 0; i < nPoints; i++)
    {
        const double dfAngle = (psArc->startang + dfAngleStep * i) * DEG_TO_RAD;

        ComputePointOnArc2D(psArc->primary_axis, psArc->secondary_axis,
                            psArc->rotation * DEG_TO_RAD, dfAngle,
                            &(pasPoints[i].x), &(pasPoints[i].y));
        pasPoints[i].x += psArc->origin.x;
        pasPoints[i].y += psArc->origin.y;
        pasPoints[i].z = psArc->origin.z;
    }

    return TRUE;
}

/************************************************************************/
/*                           DGNStrokeCurve()                           */
/************************************************************************/

/**
 * Generate a polyline approximation of an curve.
 *
 * Produce a series of equidistant points along a microstation curve element.
 * Currently this only works for 2D.
 *
 * @param hFile the DGN file to which the arc belongs (currently not used).
 * @param psCurve the curve to be approximated.
 * @param nPoints the number of points to use to approximate the curve.
 * @param pasPoints the array of points into which to put the results.
 * There must be room for at least nPoints points.
 *
 * @return TRUE on success or FALSE on failure.
 */

int DGNStrokeCurve(CPL_UNUSED DGNHandle hFile, DGNElemMultiPoint *psCurve,
                   int nPoints, DGNPoint *pasPoints)
{
    const int nDGNPoints = psCurve->num_vertices;

    if (nDGNPoints < 6)
        return FALSE;

    if (nPoints < nDGNPoints - 4)
        return FALSE;

    /* -------------------------------------------------------------------- */
    /*      Compute the Compute the slopes/distances of the segments.       */
    /* -------------------------------------------------------------------- */
    double *padfMx =
        static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
    double *padfMy =
        static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
    double *padfD =
        static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
    double *padfTx =
        static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));
    double *padfTy =
        static_cast<double *>(CPLMalloc(sizeof(double) * nDGNPoints));

    double dfTotalD = 0.0;

    DGNPoint *pasDGNPoints = psCurve->vertices;

    for (int k = 0; k < nDGNPoints - 1; k++)
    {
        /* coverity[overrun-local] */
        padfD[k] = sqrt((pasDGNPoints[k + 1].x - pasDGNPoints[k].x) *
                            (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) +
                        (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) *
                            (pasDGNPoints[k + 1].y - pasDGNPoints[k].y));
        if (padfD[k] == 0.0)
        {
            padfD[k] = 0.0001;
            padfMx[k] = 0.0;
            padfMy[k] = 0.0;
        }
        else
        {
            padfMx[k] = (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k];
            padfMy[k] = (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k];
        }

        if (k > 1 && k < nDGNPoints - 3)
            dfTotalD += padfD[k];
    }

    /* -------------------------------------------------------------------- */
    /*      Compute the Tx, and Ty coefficients for each segment.           */
    /* -------------------------------------------------------------------- */
    for (int k = 2; k < nDGNPoints - 2; k++)
    {
        if (fabs(padfMx[k + 1] - padfMx[k]) == 0.0 &&
            fabs(padfMx[k - 1] - padfMx[k - 2]) == 0.0)
        {
            padfTx[k] = (padfMx[k] + padfMx[k - 1]) / 2;
        }
        else
        {
            padfTx[k] = (padfMx[k - 1] * fabs(padfMx[k + 1] - padfMx[k]) +
                         padfMx[k] * fabs(padfMx[k - 1] - padfMx[k - 2])) /
                        (std::abs(padfMx[k + 1] - padfMx[k]) +
                         std::abs(padfMx[k - 1] - padfMx[k - 2]));
        }

        if (fabs(padfMy[k + 1] - padfMy[k]) == 0.0 &&
            fabs(padfMy[k - 1] - padfMy[k - 2]) == 0.0)
        {
            padfTy[k] = (padfMy[k] + padfMy[k - 1]) / 2;
        }
        else
        {
            padfTy[k] = (padfMy[k - 1] * fabs(padfMy[k + 1] - padfMy[k]) +
                         padfMy[k] * fabs(padfMy[k - 1] - padfMy[k - 2])) /
                        (std::abs(padfMy[k + 1] - padfMy[k]) +
                         std::abs(padfMy[k - 1] - padfMy[k - 2]));
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Determine a step size in D.  We scale things so that we have    */
    /*      roughly equidistant steps in D, but assume we also want to      */
    /*      include every node along the way.                               */
    /* -------------------------------------------------------------------- */
    double dfStepSize = dfTotalD / (nPoints - (nDGNPoints - 4) - 1);

    /* ==================================================================== */
    /*      Process each of the segments.                                   */
    /* ==================================================================== */
    double dfD = dfStepSize;
    int iOutPoint = 0;

    for (int k = 2; k < nDGNPoints - 3; k++)
    {
        /* --------------------------------------------------------------------
         */
        /*      Compute the "x" coefficients for this segment. */
        /* --------------------------------------------------------------------
         */
        const double dfCx = padfTx[k];
        const double dfBx =
            (3.0 * (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k] -
             2.0 * padfTx[k] - padfTx[k + 1]) /
            padfD[k];
        const double dfAx =
            (padfTx[k] + padfTx[k + 1] -
             2 * (pasDGNPoints[k + 1].x - pasDGNPoints[k].x) / padfD[k]) /
            (padfD[k] * padfD[k]);

        /* --------------------------------------------------------------------
         */
        /*      Compute the Y coefficients for this segment. */
        /* --------------------------------------------------------------------
         */
        const double dfCy = padfTy[k];
        const double dfBy =
            (3.0 * (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k] -
             2.0 * padfTy[k] - padfTy[k + 1]) /
            padfD[k];
        const double dfAy =
            (padfTy[k] + padfTy[k + 1] -
             2 * (pasDGNPoints[k + 1].y - pasDGNPoints[k].y) / padfD[k]) /
            (padfD[k] * padfD[k]);

        /* --------------------------------------------------------------------
         */
        /*      Add the start point for this segment. */
        /* --------------------------------------------------------------------
         */
        pasPoints[iOutPoint].x = pasDGNPoints[k].x;
        pasPoints[iOutPoint].y = pasDGNPoints[k].y;
        pasPoints[iOutPoint].z = 0.0;
        iOutPoint++;

        /* --------------------------------------------------------------------
         */
        /*      Step along, adding intermediate points. */
        /* --------------------------------------------------------------------
         */
        while (dfD < padfD[k] && iOutPoint < nPoints - (nDGNPoints - k - 1))
        {
            pasPoints[iOutPoint].x = dfAx * dfD * dfD * dfD + dfBx * dfD * dfD +
                                     dfCx * dfD + pasDGNPoints[k].x;
            pasPoints[iOutPoint].y = dfAy * dfD * dfD * dfD + dfBy * dfD * dfD +
                                     dfCy * dfD + pasDGNPoints[k].y;
            pasPoints[iOutPoint].z = 0.0;
            iOutPoint++;

            dfD += dfStepSize;
        }

        dfD -= padfD[k];
    }

    /* -------------------------------------------------------------------- */
    /*      Add the start point for this segment.                           */
    /* -------------------------------------------------------------------- */
    while (iOutPoint < nPoints)
    {
        pasPoints[iOutPoint].x = pasDGNPoints[nDGNPoints - 3].x;
        pasPoints[iOutPoint].y = pasDGNPoints[nDGNPoints - 3].y;
        pasPoints[iOutPoint].z = 0.0;
        iOutPoint++;
    }

    /* -------------------------------------------------------------------- */
    /*      Cleanup.                                                        */
    /* -------------------------------------------------------------------- */
    CPLFree(padfMx);
    CPLFree(padfMy);
    CPLFree(padfD);
    CPLFree(padfTx);
    CPLFree(padfTy);

    return TRUE;
}

/************************************************************************/
/*                                main()                                */
/*                                                                      */
/*      test mainline                                                   */
/************************************************************************/
#ifdef notdef
int main(int argc, char **argv)

{
    if (argc != 5)
    {
        printf(  // ok
            "Usage: stroke primary_axis secondary_axis axis_rotation angle\n");
        exit(1);
    }

    const double dfPrimary = CPLAtof(argv[1]);
    const double dfSecondary = CPLAtof(argv[2]);
    const double dfAxisRotation = CPLAtof(argv[3]) / 180 * M_PI;
    const double dfAngle = CPLAtof(argv[4]) / 180 * M_PI;

    double dfX = 0.0;
    double dfY = 0.0;
    ComputePointOnArc2D(dfPrimary, dfSecondary, dfAxisRotation, dfAngle, &dfX,
                        &dfY);

    printf("X=%.2f, Y=%.2f\n", dfX, dfY);  // ok

    return 0;
}

#endif