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 : * @since GDAL 2.1
83 : */
84 4 : int GDALHasTriangulation()
85 : {
86 : #if HAVE_INTERNAL_OR_EXTERNAL_QHULL
87 4 : return TRUE;
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 5 : GDALTriangulation *GDALTriangulationCreateDelaunay(int nPoints,
108 : const double *padfX,
109 : const double *padfY)
110 : {
111 : #if HAVE_INTERNAL_OR_EXTERNAL_QHULL
112 : coordT *points;
113 : int i, j;
114 5 : GDALTriangulation *psDT = NULL;
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 5 : char *pszTempFilename = NULL;
121 5 : FILE *fpTemp = NULL;
122 :
123 : qhT qh_qh;
124 5 : qhT *qh = &qh_qh;
125 :
126 5 : QHULL_LIB_CHECK /* Check for compatible library */
127 :
128 5 : points = (coordT *)VSI_MALLOC2_VERBOSE(sizeof(double) * 2, nPoints);
129 5 : if (points == NULL)
130 : {
131 0 : return NULL;
132 : }
133 14681 : for (i = 0; i < nPoints; i++)
134 : {
135 14676 : points[2 * i] = padfX[i];
136 14676 : points[2 * i + 1] = padfY[i];
137 : }
138 :
139 5 : qh_meminit(qh, NULL);
140 :
141 5 : if (CPLTestBoolean(CPLGetConfigOption("QHULL_LOG_TO_TEMP_FILE", "NO")))
142 : {
143 2 : pszTempFilename = CPLStrdup(CPLGenerateTempFilename(NULL));
144 2 : fpTemp = fopen(pszTempFilename, "wb");
145 : }
146 5 : if (fpTemp == NULL)
147 3 : fpTemp = stderr;
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 5 : int ret = qh_new_qhull(qh, 2, nPoints, points, FALSE /* ismalloc */,
155 : "qhull d Qbb Qc Qz Qt", NULL, fpTemp);
156 5 : if (fpTemp != stderr)
157 : {
158 2 : fclose(fpTemp);
159 : }
160 5 : if (pszTempFilename != NULL)
161 : {
162 2 : VSIUnlink(pszTempFilename);
163 2 : VSIFree(pszTempFilename);
164 : }
165 :
166 5 : VSIFree(points);
167 5 : points = NULL;
168 :
169 5 : if (ret != 0)
170 : {
171 2 : CPLError(CE_Failure, CPLE_AppDefined, "Delaunay triangulation failed");
172 2 : goto end;
173 : }
174 :
175 : /* Establish a map from QHull facet id to the index in our array of
176 : * sequential facets */
177 : panMapQHFacetIdToFacetIdx =
178 3 : (int *)VSI_MALLOC2_VERBOSE(sizeof(int), qh->facet_id);
179 3 : if (panMapQHFacetIdToFacetIdx == NULL)
180 : {
181 0 : goto end;
182 : }
183 3 : memset(panMapQHFacetIdToFacetIdx, 0xFF, sizeof(int) * qh->facet_id);
184 :
185 29335 : for (j = 0, facet = qh->facet_list; facet != NULL && facet->next != NULL;
186 29332 : facet = facet->next)
187 : {
188 29332 : if (facet->upperdelaunay != qh->UPPERdelaunay)
189 500 : continue;
190 :
191 57664 : if (qh_setsize(qh, facet->vertices) != 3 ||
192 28832 : qh_setsize(qh, facet->neighbors) != 3)
193 : {
194 0 : 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 0 : VSIFree(panMapQHFacetIdToFacetIdx);
199 0 : goto end;
200 : }
201 :
202 28832 : CPLAssert(facet->id < qh->facet_id);
203 28832 : panMapQHFacetIdToFacetIdx[facet->id] = j++;
204 : }
205 :
206 3 : pasFacets = (GDALTriFacet *)VSI_MALLOC2_VERBOSE(j, sizeof(GDALTriFacet));
207 3 : if (pasFacets == NULL)
208 : {
209 0 : VSIFree(panMapQHFacetIdToFacetIdx);
210 0 : goto end;
211 : }
212 :
213 3 : psDT = (GDALTriangulation *)CPLCalloc(1, sizeof(GDALTriangulation));
214 3 : psDT->nFacets = j;
215 3 : psDT->pasFacets = pasFacets;
216 :
217 : // Store vertex and neighbor information for each triangle.
218 29335 : for (facet = qh->facet_list; facet != NULL && facet->next != NULL;
219 29332 : facet = facet->next)
220 : {
221 : int k;
222 29332 : if (facet->upperdelaunay != qh->UPPERdelaunay)
223 500 : continue;
224 28832 : k = panMapQHFacetIdToFacetIdx[facet->id];
225 57664 : pasFacets[k].anVertexIdx[0] =
226 28832 : qh_pointid(qh, ((vertexT *)facet->vertices->e[0].p)->point);
227 57664 : pasFacets[k].anVertexIdx[1] =
228 28832 : qh_pointid(qh, ((vertexT *)facet->vertices->e[1].p)->point);
229 57664 : pasFacets[k].anVertexIdx[2] =
230 28832 : qh_pointid(qh, ((vertexT *)facet->vertices->e[2].p)->point);
231 28832 : pasFacets[k].anNeighborIdx[0] =
232 28832 : panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[0].p)->id];
233 28832 : pasFacets[k].anNeighborIdx[1] =
234 28832 : panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[1].p)->id];
235 28832 : pasFacets[k].anNeighborIdx[2] =
236 28832 : panMapQHFacetIdToFacetIdx[((facetT *)facet->neighbors->e[2].p)->id];
237 : }
238 :
239 3 : VSIFree(panMapQHFacetIdToFacetIdx);
240 :
241 5 : end:
242 5 : qh_freeqhull(qh, !qh_ALL);
243 5 : qh_memfreeshort(qh, &curlong, &totlong);
244 :
245 5 : return psDT;
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 5 : void GDALTriangulationFree(GDALTriangulation *psDT)
270 : {
271 5 : if (psDT)
272 : {
273 3 : VSIFree(psDT->pasFacets);
274 3 : VSIFree(psDT->pasFacetCoefficients);
275 3 : VSIFree(psDT);
276 : }
277 5 : }
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 4 : int GDALTriangulationComputeBarycentricCoefficients(GDALTriangulation *psDT,
296 : const double *padfX,
297 : const double *padfY)
298 : {
299 : int i;
300 :
301 4 : if (psDT->pasFacetCoefficients != NULL)
302 : {
303 1 : return TRUE;
304 : }
305 3 : psDT->pasFacetCoefficients =
306 3 : (GDALTriBarycentricCoefficients *)VSI_MALLOC2_VERBOSE(
307 : sizeof(GDALTriBarycentricCoefficients), psDT->nFacets);
308 3 : if (psDT->pasFacetCoefficients == NULL)
309 : {
310 0 : return FALSE;
311 : }
312 :
313 28835 : for (i = 0; i < psDT->nFacets; i++)
314 : {
315 28832 : GDALTriFacet *psFacet = &(psDT->pasFacets[i]);
316 28832 : GDALTriBarycentricCoefficients *psCoeffs =
317 28832 : &(psDT->pasFacetCoefficients[i]);
318 28832 : double dfX1 = padfX[psFacet->anVertexIdx[0]];
319 28832 : double dfY1 = padfY[psFacet->anVertexIdx[0]];
320 28832 : double dfX2 = padfX[psFacet->anVertexIdx[1]];
321 28832 : double dfY2 = padfY[psFacet->anVertexIdx[1]];
322 28832 : double dfX3 = padfX[psFacet->anVertexIdx[2]];
323 28832 : double dfY3 = padfY[psFacet->anVertexIdx[2]];
324 : /* See https://en.wikipedia.org/wiki/Barycentric_coordinate_system */
325 28832 : double dfDenom =
326 28832 : (dfY2 - dfY3) * (dfX1 - dfX3) + (dfX3 - dfX2) * (dfY1 - dfY3);
327 28832 : if (fabs(dfDenom) < 1e-5)
328 : {
329 : // Degenerate triangle
330 28 : psCoeffs->dfMul1X = 0.0;
331 28 : psCoeffs->dfMul1Y = 0.0;
332 28 : psCoeffs->dfMul2X = 0.0;
333 28 : psCoeffs->dfMul2Y = 0.0;
334 28 : psCoeffs->dfCstX = 0.0;
335 28 : psCoeffs->dfCstY = 0.0;
336 : }
337 : else
338 : {
339 28804 : psCoeffs->dfMul1X = (dfY2 - dfY3) / dfDenom;
340 28804 : psCoeffs->dfMul1Y = (dfX3 - dfX2) / dfDenom;
341 28804 : psCoeffs->dfMul2X = (dfY3 - dfY1) / dfDenom;
342 28804 : psCoeffs->dfMul2Y = (dfX1 - dfX3) / dfDenom;
343 28804 : psCoeffs->dfCstX = dfX3;
344 28804 : psCoeffs->dfCstY = dfY3;
345 : }
346 : }
347 3 : return TRUE;
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 13373 : int GDALTriangulationComputeBarycentricCoordinates(
378 : const GDALTriangulation *psDT, int nFacetIdx, double dfX, double dfY,
379 : double *pdfL1, double *pdfL2, double *pdfL3)
380 : {
381 : const GDALTriBarycentricCoefficients *psCoeffs;
382 13373 : if (psDT->pasFacetCoefficients == NULL)
383 : {
384 1 : CPLError(CE_Failure, CPLE_AppDefined,
385 : "GDALTriangulationComputeBarycentricCoefficients() should be "
386 : "called before");
387 1 : return FALSE;
388 : }
389 13372 : CPLAssert(nFacetIdx >= 0 && nFacetIdx < psDT->nFacets);
390 13471 : psCoeffs = &(psDT->pasFacetCoefficients[nFacetIdx]);
391 13471 : *pdfL1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);
392 13471 : *pdfL2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);
393 13471 : *pdfL3 = BARYC_COORD_L3(*pdfL1, *pdfL2);
394 :
395 13471 : return TRUE;
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 10683 : int GDALTriangulationFindFacetBruteForce(const GDALTriangulation *psDT,
423 : double dfX, double dfY,
424 : int *panOutputFacetIdx)
425 : {
426 : int nFacetIdx;
427 10683 : *panOutputFacetIdx = -1;
428 10683 : if (psDT->pasFacetCoefficients == NULL)
429 : {
430 1 : CPLError(CE_Failure, CPLE_AppDefined,
431 : "GDALTriangulationComputeBarycentricCoefficients() should be "
432 : "called before");
433 1 : return FALSE;
434 : }
435 284815 : for (nFacetIdx = 0; nFacetIdx < psDT->nFacets; nFacetIdx++)
436 : {
437 : double l1, l2, l3;
438 274141 : const GDALTriBarycentricCoefficients *psCoeffs =
439 274141 : &(psDT->pasFacetCoefficients[nFacetIdx]);
440 274141 : if (psCoeffs->dfMul1X == 0.0 && psCoeffs->dfMul2X == 0.0 &&
441 270806 : psCoeffs->dfMul1Y == 0.0 && psCoeffs->dfMul2Y == 0.0)
442 : {
443 : // Degenerate triangle
444 274121 : continue;
445 : }
446 20 : l1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);
447 20 : if (l1 < -EPS)
448 : {
449 2 : int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[0];
450 2 : if (neighbor < 0)
451 : {
452 2 : *panOutputFacetIdx = nFacetIdx;
453 2 : return FALSE;
454 : }
455 0 : continue;
456 : }
457 18 : if (l1 > 1 + EPS)
458 7 : continue;
459 11 : l2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);
460 11 : if (l2 < -EPS)
461 : {
462 3 : int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[1];
463 3 : if (neighbor < 0)
464 : {
465 2 : *panOutputFacetIdx = nFacetIdx;
466 2 : return FALSE;
467 : }
468 1 : continue;
469 : }
470 8 : if (l2 > 1 + EPS)
471 0 : continue;
472 8 : l3 = BARYC_COORD_L3(l1, l2);
473 8 : if (l3 < -EPS)
474 : {
475 4 : int neighbor = psDT->pasFacets[nFacetIdx].anNeighborIdx[2];
476 4 : if (neighbor < 0)
477 : {
478 0 : *panOutputFacetIdx = nFacetIdx;
479 0 : return FALSE;
480 : }
481 4 : continue;
482 : }
483 4 : if (l3 > 1 + EPS)
484 0 : continue;
485 4 : *panOutputFacetIdx = nFacetIdx;
486 4 : return TRUE;
487 : }
488 10674 : return FALSE;
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 24803 : int GDALTriangulationFindFacetDirected(const GDALTriangulation *psDT,
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 24803 : *panOutputFacetIdx = -1;
526 24803 : if (psDT->pasFacetCoefficients == NULL)
527 : {
528 1 : CPLError(CE_Failure, CPLE_AppDefined,
529 : "GDALTriangulationComputeBarycentricCoefficients() should be "
530 : "called before");
531 1 : return FALSE;
532 : }
533 24802 : CPLAssert(nFacetIdx >= 0 && nFacetIdx < psDT->nFacets);
534 :
535 24904 : nIterMax = 2 + psDT->nFacets / 4;
536 80204 : for (k = 0; k < nIterMax; k++)
537 : {
538 : double l1, l2, l3;
539 75862 : int bMatch = TRUE;
540 75862 : const GDALTriFacet *psFacet = &(psDT->pasFacets[nFacetIdx]);
541 75862 : const GDALTriBarycentricCoefficients *psCoeffs =
542 75862 : &(psDT->pasFacetCoefficients[nFacetIdx]);
543 75862 : if (psCoeffs->dfMul1X == 0.0 && psCoeffs->dfMul2X == 0.0 &&
544 10670 : psCoeffs->dfMul1Y == 0.0 && psCoeffs->dfMul2Y == 0.0)
545 : {
546 : // Degenerate triangle
547 10668 : break;
548 : }
549 65194 : l1 = BARYC_COORD_L1(psCoeffs, dfX, dfY);
550 65194 : if (l1 < -EPS)
551 : {
552 14904 : int neighbor = psFacet->anNeighborIdx[0];
553 14904 : if (neighbor < 0)
554 : {
555 : #ifdef DEBUG_VERBOSE
556 : CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",
557 : nFacetIdx, k, nFacetIdxInitial);
558 : #endif
559 432 : *panOutputFacetIdx = nFacetIdx;
560 432 : return FALSE;
561 : }
562 14472 : nFacetIdx = neighbor;
563 14472 : continue;
564 : }
565 50290 : else if (l1 > 1 + EPS)
566 12548 : bMatch = FALSE; // outside or degenerate
567 :
568 50290 : l2 = BARYC_COORD_L2(psCoeffs, dfX, dfY);
569 50290 : if (l2 < -EPS)
570 : {
571 20164 : int neighbor = psFacet->anNeighborIdx[1];
572 20164 : if (neighbor < 0)
573 : {
574 : #ifdef DEBUG_VERBOSE
575 : CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",
576 : nFacetIdx, k, nFacetIdxInitial);
577 : #endif
578 27 : *panOutputFacetIdx = nFacetIdx;
579 27 : return FALSE;
580 : }
581 20137 : nFacetIdx = neighbor;
582 20137 : continue;
583 : }
584 30126 : else if (l2 > 1 + EPS)
585 11020 : bMatch = FALSE; // outside or degenerate
586 :
587 30126 : l3 = BARYC_COORD_L3(l1, l2);
588 30126 : if (l3 < -EPS)
589 : {
590 20729 : int neighbor = psFacet->anNeighborIdx[2];
591 20729 : if (neighbor < 0)
592 : {
593 : #ifdef DEBUG_VERBOSE
594 : CPLDebug("GDAL", "Outside %d in %d iters (initial = %d)",
595 : nFacetIdx, k, nFacetIdxInitial);
596 : #endif
597 38 : *panOutputFacetIdx = nFacetIdx;
598 38 : return FALSE;
599 : }
600 20691 : nFacetIdx = neighbor;
601 20691 : continue;
602 : }
603 9397 : else if (l3 > 1 + EPS)
604 0 : bMatch = FALSE; // outside or degenerate
605 :
606 9397 : if (bMatch)
607 : {
608 : #ifdef DEBUG_VERBOSE
609 : CPLDebug("GDAL", "Inside %d in %d iters (initial = %d)", nFacetIdx,
610 : k, nFacetIdxInitial);
611 : #endif
612 13547 : *panOutputFacetIdx = nFacetIdx;
613 13547 : return TRUE;
614 : }
615 : else
616 : {
617 0 : break;
618 : }
619 : }
620 :
621 : static int nDebugMsgCount = 0;
622 10860 : if (nDebugMsgCount <= 20)
623 : {
624 22 : CPLDebug("GDAL", "Using brute force lookup%s",
625 22 : (nDebugMsgCount == 20)
626 : ? " (this debug message will no longer be emitted)"
627 : : "");
628 22 : nDebugMsgCount++;
629 : }
630 :
631 10860 : return GDALTriangulationFindFacetBruteForce(psDT, dfX, dfY,
632 : panOutputFacetIdx);
633 : }
|
