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 : }
|
