LCOV - code coverage report
Current view: top level - alg - delaunay.c (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 204 219 93.2 %
Date: 2025-10-22 13:51:22 Functions: 7 7 100.0 %

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

Generated by: LCOV version 1.14