15899162844

Committed 26 Jun 2025 10:14AM UTC coverage: 71.088% (+0.004%) from 71.084%

Build # 15899162844

Build Type

Pull #12623

github

Committed by

web-flow

Commit Message

Merge c704a8392 into f5cb024d4

Pull Request Pull Request #12623: gdal raster overview add: add a --overview-src option

Run Details

209 of 244 new or added lines in 5 files covered. (85.66%)

96 existing lines in 44 files now uncovered.

574014 of 807474 relevant lines covered (71.09%)

250815.03 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.24

/alg/delaunay.c

/******************************************************************************
 *
 * 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);
}

1	/******************************************************************************
2	*
3	* Project: GDAL algorithms
4	* Purpose: Delaunay triangulation
5	* Author: Even Rouault, even.rouault at spatialys.com
6	*
7	******************************************************************************
8	* Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>
9	*
10	* SPDX-License-Identifier: MIT
11	****************************************************************************/
12
13	#if defined(__MINGW32__) \|\| defined(__MINGW64__)
14	/* This avoids i586-mingw32msvc/include/direct.h from including libqhull/io.h
15	* ... */
16	#define _DIRECT_H_
17	/* For __MINGW64__ */
18	#define _INC_DIRECT
19	#define _INC_STAT
20	#endif
21
22	#if defined(INTERNAL_QHULL) && !defined(DONT_DEPRECATE_SPRINTF)
23	#define DONT_DEPRECATE_SPRINTF
24	#endif
25
26	#include "cpl_error.h"
27	#include "cpl_conv.h"
28	#include "cpl_string.h"
29	#include "gdal_alg.h"
30
31	#include <stdio.h>
32	#include <stdlib.h>
33	#include <string.h>
34	#include <ctype.h>
35	#include <math.h>
36
37	#if defined(INTERNAL_QHULL) \|\| defined(EXTERNAL_QHULL)
38	#define HAVE_INTERNAL_OR_EXTERNAL_QHULL 1
39	#endif
40
41	#if HAVE_INTERNAL_OR_EXTERNAL_QHULL
42	#ifdef INTERNAL_QHULL
43
44	#include "internal_qhull_headers.h"
45
46	#else /* INTERNAL_QHULL */
47
48	#ifdef _MSC_VER
49	#pragma warning(push)
50	#pragma warning( \
51	disable : 4324) // 'qhT': structure was padded due to alignment specifier
52	#endif
53
54	#if defined(__clang__)
55	#pragma clang diagnostic push
56	#pragma clang diagnostic ignored "-Wdocumentation"
57	#endif
58
59	#include "libqhull_r/libqhull_r.h"
60	#include "libqhull_r/qset_r.h"
61
62	#if defined(__clang__)
63	#pragma clang diagnostic pop
64	#endif
65
66	#ifdef _MSC_VER
67	#pragma warning(pop)
68	#endif
69
70	#endif /* INTERNAL_QHULL */
71
72	#endif /* HAVE_INTERNAL_OR_EXTERNAL_QHULL*/
73
74	/************************************************************************/
75	/* GDALHasTriangulation() */
76	/************************************************************************/
77
78	/** Returns if GDAL is built with Delaunay triangulation support.
79	*
80	* @return TRUE if GDAL is built with Delaunay triangulation support.
81	*
82	* @since GDAL 2.1
83	*/
84	int GDALHasTriangulation()	9✔
85	{
86	#if HAVE_INTERNAL_OR_EXTERNAL_QHULL
87	return TRUE;	9✔
88	#else
89	return FALSE;
90	#endif
91	}
92
93	/************************************************************************/
94	/* GDALTriangulationCreateDelaunay() */
95	/************************************************************************/
96
97	/** Computes a Delaunay triangulation of the passed points
98	*
99	* @param nPoints number of points
100	* @param padfX x coordinates of the points.
101	* @param padfY y coordinates of the points.
102	* @return triangulation that must be freed with GDALTriangulationFree(), or
103	* NULL in case of error.
104	*
105	* @since GDAL 2.1
106	*/
107	GDALTriangulation *GDALTriangulationCreateDelaunay(int nPoints,	10✔
108	const double *padfX,
109	const double *padfY)
110	{
111	#if HAVE_INTERNAL_OR_EXTERNAL_QHULL
112	coordT *points;
113	int i, j;
114	GDALTriangulation *psDT = NULL;	10✔
115	facetT *facet;
116	GDALTriFacet *pasFacets;
117	int panMapQHFacetIdToFacetIdx; / map from QHull facet ID to the index of
118	our GDALTriFacet* array */
119	int curlong, totlong; /* memory remaining after qh_memfreeshort */
120	char *pszTempFilename = NULL;	10✔
121	FILE *fpTemp = NULL;	10✔
122
123	qhT qh_qh;
124	qhT *qh = &qh_qh;	10✔
125
126	QHULL_LIB_CHECK /* Check for compatible library */	10✔
127
128	points = (coordT )VSI_MALLOC2_VERBOSE(sizeof(double) 2, nPoints);	10✔
129	if (points == NULL)	10✔
130	{
131	return NULL;	×
132	}
133	for (i = 0; i < nPoints; i++)	14,711✔
134	{
135	points[2 * i] = padfX[i];	14,701✔
136	points[2 * i + 1] = padfY[i];	14,701✔
137	}
138
139	qh_meminit(qh, NULL);	10✔
140
141	if (CPLTestBoolean(CPLGetConfigOption("QHULL_LOG_TO_TEMP_FILE", "NO")))	10✔
142	{
143	pszTempFilename = CPLStrdup(CPLGenerateTempFilename(NULL));	2✔
144	fpTemp = fopen(pszTempFilename, "wb");	2✔
145	}
146	if (fpTemp == NULL)	10✔
147	fpTemp = stderr;	8✔
148
149	/* d: Delaunay */
150	/* Qbb: scale last coordinate to [0,m] for Delaunay */
151	/* Qc: keep coplanar points with nearest facet */
152	/* Qz: add a point-at-infinity for Delaunay triangulation */
153	/* Qt: triangulated output */
154	int ret = qh_new_qhull(qh, 2, nPoints, points, FALSE /* ismalloc */,	10✔
155	"qhull d Qbb Qc Qz Qt", NULL, fpTemp);
156	if (fpTemp != stderr)	10✔
157	{
158	fclose(fpTemp);	2✔
159	}
160	if (pszTempFilename != NULL)	10✔
161	{
162	VSIUnlink(pszTempFilename);	2✔
163	VSIFree(pszTempFilename);	2✔
164	}
165
166	VSIFree(points);	10✔
167	points = NULL;	10✔
168
169	if (ret != 0)	10✔
170	{
171	CPLError(CE_Failure, CPLE_AppDefined, "Delaunay triangulation failed");	2✔
172	goto end;	2✔
173	}
174
175	/* Establish a map from QHull facet id to the index in our array of
176	* sequential facets */
177	panMapQHFacetIdToFacetIdx =
178	(int *)VSI_MALLOC2_VERBOSE(sizeof(int), qh->facet_id);	8✔
179	if (panMapQHFacetIdToFacetIdx == NULL)	8✔
180	{
181	goto end;	×
182	}
183	memset(panMapQHFacetIdToFacetIdx, 0xFF, sizeof(int) * qh->facet_id);	8✔
184
185	for (j = 0, facet = qh->facet_list; facet != NULL && facet->next != NULL;	29,380✔
186	facet = facet->next)	29,372✔
187	{
188	if (facet->upperdelaunay != qh->UPPERdelaunay)	29,372✔
189	continue;	520✔
190
191	if (qh_setsize(qh, facet->vertices) != 3 \|\|	57,704✔
192	qh_setsize(qh, facet->neighbors) != 3)	28,852✔
193	{
194	CPLError(CE_Failure, CPLE_AppDefined,	×
195	"Triangulation resulted in non triangular facet %d: "
196	"vertices=%d",
197	facet->id, qh_setsize(qh, facet->vertices));
198	VSIFree(panMapQHFacetIdToFacetIdx);	×
199	goto end;	×
200	}
201
202	CPLAssert(facet->id < qh->facet_id);	28,852✔
203	panMapQHFacetIdToFacetIdx[facet->id] = j++;	28,852✔
204	}
205
206	pasFacets = (GDALTriFacet *)VSI_MALLOC2_VERBOSE(j, sizeof(GDALTriFacet));	8✔
207	if (pasFacets == NULL)	8✔
208	{
209	VSIFree(panMapQHFacetIdToFacetIdx);	×
210	goto end;	×
211	}
212
213	psDT = (GDALTriangulation *)CPLCalloc(1, sizeof(GDALTriangulation));	8✔
214	psDT->nFacets = j;	8✔
215	psDT->pasFacets = pasFacets;	8✔
216
217	// Store vertex and neighbor information for each triangle.
218	for (facet = qh->facet_list; facet != NULL && facet->next != NULL;	29,380✔
219	facet = facet->next)	29,372✔
220	{
221	int k;
222	if (facet->upperdelaunay != qh->UPPERdelaunay)	29,372✔
223	continue;	520✔
224	k = panMapQHFacetIdToFacetIdx[facet->id];	28,852✔
225	pasFacets[k].anVertexIdx[0] =	57,704✔
226	qh_pointid(qh, ((vertexT *)facet->vertices->e[0].p)->point);	28,852✔
227	pasFacets[k].anVertexIdx[1] =	57,704✔
228	qh_pointid(qh, ((vertexT *)facet->vertices->e[1].p)->point);	28,852✔
229	pasFacets[k].anVertexIdx[2] =	57,704✔
230	qh_pointid(qh, ((vertexT *)facet->vertices->e[2].p)->point);	28,852✔
231	pasFacets[k].anNeighborIdx[0] =	28,852✔
232	panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[0].p)->id];	28,852✔
233	pasFacets[k].anNeighborIdx[1] =	28,852✔
234	panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[1].p)->id];	28,852✔
235	pasFacets[k].anNeighborIdx[2] =	28,852✔
236	panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[2].p)->id];	28,852✔
237	}
238
239	VSIFree(panMapQHFacetIdToFacetIdx);	8✔
240
241	end:	10✔
242	qh_freeqhull(qh, !qh_ALL);	10✔
243	qh_memfreeshort(qh, &curlong, &totlong);	10✔
244
245	return psDT;	10✔
246	#else /* HAVE_INTERNAL_OR_EXTERNAL_QHULL */
247
248	/* Suppress unused argument warnings. */
249	(void)nPoints;
250	(void)padfX;
251	(void)padfY;
252
253	CPLError(CE_Failure, CPLE_NotSupported,
254	"GDALTriangulationCreateDelaunay() unavailable since GDAL built "
255	"without QHull support");
256	return NULL;
257	#endif /* HAVE_INTERNAL_OR_EXTERNAL_QHULL */
258	}
259
260	/************************************************************************/
261	/* GDALTriangulationFree() */
262	/************************************************************************/
263
264	/** Free a triangulation.
265	*
266	* @param psDT triangulation.
267	* @since GDAL 2.1
268	*/
269	void GDALTriangulationFree(GDALTriangulation *psDT)	10✔
270	{
271	if (psDT)	10✔
272	{
273	VSIFree(psDT->pasFacets);	8✔
274	VSIFree(psDT->pasFacetCoefficients);	8✔
275	VSIFree(psDT);	8✔
276	}
277	}	10✔
278
279	/************************************************************************/
280	/* GDALTriangulationComputeBarycentricCoefficients() */
281	/************************************************************************/
282
283	/** Computes barycentric coefficients for each triangles of the triangulation.
284	*
285	* @param psDT triangulation.
286	* @param padfX x coordinates of the points. Must be identical to the one passed
287	* to GDALTriangulationCreateDelaunay().
288	* @param padfY y coordinates of the points. Must be identical to the one passed
289	* to GDALTriangulationCreateDelaunay().
290	*
291	* @return TRUE in case of success.
292	*
293	* @since GDAL 2.1
294	*/
295	int GDALTriangulationComputeBarycentricCoefficients(GDALTriangulation *psDT,	9✔
296	const double *padfX,
297	const double *padfY)
298	{
299	int i;
300
301	if (psDT->pasFacetCoefficients != NULL)	9✔
302	{
303	return TRUE;	1✔
304	}
305	psDT->pasFacetCoefficients =	8✔
306	(GDALTriBarycentricCoefficients *)VSI_MALLOC2_VERBOSE(	8✔
307	sizeof(GDALTriBarycentricCoefficients), psDT->nFacets);
308	if (psDT->pasFacetCoefficients == NULL)	8✔
309	{
310	return FALSE;	×
311	}
312
313	for (i = 0; i < psDT->nFacets; i++)	28,860✔
314	{
315	GDALTriFacet *psFacet = &(psDT->pasFacets[i]);	28,852✔
316	GDALTriBarycentricCoefficients *psCoeffs =	28,852✔
317	&(psDT->pasFacetCoefficients[i]);	28,852✔
318	double dfX1 = padfX[psFacet->anVertexIdx[0]];	28,852✔
319	double dfY1 = padfY[psFacet->anVertexIdx[0]];	28,852✔
320	double dfX2 = padfX[psFacet->anVertexIdx[1]];	28,852✔
321	double dfY2 = padfY[psFacet->anVertexIdx[1]];	28,852✔
322	double dfX3 = padfX[psFacet->anVertexIdx[2]];	28,852✔
323	double dfY3 = padfY[psFacet->anVertexIdx[2]];	28,852✔
324	/* See https://en.wikipedia.org/wiki/Barycentric_coordinate_system */
325	double dfDenom =	28,852✔
326	(dfY2 - dfY3) * (dfX1 - dfX3) + (dfX3 - dfX2) * (dfY1 - dfY3);	28,852✔
327	if (fabs(dfDenom) < 1e-5)	28,852✔
328	{
329	// Degenerate triangle
330	psCoeffs->dfMul1X = 0.0;	28✔
331	psCoeffs->dfMul1Y = 0.0;	28✔
332	psCoeffs->dfMul2X = 0.0;	28✔
333	psCoeffs->dfMul2Y = 0.0;	28✔
334	psCoeffs->dfCstX = 0.0;	28✔
335	psCoeffs->dfCstY = 0.0;	28✔
336	}
337	else
338	{
339	psCoeffs->dfMul1X = (dfY2 - dfY3) / dfDenom;	28,824✔
340	psCoeffs->dfMul1Y = (dfX3 - dfX2) / dfDenom;	28,824✔
341	psCoeffs->dfMul2X = (dfY3 - dfY1) / dfDenom;	28,824✔
342	psCoeffs->dfMul2Y = (dfX1 - dfX3) / dfDenom;	28,824✔
343	psCoeffs->dfCstX = dfX3;	28,824✔
344	psCoeffs->dfCstY = dfY3;	28,824✔
345	}
346	}
347	return TRUE;	8✔
348	}
349
350	/************************************************************************/
351	/* GDALTriangulationComputeBarycentricCoordinates() */
352	/************************************************************************/
353
354	#define BARYC_COORD_L1(psCoeffs, dfX, dfY) \
355	(psCoeffs->dfMul1X * ((dfX)-psCoeffs->dfCstX) + \
356	psCoeffs->dfMul1Y * ((dfY)-psCoeffs->dfCstY))
357	#define BARYC_COORD_L2(psCoeffs, dfX, dfY) \
358	(psCoeffs->dfMul2X * ((dfX)-psCoeffs->dfCstX) + \
359	psCoeffs->dfMul2Y * ((dfY)-psCoeffs->dfCstY))
360	#define BARYC_COORD_L3(l1, l2) (1 - (l1) - (l2))
361
362	/** Computes the barycentric coordinates of a point.
363	*
364	* @param psDT triangulation.
365	* @param nFacetIdx index of the triangle in the triangulation
366	* @param dfX x coordinate of the point.
367	* @param dfY y coordinate of the point.
368	* @param pdfL1 (output) pointer to the 1st barycentric coordinate.
369	* @param pdfL2 (output) pointer to the 2nd barycentric coordinate.
370	* @param pdfL3 (output) pointer to the 2nd barycentric coordinate.
371	*
372	* @return TRUE in case of success.
373	*
374	* @since GDAL 2.1
375	*/
376
377	int GDALTriangulationComputeBarycentricCoordinates(	63,099✔
378	const GDALTriangulation *psDT, int nFacetIdx, double dfX, double dfY,
379	double pdfL1, double pdfL2, double *pdfL3)
380	{
381	const GDALTriBarycentricCoefficients *psCoeffs;
382	if (psDT->pasFacetCoefficients == NULL)	63,099✔
383	{
384	CPLError(CE_Failure, CPLE_AppDefined,	1✔
385	"GDALTriangulationComputeBarycentricCoefficients() should be "
386	"called before");
387	return FALSE;	1✔
388	}
389	CPLAssert(nFacetIdx >= 0 && nFacetIdx < psDT->nFacets);	63,098✔
390	psCoeffs = &(psDT->pasFacetCoefficients[nFacetIdx]);	73,464✔
391	*pdfL1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);	73,464✔
392	*pdfL2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);	73,464✔
393	pdfL3 = BARYC_COORD_L3(pdfL1, *pdfL2);	73,464✔
394
395	return TRUE;	73,464✔
396	}
397
398	/************************************************************************/
399	/* GDALTriangulationFindFacetBruteForce() */
400	/************************************************************************/
401
402	#define EPS 1e-10
403
404	/** Returns the index of the triangle that contains the point by iterating
405	* over all triangles.
406	*
407	* If the function returns FALSE and *panOutputFacetIdx >= 0, then it means
408	* the point is outside the hull of the triangulation, and *panOutputFacetIdx
409	* is the closest triangle to the point.
410	*
411	* @param psDT triangulation.
412	* @param dfX x coordinate of the point.
413	* @param dfY y coordinate of the point.
414	* @param panOutputFacetIdx (output) pointer to the index of the triangle,
415	* or -1 in case of failure.
416	*
417	* @return index >= 0 of the triangle in case of success, -1 otherwise.
418	*
419	* @since GDAL 2.1
420	*/
421
422	int GDALTriangulationFindFacetBruteForce(const GDALTriangulation *psDT,	10,635✔
423	double dfX, double dfY,
424	int *panOutputFacetIdx)
425	{
426	int nFacetIdx;
427	*panOutputFacetIdx = -1;	10,635✔
428	if (psDT->pasFacetCoefficients == NULL)	10,635✔
429	{
430	CPLError(CE_Failure, CPLE_AppDefined,	1✔
431	"GDALTriangulationComputeBarycentricCoefficients() should be "
432	"called before");
433	return FALSE;	1✔
434	}
435	for (nFacetIdx = 0; nFacetIdx < psDT->nFacets; nFacetIdx++)	261,574✔
436	{
437	double l1, l2, l3;
438	const GDALTriBarycentricCoefficients *psCoeffs =	250,948✔
439	&(psDT->pasFacetCoefficients[nFacetIdx]);	250,948✔
440	if (psCoeffs->dfMul1X == 0.0 && psCoeffs->dfMul2X == 0.0 &&	250,948✔
441	psCoeffs->dfMul1Y == 0.0 && psCoeffs->dfMul2Y == 0.0)	245,718✔
442	{
443	// Degenerate triangle
444	continue;	250,928✔
445	}
446	l1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);	20✔
447	if (l1 < -EPS)	20✔
448	{
449	int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[0];	2✔
450	if (neighbor < 0)	2✔
451	{
452	*panOutputFacetIdx = nFacetIdx;	2✔
453	return FALSE;	2✔
454	}
455	continue;	×
456	}
457	if (l1 > 1 + EPS)	18✔
458	continue;	7✔
459	l2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);	11✔
460	if (l2 < -EPS)	11✔
461	{
462	int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[1];	3✔
463	if (neighbor < 0)	3✔
464	{
465	*panOutputFacetIdx = nFacetIdx;	2✔
466	return FALSE;	2✔
467	}
468	continue;	1✔
469	}
470	if (l2 > 1 + EPS)	8✔
471	continue;	×
472	l3 = BARYC_COORD_L3(l1, l2);	8✔
473	if (l3 < -EPS)	8✔
474	{
475	int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[2];	4✔
476	if (neighbor < 0)	4✔
477	{
478	*panOutputFacetIdx = nFacetIdx;	×
479	return FALSE;	×
480	}
481	continue;	4✔
482	}
483	if (l3 > 1 + EPS)	4✔
484	continue;	×
485	*panOutputFacetIdx = nFacetIdx;	4✔
486	return TRUE;	4✔
487	}
488	return FALSE;	10,626✔
489	}
490
491	/************************************************************************/
492	/* GDALTriangulationFindFacetDirected() */
493	/************************************************************************/
494
495	#define EPS 1e-10
496
497	/** Returns the index of the triangle that contains the point by walking in
498	* the triangulation.
499	*
500	* If the function returns FALSE and *panOutputFacetIdx >= 0, then it means
501	* the point is outside the hull of the triangulation, and *panOutputFacetIdx
502	* is the closest triangle to the point.
503	*
504	* @param psDT triangulation.
505	* @param nFacetIdx index of first triangle to start with.
506	* Must be >= 0 && < psDT->nFacets
507	* @param dfX x coordinate of the point.
508	* @param dfY y coordinate of the point.
509	* @param panOutputFacetIdx (output) pointer to the index of the triangle,
510	* or -1 in case of failure.
511	*
512	* @return TRUE in case of success, FALSE otherwise.
513	*
514	* @since GDAL 2.1
515	*/
516
517	int GDALTriangulationFindFacetDirected(const GDALTriangulation *psDT,	311,838✔
518	int nFacetIdx, double dfX, double dfY,
519	int *panOutputFacetIdx)
520	{
521	#ifdef DEBUG_VERBOSE
522	const int nFacetIdxInitial = nFacetIdx;
523	#endif
524	int k, nIterMax;
525	*panOutputFacetIdx = -1;	311,838✔
526	if (psDT->pasFacetCoefficients == NULL)	311,838✔
527	{
528	CPLError(CE_Failure, CPLE_AppDefined,	1✔
529	"GDALTriangulationComputeBarycentricCoefficients() should be "
530	"called before");
531	return FALSE;	1✔
532	}
533	CPLAssert(nFacetIdx >= 0 && nFacetIdx < psDT->nFacets);	311,837✔
534
535	nIterMax = 2 + psDT->nFacets / 4;	316,893✔
536	for (k = 0; k < nIterMax; k++)	354,305✔
537	{
538	double l1, l2, l3;
539	int bMatch = TRUE;	359,462✔
540	const GDALTriFacet *psFacet = &(psDT->pasFacets[nFacetIdx]);	359,462✔
541	const GDALTriBarycentricCoefficients *psCoeffs =	359,462✔
542	&(psDT->pasFacetCoefficients[nFacetIdx]);	359,462✔
543	if (psCoeffs->dfMul1X == 0.0 && psCoeffs->dfMul2X == 0.0 &&	359,462✔
544	psCoeffs->dfMul1Y == 0.0 && psCoeffs->dfMul2Y == 0.0)	10,547✔
545	{
546	// Degenerate triangle
547	break;	10,649✔
548	}
549	l1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);	348,813✔
550	if (l1 < -EPS)	348,813✔
551	{
552	int neighbor = psFacet->anNeighborIdx[0];	230,664✔
553	if (neighbor < 0)	230,664✔
554	{
555	#ifdef DEBUG_VERBOSE
556	CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",
557	nFacetIdx, k, nFacetIdxInitial);
558	#endif
559	*panOutputFacetIdx = nFacetIdx;	236,798✔
560	return FALSE;	236,798✔
561	}
UNCOV 562	nFacetIdx = neighbor;	×
UNCOV 563	continue;	×
564	}
565	else if (l1 > 1 + EPS)	118,149✔
566	bMatch = FALSE; // outside or degenerate	13,180✔
567
568	l2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);	118,149✔
569	if (l2 < -EPS)	118,149✔
570	{
571	int neighbor = psFacet->anNeighborIdx[1];	21,699✔
572	if (neighbor < 0)	21,699✔
573	{
574	#ifdef DEBUG_VERBOSE
575	CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",
576	nFacetIdx, k, nFacetIdxInitial);
577	#endif
578	*panOutputFacetIdx = nFacetIdx;	27✔
579	return FALSE;	27✔
580	}
581	nFacetIdx = neighbor;	21,672✔
582	continue;	21,672✔
583	}
584	else if (l2 > 1 + EPS)	96,450✔
585	bMatch = FALSE; // outside or degenerate	11,307✔
586
587	l3 = BARYC_COORD_L3(l1, l2);	96,450✔
588	if (l3 < -EPS)	96,450✔
589	{
590	int neighbor = psFacet->anNeighborIdx[2];	21,912✔
591	if (neighbor < 0)	21,912✔
592	{
593	#ifdef DEBUG_VERBOSE
594	CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",
595	nFacetIdx, k, nFacetIdxInitial);
596	#endif
597	*panOutputFacetIdx = nFacetIdx;	38✔
598	return FALSE;	38✔
599	}
600	nFacetIdx = neighbor;	21,874✔
601	continue;	21,874✔
602	}
603	else if (l3 > 1 + EPS)	74,538✔
604	bMatch = FALSE; // outside or degenerate	×
605
606	if (bMatch)	74,538✔
607	{
608	#ifdef DEBUG_VERBOSE
609	CPLDebug("GDAL", "Inside %d in %d iters (initial = %d)", nFacetIdx,
610	k, nFacetIdxInitial);
611	#endif
612	*panOutputFacetIdx = nFacetIdx;	77,808✔
613	return TRUE;	77,808✔
614	}
615	else
616	{
617	break;	×
618	}
619	}
620
621	static int nDebugMsgCount = 0;
622	if (nDebugMsgCount <= 20)	2,222✔
623	{
624	CPLDebug("GDAL", "Using brute force lookup%s",	20✔
625	(nDebugMsgCount == 20)	20✔
626	? " (this debug message will no longer be emitted)"
627	: "");
628	nDebugMsgCount++;	21✔
629	}
630
631	return GDALTriangulationFindFacetBruteForce(psDT, dfX, dfY,	2,223✔
632	panOutputFacetIdx);
633	}

OSGeo / gdal / 15899162844

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous