/******************************************************************************
 *
 * 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>

CPL_CVSID("$Id: dgnstroke.cpp 1dcdb5c8f617002a6e034d4d98dfdf080f709325 2018-01-14 01:50:09Z Kurt Schwehr $")

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