1 /* -*- coding: utf-8 -*- */
2 /**
3 \ingroup geom
4 @{
5 \file dpeuckerp.c
6 \brief Definición de funciones para el aligerado de polilíneas en el plano
7 basadas en el algoritmo de Douglas-Peucker.
8 \author José Luis García Pallero, jgpallero@gmail.com
9 \note Este fichero contiene funciones paralelizadas con OpenMP.
10 \date 04 de julio de 2011
11 \copyright
12 Copyright (c) 2011-2020, José Luis García Pallero. All rights reserved.
13 \par
14 Redistribution and use in source and binary forms, with or without modification,
15 are permitted provided that the following conditions are met:
16 \par
17 - Redistributions of source code must retain the above copyright notice, this
18 list of conditions and the following disclaimer.
19 - Redistributions in binary form must reproduce the above copyright notice, this
20 list of conditions and the following disclaimer in the documentation and/or
21 other materials provided with the distribution.
22 - Neither the name of the copyright holders nor the names of its contributors
23 may be used to endorse or promote products derived from this software without
24 specific prior written permission.
25 \par
26 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
27 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
28 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
29 DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDER BE LIABLE FOR ANY DIRECT,
30 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
31 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
34 OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
35 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 */
37 /******************************************************************************/
38 /******************************************************************************/
39 #include"libgeoc/dpeuckerp.h"
40 /******************************************************************************/
41 /******************************************************************************/
GeocParOmpDpeuckerp(char version[])42 int GeocParOmpDpeuckerp(char version[])
43 {
44 //comprobamos si hay paralelización
45 #if defined(_OPENMP)
46 //comprobamos si hay que extraer versión
47 if(version!=NULL)
48 {
49 //calculamos la versión
50 VersionOpenMP(_OPENMP,version);
51 }
52 //salimos de la función
53 return 1;
54 #else
55 if(version!=NULL)
56 {
57 //utilizamos la variable version para que no dé warming al compilar
58 strcpy(version,"");
59 }
60 //salimos de la función
61 return 0;
62 #endif
63 }
64 /******************************************************************************/
65 /******************************************************************************/
DouglasPeuckerOriginalPlano(const double * x,const double * y,const size_t nPtos,const size_t incX,const size_t incY,const double tol,const size_t posIni,const size_t posFin,char * usados)66 void DouglasPeuckerOriginalPlano(const double* x,
67 const double* y,
68 const size_t nPtos,
69 const size_t incX,
70 const size_t incY,
71 const double tol,
72 const size_t posIni,
73 const size_t posFin,
74 char* usados)
75 {
76 //índice para recorrer bucles
77 size_t i=0;
78 //distancia de los puntos al segmento base
79 double dist=0.0;
80 //valor absoluto de la tolerancia
81 double atol=fabs(tol);
82 //posición de la distancia máxima
83 size_t pos=0;
84 ////////////////////////////////////////////////////////////////////////////
85 ////////////////////////////////////////////////////////////////////////////
86 //comprobamos casos especiales
87 if((nPtos<=2)||(atol==0.0))
88 {
89 //se usan todos los puntos
90 for(i=0;i<nPtos;i++)
91 {
92 usados[i] = 1;
93 }
94 //salimos de la función
95 return;
96 }
97 ////////////////////////////////////////////////////////////////////////////
98 ////////////////////////////////////////////////////////////////////////////
99 //sólo continuamos si los puntos extremos no están seguidos
100 if(posIni!=(posFin-1))
101 {
102 //calculamos la distancia máxima de los puntos de trabajo al segmento
103 //formado por los extremos
104 dist = DouglasPeuckerDistMaxPlano(x,y,incX,incY,posIni,posFin,&pos);
105 //comprobamos si esa distancia está fuera de tolerancia
106 if(dist>atol)
107 {
108 //indicamos que el punto se usa
109 usados[pos] = 1;
110 //paralelización con OpenMP
111 #if defined(_OPENMP)
112 #pragma omp parallel sections default(none) \
113 shared(nPtos,incX,x,incY,y,posIni,posFin,pos,tol,usados)
114 #endif
115 {
116 #if defined(_OPENMP)
117 #pragma omp section
118 #endif
119 //aplicamos el algoritmo a la parte anterior al punto de trabajo
120 DouglasPeuckerOriginalPlano(x,y,nPtos,incX,incY,tol,posIni,pos,
121 usados);
122 #if defined(_OPENMP)
123 #pragma omp section
124 #endif
125 //aplicamos el algoritmo a la parte posterior al punto de trabajo
126 DouglasPeuckerOriginalPlano(x,y,nPtos,incX,incY,tol,pos,posFin,
127 usados);
128 } // --> fin del #pragma omp parallel sections
129 }
130 }
131 ////////////////////////////////////////////////////////////////////////////
132 ////////////////////////////////////////////////////////////////////////////
133 //salimos de la función
134 return;
135 }
136 /******************************************************************************/
137 /******************************************************************************/
DouglasPeuckerRobustoPlano(const double * x,const double * y,const size_t nPtos,const size_t incX,const size_t incY,const double tol,const int paralelizaTol,const enum GEOC_DPEUCKER_ROBUSTO robusto,const size_t nSegRobOrig,const size_t nSegRobAuto,size_t * nPtosSal)138 size_t* DouglasPeuckerRobustoPlano(const double* x,
139 const double* y,
140 const size_t nPtos,
141 const size_t incX,
142 const size_t incY,
143 const double tol,
144 const int paralelizaTol,
145 const enum GEOC_DPEUCKER_ROBUSTO robusto,
146 const size_t nSegRobOrig,
147 const size_t nSegRobAuto,
148 size_t* nPtosSal)
149 {
150 //índices para recorrer bucles
151 size_t i=0,j=0;
152 //valor absoluto de la tolerancia
153 double atol=fabs(tol);
154 //variable indicadora de punto en tolerancia
155 int entol=0;
156 //identificador de caso especial
157 int hayCasoEspecial=0;
158 //identificadores de utilización de algoritmos semi robustos
159 int robOrig=0,robAuto=0;
160 //número de elementos de trabajo internos del vector de salida
161 size_t nElem=0;
162 //identificador de paralelización
163 int paraleliza=0;
164 //vector de salida
165 size_t* sal=NULL;
166 ////////////////////////////////////////////////////////////////////////////
167 ////////////////////////////////////////////////////////////////////////////
168 //comprobamos casos especiales
169 sal = CasosEspecialesAligeraPolilinea(x,y,nPtos,incX,incY,atol,nPtosSal,
170 &hayCasoEspecial);
171 //comprobamos si ha habido algún caso especial
172 if(hayCasoEspecial)
173 {
174 //comprobamos si ha ocurrido algún error de asignación de memoria
175 if(nPtos&&(sal==NULL))
176 {
177 //mensaje de error
178 GEOC_ERROR("Error de asignación de memoria");
179 //salimos de la función
180 return NULL;
181 }
182 else
183 {
184 //salimos de la función
185 return sal;
186 }
187 }
188 ////////////////////////////////////////////////////////////////////////////
189 ////////////////////////////////////////////////////////////////////////////
190 //comprobamos si se utiliza algoritmo de intersección con línea original
191 if((robusto==GeocDPeuckerRobSi)||(robusto==GeocDPeuckerRobOrig))
192 {
193 robOrig = 1;
194 }
195 //comprobamos si se utiliza algoritmo de intersección con línea generada
196 if((robusto==GeocDPeuckerRobSi)||(robusto==GeocDPeuckerRobAuto))
197 {
198 robAuto = 1;
199 }
200 ////////////////////////////////////////////////////////////////////////////
201 ////////////////////////////////////////////////////////////////////////////
202 #if defined(_OPENMP)
203 //comprobamos si hay más de un procesador
204 if(omp_get_num_procs()>1)
205 {
206 //indicamos que hay paralelización
207 paraleliza = 1;
208 }
209 #endif
210 ////////////////////////////////////////////////////////////////////////////
211 ////////////////////////////////////////////////////////////////////////////
212 //inicializamos el indicador interno de tamaño del vector
213 nElem = GEOC_DPEUCKER_BUFFER_PTOS;
214 //asignamos memoria para el vector de salida
215 sal = (size_t*)malloc(nElem*sizeof(size_t));
216 //comprobamos los posibles errores
217 if(sal==NULL)
218 {
219 //mensaje de error
220 GEOC_ERROR("Error de asignación de memoria");
221 //salimos de la función
222 return NULL;
223 }
224 //indicamos que el primer punto siempre se usa
225 *nPtosSal = 1;
226 sal[0] = 0;
227 //puntos de trabajo para iniciar los cálculos
228 i = 0;
229 j = 2;
230 ////////////////////////////////////////////////////////////////////////////
231 ////////////////////////////////////////////////////////////////////////////
232 //entramos en un bucle mientras no hayamos llegado hasta el último punto
233 while(j<nPtos)
234 {
235 //comprobamos si los puntos intermedios están en tolerancia
236 //sólo paralelizamos si el número de puntos intermedios es mayor que uno
237 if(paraleliza&¶lelizaTol&&((j-i-1)>1))
238 {
239 //aplicamos el algoritmo en paralelo
240 entol = DouglasPeuckerPuntosEnTolPlanoOMP(x,y,incX,incY,atol,i,j,
241 i+1,j-1);
242 }
243 else
244 {
245 //aplicamos el algoritmo en serie
246 entol = DouglasPeuckerPuntosEnTolPlanoSerie(x,y,incX,incY,atol,i,j,
247 i+1,j-1);
248 }
249 //comprobamos si todos los puntos están en tolerancia
250 if(entol)
251 {
252 //pasamos al siguiente punto como extremo del segmento
253 j++;
254 }
255 else
256 {
257 //el punto final será el anterior al actual, ya que con el actual
258 //hay al menos un vértice fuera de tolerancia y con el anterior se
259 //comprobó en el paso previo del bucle que no había ningún vértice
260 //fuera
261 j--;
262 //aplicación del algoritmo de intersección con puntos originales
263 if(robOrig)
264 {
265 //aplicamos el algoritmo
266 DouglasPeuckerRobIntersecOrigPlano(x,y,nPtos,incX,incY,
267 nSegRobOrig,i,&j);
268 }
269 //aplicación del algoritmo de auto intersección
270 if(robAuto)
271 {
272 //aplicamos el algoritmo
273 DouglasPeuckerRobAutoIntersecPlano(x,y,incX,incY,i,&j,sal,
274 *nPtosSal,nSegRobAuto);
275 }
276 //añadimos al contador el nuevo punto
277 (*nPtosSal)++;
278 //comprobamos si hay que reasignar memoria
279 if((*nPtosSal)>nElem)
280 {
281 //añadimos otro grupo de puntos
282 nElem += GEOC_DPEUCKER_BUFFER_PTOS;
283 //asignamos memoria para el vector de salida
284 sal = (size_t*)realloc(sal,nElem*sizeof(size_t));
285 //comprobamos los posibles errores
286 if(sal==NULL)
287 {
288 //mensaje de error
289 GEOC_ERROR("Error de asignación de memoria");
290 //salimos de la función
291 return NULL;
292 }
293 }
294 //añadimos el punto al vector de salida
295 sal[(*nPtosSal)-1] = j;
296 //actualizamos los índices de los puntos de trabajo
297 i = j;
298 j = i+2;
299 }
300 }
301 ////////////////////////////////////////////////////////////////////////////
302 ////////////////////////////////////////////////////////////////////////////
303 //comprobamos si hay que añadir el último punto
304 if((sal[(*nPtosSal)-1]!=(nPtos-1))&&
305 ((x[sal[(*nPtosSal)-1]*incX]!=x[(nPtos-1)*incX])||
306 (y[sal[(*nPtosSal)-1]*incY]!=y[(nPtos-1)*incY])))
307 {
308 //añadimos al contador el último punto
309 (*nPtosSal)++;
310 //comprobamos si hay que reasignar memoria
311 if((*nPtosSal)>nElem)
312 {
313 //añadimos otro grupo de puntos
314 nElem += GEOC_DPEUCKER_BUFFER_PTOS;
315 //asignamos memoria para el vector de salida
316 sal = (size_t*)realloc(sal,nElem*sizeof(size_t));
317 //comprobamos los posibles errores
318 if(sal==NULL)
319 {
320 //mensaje de error
321 GEOC_ERROR("Error de asignación de memoria");
322 //salimos de la función
323 return NULL;
324 }
325 }
326 //asignamos el último punto
327 sal[(*nPtosSal)-1] = nPtos-1;
328 }
329 ////////////////////////////////////////////////////////////////////////////
330 ////////////////////////////////////////////////////////////////////////////
331 //comprobamos si el vector de salida tiene demasiada memoria asignada
332 if(nElem>(*nPtosSal))
333 {
334 //ajustamos el tamaño del vector de salida
335 sal = (size_t*)realloc(sal,(*nPtosSal)*sizeof(size_t));
336 //comprobamos los posibles errores
337 if(sal==NULL)
338 {
339 //mensaje de error
340 GEOC_ERROR("Error de asignación de memoria");
341 //salimos de la función
342 return NULL;
343 }
344 }
345 ////////////////////////////////////////////////////////////////////////////
346 ////////////////////////////////////////////////////////////////////////////
347 //salimos de la función
348 return sal;
349 }
350 /******************************************************************************/
351 /******************************************************************************/
DouglasPeuckerPuntosEnTolPlanoOMP(const double * x,const double * y,const size_t incX,const size_t incY,const double tol,const size_t posBaseIni,const size_t posBaseFin,const size_t posPtoIni,const size_t posPtoFin)352 int DouglasPeuckerPuntosEnTolPlanoOMP(const double* x,
353 const double* y,
354 const size_t incX,
355 const size_t incY,
356 const double tol,
357 const size_t posBaseIni,
358 const size_t posBaseFin,
359 const size_t posPtoIni,
360 const size_t posPtoFin)
361 {
362 //índice para recorrer bucles
363 size_t i=0;
364 //valor absoluto de la tolerancia
365 double atol=fabs(tol);
366 //identificador de punto fuera de tolerancia
367 int ftol=0;
368 //coordenadas de los vértices de trabajo
369 double xIni=0.0,yIni=0.0,xFin=0.0,yFin=0.0,xTrab=0.0,yTrab=0.0;
370 //longitud de la base y parámetros de rotación para el plano
371 double dx=0.0,sA=0.0,cA=0.0;
372 //distancia calculada
373 double dist=0.0;
374 //variable de salida
375 int entol=1;
376 ////////////////////////////////////////////////////////////////////////////
377 ////////////////////////////////////////////////////////////////////////////
378 //comprobamos salida rápida
379 if((posBaseIni+1)>=posBaseFin)
380 {
381 //salimos de la función
382 return entol;
383 }
384 ////////////////////////////////////////////////////////////////////////////
385 ////////////////////////////////////////////////////////////////////////////
386 //coordenadas del primer punto de la base
387 xIni = x[posBaseIni*incX];
388 yIni = y[posBaseIni*incY];
389 //coordenadas del segundo punto de la base, referidas al primero
390 xFin = x[posBaseFin*incX]-xIni;
391 yFin = y[posBaseFin*incY]-yIni;
392 //calculamos la longitud de la base y la rotación
393 DouglasPeuckerParamRotaBase(xFin,yFin,&sA,&cA,&dx);
394 ////////////////////////////////////////////////////////////////////////////
395 ////////////////////////////////////////////////////////////////////////////
396 //paralelización con OpenMP
397 #if defined(_OPENMP)
398 #pragma omp parallel for default(none) \
399 shared(posPtoIni,posPtoFin,incX,x,incY,y,xIni,yIni,dx,sA,cA,atol) \
400 private(i,xTrab,yTrab,dist) \
401 reduction(+:ftol)
402 #endif
403 //recorremos los puntos de trabajo
404 for(i=posPtoIni;i<=posPtoFin;i++)
405 {
406 //sólo calculo si no se ha encontrado ningún punto fuera de tolerancia
407 //en este hilo
408 if(!ftol)
409 {
410 //extraemos las coordenadas del vértice de trabajo y las referimos
411 //al punto inicial de la base
412 xTrab = x[i*incX]-xIni;
413 yTrab = y[i*incY]-yIni;
414 //calculamos la distancia del punto a la base
415 dist = DouglasPeuckerDistMaxPlanoAux(dx,sA,cA,xTrab,yTrab);
416 //comprobamos si está fuera de tolerancia
417 if(dist>atol)
418 {
419 //aumentamos el indicador de fuera de tolerancia
420 ftol++;
421 }
422 }
423 } // --> fin del #pragma omp parallel for
424 ////////////////////////////////////////////////////////////////////////////
425 ////////////////////////////////////////////////////////////////////////////
426 //comprobamos si hay algún punto fuera de tolerancia
427 if(ftol)
428 {
429 //indicamos que hay algún punto que no está en tolerancia
430 entol = 0;
431 }
432 ////////////////////////////////////////////////////////////////////////////
433 ////////////////////////////////////////////////////////////////////////////
434 //salimos de la función
435 return entol;
436 }
437 /******************************************************************************/
438 /******************************************************************************/
DouglasPeuckerPuntosEnTolPlanoSerie(const double * x,const double * y,const size_t incX,const size_t incY,const double tol,const size_t posBaseIni,const size_t posBaseFin,const size_t posPtoIni,const size_t posPtoFin)439 int DouglasPeuckerPuntosEnTolPlanoSerie(const double* x,
440 const double* y,
441 const size_t incX,
442 const size_t incY,
443 const double tol,
444 const size_t posBaseIni,
445 const size_t posBaseFin,
446 const size_t posPtoIni,
447 const size_t posPtoFin)
448 {
449 //índice para recorrer bucles
450 size_t i=0;
451 //valor absoluto de la tolerancia
452 double atol=fabs(tol);
453 //coordenadas de los vértices
454 double xIni=0.0,yIni=0.0,xFin=0.0,yFin=0.0,xTrab=0.0,yTrab=0.0;
455 //longitud de la base y parámetros de rotación para el plano
456 double dx=0.0,sA=0.0,cA=0.0;
457 //distancia calculada
458 double dist=0.0;
459 //variable de salida
460 int entol=1;
461 ////////////////////////////////////////////////////////////////////////////
462 ////////////////////////////////////////////////////////////////////////////
463 //comprobamos salida rápida
464 if((posBaseIni+1)>=posBaseFin)
465 {
466 //salimos de la función
467 return entol;
468 }
469 ////////////////////////////////////////////////////////////////////////////
470 ////////////////////////////////////////////////////////////////////////////
471 //coordenadas del primer punto de la base
472 xIni = x[posBaseIni*incX];
473 yIni = y[posBaseIni*incY];
474 //coordenadas del segundo punto de la base, referidas al primero
475 xFin = x[posBaseFin*incX]-xIni;
476 yFin = y[posBaseFin*incY]-yIni;
477 //calculamos la longitud de la base y la rotación
478 DouglasPeuckerParamRotaBase(xFin,yFin,&sA,&cA,&dx);
479 ////////////////////////////////////////////////////////////////////////////
480 ////////////////////////////////////////////////////////////////////////////
481 //recorremos los puntos a chequear
482 for(i=posPtoIni;i<=posPtoFin;i++)
483 {
484 //extraemos las coordenadas del vértice de trabajo y las referimos al
485 //punto inicial de la base
486 xTrab = x[i*incX]-xIni;
487 yTrab = y[i*incY]-yIni;
488 //calculamos la distancia del punto a la base
489 dist = DouglasPeuckerDistMaxPlanoAux(dx,sA,cA,xTrab,yTrab);
490 //comprobamos si está fuera de tolerancia
491 if(dist>atol)
492 {
493 //indicamos que estamos fuera de tolerancia
494 entol = 0;
495 //salimos del bucle
496 break;
497 }
498 }
499 ////////////////////////////////////////////////////////////////////////////
500 ////////////////////////////////////////////////////////////////////////////
501 //salimos de la función
502 return entol;
503 }
504 /******************************************************************************/
505 /******************************************************************************/
DouglasPeuckerRobIntersecOrigPlano(const double * x,const double * y,const size_t nPtos,const size_t incX,const size_t incY,const size_t segAUsar,const size_t posIni,size_t * posFin)506 void DouglasPeuckerRobIntersecOrigPlano(const double* x,
507 const double* y,
508 const size_t nPtos,
509 const size_t incX,
510 const size_t incY,
511 const size_t segAUsar,
512 const size_t posIni,
513 size_t* posFin)
514 {
515 //índice para recorrer bucles
516 size_t i=0;
517 //coordenadas del segmento base
518 double xA=0.0,yA=0.0,xB=0.0,yB=0.0;
519 //posición de parada para comprobar la intersección de segmentos/arcos
520 size_t posParada=0;
521 //identificación de paralelización
522 int paraleliza=0;
523 //variable identificadora de existencia de corte de segmentos
524 int corte=0;
525 ////////////////////////////////////////////////////////////////////////////
526 ////////////////////////////////////////////////////////////////////////////
527 //si los puntos de inicio y fin son contiguos, hay salida rápida
528 if((posIni+1)>=(*posFin))
529 {
530 //salimos de la función
531 return;
532 }
533 ////////////////////////////////////////////////////////////////////////////
534 ////////////////////////////////////////////////////////////////////////////
535 #if defined(_OPENMP)
536 //comprobamos si hay más de un procesador
537 if(omp_get_num_procs()>1)
538 {
539 //indicamos que hay paralelización
540 paraleliza = 1;
541 }
542 #endif
543 ////////////////////////////////////////////////////////////////////////////
544 ////////////////////////////////////////////////////////////////////////////
545 //posición de parada para el chequeo de segmentos/arcos
546 posParada = ((segAUsar==0)||(segAUsar>=(nPtos-(*posFin))))
547 ? nPtos-1
548 : (*posFin)+segAUsar;
549 ////////////////////////////////////////////////////////////////////////////
550 ////////////////////////////////////////////////////////////////////////////
551 //coordenadas del punto inicial del segmento/arco base (no cambian)
552 xA = x[posIni*incX];
553 yA = y[posIni*incY];
554 //construimos todos los segmentos/arcos base posibles
555 for(i=(*posFin);i>posIni;i--)
556 {
557 //comprobamos si estamos ante el punto posterior al inicial
558 if(i==(posIni+1))
559 {
560 //este punto es el siguiente al punto inicial del segmento/arco base
561 *posFin = i;
562 //salimos del bucle
563 break;
564 }
565 //coordenadas del punto final del segmento base
566 xB = x[i*incX];
567 yB = y[i*incY];
568 //comprobamos si hay que paralelizar
569 if(paraleliza)
570 {
571 //calculamos en paralelo
572 corte = DouglasPeuckerRobIntersecOrigPlanoOMP(xA,yA,xB,yB,x,y,incX,
573 incY,i,posParada);
574 }
575 else
576 {
577 //calculamos en serie
578 corte = DouglasPeuckerRobIntersecOrigPlanoSerie(xA,yA,xB,yB,x,y,
579 incX,incY,i,
580 posParada);
581 }
582 //comprobamos si no ha habido ninguna intersección
583 if(!corte)
584 {
585 //indicamos el índice del vértice final
586 *posFin = i;
587 //salimos del bucle de segmentos base
588 break;
589 }
590 }
591 ////////////////////////////////////////////////////////////////////////////
592 ////////////////////////////////////////////////////////////////////////////
593 //salimos de la función
594 return;
595 }
596 /******************************************************************************/
597 /******************************************************************************/
DouglasPeuckerRobIntersecOrigPlanoOMP(const double xA,const double yA,const double xB,const double yB,const double * x,const double * y,const size_t incX,const size_t incY,const size_t posIni,const size_t posFin)598 int DouglasPeuckerRobIntersecOrigPlanoOMP(const double xA,
599 const double yA,
600 const double xB,
601 const double yB,
602 const double* x,
603 const double* y,
604 const size_t incX,
605 const size_t incY,
606 const size_t posIni,
607 const size_t posFin)
608 {
609 //índice para recorrer bucles
610 size_t i=0;
611 //coordenadas de los segmentos de trabajo
612 double xC=0.0,yC=0.0,xD=0.0,yD=0.0;
613 //variable de salida
614 int corte=0;
615 ////////////////////////////////////////////////////////////////////////////
616 ////////////////////////////////////////////////////////////////////////////
617 //paralelización con OpenMP
618 #if defined(_OPENMP)
619 #pragma omp parallel for default(none) \
620 shared(posIni,posFin,incX,x,incY,y,xA,xB,yA,yB) \
621 private(i,xC,yC,xD,yD) \
622 reduction(+:corte)
623 #endif
624 //recorremos los puntos de trabajo
625 for(i=posIni;i<posFin;i++)
626 {
627 //sólo realizo cálculos si no se ha encontrado ningún corte en este hilo
628 if(!corte)
629 {
630 //puntos inicial y final del siguiente segmento/arco de trabajo
631 xC = x[i*incX];
632 yC = y[i*incY];
633 xD = x[(i+1)*incX];
634 yD = y[(i+1)*incY];
635 //sigo si los rectángulos que encierran a los segmentos se cortan
636 if(!GEOC_RECT_DISJUNTOS(GEOC_MIN(xA,xB),GEOC_MAX(xA,xB),
637 GEOC_MIN(yA,yB),GEOC_MAX(yA,yB),
638 GEOC_MIN(xC,xD),GEOC_MAX(xC,xD),
639 GEOC_MIN(yC,yD),GEOC_MAX(yC,yD)))
640 {
641 //comprobamos si hay intersección
642 corte += DouglasPeuckerRobIntersecPlano(xA,yA,xB,yB,xC,yC,xD,yD,
643 posIni,i);
644 }
645 }
646 } // --> fin del #pragma omp parallel for
647 ////////////////////////////////////////////////////////////////////////////
648 ////////////////////////////////////////////////////////////////////////////
649 //salimos de la función
650 return corte;
651 }
652 /******************************************************************************/
653 /******************************************************************************/
DouglasPeuckerRobIntersecOrigPlanoSerie(const double xA,const double yA,const double xB,const double yB,const double * x,const double * y,const size_t incX,const size_t incY,const size_t posIni,const size_t posFin)654 int DouglasPeuckerRobIntersecOrigPlanoSerie(const double xA,
655 const double yA,
656 const double xB,
657 const double yB,
658 const double* x,
659 const double* y,
660 const size_t incX,
661 const size_t incY,
662 const size_t posIni,
663 const size_t posFin)
664 {
665 //índice para recorrer bucles
666 size_t i=0;
667 //coordenadas de los segmentos de trabajo
668 double xC=0.0,yC=0.0,xD=0.0,yD=0.0;
669 //variable de salida
670 int corte=0;
671 ////////////////////////////////////////////////////////////////////////////
672 ////////////////////////////////////////////////////////////////////////////
673 //recorremos los puntos de trabajo
674 for(i=posIni;i<posFin;i++)
675 {
676 //puntos inicial y final del siguiente segmento/arco de trabajo
677 xC = x[i*incX];
678 yC = y[i*incY];
679 xD = x[(i+1)*incX];
680 yD = y[(i+1)*incY];
681 //seguimos si los rectángulos que encierran a los segmentos se cortan
682 if(!GEOC_RECT_DISJUNTOS(GEOC_MIN(xA,xB),GEOC_MAX(xA,xB),
683 GEOC_MIN(yA,yB),GEOC_MAX(yA,yB),
684 GEOC_MIN(xC,xD),GEOC_MAX(xC,xD),
685 GEOC_MIN(yC,yD),GEOC_MAX(yC,yD)))
686 {
687 //comprobamos si hay intersección
688 corte = DouglasPeuckerRobIntersecPlano(xA,yA,xB,yB,xC,yC,xD,yD,
689 posIni,i);
690 }
691 //si ha habido intersección de segmentos/arcos, salimos del bucle
692 if(corte)
693 {
694 //salimos del bucle
695 break;
696 }
697 }
698 ////////////////////////////////////////////////////////////////////////////
699 ////////////////////////////////////////////////////////////////////////////
700 //salimos de la función
701 return corte;
702 }
703 /******************************************************************************/
704 /******************************************************************************/
DouglasPeuckerRobAutoIntersecPlano(const double * x,const double * y,const size_t incX,const size_t incY,const size_t posIni,size_t * posFin,const size_t * posAlig,const size_t nPosAlig,const size_t segAUsar)705 void DouglasPeuckerRobAutoIntersecPlano(const double* x,
706 const double* y,
707 const size_t incX,
708 const size_t incY,
709 const size_t posIni,
710 size_t* posFin,
711 const size_t* posAlig,
712 const size_t nPosAlig,
713 const size_t segAUsar)
714 {
715 //índice para recorrer bucles
716 size_t i=0;
717 //coordenadas del segmento base
718 double xA=0.0,yA=0.0,xB=0.0,yB=0.0;
719 //posición de parada para comprobar la intersección de segmentos/arcos
720 size_t posParada=0;
721 //identificación de paralelización
722 int paraleliza=0;
723 //variable identificadora de existencia de corte de segmentos
724 int corte=0;
725 ////////////////////////////////////////////////////////////////////////////
726 ////////////////////////////////////////////////////////////////////////////
727 //si los puntos de inicio y fin son contiguos, hay salida rápida
728 if((posIni+1)>=(*posFin))
729 {
730 //salimos de la función
731 return;
732 }
733 ////////////////////////////////////////////////////////////////////////////
734 ////////////////////////////////////////////////////////////////////////////
735 #if defined(_OPENMP)
736 //comprobamos si hay más de un procesador
737 if(omp_get_num_procs()>1)
738 {
739 //indicamos que hay paralelización
740 paraleliza = 1;
741 }
742 #endif
743 ////////////////////////////////////////////////////////////////////////////
744 ////////////////////////////////////////////////////////////////////////////
745 //posición de parada en el vector posAlig para el chequeo de segmentos/arcos
746 posParada = ((segAUsar==0)||(segAUsar>=(nPosAlig-1)))
747 ? 0
748 : nPosAlig-1-segAUsar;
749 //coordenadas del punto inicial del segmento base (no cambian)
750 xA = x[posIni*incX];
751 yA = y[posIni*incY];
752 //construimos todos los segmentos base posibles
753 for(i=(*posFin);i>posIni;i--)
754 {
755 //comprobamos si estamos ante el punto posterior al inicial
756 if(i==(posIni+1))
757 {
758 //este punto es el siguiente al punto inicial del segmento base
759 *posFin = i;
760 //salimos del bucle
761 break;
762 }
763 //coordenadas del punto final del segmento base
764 xB = x[i*incX];
765 yB = y[i*incY];
766 //comprobamos si hay que paralelizar
767 if(paraleliza)
768 {
769 //calculamos en paralelo
770 corte = DouglasPeuckerRobAutoIntersecPlanoOMP(xA,yA,xB,yB,x,y,incX,
771 incY,posAlig,nPosAlig,
772 nPosAlig-1,posParada);
773 }
774 else
775 {
776 //calculamos en serie
777 corte = DouglasPeuckerRobAutoIntersecPlanoSerie(xA,yA,xB,yB,x,y,
778 incX,incY,posAlig,
779 nPosAlig,nPosAlig-1,
780 posParada);
781 }
782 //comprobamos si no ha habido ninguna intersección
783 if(!corte)
784 {
785 //indicamos el índice del vértice final
786 *posFin = i;
787 //salimos del bucle de segmentos base
788 break;
789 }
790 }
791 ////////////////////////////////////////////////////////////////////////////
792 ////////////////////////////////////////////////////////////////////////////
793 //salimos de la función
794 return;
795 }
796 /******************************************************************************/
797 /******************************************************************************/
DouglasPeuckerRobAutoIntersecPlanoOMP(const double xA,const double yA,const double xB,const double yB,const double * x,const double * y,const size_t incX,const size_t incY,const size_t * posAlig,const size_t nPosAlig,const size_t posIni,const size_t posFin)798 int DouglasPeuckerRobAutoIntersecPlanoOMP(const double xA,
799 const double yA,
800 const double xB,
801 const double yB,
802 const double* x,
803 const double* y,
804 const size_t incX,
805 const size_t incY,
806 const size_t* posAlig,
807 const size_t nPosAlig,
808 const size_t posIni,
809 const size_t posFin)
810 {
811 //índice para recorrer bucles
812 size_t i=0;
813 //coordenadas de los segmentos de trabajo
814 double xC=0.0,yC=0.0,xD=0.0,yD=0.0;
815 //variable de salida
816 int corte=0;
817 ////////////////////////////////////////////////////////////////////////////
818 ////////////////////////////////////////////////////////////////////////////
819 //paralelización con OpenMP
820 #if defined(_OPENMP)
821 #pragma omp parallel for default(none) \
822 shared(posIni,posFin,incX,x,posAlig,incY,y,xA,xB,yA,yB,nPosAlig) \
823 private(i,xC,yC,xD,yD) \
824 reduction(+:corte)
825 #endif
826 //recorremos los puntos de trabajo
827 for(i=posIni;i>posFin;i--)
828 {
829 //sólo realizo cálculos si no se ha encontrado ningún corte en este hilo
830 if(!corte)
831 {
832 //puntos inicial y final del siguiente segmento/arco de trabajo
833 xC = x[posAlig[i]*incX];
834 yC = y[posAlig[i]*incY];
835 xD = x[posAlig[i-1]*incX];
836 yD = y[posAlig[i-1]*incY];
837 //seguimos si los rectángulos que encierran a los segmentos se cortan
838 if(!GEOC_RECT_DISJUNTOS(GEOC_MIN(xA,xB),GEOC_MAX(xA,xB),
839 GEOC_MIN(yA,yB),GEOC_MAX(yA,yB),
840 GEOC_MIN(xC,xD),GEOC_MAX(xC,xD),
841 GEOC_MIN(yC,yD),GEOC_MAX(yC,yD)))
842 {
843 //comprobamos si hay intersección
844 corte += DouglasPeuckerRobIntersecPlano(xB,yB,xA,yA,xC,yC,xD,yD,
845 posAlig[nPosAlig-1],
846 posAlig[i]);
847 }
848 }
849 } // --> fin del #pragma omp parallel for
850 ////////////////////////////////////////////////////////////////////////////
851 ////////////////////////////////////////////////////////////////////////////
852 //salimos de la función
853 return corte;
854 }
855 /******************************************************************************/
856 /******************************************************************************/
DouglasPeuckerRobAutoIntersecPlanoSerie(const double xA,const double yA,const double xB,const double yB,const double * x,const double * y,const size_t incX,const size_t incY,const size_t * posAlig,const size_t nPosAlig,const size_t posIni,const size_t posFin)857 int DouglasPeuckerRobAutoIntersecPlanoSerie(const double xA,
858 const double yA,
859 const double xB,
860 const double yB,
861 const double* x,
862 const double* y,
863 const size_t incX,
864 const size_t incY,
865 const size_t* posAlig,
866 const size_t nPosAlig,
867 const size_t posIni,
868 const size_t posFin)
869 {
870 //índice para recorrer bucles
871 size_t i=0;
872 //coordenadas de los segmentos de trabajo
873 double xC=0.0,yC=0.0,xD=0.0,yD=0.0;
874 //variable de salida
875 int corte=0;
876 ////////////////////////////////////////////////////////////////////////////
877 ////////////////////////////////////////////////////////////////////////////
878 //recorremos los puntos de trabajo
879 for(i=posIni;i>posFin;i--)
880 {
881 //puntos inicial y final del siguiente segmento/arco de trabajo
882 xC = x[posAlig[i]*incX];
883 yC = y[posAlig[i]*incY];
884 xD = x[posAlig[i-1]*incX];
885 yD = y[posAlig[i-1]*incY];
886 //seguimos si los rectángulos que encierran a los segmentos se cortan
887 if(!GEOC_RECT_DISJUNTOS(GEOC_MIN(xA,xB),GEOC_MAX(xA,xB),
888 GEOC_MIN(yA,yB),GEOC_MAX(yA,yB),
889 GEOC_MIN(xC,xD),GEOC_MAX(xC,xD),
890 GEOC_MIN(yC,yD),GEOC_MAX(yC,yD)))
891 {
892 //comprobamos si hay intersección
893 corte = DouglasPeuckerRobIntersecPlano(xB,yB,xA,yA,xC,yC,xD,yD,
894 posAlig[nPosAlig-1],
895 posAlig[i]);
896 }
897 //si ha habido intersección de segmentos/arcos, salimos del bucle
898 if(corte)
899 {
900 //salimos del bucle
901 break;
902 }
903 }
904 ////////////////////////////////////////////////////////////////////////////
905 ////////////////////////////////////////////////////////////////////////////
906 //salimos de la función
907 return corte;
908 }
909 /******************************************************************************/
910 /******************************************************************************/
DouglasPeuckerRobIntersecPlano(const double xA,const double yA,const double xB,const double yB,const double xC,const double yC,const double xD,const double yD,const size_t posFinAB,const size_t posIniCD)911 int DouglasPeuckerRobIntersecPlano(const double xA,
912 const double yA,
913 const double xB,
914 const double yB,
915 const double xC,
916 const double yC,
917 const double xD,
918 const double yD,
919 const size_t posFinAB,
920 const size_t posIniCD)
921 {
922 //variables auxiliares
923 double xAux=0.0,yAux=0.0;
924 //identificador de intersección
925 int inter=GEOC_DPEUCKER_NO_INTERSEC;
926 //variable de salida
927 int corte=0;
928 ////////////////////////////////////////////////////////////////////////////
929 ////////////////////////////////////////////////////////////////////////////
930 //seguimos si los rectángulos que encierran a los segmentos se cortan
931 if(!GEOC_RECT_DISJUNTOS(GEOC_MIN(xA,xB),GEOC_MAX(xA,xB),
932 GEOC_MIN(yA,yB),GEOC_MAX(yA,yB),
933 GEOC_MIN(xC,xD),GEOC_MAX(xC,xD),
934 GEOC_MIN(yC,yD),GEOC_MAX(yC,yD)))
935 {
936 //compruebo si los dos segmentos son contiguos
937 if(posFinAB==posIniCD)
938 {
939 //compruebo intersección con la función completa (lenta)
940 inter = IntersecSegmentos2D(xA,yA,xB,yB,xC,yC,xD,yD,
941 &xAux,&yAux);
942 //compruebo si es la sucesión de segmento inicial+final
943 if((inter!=GEOC_SEG_INTERSEC_MISMO_SEG)&&
944 (inter!=GEOC_SEG_INTERSEC_COLIN))
945 {
946 //en este caso, no hay intersección
947 inter = GEOC_DPEUCKER_NO_INTERSEC;
948 }
949 }
950 else
951 {
952 //compruebo intersección con la función simple (rápida)
953 inter = IntersecSegmentos2DSimple(xA,yA,xB,yB,xC,yC,xD,yD);
954 }
955 //unificamos los identificadores de intersección
956 if(!((inter==GEOC_SEG_NO_INTERSEC)||(inter==GEOC_DPEUCKER_NO_INTERSEC)))
957 {
958 //hay intersección
959 corte = 1;
960 }
961 }
962 ////////////////////////////////////////////////////////////////////////////
963 ////////////////////////////////////////////////////////////////////////////
964 //salimos de la función
965 return corte;
966 }
967 /******************************************************************************/
968 /******************************************************************************/
DouglasPeuckerDistMaxPlano(const double * x,const double * y,const size_t incX,const size_t incY,const size_t posIni,const size_t posFin,size_t * pos)969 double DouglasPeuckerDistMaxPlano(const double* x,
970 const double* y,
971 const size_t incX,
972 const size_t incY,
973 const size_t posIni,
974 const size_t posFin,
975 size_t* pos)
976 {
977 //índice para recorrer bucles
978 size_t i=0;
979 //coordenadas de los extremos del segmento base y del punto de trabajo
980 double xIni=0.0,yIni=0.0,xFin=0.0,yFin=0.0,xP=0.0,yP=0.0;
981 //razones trigonométricas del ángulo de giro del sistema de coordenadas
982 double sA=0.0,cA=0.0;
983 //longitudes auxiliares
984 double dx=0.0,lonAux=0.0;
985 //variable de salida
986 double lon=-1.0;
987 ////////////////////////////////////////////////////////////////////////////
988 ////////////////////////////////////////////////////////////////////////////
989 //comprobamos si los dos puntos están seguidos
990 if((posIni+1)==posFin)
991 {
992 //la posición que devolvemos es la del punto inicial
993 *pos = posIni;
994 //la distancia devuelta es -1.0
995 return lon;
996 }
997 ////////////////////////////////////////////////////////////////////////////
998 ////////////////////////////////////////////////////////////////////////////
999 //coordenadas del primer punto de la base
1000 xIni = x[posIni*incX];
1001 yIni = y[posIni*incY];
1002 //coordenadas del segundo punto de la base, referidas al primero
1003 xFin = x[posFin*incX]-xIni;
1004 yFin = y[posFin*incY]-yIni;
1005 //calculamos la longitud de la base y la rotación
1006 DouglasPeuckerParamRotaBase(xFin,yFin,&sA,&cA,&dx);
1007 //recorremos los puntos entre los extremos
1008 for(i=posIni+1;i<posFin;i++)
1009 {
1010 //coordenadas del punto de trabajo referidas al punto inicial de la base
1011 xP = x[i*incX]-xIni;
1012 yP = y[i*incY]-yIni;
1013 //calculamos la distancia del punto de trabajo al segmento base
1014 lonAux = DouglasPeuckerDistMaxPlanoAux(dx,sA,cA,xP,yP);
1015 //comprobamos si la distancia calculada es mayor que la anterior
1016 if(lonAux>lon)
1017 {
1018 //actualizamos la distancia máxima
1019 lon = lonAux;
1020 //guardamos la posición del punto
1021 *pos = i;
1022 }
1023 }
1024 ////////////////////////////////////////////////////////////////////////////
1025 ////////////////////////////////////////////////////////////////////////////
1026 //salimos de la función
1027 return lon;
1028 }
1029 /******************************************************************************/
1030 /******************************************************************************/
DouglasPeuckerParamRotaBase(const double xBase2RB1,const double yBase2RB1,double * sAlfa,double * cAlfa,double * lonBase)1031 void DouglasPeuckerParamRotaBase(const double xBase2RB1,
1032 const double yBase2RB1,
1033 double* sAlfa,
1034 double* cAlfa,
1035 double* lonBase)
1036 {
1037 //álgulo de rotación
1038 double alfa=0.0;
1039 ////////////////////////////////////////////////////////////////////////////
1040 ////////////////////////////////////////////////////////////////////////////
1041 //ángulo a rotar el sistema de coordenadas para llevar el eje X a coincidir
1042 //con el segmento base
1043 alfa = atan2(yBase2RB1,xBase2RB1);
1044 //calculamos las razones trigonométricas del ángulo de rotación
1045 *sAlfa = sin(alfa);
1046 *cAlfa = cos(alfa);
1047 //la longitud del segmento base será el valor absoluto de la coordenada X
1048 //del extremo final del segmento base en el sistema de coordenadas rotado
1049 *lonBase = fabs((*cAlfa)*xBase2RB1+(*sAlfa)*yBase2RB1);
1050 ////////////////////////////////////////////////////////////////////////////
1051 ////////////////////////////////////////////////////////////////////////////
1052 //salimos de la función
1053 return;
1054 }
1055 /******************************************************************************/
1056 /******************************************************************************/
DouglasPeuckerDistMaxPlanoAux(const double lonBase,const double sAlfa,const double cAlfa,const double xVertRB1,const double yVertRB1)1057 double DouglasPeuckerDistMaxPlanoAux(const double lonBase,
1058 const double sAlfa,
1059 const double cAlfa,
1060 const double xVertRB1,
1061 const double yVertRB1)
1062 {
1063 //coordenadas del vértice de trabajo en el sistema rotado
1064 double xVert1=0.0,yVert1=0.0;
1065 //variable de salida
1066 double lon=0.0;
1067 ////////////////////////////////////////////////////////////////////////////
1068 ////////////////////////////////////////////////////////////////////////////
1069 //transformamos el vértice de trabajo al nuevo sistema de coordenadas
1070 xVert1 = cAlfa*xVertRB1+sAlfa*yVertRB1;
1071 yVert1 = -sAlfa*xVertRB1+cAlfa*yVertRB1;
1072 //comprobamos la posición del punto con respecto al segmento base
1073 if((xVert1>=0.0)&&(xVert1<=lonBase))
1074 {
1075 //el punto de trabajo está entre los extremos del segmento base
1076 //su distancia hasta él es el valor absoluto de su coordenada Y en
1077 //el sistema rotado
1078 lon = fabs(yVert1);
1079 }
1080 else if(xVert1<0.0)
1081 {
1082 //el punto de trabajo está a la izquierda del punto inicial de la
1083 //base, luego su distancia hasta el segmento será la distancia
1084 //hasta dicho punto inicial
1085 //Konrad Ebisch (2002), A correction to the Douglas-Peucker line
1086 //generalization algorithm, Computers and Geosciences, vol. 28,
1087 //págs. 995 a 997
1088 lon = sqrt(xVert1*xVert1+yVert1*yVert1);
1089 }
1090 else
1091 {
1092 //el punto de trabajo está a la derecha del punto final de la base,
1093 //luego su distancia hasta el segmento será la distancia hasta dicho
1094 //punto final
1095 //Konrad Ebisch (2002), A correction to the Douglas-Peucker line
1096 //generalization algorithm, Computers and Geosciences, vol. 28,
1097 //págs. 995 a 997
1098 lon = sqrt((xVert1-lonBase)*(xVert1-lonBase)+yVert1*yVert1);
1099 }
1100 ////////////////////////////////////////////////////////////////////////////
1101 ////////////////////////////////////////////////////////////////////////////
1102 //salimos de la función
1103 return lon;
1104 }
1105 /******************************************************************************/
1106 /******************************************************************************/
1107 /** @} */
1108 /******************************************************************************/
1109 /******************************************************************************/
1110 /* kate: encoding utf-8; end-of-line unix; syntax c; indent-mode cstyle; */
1111 /* kate: replace-tabs on; space-indent on; tab-indents off; indent-width 4; */
1112 /* kate: line-numbers on; folding-markers on; remove-trailing-space on; */
1113 /* kate: backspace-indents on; show-tabs on; */
1114 /* kate: word-wrap-column 80; word-wrap-marker-color #D2D2D2; word-wrap off; */
1115