/******************************************************************************

 *

 * Project:  GDAL algorithms

 * Purpose:  Delaunay triangulation

 * Author:   Even Rouault, even.rouault at spatialys.com

 *

 ******************************************************************************

 * Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#if defined(__MINGW32__) || defined(__MINGW64__)

/*  This avoids i586-mingw32msvc/include/direct.h from including libqhull/io.h

 * ... */

#define _DIRECT_H_

/* For __MINGW64__ */

#define _INC_DIRECT

#define _INC_STAT

#endif

#if defined(INTERNAL_QHULL) && !defined(DONT_DEPRECATE_SPRINTF)

#define DONT_DEPRECATE_SPRINTF

#endif

#include "cpl_error.h"

#include "cpl_conv.h"

#include "cpl_string.h"

#include "gdal_alg.h"

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <ctype.h>

#include <math.h>

#if defined(INTERNAL_QHULL) || defined(EXTERNAL_QHULL)

#define HAVE_INTERNAL_OR_EXTERNAL_QHULL 1

#endif

#if HAVE_INTERNAL_OR_EXTERNAL_QHULL

#ifdef INTERNAL_QHULL

#include "internal_qhull_headers.h"

#else /* INTERNAL_QHULL */

#ifdef _MSC_VER

#pragma warning(push)

#pragma warning(                                                               \

    disable : 4324)  // 'qhT': structure was padded due to alignment specifier

#endif

#if defined(__clang__)

#pragma clang diagnostic push

#pragma clang diagnostic ignored "-Wdocumentation"

#endif

#include "libqhull_r/libqhull_r.h"

#include "libqhull_r/qset_r.h"

#if defined(__clang__)

#pragma clang diagnostic pop

#endif

#ifdef _MSC_VER

#pragma warning(pop)

#endif

#endif /* INTERNAL_QHULL */

#endif /* HAVE_INTERNAL_OR_EXTERNAL_QHULL*/

/************************************************************************/

/*                       GDALHasTriangulation()                         */

/************************************************************************/

/** Returns if GDAL is built with Delaunay triangulation support.

 *

 * @return TRUE if GDAL is built with Delaunay triangulation support.

 *

 * @since GDAL 2.1

 */

int GDALHasTriangulation()

{

#if HAVE_INTERNAL_OR_EXTERNAL_QHULL

    return TRUE;

#else

    return FALSE;

#endif

}

/************************************************************************/

/*                   GDALTriangulationCreateDelaunay()                  */

/************************************************************************/

/** Computes a Delaunay triangulation of the passed points

 *

 * @param nPoints number of points

 * @param padfX x coordinates of the points.

 * @param padfY y coordinates of the points.

 * @return triangulation that must be freed with GDALTriangulationFree(), or

 *         NULL in case of error.

 *

 * @since GDAL 2.1

 */

GDALTriangulation *GDALTriangulationCreateDelaunay(int nPoints,

                                                   const double *padfX,

                                                   const double *padfY)

