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