{

#if HAVE_INTERNAL_OR_EXTERNAL_QHULL

    coordT *points;

    int i, j;

    GDALTriangulation *psDT = NULL;

    facetT *facet;

    GDALTriFacet *pasFacets;

    int *panMapQHFacetIdToFacetIdx; /* map from QHull facet ID to the index of

                                       our GDALTriFacet* array */

    int curlong, totlong;           /* memory remaining after qh_memfreeshort */

    char *pszTempFilename = NULL;

    FILE *fpTemp = NULL;

    qhT qh_qh;

    qhT *qh = &qh_qh;

    QHULL_LIB_CHECK /* Check for compatible library */

        points = (coordT *)VSI_MALLOC2_VERBOSE(sizeof(double) * 2, nPoints);

    if (points == NULL)

    {

        return NULL;

    }

    for (i = 0; i < nPoints; i++)

    {

        points[2 * i] = padfX[i];

        points[2 * i + 1] = padfY[i];

    }

    qh_meminit(qh, NULL);

    if (CPLTestBoolean(CPLGetConfigOption("QHULL_LOG_TO_TEMP_FILE", "NO")))

    {

        pszTempFilename = CPLStrdup(CPLGenerateTempFilename(NULL));

        fpTemp = fopen(pszTempFilename, "wb");

    }

    if (fpTemp == NULL)

        fpTemp = stderr;

    /* d: Delaunay */

    /* Qbb: scale last coordinate to [0,m] for Delaunay */

    /* Qc: keep coplanar points with nearest facet */

    /* Qz: add a point-at-infinity for Delaunay triangulation */

    /* Qt: triangulated output */

    int ret = qh_new_qhull(qh, 2, nPoints, points, FALSE /* ismalloc */,

                           "qhull d Qbb Qc Qz Qt", NULL, fpTemp);

    if (fpTemp != stderr)

    {

        fclose(fpTemp);

    }

    if (pszTempFilename != NULL)

    {

        VSIUnlink(pszTempFilename);

        VSIFree(pszTempFilename);

    }

    VSIFree(points);

    points = NULL;

    if (ret != 0)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "Delaunay triangulation failed");

        goto end;

    }

    /* Establish a map from QHull facet id to the index in our array of

     * sequential facets */

    panMapQHFacetIdToFacetIdx =

        (int *)VSI_MALLOC2_VERBOSE(sizeof(int), qh->facet_id);

    if (panMapQHFacetIdToFacetIdx == NULL)

    {

        goto end;

    }

    memset(panMapQHFacetIdToFacetIdx, 0xFF, sizeof(int) * qh->facet_id);

    for (j = 0, facet = qh->facet_list; facet != NULL && facet->next != NULL;

         facet = facet->next)

    {

        if (facet->upperdelaunay != qh->UPPERdelaunay)

            continue;

        if (qh_setsize(qh, facet->vertices) != 3 ||

            qh_setsize(qh, facet->neighbors) != 3)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Triangulation resulted in non triangular facet %d: "

                     "vertices=%d",

                     facet->id, qh_setsize(qh, facet->vertices));

            VSIFree(panMapQHFacetIdToFacetIdx);

            goto end;

        }

        CPLAssert(facet->id < qh->facet_id);

        panMapQHFacetIdToFacetIdx[facet->id] = j++;

    }

    pasFacets = (GDALTriFacet *)VSI_MALLOC2_VERBOSE(j, sizeof(GDALTriFacet));

    if (pasFacets == NULL)

    {

        VSIFree(panMapQHFacetIdToFacetIdx);

        goto end;

    }

    psDT = (GDALTriangulation *)CPLCalloc(1, sizeof(GDALTriangulation));

    psDT->nFacets = j;

    psDT->pasFacets = pasFacets;

    // Store vertex and neighbor information for each triangle.

    for (facet = qh->facet_list; facet != NULL && facet->next != NULL;

         facet = facet->next)

    {

        int k;

        if (facet->upperdelaunay != qh->UPPERdelaunay)

            continue;

        k = panMapQHFacetIdToFacetIdx[facet->id];

        pasFacets[k].anVertexIdx[0] =

            qh_pointid(qh, ((vertexT *)facet->vertices->e[0].p)->point);

        pasFacets[k].anVertexIdx[1] =

            qh_pointid(qh, ((vertexT *)facet->vertices->e[1].p)->point);

        pasFacets[k].anVertexIdx[2] =

            qh_pointid(qh, ((vertexT *)facet->vertices->e[2].p)->point);

        pasFacets[k].anNeighborIdx[0] =

            panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[0].p)->id];

        pasFacets[k].anNeighborIdx[1] =

            panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[1].p)->id];

        pasFacets[k].anNeighborIdx[2] =

            panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[2].p)->id];

    }

    VSIFree(panMapQHFacetIdToFacetIdx);

end:

    qh_freeqhull(qh, !qh_ALL);

    qh_memfreeshort(qh, &curlong, &totlong);

    return psDT;

#else  /* HAVE_INTERNAL_OR_EXTERNAL_QHULL */

    /* Suppress unused argument warnings. */

    (void)nPoints;

    (void)padfX;

    (void)padfY;

    CPLError(CE_Failure, CPLE_NotSupported,

             "GDALTriangulationCreateDelaunay() unavailable since GDAL built "

             "without QHull support");

    return NULL;

#endif /* HAVE_INTERNAL_OR_EXTERNAL_QHULL */

}

/************************************************************************/

/*                       GDALTriangulationFree()                        */

/************************************************************************/

/** Free a triangulation.

 *

 * @param psDT triangulation.

 * @since GDAL 2.1

 */

void GDALTriangulationFree(GDALTriangulation *psDT)

{

    if (psDT)

    {

        VSIFree(psDT->pasFacets);

        VSIFree(psDT->pasFacetCoefficients);

        VSIFree(psDT);

    }

}

/************************************************************************/

/*               GDALTriangulationComputeBarycentricCoefficients()      */

/************************************************************************/

/** Computes barycentric coefficients for each triangles of the triangulation.

 *

 * @param psDT triangulation.

 * @param padfX x coordinates of the points. Must be identical to the one passed

 *              to GDALTriangulationCreateDelaunay().

 * @param padfY y coordinates of the points. Must be identical to the one passed

 *              to GDALTriangulationCreateDelaunay().

 *

 * @return TRUE in case of success.

 *

 * @since GDAL 2.1

 */

int GDALTriangulationComputeBarycentricCoefficients(GDALTriangulation *psDT,

                                                    const double *padfX,

                                                    const double *padfY)

{

    int i;

    if (psDT->pasFacetCoefficients != NULL)

    {

        return TRUE;

    }

    psDT->pasFacetCoefficients =

        (GDALTriBarycentricCoefficients *)VSI_MALLOC2_VERBOSE(

            sizeof(GDALTriBarycentricCoefficients), psDT->nFacets);

    if (psDT->pasFacetCoefficients == NULL)

    {

        return FALSE;

    }

    for (i = 0; i < psDT->nFacets; i++)

    {

        GDALTriFacet *psFacet = &(psDT->pasFacets[i]);

        GDALTriBarycentricCoefficients *psCoeffs =

            &(psDT->pasFacetCoefficients[i]);

        double dfX1 = padfX[psFacet->anVertexIdx[0]];

        double dfY1 = padfY[psFacet->anVertexIdx[0]];

        double dfX2 = padfX[psFacet->anVertexIdx[1]];

        double dfY2 = padfY[psFacet->anVertexIdx[1]];

        double dfX3 = padfX[psFacet->anVertexIdx[2]];

        double dfY3 = padfY[psFacet->anVertexIdx[2]];

        /* See https://en.wikipedia.org/wiki/Barycentric_coordinate_system */

        double dfDenom =

            (dfY2 - dfY3) * (dfX1 - dfX3) + (dfX3 - dfX2) * (dfY1 - dfY3);

        if (fabs(dfDenom) < 1e-5)

        {

            // Degenerate triangle

            psCoeffs->dfMul1X = 0.0;

            psCoeffs->dfMul1Y = 0.0;

            psCoeffs->dfMul2X = 0.0;

            psCoeffs->dfMul2Y = 0.0;

            psCoeffs->dfCstX = 0.0;

            psCoeffs->dfCstY = 0.0;

        }

        else

        {

            psCoeffs->dfMul1X = (dfY2 - dfY3) / dfDenom;

            psCoeffs->dfMul1Y = (dfX3 - dfX2) / dfDenom;

            psCoeffs->dfMul2X = (dfY3 - dfY1) / dfDenom;

            psCoeffs->dfMul2Y = (dfX1 - dfX3) / dfDenom;

            psCoeffs->dfCstX = dfX3;

            psCoeffs->dfCstY = dfY3;

        }

    }

    return TRUE;

}

/************************************************************************/

/*               GDALTriangulationComputeBarycentricCoordinates()       */

/************************************************************************/

#define BARYC_COORD_L1(psCoeffs, dfX, dfY)                                     \

    (psCoeffs->dfMul1X * ((dfX)-psCoeffs->dfCstX) +                            \

     psCoeffs->dfMul1Y * ((dfY)-psCoeffs->dfCstY))

#define BARYC_COORD_L2(psCoeffs, dfX, dfY)                                     \

    (psCoeffs->dfMul2X * ((dfX)-psCoeffs->dfCstX) +                            \

     psCoeffs->dfMul2Y * ((dfY)-psCoeffs->dfCstY))

#define BARYC_COORD_L3(l1, l2) (1 - (l1) - (l2))

/** Computes the barycentric coordinates of a point.

 *

 * @param psDT triangulation.

 * @param nFacetIdx index of the triangle in the triangulation

 * @param dfX x coordinate of the point.

 * @param dfY y coordinate of the point.

 * @param pdfL1 (output) pointer to the 1st barycentric coordinate.

 * @param pdfL2 (output) pointer to the 2nd barycentric coordinate.

 * @param pdfL3 (output) pointer to the 2nd barycentric coordinate.

 *

 * @return TRUE in case of success.

 *

 * @since GDAL 2.1

 */

int GDALTriangulationComputeBarycentricCoordinates(

    const GDALTriangulation *psDT, int nFacetIdx, double dfX, double dfY,

    double *pdfL1, double *pdfL2, double *pdfL3)

{

    const GDALTriBarycentricCoefficients *psCoeffs;

    if (psDT->pasFacetCoefficients == NULL)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "GDALTriangulationComputeBarycentricCoefficients() should be "

                 "called before");

        return FALSE;

    }

    CPLAssert(nFacetIdx >= 0 && nFacetIdx < psDT->nFacets);

    psCoeffs = &(psDT->pasFacetCoefficients[nFacetIdx]);

    *pdfL1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);

    *pdfL2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);

    *pdfL3 = BARYC_COORD_L3(*pdfL1, *pdfL2);

    return TRUE;

}

/************************************************************************/

/*               GDALTriangulationFindFacetBruteForce()                 */

/************************************************************************/

#define EPS 1e-10

/** Returns the index of the triangle that contains the point by iterating

 * over all triangles.

 *

 * If the function returns FALSE and *panOutputFacetIdx >= 0, then it means

 * the point is outside the hull of the triangulation, and *panOutputFacetIdx

 * is the closest triangle to the point.

 *

 * @param psDT triangulation.

 * @param dfX x coordinate of the point.

 * @param dfY y coordinate of the point.

 * @param panOutputFacetIdx (output) pointer to the index of the triangle,

 *                          or -1 in case of failure.

 *

 * @return index >= 0 of the triangle in case of success, -1 otherwise.

 *

 * @since GDAL 2.1

 */

int GDALTriangulationFindFacetBruteForce(const GDALTriangulation *psDT,

                                         double dfX, double dfY,

                                         int *panOutputFacetIdx)

{

    int nFacetIdx;

    *panOutputFacetIdx = -1;

    if (psDT->pasFacetCoefficients == NULL)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "GDALTriangulationComputeBarycentricCoefficients() should be "

                 "called before");

        return FALSE;

    }

    for (nFacetIdx = 0; nFacetIdx < psDT->nFacets; nFacetIdx++)

    {

        double l1, l2, l3;

        const GDALTriBarycentricCoefficients *psCoeffs =

            &(psDT->pasFacetCoefficients[nFacetIdx]);

        if (psCoeffs->dfMul1X == 0.0 && psCoeffs->dfMul2X == 0.0 &&

            psCoeffs->dfMul1Y == 0.0 && psCoeffs->dfMul2Y == 0.0)

        {

            // Degenerate triangle

            continue;

        }

        l1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);

        if (l1 < -EPS)

        {

            int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[0];

            if (neighbor < 0)

            {

                *panOutputFacetIdx = nFacetIdx;

                return FALSE;

            }

            continue;

        }

        if (l1 > 1 + EPS)

            continue;

        l2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);

        if (l2 < -EPS)

        {

            int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[1];

            if (neighbor < 0)

            {

                *panOutputFacetIdx = nFacetIdx;

                return FALSE;

            }

            continue;

        }

        if (l2 > 1 + EPS)

            continue;

        l3 = BARYC_COORD_L3(l1, l2);

        if (l3 < -EPS)

        {

            int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[2];

            if (neighbor < 0)

            {

                *panOutputFacetIdx = nFacetIdx;

                return FALSE;

            }

            continue;

        }

        if (l3 > 1 + EPS)

            continue;

        *panOutputFacetIdx = nFacetIdx;

        return TRUE;

    }

    return FALSE;

}

/************************************************************************/

/*               GDALTriangulationFindFacetDirected()                   */

/************************************************************************/

#define EPS 1e-10

/** Returns the index of the triangle that contains the point by walking in

 * the triangulation.

 *

 * If the function returns FALSE and *panOutputFacetIdx >= 0, then it means

 * the point is outside the hull of the triangulation, and *panOutputFacetIdx

 * is the closest triangle to the point.

 *

 * @param psDT triangulation.

 * @param nFacetIdx index of first triangle to start with.

 *                  Must be >= 0 && < psDT->nFacets

 * @param dfX x coordinate of the point.

 * @param dfY y coordinate of the point.

 * @param panOutputFacetIdx (output) pointer to the index of the triangle,

 *                          or -1 in case of failure.

 *

 * @return TRUE in case of success, FALSE otherwise.

 *

 * @since GDAL 2.1

 */

int GDALTriangulationFindFacetDirected(const GDALTriangulation *psDT,

                                       int nFacetIdx, double dfX, double dfY,

                                       int *panOutputFacetIdx)

{

#ifdef DEBUG_VERBOSE

    const int nFacetIdxInitial = nFacetIdx;

#endif

    int k, nIterMax;

    *panOutputFacetIdx = -1;

    if (psDT->pasFacetCoefficients == NULL)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "GDALTriangulationComputeBarycentricCoefficients() should be "

                 "called before");

        return FALSE;

    }

    CPLAssert(nFacetIdx >= 0 && nFacetIdx < psDT->nFacets);

    nIterMax = 2 + psDT->nFacets / 4;

    for (k = 0; k < nIterMax; k++)

    {

        double l1, l2, l3;

        int bMatch = TRUE;

        const GDALTriFacet *psFacet = &(psDT->pasFacets[nFacetIdx]);

        const GDALTriBarycentricCoefficients *psCoeffs =

            &(psDT->pasFacetCoefficients[nFacetIdx]);

        if (psCoeffs->dfMul1X == 0.0 && psCoeffs->dfMul2X == 0.0 &&

            psCoeffs->dfMul1Y == 0.0 && psCoeffs->dfMul2Y == 0.0)

        {

            // Degenerate triangle

            break;

        }

        l1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);

        if (l1 < -EPS)

        {

            int neighbor = psFacet->anNeighborIdx[0];

            if (neighbor < 0)

            {

#ifdef DEBUG_VERBOSE

                CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",

                         nFacetIdx, k, nFacetIdxInitial);

#endif

                *panOutputFacetIdx = nFacetIdx;

                return FALSE;

            }

            nFacetIdx = neighbor;

            continue;

        }

        else if (l1 > 1 + EPS)

            bMatch = FALSE;  // outside or degenerate

        l2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);

        if (l2 < -EPS)

        {

            int neighbor = psFacet->anNeighborIdx[1];

            if (neighbor < 0)

            {

#ifdef DEBUG_VERBOSE

                CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",

                         nFacetIdx, k, nFacetIdxInitial);

#endif

                *panOutputFacetIdx = nFacetIdx;

                return FALSE;

            }

            nFacetIdx = neighbor;

            continue;

        }

        else if (l2 > 1 + EPS)

            bMatch = FALSE;  // outside or degenerate

        l3 = BARYC_COORD_L3(l1, l2);

        if (l3 < -EPS)

        {

            int neighbor = psFacet->anNeighborIdx[2];

            if (neighbor < 0)

            {

#ifdef DEBUG_VERBOSE

                CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",

                         nFacetIdx, k, nFacetIdxInitial);

#endif

                *panOutputFacetIdx = nFacetIdx;

                return FALSE;

            }

            nFacetIdx = neighbor;

            continue;

        }

        else if (l3 > 1 + EPS)

            bMatch = FALSE;  // outside or degenerate

        if (bMatch)

        {

#ifdef DEBUG_VERBOSE

            CPLDebug("GDAL", "Inside %d in %d iters (initial = %d)", nFacetIdx,

                     k, nFacetIdxInitial);

#endif

            *panOutputFacetIdx = nFacetIdx;

            return TRUE;

        }

        else

        {

            break;

        }

    }

    static int nDebugMsgCount = 0;

    if (nDebugMsgCount <= 20)

    {

        CPLDebug("GDAL", "Using brute force lookup%s",

                 (nDebugMsgCount == 20)

                     ? " (this debug message will no longer be emitted)"

                     : "");

        nDebugMsgCount++;

    }

    return GDALTriangulationFindFacetBruteForce(psDT, dfX, dfY,

                                                panOutputFacetIdx);

}