← Back to branch summary

~ubuntu-branches/ubuntu/wily/spatialite/wily-proposed

« back to all changes in this revision

Viewing changes to src/gaiageo/gg_relations.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-07-14 11:57:46 UTC
  • mfrom: (16.1.1 sid)
  • Revision ID: package-import@ubuntu.com-20150714115746-e2iljfmb5sq7o5hh
Tags: 4.3.0-1
Move from experimental to unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gaiageo/gg_relations.c
2
2
 
3
3
 gg_relations.c -- Gaia spatial relations
4
4
    
5
 
 version 4.1, 2013 May 8
 
5
 version 4.3, 2015 June 29
6
6
 
7
7
 Author: Sandro Furieri a.furieri@lqt.it
8
8
 
24
24
 
25
25
The Initial Developer of the Original Code is Alessandro Furieri
26
26
 
27
 
Portions created by the Initial Developer are Copyright (C) 2008-2013
 
27
Portions created by the Initial Developer are Copyright (C) 2008-2015
28
28
the Initial Developer. All Rights Reserved.
29
29
 
30
30
Contributor(s):
82
82
    p->preparedGeosGeom = NULL;
83
83
}
84
84
 
 
85
SPATIALITE_PRIVATE void
 
86
splite_free_geos_cache_item_r (const void *p_cache,
 
87
                               struct splite_geos_cache_item *p)
 
88
{
 
89
#ifndef OMIT_GEOS               /* including GEOS */
 
90
    struct splite_internal_cache *cache =
 
91
        (struct splite_internal_cache *) p_cache;
 
92
    GEOSContextHandle_t handle = NULL;
 
93
    if (cache == NULL)
 
94
      {
 
95
          splite_free_geos_cache_item (p);
 
96
          return;
 
97
      }
 
98
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
99
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
100
      {
 
101
          splite_free_geos_cache_item (p);
 
102
          return;
 
103
      }
 
104
    handle = cache->GEOS_handle;
 
105
    if (handle == NULL)
 
106
      {
 
107
          splite_free_geos_cache_item (p);
 
108
          return;
 
109
      }
 
110
    if (p->preparedGeosGeom)
 
111
        GEOSPreparedGeom_destroy_r (handle, p->preparedGeosGeom);
 
112
    if (p->geosGeom)
 
113
        GEOSGeom_destroy_r (handle, p->geosGeom);
 
114
#endif
 
115
    p->geosGeom = NULL;
 
116
    p->preparedGeosGeom = NULL;
 
117
}
 
118
 
85
119
GAIAGEO_DECLARE void
86
120
gaiaResetGeosMsg ()
87
121
{
121
155
GAIAGEO_DECLARE void
122
156
gaiaSetGeosErrorMsg (const char *msg)
123
157
{
124
 
/* return the latest GEOS error message */
 
158
/* setting the latest GEOS error message */
125
159
    int len;
126
160
    if (gaia_geos_error_msg != NULL)
127
161
        free (gaia_geos_error_msg);
136
170
GAIAGEO_DECLARE void
137
171
gaiaSetGeosWarningMsg (const char *msg)
138
172
{
139
 
/* return the latest GEOS error message */
 
173
/* setting the latest GEOS error message */
140
174
    int len;
141
175
    if (gaia_geos_warning_msg != NULL)
142
176
        free (gaia_geos_warning_msg);
151
185
GAIAGEO_DECLARE void
152
186
gaiaSetGeosAuxErrorMsg (const char *msg)
153
187
{
154
 
/* return the latest GEOS (auxiliary) error message */
 
188
/* setting the latest GEOS (auxiliary) error message */
155
189
    int len;
156
190
    if (gaia_geosaux_error_msg != NULL)
157
191
        free (gaia_geosaux_error_msg);
345
379
               gaiaGeomCollPtr * geom)
346
380
{
347
381
/* handling the internal GEOS cache */
348
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
349
382
    struct splite_geos_cache_item *p1 = &(cache->cacheItem1);
350
383
    struct splite_geos_cache_item *p2 = &(cache->cacheItem2);
351
384
    uLong crc1 = crc32 (0L, blob1, size1);
352
385
    uLong crc2 = crc32 (0L, blob2, size2);
 
386
    GEOSContextHandle_t handle = NULL;
 
387
    if (cache == NULL)
 
388
        return 0;
 
389
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
390
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
391
        return 0;
 
392
    handle = cache->GEOS_handle;
 
393
    if (handle == NULL)
 
394
        return 0;
353
395
 
354
396
/* checking the first cache item */
355
397
    if (evalGeosCacheItem (blob1, size1, crc1, p1))
358
400
          if (p1->preparedGeosGeom == NULL)
359
401
            {
360
402
                /* preparing the GeosGeometries */
361
 
                p1->geosGeom = gaiaToGeos (geom1);
 
403
                p1->geosGeom = gaiaToGeos_r (cache, geom1);
362
404
                if (p1->geosGeom)
363
405
                  {
364
406
                      p1->preparedGeosGeom =
365
 
                          (void *) GEOSPrepare (p1->geosGeom);
 
407
                          (void *) GEOSPrepare_r (handle, p1->geosGeom);
366
408
                      if (p1->preparedGeosGeom == NULL)
367
409
                        {
368
410
                            /* unexpected failure */
369
 
                            GEOSGeom_destroy (p1->geosGeom);
 
411
                            GEOSGeom_destroy_r (handle, p1->geosGeom);
370
412
                            p1->geosGeom = NULL;
371
413
                        }
372
414
                  }
388
430
          if (p2->preparedGeosGeom == NULL)
389
431
            {
390
432
                /* preparing the GeosGeometries */
391
 
                p2->geosGeom = gaiaToGeos (geom2);
 
433
                p2->geosGeom = gaiaToGeos_r (cache, geom2);
392
434
                if (p2->geosGeom)
393
435
                  {
394
436
                      p2->preparedGeosGeom =
395
 
                          (void *) GEOSPrepare (p2->geosGeom);
 
437
                          (void *) GEOSPrepare_r (handle, p2->geosGeom);
396
438
                      if (p2->preparedGeosGeom == NULL)
397
439
                        {
398
440
                            /* unexpected failure */
399
 
                            GEOSGeom_destroy (p2->geosGeom);
 
441
                            GEOSGeom_destroy_r (handle, p2->geosGeom);
400
442
                            p2->geosGeom = NULL;
401
443
                        }
402
444
                  }
416
458
    p1->gaiaBlobSize = size1;
417
459
    p1->crc32 = crc1;
418
460
    if (p1->preparedGeosGeom)
419
 
        GEOSPreparedGeom_destroy (p1->preparedGeosGeom);
 
461
        GEOSPreparedGeom_destroy_r (handle, p1->preparedGeosGeom);
420
462
    if (p1->geosGeom)
421
 
        GEOSGeom_destroy (p1->geosGeom);
 
463
        GEOSGeom_destroy_r (handle, p1->geosGeom);
422
464
    p1->geosGeom = NULL;
423
465
    p1->preparedGeosGeom = NULL;
424
466
 
427
469
    p2->gaiaBlobSize = size2;
428
470
    p2->crc32 = crc2;
429
471
    if (p2->preparedGeosGeom)
430
 
        GEOSPreparedGeom_destroy (p2->preparedGeosGeom);
 
472
        GEOSPreparedGeom_destroy_r (handle, p2->preparedGeosGeom);
431
473
    if (p2->geosGeom)
432
 
        GEOSGeom_destroy (p2->geosGeom);
 
474
        GEOSGeom_destroy_r (handle, p2->geosGeom);
433
475
    p2->geosGeom = NULL;
434
476
    p2->preparedGeosGeom = NULL;
435
 
#endif /* end GEOS_ADVANCED */
436
477
 
437
478
    return 0;
438
479
}
444
485
    int ret;
445
486
    GEOSGeometry *g1;
446
487
    GEOSGeometry *g2;
 
488
    gaiaResetGeosMsg ();
447
489
    if (!geom1 || !geom2)
448
490
        return -1;
449
491
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
462
504
}
463
505
 
464
506
GAIAGEO_DECLARE int
 
507
gaiaGeomCollEquals_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
508
                      gaiaGeomCollPtr geom2)
 
509
{
 
510
/* checks if two Geometries are "spatially equal" */
 
511
    int ret;
 
512
    GEOSGeometry *g1;
 
513
    GEOSGeometry *g2;
 
514
    struct splite_internal_cache *cache =
 
515
        (struct splite_internal_cache *) p_cache;
 
516
    GEOSContextHandle_t handle = NULL;
 
517
    if (cache == NULL)
 
518
        return -1;
 
519
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
520
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
521
        return -1;
 
522
    handle = cache->GEOS_handle;
 
523
    if (handle == NULL)
 
524
        return -1;
 
525
    gaiaResetGeosMsg_r (cache);
 
526
    if (!geom1 || !geom2)
 
527
        return -1;
 
528
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
529
        return -1;
 
530
 
 
531
/* quick check based on MBRs comparison */
 
532
    if (!splite_mbr_equals (geom1, geom2))
 
533
        return 0;
 
534
 
 
535
    g1 = gaiaToGeos_r (cache, geom1);
 
536
    g2 = gaiaToGeos_r (cache, geom2);
 
537
    ret = GEOSEquals_r (handle, g1, g2);
 
538
    GEOSGeom_destroy_r (handle, g1);
 
539
    GEOSGeom_destroy_r (handle, g2);
 
540
    return ret;
 
541
}
 
542
 
 
543
GAIAGEO_DECLARE int
465
544
gaiaGeomCollIntersects (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
466
545
{
467
546
/* checks if two Geometries do "spatially intersects" */
468
547
    int ret;
469
548
    GEOSGeometry *g1;
470
549
    GEOSGeometry *g2;
 
550
    gaiaResetGeosMsg ();
471
551
    if (!geom1 || !geom2)
472
552
        return -1;
473
553
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
485
565
    return ret;
486
566
}
487
567
 
488
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
489
 
 
490
 
GAIAGEO_DECLARE int
491
 
gaiaGeomCollPreparedIntersects (void *p_cache, gaiaGeomCollPtr geom1,
 
568
GAIAGEO_DECLARE int
 
569
gaiaGeomCollIntersects_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
570
                          gaiaGeomCollPtr geom2)
 
571
{
 
572
/* checks if two Geometries do "spatially intersects" */
 
573
    int ret;
 
574
    GEOSGeometry *g1;
 
575
    GEOSGeometry *g2;
 
576
    struct splite_internal_cache *cache =
 
577
        (struct splite_internal_cache *) p_cache;
 
578
    GEOSContextHandle_t handle = NULL;
 
579
    if (cache == NULL)
 
580
        return -1;
 
581
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
582
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
583
        return -1;
 
584
    handle = cache->GEOS_handle;
 
585
    if (handle == NULL)
 
586
        return -1;
 
587
    gaiaResetGeosMsg_r (cache);
 
588
    if (!geom1 || !geom2)
 
589
        return -1;
 
590
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
591
        return -1;
 
592
 
 
593
/* quick check based on MBRs comparison */
 
594
    if (!splite_mbr_overlaps (geom1, geom2))
 
595
        return 0;
 
596
 
 
597
    g1 = gaiaToGeos_r (cache, geom1);
 
598
    g2 = gaiaToGeos_r (cache, geom2);
 
599
    ret = GEOSIntersects_r (handle, g1, g2);
 
600
    GEOSGeom_destroy_r (handle, g1);
 
601
    GEOSGeom_destroy_r (handle, g2);
 
602
    return ret;
 
603
}
 
604
 
 
605
GAIAGEO_DECLARE int
 
606
gaiaGeomCollPreparedIntersects (const void *p_cache, gaiaGeomCollPtr geom1,
492
607
                                unsigned char *blob1, int size1,
493
608
                                gaiaGeomCollPtr geom2, unsigned char *blob2,
494
609
                                int size2)
501
616
    GEOSGeometry *g1;
502
617
    GEOSGeometry *g2;
503
618
    gaiaGeomCollPtr geom;
 
619
    GEOSContextHandle_t handle = NULL;
 
620
    if (cache == NULL)
 
621
        return -1;
 
622
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
623
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
624
        return -1;
 
625
    handle = cache->GEOS_handle;
 
626
    if (handle == NULL)
 
627
        return -1;
 
628
    gaiaResetGeosMsg_r (cache);
504
629
    if (!geom1 || !geom2)
505
630
        return -1;
506
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
631
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
507
632
        return -1;
508
633
 
509
634
/* quick check based on MBRs comparison */
514
639
    if (evalGeosCache
515
640
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
516
641
      {
517
 
          g2 = gaiaToGeos (geom);
518
 
          ret = GEOSPreparedIntersects (gPrep, g2);
519
 
          GEOSGeom_destroy (g2);
 
642
          g2 = gaiaToGeos_r (cache, geom);
 
643
          ret = GEOSPreparedIntersects_r (handle, gPrep, g2);
 
644
          GEOSGeom_destroy_r (handle, g2);
520
645
          return ret;
521
646
      }
522
 
    g1 = gaiaToGeos (geom1);
523
 
    g2 = gaiaToGeos (geom2);
524
 
    ret = GEOSIntersects (g1, g2);
525
 
    GEOSGeom_destroy (g1);
526
 
    GEOSGeom_destroy (g2);
 
647
    g1 = gaiaToGeos_r (cache, geom1);
 
648
    g2 = gaiaToGeos_r (cache, geom2);
 
649
    ret = GEOSIntersects_r (handle, g1, g2);
 
650
    GEOSGeom_destroy_r (handle, g1);
 
651
    GEOSGeom_destroy_r (handle, g2);
527
652
    return ret;
528
653
}
529
654
 
530
 
#endif /* end GEOS_ADVANCED */
531
 
 
532
655
GAIAGEO_DECLARE int
533
656
gaiaGeomCollDisjoint (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
534
657
{
536
659
    int ret;
537
660
    GEOSGeometry *g1;
538
661
    GEOSGeometry *g2;
 
662
    gaiaResetGeosMsg ();
539
663
    if (!geom1 || !geom2)
540
664
        return -1;
541
665
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
553
677
    return ret;
554
678
}
555
679
 
556
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
557
 
 
558
 
GAIAGEO_DECLARE int
559
 
gaiaGeomCollPreparedDisjoint (void *p_cache, gaiaGeomCollPtr geom1,
 
680
GAIAGEO_DECLARE int
 
681
gaiaGeomCollDisjoint_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
682
                        gaiaGeomCollPtr geom2)
 
683
{
 
684
/* checks if two Geometries are "spatially disjoint" */
 
685
    int ret;
 
686
    GEOSGeometry *g1;
 
687
    GEOSGeometry *g2;
 
688
    struct splite_internal_cache *cache =
 
689
        (struct splite_internal_cache *) p_cache;
 
690
    GEOSContextHandle_t handle = NULL;
 
691
    if (cache == NULL)
 
692
        return -1;
 
693
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
694
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
695
        return -1;
 
696
    handle = cache->GEOS_handle;
 
697
    if (handle == NULL)
 
698
        return -1;
 
699
    gaiaResetGeosMsg_r (cache);
 
700
    if (!geom1 || !geom2)
 
701
        return -1;
 
702
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
703
        return -1;
 
704
 
 
705
/* quick check based on MBRs comparison */
 
706
    if (!splite_mbr_overlaps (geom1, geom2))
 
707
        return 1;
 
708
 
 
709
    g1 = gaiaToGeos_r (cache, geom1);
 
710
    g2 = gaiaToGeos_r (cache, geom2);
 
711
    ret = GEOSDisjoint_r (handle, g1, g2);
 
712
    GEOSGeom_destroy_r (handle, g1);
 
713
    GEOSGeom_destroy_r (handle, g2);
 
714
    return ret;
 
715
}
 
716
 
 
717
GAIAGEO_DECLARE int
 
718
gaiaGeomCollPreparedDisjoint (const void *p_cache, gaiaGeomCollPtr geom1,
560
719
                              unsigned char *blob1, int size1,
561
720
                              gaiaGeomCollPtr geom2, unsigned char *blob2,
562
721
                              int size2)
563
722
{
564
723
/* checks if two Geometries are "spatially disjoint" */
565
724
    int ret;
566
 
    struct splite_internal_cache *cache =
567
 
        (struct splite_internal_cache *) p_cache;
568
725
    GEOSPreparedGeometry *gPrep;
569
726
    GEOSGeometry *g1;
570
727
    GEOSGeometry *g2;
571
728
    gaiaGeomCollPtr geom;
 
729
    GEOSContextHandle_t handle = NULL;
 
730
    struct splite_internal_cache *cache =
 
731
        (struct splite_internal_cache *) p_cache;
 
732
    if (cache == NULL)
 
733
        return -1;
 
734
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
735
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
736
        return -1;
 
737
    handle = cache->GEOS_handle;
 
738
    if (handle == NULL)
 
739
        return -1;
 
740
    gaiaResetGeosMsg_r (cache);
572
741
    if (!geom1 || !geom2)
573
742
        return -1;
574
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
743
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
575
744
        return -1;
576
745
 
577
746
/* quick check based on MBRs comparison */
582
751
    if (evalGeosCache
583
752
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
584
753
      {
585
 
          g2 = gaiaToGeos (geom);
586
 
          ret = GEOSPreparedDisjoint (gPrep, g2);
587
 
          GEOSGeom_destroy (g2);
 
754
          g2 = gaiaToGeos_r (cache, geom);
 
755
          ret = GEOSPreparedDisjoint_r (handle, gPrep, g2);
 
756
          GEOSGeom_destroy_r (handle, g2);
588
757
          return ret;
589
758
      }
590
759
 
591
 
    g1 = gaiaToGeos (geom1);
592
 
    g2 = gaiaToGeos (geom2);
593
 
    ret = GEOSDisjoint (g1, g2);
594
 
    GEOSGeom_destroy (g1);
595
 
    GEOSGeom_destroy (g2);
 
760
    g1 = gaiaToGeos_r (cache, geom1);
 
761
    g2 = gaiaToGeos_r (cache, geom2);
 
762
    ret = GEOSDisjoint_r (handle, g1, g2);
 
763
    GEOSGeom_destroy_r (handle, g1);
 
764
    GEOSGeom_destroy_r (handle, g2);
596
765
    return ret;
597
766
}
598
767
 
599
 
#endif /* end GEOS_ADVANCED */
600
 
 
601
768
GAIAGEO_DECLARE int
602
769
gaiaGeomCollOverlaps (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
603
770
{
605
772
    int ret;
606
773
    GEOSGeometry *g1;
607
774
    GEOSGeometry *g2;
 
775
    gaiaResetGeosMsg ();
608
776
    if (!geom1 || !geom2)
609
777
        return -1;
610
778
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
622
790
    return ret;
623
791
}
624
792
 
625
 
 
626
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
627
 
 
628
 
GAIAGEO_DECLARE int
629
 
gaiaGeomCollPreparedOverlaps (void *p_cache, gaiaGeomCollPtr geom1,
 
793
GAIAGEO_DECLARE int
 
794
gaiaGeomCollOverlaps_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
795
                        gaiaGeomCollPtr geom2)
 
796
{
 
797
/* checks if two Geometries do "spatially overlaps" */
 
798
    int ret;
 
799
    GEOSGeometry *g1;
 
800
    GEOSGeometry *g2;
 
801
    struct splite_internal_cache *cache =
 
802
        (struct splite_internal_cache *) p_cache;
 
803
    GEOSContextHandle_t handle = NULL;
 
804
    if (cache == NULL)
 
805
        return -1;
 
806
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
807
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
808
        return -1;
 
809
    handle = cache->GEOS_handle;
 
810
    if (handle == NULL)
 
811
        return -1;
 
812
    gaiaResetGeosMsg_r (cache);
 
813
    if (!geom1 || !geom2)
 
814
        return -1;
 
815
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
816
        return -1;
 
817
 
 
818
/* quick check based on MBRs comparison */
 
819
    if (!splite_mbr_overlaps (geom1, geom2))
 
820
        return 0;
 
821
 
 
822
    g1 = gaiaToGeos_r (cache, geom1);
 
823
    g2 = gaiaToGeos_r (cache, geom2);
 
824
    ret = GEOSOverlaps_r (handle, g1, g2);
 
825
    GEOSGeom_destroy_r (handle, g1);
 
826
    GEOSGeom_destroy_r (handle, g2);
 
827
    return ret;
 
828
}
 
829
 
 
830
GAIAGEO_DECLARE int
 
831
gaiaGeomCollPreparedOverlaps (const void *p_cache, gaiaGeomCollPtr geom1,
630
832
                              unsigned char *blob1, int size1,
631
833
                              gaiaGeomCollPtr geom2, unsigned char *blob2,
632
834
                              int size2)
639
841
    GEOSGeometry *g1;
640
842
    GEOSGeometry *g2;
641
843
    gaiaGeomCollPtr geom;
 
844
    GEOSContextHandle_t handle = NULL;
 
845
    if (cache == NULL)
 
846
        return -1;
 
847
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
848
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
849
        return -1;
 
850
    handle = cache->GEOS_handle;
 
851
    if (handle == NULL)
 
852
        return -1;
 
853
    gaiaResetGeosMsg_r (cache);
642
854
    if (!geom1 || !geom2)
643
855
        return -1;
644
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
856
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
645
857
        return -1;
646
858
 
647
859
/* quick check based on MBRs comparison */
652
864
    if (evalGeosCache
653
865
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
654
866
      {
655
 
          g2 = gaiaToGeos (geom);
656
 
          ret = GEOSPreparedOverlaps (gPrep, g2);
657
 
          GEOSGeom_destroy (g2);
 
867
          g2 = gaiaToGeos_r (cache, geom);
 
868
          ret = GEOSPreparedOverlaps_r (handle, gPrep, g2);
 
869
          GEOSGeom_destroy_r (handle, g2);
658
870
          return ret;
659
871
      }
660
872
 
661
 
    g1 = gaiaToGeos (geom1);
662
 
    g2 = gaiaToGeos (geom2);
663
 
    ret = GEOSOverlaps (g1, g2);
664
 
    GEOSGeom_destroy (g1);
665
 
    GEOSGeom_destroy (g2);
 
873
    g1 = gaiaToGeos_r (cache, geom1);
 
874
    g2 = gaiaToGeos_r (cache, geom2);
 
875
    ret = GEOSOverlaps_r (handle, g1, g2);
 
876
    GEOSGeom_destroy_r (handle, g1);
 
877
    GEOSGeom_destroy_r (handle, g2);
666
878
    return ret;
667
879
}
668
880
 
669
 
#endif /* end GEOS_ADVANCED */
670
 
 
671
881
GAIAGEO_DECLARE int
672
882
gaiaGeomCollCrosses (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
673
883
{
675
885
    int ret;
676
886
    GEOSGeometry *g1;
677
887
    GEOSGeometry *g2;
 
888
    gaiaResetGeosMsg ();
678
889
    if (!geom1 || !geom2)
679
890
        return -1;
680
891
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
692
903
    return ret;
693
904
}
694
905
 
695
 
 
696
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
697
 
 
698
 
GAIAGEO_DECLARE int
699
 
gaiaGeomCollPreparedCrosses (void *p_cache, gaiaGeomCollPtr geom1,
 
906
GAIAGEO_DECLARE int
 
907
gaiaGeomCollCrosses_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
908
                       gaiaGeomCollPtr geom2)
 
909
{
 
910
/* checks if two Geometries do "spatially crosses" */
 
911
    int ret;
 
912
    GEOSGeometry *g1;
 
913
    GEOSGeometry *g2;
 
914
    struct splite_internal_cache *cache =
 
915
        (struct splite_internal_cache *) p_cache;
 
916
    GEOSContextHandle_t handle = NULL;
 
917
    if (cache == NULL)
 
918
        return -1;
 
919
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
920
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
921
        return -1;
 
922
    handle = cache->GEOS_handle;
 
923
    if (handle == NULL)
 
924
        return -1;
 
925
    gaiaResetGeosMsg_r (cache);
 
926
    if (!geom1 || !geom2)
 
927
        return -1;
 
928
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
929
        return -1;
 
930
 
 
931
/* quick check based on MBRs comparison */
 
932
    if (!splite_mbr_overlaps (geom1, geom2))
 
933
        return 0;
 
934
 
 
935
    g1 = gaiaToGeos_r (cache, geom1);
 
936
    g2 = gaiaToGeos_r (cache, geom2);
 
937
    ret = GEOSCrosses_r (handle, g1, g2);
 
938
    GEOSGeom_destroy_r (handle, g1);
 
939
    GEOSGeom_destroy_r (handle, g2);
 
940
    return ret;
 
941
}
 
942
 
 
943
GAIAGEO_DECLARE int
 
944
gaiaGeomCollPreparedCrosses (const void *p_cache, gaiaGeomCollPtr geom1,
700
945
                             unsigned char *blob1, int size1,
701
946
                             gaiaGeomCollPtr geom2, unsigned char *blob2,
702
947
                             int size2)
709
954
    GEOSGeometry *g1;
710
955
    GEOSGeometry *g2;
711
956
    gaiaGeomCollPtr geom;
 
957
    GEOSContextHandle_t handle = NULL;
 
958
    if (cache == NULL)
 
959
        return -1;
 
960
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
961
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
962
        return -1;
 
963
    handle = cache->GEOS_handle;
 
964
    if (handle == NULL)
 
965
        return -1;
 
966
    gaiaResetGeosMsg_r (cache);
712
967
    if (!geom1 || !geom2)
713
968
        return -1;
714
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
969
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
715
970
        return -1;
716
971
 
717
972
/* quick check based on MBRs comparison */
722
977
    if (evalGeosCache
723
978
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
724
979
      {
725
 
          g2 = gaiaToGeos (geom);
726
 
          ret = GEOSPreparedCrosses (gPrep, g2);
727
 
          GEOSGeom_destroy (g2);
 
980
          g2 = gaiaToGeos_r (cache, geom);
 
981
          ret = GEOSPreparedCrosses_r (handle, gPrep, g2);
 
982
          GEOSGeom_destroy_r (handle, g2);
728
983
          return ret;
729
984
      }
730
985
 
731
 
    g1 = gaiaToGeos (geom1);
732
 
    g2 = gaiaToGeos (geom2);
733
 
    ret = GEOSCrosses (g1, g2);
734
 
    GEOSGeom_destroy (g1);
735
 
    GEOSGeom_destroy (g2);
 
986
    g1 = gaiaToGeos_r (cache, geom1);
 
987
    g2 = gaiaToGeos_r (cache, geom2);
 
988
    ret = GEOSCrosses_r (handle, g1, g2);
 
989
    GEOSGeom_destroy_r (handle, g1);
 
990
    GEOSGeom_destroy_r (handle, g2);
736
991
    return ret;
737
992
}
738
993
 
739
 
#endif /* end GEOS_ADVANCED */
740
 
 
741
994
GAIAGEO_DECLARE int
742
995
gaiaGeomCollTouches (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
743
996
{
745
998
    int ret;
746
999
    GEOSGeometry *g1;
747
1000
    GEOSGeometry *g2;
 
1001
    gaiaResetGeosMsg ();
748
1002
    if (!geom1 || !geom2)
749
1003
        return -1;
750
1004
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
762
1016
    return ret;
763
1017
}
764
1018
 
765
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
766
 
 
767
 
GAIAGEO_DECLARE int
768
 
gaiaGeomCollPreparedTouches (void *p_cache, gaiaGeomCollPtr geom1,
 
1019
GAIAGEO_DECLARE int
 
1020
gaiaGeomCollTouches_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1021
                       gaiaGeomCollPtr geom2)
 
1022
{
 
1023
/* checks if two Geometries do "spatially touches" */
 
1024
    int ret;
 
1025
    GEOSGeometry *g1;
 
1026
    GEOSGeometry *g2;
 
1027
    struct splite_internal_cache *cache =
 
1028
        (struct splite_internal_cache *) p_cache;
 
1029
    GEOSContextHandle_t handle = NULL;
 
1030
    if (cache == NULL)
 
1031
        return -1;
 
1032
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1033
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1034
        return -1;
 
1035
    handle = cache->GEOS_handle;
 
1036
    if (handle == NULL)
 
1037
        return -1;
 
1038
    gaiaResetGeosMsg_r (cache);
 
1039
    if (!geom1 || !geom2)
 
1040
        return -1;
 
1041
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1042
        return -1;
 
1043
 
 
1044
/* quick check based on MBRs comparison */
 
1045
    if (!splite_mbr_overlaps (geom1, geom2))
 
1046
        return 0;
 
1047
 
 
1048
    g1 = gaiaToGeos_r (cache, geom1);
 
1049
    g2 = gaiaToGeos_r (cache, geom2);
 
1050
    ret = GEOSTouches_r (handle, g1, g2);
 
1051
    GEOSGeom_destroy_r (handle, g1);
 
1052
    GEOSGeom_destroy_r (handle, g2);
 
1053
    return ret;
 
1054
}
 
1055
 
 
1056
GAIAGEO_DECLARE int
 
1057
gaiaGeomCollPreparedTouches (const void *p_cache, gaiaGeomCollPtr geom1,
769
1058
                             unsigned char *blob1, int size1,
770
1059
                             gaiaGeomCollPtr geom2, unsigned char *blob2,
771
1060
                             int size2)
778
1067
    GEOSGeometry *g2;
779
1068
    GEOSPreparedGeometry *gPrep;
780
1069
    gaiaGeomCollPtr geom;
 
1070
    GEOSContextHandle_t handle = NULL;
 
1071
    if (cache == NULL)
 
1072
        return -1;
 
1073
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1074
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1075
        return -1;
 
1076
    handle = cache->GEOS_handle;
 
1077
    if (handle == NULL)
 
1078
        return -1;
 
1079
    gaiaResetGeosMsg_r (cache);
781
1080
    if (!geom1 || !geom2)
782
1081
        return -1;
783
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
1082
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
784
1083
        return -1;
785
1084
 
786
1085
/* quick check based on MBRs comparison */
791
1090
    if (evalGeosCache
792
1091
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
793
1092
      {
794
 
          g2 = gaiaToGeos (geom);
795
 
          ret = GEOSPreparedTouches (gPrep, g2);
796
 
          GEOSGeom_destroy (g2);
 
1093
          g2 = gaiaToGeos_r (cache, geom);
 
1094
          ret = GEOSPreparedTouches_r (handle, gPrep, g2);
 
1095
          GEOSGeom_destroy_r (handle, g2);
797
1096
          return ret;
798
1097
      }
799
1098
 
800
 
    g1 = gaiaToGeos (geom1);
801
 
    g2 = gaiaToGeos (geom2);
802
 
    ret = GEOSTouches (g1, g2);
803
 
    GEOSGeom_destroy (g1);
804
 
    GEOSGeom_destroy (g2);
 
1099
    g1 = gaiaToGeos_r (cache, geom1);
 
1100
    g2 = gaiaToGeos_r (cache, geom2);
 
1101
    ret = GEOSTouches_r (handle, g1, g2);
 
1102
    GEOSGeom_destroy_r (handle, g1);
 
1103
    GEOSGeom_destroy_r (handle, g2);
805
1104
    return ret;
806
1105
}
807
1106
 
808
 
#endif /* end GEOS_ADVANCED */
809
 
 
810
1107
GAIAGEO_DECLARE int
811
1108
gaiaGeomCollWithin (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
812
1109
{
814
1111
    int ret;
815
1112
    GEOSGeometry *g1;
816
1113
    GEOSGeometry *g2;
 
1114
    gaiaResetGeosMsg ();
817
1115
    if (!geom1 || !geom2)
818
1116
        return -1;
819
1117
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
831
1129
    return ret;
832
1130
}
833
1131
 
834
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
835
 
 
836
 
GAIAGEO_DECLARE int
837
 
gaiaGeomCollPreparedWithin (void *p_cache, gaiaGeomCollPtr geom1,
 
1132
GAIAGEO_DECLARE int
 
1133
gaiaGeomCollWithin_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1134
                      gaiaGeomCollPtr geom2)
 
1135
{
 
1136
/* checks if GEOM-1 is completely contained within GEOM-2 */
 
1137
    int ret;
 
1138
    GEOSGeometry *g1;
 
1139
    GEOSGeometry *g2;
 
1140
    struct splite_internal_cache *cache =
 
1141
        (struct splite_internal_cache *) p_cache;
 
1142
    GEOSContextHandle_t handle = NULL;
 
1143
    if (cache == NULL)
 
1144
        return -1;
 
1145
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1146
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1147
        return -1;
 
1148
    handle = cache->GEOS_handle;
 
1149
    if (handle == NULL)
 
1150
        return -1;
 
1151
    gaiaResetGeosMsg_r (cache);
 
1152
    if (!geom1 || !geom2)
 
1153
        return -1;
 
1154
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1155
        return -1;
 
1156
 
 
1157
/* quick check based on MBRs comparison */
 
1158
    if (!splite_mbr_within (geom1, geom2))
 
1159
        return 0;
 
1160
 
 
1161
    g1 = gaiaToGeos_r (cache, geom1);
 
1162
    g2 = gaiaToGeos_r (cache, geom2);
 
1163
    ret = GEOSWithin_r (handle, g1, g2);
 
1164
    GEOSGeom_destroy_r (handle, g1);
 
1165
    GEOSGeom_destroy_r (handle, g2);
 
1166
    return ret;
 
1167
}
 
1168
 
 
1169
GAIAGEO_DECLARE int
 
1170
gaiaGeomCollPreparedWithin (const void *p_cache, gaiaGeomCollPtr geom1,
838
1171
                            unsigned char *blob1, int size1,
839
1172
                            gaiaGeomCollPtr geom2, unsigned char *blob2,
840
1173
                            int size2)
847
1180
    GEOSGeometry *g1;
848
1181
    GEOSGeometry *g2;
849
1182
    gaiaGeomCollPtr geom;
 
1183
    GEOSContextHandle_t handle = NULL;
 
1184
    if (cache == NULL)
 
1185
        return -1;
 
1186
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1187
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1188
        return -1;
 
1189
    handle = cache->GEOS_handle;
 
1190
    if (handle == NULL)
 
1191
        return -1;
 
1192
    gaiaResetGeosMsg_r (cache);
850
1193
    if (!geom1 || !geom2)
851
1194
        return -1;
852
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
1195
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
853
1196
        return -1;
854
1197
 
855
1198
/* quick check based on MBRs comparison */
860
1203
    if (evalGeosCache
861
1204
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
862
1205
      {
863
 
          g2 = gaiaToGeos (geom);
 
1206
          g2 = gaiaToGeos_r (cache, geom);
864
1207
          if (geom == geom2)
865
 
              ret = GEOSPreparedWithin (gPrep, g2);
 
1208
              ret = GEOSPreparedWithin_r (handle, gPrep, g2);
866
1209
          else
867
 
              ret = GEOSPreparedContains (gPrep, g2);
868
 
          GEOSGeom_destroy (g2);
 
1210
              ret = GEOSPreparedContains_r (handle, gPrep, g2);
 
1211
          GEOSGeom_destroy_r (handle, g2);
869
1212
          return ret;
870
1213
      }
871
1214
 
872
 
    g1 = gaiaToGeos (geom1);
873
 
    g2 = gaiaToGeos (geom2);
874
 
    ret = GEOSWithin (g1, g2);
875
 
    GEOSGeom_destroy (g1);
876
 
    GEOSGeom_destroy (g2);
 
1215
    g1 = gaiaToGeos_r (cache, geom1);
 
1216
    g2 = gaiaToGeos_r (cache, geom2);
 
1217
    ret = GEOSWithin_r (handle, g1, g2);
 
1218
    GEOSGeom_destroy_r (handle, g1);
 
1219
    GEOSGeom_destroy_r (handle, g2);
877
1220
    return ret;
878
1221
}
879
1222
 
880
 
#endif /* end GEOS_ADVANCED */
881
 
 
882
1223
GAIAGEO_DECLARE int
883
1224
gaiaGeomCollContains (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
884
1225
{
886
1227
    int ret;
887
1228
    GEOSGeometry *g1;
888
1229
    GEOSGeometry *g2;
 
1230
    gaiaResetGeosMsg ();
889
1231
    if (!geom1 || !geom2)
890
1232
        return -1;
891
1233
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
903
1245
    return ret;
904
1246
}
905
1247
 
906
 
#ifdef GEOS_ADVANCED            /* only if GEOS advanced features are enable */
907
 
 
908
 
GAIAGEO_DECLARE int
909
 
gaiaGeomCollPreparedContains (void *p_cache, gaiaGeomCollPtr geom1,
 
1248
GAIAGEO_DECLARE int
 
1249
gaiaGeomCollContains_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1250
                        gaiaGeomCollPtr geom2)
 
1251
{
 
1252
/* checks if GEOM-1 completely contains GEOM-2 */
 
1253
    int ret;
 
1254
    GEOSGeometry *g1;
 
1255
    GEOSGeometry *g2;
 
1256
    struct splite_internal_cache *cache =
 
1257
        (struct splite_internal_cache *) p_cache;
 
1258
    GEOSContextHandle_t handle = NULL;
 
1259
    if (cache == NULL)
 
1260
        return -1;
 
1261
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1262
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1263
        return -1;
 
1264
    handle = cache->GEOS_handle;
 
1265
    if (handle == NULL)
 
1266
        return -1;
 
1267
    gaiaResetGeosMsg_r (cache);
 
1268
    if (!geom1 || !geom2)
 
1269
        return -1;
 
1270
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1271
        return -1;
 
1272
 
 
1273
/* quick check based on MBRs comparison */
 
1274
    if (!splite_mbr_contains (geom1, geom2))
 
1275
        return 0;
 
1276
 
 
1277
    g1 = gaiaToGeos_r (cache, geom1);
 
1278
    g2 = gaiaToGeos_r (cache, geom2);
 
1279
    ret = GEOSContains_r (handle, g1, g2);
 
1280
    GEOSGeom_destroy_r (handle, g1);
 
1281
    GEOSGeom_destroy_r (handle, g2);
 
1282
    return ret;
 
1283
}
 
1284
 
 
1285
GAIAGEO_DECLARE int
 
1286
gaiaGeomCollPreparedContains (const void *p_cache, gaiaGeomCollPtr geom1,
910
1287
                              unsigned char *blob1, int size1,
911
1288
                              gaiaGeomCollPtr geom2, unsigned char *blob2,
912
1289
                              int size2)
919
1296
    GEOSGeometry *g1;
920
1297
    GEOSGeometry *g2;
921
1298
    gaiaGeomCollPtr geom;
 
1299
    GEOSContextHandle_t handle = NULL;
 
1300
    if (cache == NULL)
 
1301
        return -1;
 
1302
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1303
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1304
        return -1;
 
1305
    handle = cache->GEOS_handle;
 
1306
    if (handle == NULL)
 
1307
        return -1;
 
1308
    gaiaResetGeosMsg_r (cache);
922
1309
    if (!geom1 || !geom2)
923
1310
        return -1;
924
 
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
 
1311
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
925
1312
        return -1;
926
1313
 
927
1314
/* quick check based on MBRs comparison */
932
1319
    if (evalGeosCache
933
1320
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
934
1321
      {
935
 
          g2 = gaiaToGeos (geom);
 
1322
          g2 = gaiaToGeos_r (cache, geom);
936
1323
          if (geom == geom2)
937
 
              ret = GEOSPreparedContains (gPrep, g2);
 
1324
              ret = GEOSPreparedContains_r (handle, gPrep, g2);
938
1325
          else
939
 
              ret = GEOSPreparedWithin (gPrep, g2);
940
 
          GEOSGeom_destroy (g2);
 
1326
              ret = GEOSPreparedWithin_r (handle, gPrep, g2);
 
1327
          GEOSGeom_destroy_r (handle, g2);
941
1328
          return ret;
942
1329
      }
943
1330
 
944
 
    g1 = gaiaToGeos (geom1);
945
 
    g2 = gaiaToGeos (geom2);
946
 
    ret = GEOSContains (g1, g2);
947
 
    GEOSGeom_destroy (g1);
948
 
    GEOSGeom_destroy (g2);
 
1331
    g1 = gaiaToGeos_r (cache, geom1);
 
1332
    g2 = gaiaToGeos_r (cache, geom2);
 
1333
    ret = GEOSContains_r (handle, g1, g2);
 
1334
    GEOSGeom_destroy_r (handle, g1);
 
1335
    GEOSGeom_destroy_r (handle, g2);
949
1336
    return ret;
950
1337
}
951
1338
 
952
 
#endif /* end GEOS_ADVANCED */
953
 
 
954
1339
GAIAGEO_DECLARE int
955
1340
gaiaGeomCollRelate (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
956
1341
                    const char *pattern)
959
1344
    int ret;
960
1345
    GEOSGeometry *g1;
961
1346
    GEOSGeometry *g2;
 
1347
    gaiaResetGeosMsg ();
962
1348
    if (!geom1 || !geom2)
963
1349
        return -1;
964
1350
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
974
1360
}
975
1361
 
976
1362
GAIAGEO_DECLARE int
 
1363
gaiaGeomCollRelate_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1364
                      gaiaGeomCollPtr geom2, const char *pattern)
 
1365
{
 
1366
/* checks if if GEOM-1 and GEOM-2 have a spatial relationship as specified by the pattern Matrix */
 
1367
    int ret;
 
1368
    GEOSGeometry *g1;
 
1369
    GEOSGeometry *g2;
 
1370
    struct splite_internal_cache *cache =
 
1371
        (struct splite_internal_cache *) p_cache;
 
1372
    GEOSContextHandle_t handle = NULL;
 
1373
    if (cache == NULL)
 
1374
        return -1;
 
1375
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1376
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1377
        return -1;
 
1378
    handle = cache->GEOS_handle;
 
1379
    if (handle == NULL)
 
1380
        return -1;
 
1381
    gaiaResetGeosMsg_r (cache);
 
1382
    if (!geom1 || !geom2)
 
1383
        return -1;
 
1384
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1385
        return -1;
 
1386
    g1 = gaiaToGeos_r (cache, geom1);
 
1387
    g2 = gaiaToGeos_r (cache, geom2);
 
1388
    ret = GEOSRelatePattern_r (handle, g1, g2, pattern);
 
1389
    GEOSGeom_destroy_r (handle, g1);
 
1390
    GEOSGeom_destroy_r (handle, g2);
 
1391
    if (ret == 2)
 
1392
        return -1;
 
1393
    return ret;
 
1394
}
 
1395
 
 
1396
GAIAGEO_DECLARE int
977
1397
gaiaGeomCollLength (gaiaGeomCollPtr geom, double *xlength)
978
1398
{
979
1399
/* computes the total length for this Geometry */
980
1400
    double length;
981
1401
    int ret;
982
1402
    GEOSGeometry *g;
 
1403
    gaiaResetGeosMsg ();
983
1404
    if (!geom)
984
1405
        return 0;
985
1406
    if (gaiaIsToxic (geom))
993
1414
}
994
1415
 
995
1416
GAIAGEO_DECLARE int
 
1417
gaiaGeomCollLength_r (const void *p_cache, gaiaGeomCollPtr geom,
 
1418
                      double *xlength)
 
1419
{
 
1420
/* computes the total length for this Geometry */
 
1421
    double length;
 
1422
    int ret;
 
1423
    GEOSGeometry *g;
 
1424
    struct splite_internal_cache *cache =
 
1425
        (struct splite_internal_cache *) p_cache;
 
1426
    GEOSContextHandle_t handle = NULL;
 
1427
    if (cache == NULL)
 
1428
        return -1;
 
1429
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1430
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1431
        return -1;
 
1432
    handle = cache->GEOS_handle;
 
1433
    if (handle == NULL)
 
1434
        return -1;
 
1435
    gaiaResetGeosMsg_r (cache);
 
1436
    if (!geom)
 
1437
        return 0;
 
1438
    if (gaiaIsToxic_r (cache, geom))
 
1439
        return 0;
 
1440
    g = gaiaToGeos_r (cache, geom);
 
1441
    ret = GEOSLength_r (handle, g, &length);
 
1442
    GEOSGeom_destroy_r (handle, g);
 
1443
    if (ret)
 
1444
        *xlength = length;
 
1445
    return ret;
 
1446
}
 
1447
 
 
1448
GAIAGEO_DECLARE int
996
1449
gaiaGeomCollLengthOrPerimeter (gaiaGeomCollPtr geom, int perimeter,
997
1450
                               double *xlength)
998
1451
{
1003
1456
    int mode = GAIA2GEOS_ONLY_LINESTRINGS;
1004
1457
    if (perimeter)
1005
1458
        mode = GAIA2GEOS_ONLY_POLYGONS;
 
1459
    gaiaResetGeosMsg ();
1006
1460
    if (!geom)
1007
1461
        return 0;
1008
1462
    if (gaiaIsToxic (geom))
1021
1475
}
1022
1476
 
1023
1477
GAIAGEO_DECLARE int
 
1478
gaiaGeomCollLengthOrPerimeter_r (const void *p_cache, gaiaGeomCollPtr geom,
 
1479
                                 int perimeter, double *xlength)
 
1480
{
 
1481
/* computes the total length or perimeter for this Geometry */
 
1482
    double length;
 
1483
    int ret;
 
1484
    int mode = GAIA2GEOS_ONLY_LINESTRINGS;
 
1485
    GEOSGeometry *g;
 
1486
    struct splite_internal_cache *cache =
 
1487
        (struct splite_internal_cache *) p_cache;
 
1488
    GEOSContextHandle_t handle = NULL;
 
1489
    if (cache == NULL)
 
1490
        return -1;
 
1491
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1492
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1493
        return -1;
 
1494
    handle = cache->GEOS_handle;
 
1495
    if (handle == NULL)
 
1496
        return -1;
 
1497
    if (perimeter)
 
1498
        mode = GAIA2GEOS_ONLY_POLYGONS;
 
1499
    gaiaResetGeosMsg_r (cache);
 
1500
    if (!geom)
 
1501
        return 0;
 
1502
    if (gaiaIsToxic_r (cache, geom))
 
1503
        return 0;
 
1504
    g = gaiaToGeosSelective_r (cache, geom, mode);
 
1505
    if (g == NULL)
 
1506
      {
 
1507
          *xlength = 0.0;
 
1508
          return 1;
 
1509
      }
 
1510
    ret = GEOSLength_r (handle, g, &length);
 
1511
    GEOSGeom_destroy_r (handle, g);
 
1512
    if (ret)
 
1513
        *xlength = length;
 
1514
    return ret;
 
1515
}
 
1516
 
 
1517
GAIAGEO_DECLARE int
1024
1518
gaiaGeomCollArea (gaiaGeomCollPtr geom, double *xarea)
1025
1519
{
1026
1520
/* computes the total area for this Geometry */
1027
1521
    double area;
1028
1522
    int ret;
1029
1523
    GEOSGeometry *g;
 
1524
    gaiaResetGeosMsg ();
1030
1525
    if (!geom)
1031
1526
        return 0;
1032
1527
    if (gaiaIsToxic (geom))
1040
1535
}
1041
1536
 
1042
1537
GAIAGEO_DECLARE int
 
1538
gaiaGeomCollArea_r (const void *p_cache, gaiaGeomCollPtr geom, double *xarea)
 
1539
{
 
1540
/* computes the total area for this Geometry */
 
1541
    double area;
 
1542
    int ret;
 
1543
    GEOSGeometry *g;
 
1544
    struct splite_internal_cache *cache =
 
1545
        (struct splite_internal_cache *) p_cache;
 
1546
    GEOSContextHandle_t handle = NULL;
 
1547
    if (cache == NULL)
 
1548
        return -1;
 
1549
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1550
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1551
        return -1;
 
1552
    handle = cache->GEOS_handle;
 
1553
    if (handle == NULL)
 
1554
        return -1;
 
1555
    gaiaResetGeosMsg_r (cache);
 
1556
    if (!geom)
 
1557
        return 0;
 
1558
    if (gaiaIsToxic_r (cache, geom))
 
1559
        return 0;
 
1560
    g = gaiaToGeos_r (cache, geom);
 
1561
    ret = GEOSArea_r (handle, g, &area);
 
1562
    GEOSGeom_destroy_r (handle, g);
 
1563
    if (ret)
 
1564
        *xarea = area;
 
1565
    return ret;
 
1566
}
 
1567
 
 
1568
GAIAGEO_DECLARE int
1043
1569
gaiaGeomCollDistance (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
1044
1570
                      double *xdist)
1045
1571
{
1048
1574
    int ret;
1049
1575
    GEOSGeometry *g1;
1050
1576
    GEOSGeometry *g2;
 
1577
    gaiaResetGeosMsg ();
1051
1578
    if (!geom1 || !geom2)
1052
1579
        return 0;
1053
1580
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
1062
1589
    return ret;
1063
1590
}
1064
1591
 
 
1592
GAIAGEO_DECLARE int
 
1593
gaiaGeomCollDistance_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1594
                        gaiaGeomCollPtr geom2, double *xdist)
 
1595
{
 
1596
/* computes the minimum distance intercurring between GEOM-1 and GEOM-2 */
 
1597
    double dist;
 
1598
    int ret;
 
1599
    GEOSGeometry *g1;
 
1600
    GEOSGeometry *g2;
 
1601
    struct splite_internal_cache *cache =
 
1602
        (struct splite_internal_cache *) p_cache;
 
1603
    GEOSContextHandle_t handle = NULL;
 
1604
    if (cache == NULL)
 
1605
        return -1;
 
1606
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1607
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1608
        return -1;
 
1609
    handle = cache->GEOS_handle;
 
1610
    if (handle == NULL)
 
1611
        return -1;
 
1612
    gaiaResetGeosMsg_r (cache);
 
1613
    if (!geom1 || !geom2)
 
1614
        return 0;
 
1615
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1616
        return 0;
 
1617
    g1 = gaiaToGeos_r (cache, geom1);
 
1618
    g2 = gaiaToGeos_r (cache, geom2);
 
1619
    ret = GEOSDistance_r (handle, g1, g2, &dist);
 
1620
    GEOSGeom_destroy_r (handle, g1);
 
1621
    GEOSGeom_destroy_r (handle, g2);
 
1622
    if (ret)
 
1623
        *xdist = dist;
 
1624
    return ret;
 
1625
}
 
1626
 
1065
1627
GAIAGEO_DECLARE gaiaGeomCollPtr
1066
1628
gaiaGeometryIntersection (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
1067
1629
{
1070
1632
    GEOSGeometry *g1;
1071
1633
    GEOSGeometry *g2;
1072
1634
    GEOSGeometry *g3;
 
1635
    gaiaResetGeosMsg ();
1073
1636
    if (!geom1 || !geom2)
1074
1637
        return NULL;
1075
1638
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
1102
1665
}
1103
1666
 
1104
1667
GAIAGEO_DECLARE gaiaGeomCollPtr
 
1668
gaiaGeometryIntersection_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1669
                            gaiaGeomCollPtr geom2)
 
1670
{
 
1671
/* builds a new geometry representing the "spatial intersection" of GEOM-1 and GEOM-2 */
 
1672
    gaiaGeomCollPtr geo;
 
1673
    GEOSGeometry *g1;
 
1674
    GEOSGeometry *g2;
 
1675
    GEOSGeometry *g3;
 
1676
    struct splite_internal_cache *cache =
 
1677
        (struct splite_internal_cache *) p_cache;
 
1678
    GEOSContextHandle_t handle = NULL;
 
1679
    if (cache == NULL)
 
1680
        return NULL;
 
1681
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1682
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1683
        return NULL;
 
1684
    handle = cache->GEOS_handle;
 
1685
    if (handle == NULL)
 
1686
        return NULL;
 
1687
    gaiaResetGeosMsg_r (cache);
 
1688
    if (!geom1 || !geom2)
 
1689
        return NULL;
 
1690
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1691
        return NULL;
 
1692
 
 
1693
/* quick check based on MBRs comparison */
 
1694
    if (!splite_mbr_overlaps (geom1, geom2))
 
1695
        return NULL;
 
1696
 
 
1697
    g1 = gaiaToGeos_r (cache, geom1);
 
1698
    g2 = gaiaToGeos_r (cache, geom2);
 
1699
    g3 = GEOSIntersection_r (handle, g1, g2);
 
1700
    GEOSGeom_destroy_r (handle, g1);
 
1701
    GEOSGeom_destroy_r (handle, g2);
 
1702
    if (!g3)
 
1703
        return NULL;
 
1704
    if (geom1->DimensionModel == GAIA_XY_Z)
 
1705
        geo = gaiaFromGeos_XYZ_r (cache, g3);
 
1706
    else if (geom1->DimensionModel == GAIA_XY_M)
 
1707
        geo = gaiaFromGeos_XYM_r (cache, g3);
 
1708
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
 
1709
        geo = gaiaFromGeos_XYZM_r (cache, g3);
 
1710
    else
 
1711
        geo = gaiaFromGeos_XY_r (cache, g3);
 
1712
    GEOSGeom_destroy_r (handle, g3);
 
1713
    if (geo == NULL)
 
1714
        return NULL;
 
1715
    geo->Srid = geom1->Srid;
 
1716
    return geo;
 
1717
}
 
1718
 
 
1719
GAIAGEO_DECLARE gaiaGeomCollPtr
1105
1720
gaiaGeometryUnion (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
1106
1721
{
1107
1722
/* builds a new geometry representing the "spatial union" of GEOM-1 and GEOM-2 */
1109
1724
    GEOSGeometry *g1;
1110
1725
    GEOSGeometry *g2;
1111
1726
    GEOSGeometry *g3;
 
1727
    gaiaResetGeosMsg ();
1112
1728
    if (!geom1 || !geom2)
1113
1729
        return NULL;
1114
1730
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
1143
1759
}
1144
1760
 
1145
1761
GAIAGEO_DECLARE gaiaGeomCollPtr
 
1762
gaiaGeometryUnion_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1763
                     gaiaGeomCollPtr geom2)
 
1764
{
 
1765
/* builds a new geometry representing the "spatial union" of GEOM-1 and GEOM-2 */
 
1766
    gaiaGeomCollPtr geo;
 
1767
    GEOSGeometry *g1;
 
1768
    GEOSGeometry *g2;
 
1769
    GEOSGeometry *g3;
 
1770
    struct splite_internal_cache *cache =
 
1771
        (struct splite_internal_cache *) p_cache;
 
1772
    GEOSContextHandle_t handle = NULL;
 
1773
    if (cache == NULL)
 
1774
        return NULL;
 
1775
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1776
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1777
        return NULL;
 
1778
    handle = cache->GEOS_handle;
 
1779
    if (handle == NULL)
 
1780
        return NULL;
 
1781
    gaiaResetGeosMsg_r (cache);
 
1782
    if (!geom1 || !geom2)
 
1783
        return NULL;
 
1784
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
1785
        return NULL;
 
1786
    g1 = gaiaToGeos_r (cache, geom1);
 
1787
    g2 = gaiaToGeos_r (cache, geom2);
 
1788
    g3 = GEOSUnion_r (handle, g1, g2);
 
1789
    GEOSGeom_destroy_r (handle, g1);
 
1790
    GEOSGeom_destroy_r (handle, g2);
 
1791
    if (geom1->DimensionModel == GAIA_XY_Z)
 
1792
        geo = gaiaFromGeos_XYZ_r (cache, g3);
 
1793
    else if (geom1->DimensionModel == GAIA_XY_M)
 
1794
        geo = gaiaFromGeos_XYM_r (cache, g3);
 
1795
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
 
1796
        geo = gaiaFromGeos_XYZM_r (cache, g3);
 
1797
    else
 
1798
        geo = gaiaFromGeos_XY_r (cache, g3);
 
1799
    GEOSGeom_destroy_r (handle, g3);
 
1800
    if (geo == NULL)
 
1801
        return NULL;
 
1802
    geo->Srid = geom1->Srid;
 
1803
    if (geo->DeclaredType == GAIA_POINT &&
 
1804
        geom1->DeclaredType == GAIA_MULTIPOINT)
 
1805
        geo->DeclaredType = GAIA_MULTIPOINT;
 
1806
    if (geo->DeclaredType == GAIA_LINESTRING &&
 
1807
        geom1->DeclaredType == GAIA_MULTILINESTRING)
 
1808
        geo->DeclaredType = GAIA_MULTILINESTRING;
 
1809
    if (geo->DeclaredType == GAIA_POLYGON &&
 
1810
        geom1->DeclaredType == GAIA_MULTIPOLYGON)
 
1811
        geo->DeclaredType = GAIA_MULTIPOLYGON;
 
1812
    return geo;
 
1813
}
 
1814
 
 
1815
GAIAGEO_DECLARE gaiaGeomCollPtr
1146
1816
gaiaUnionCascaded (gaiaGeomCollPtr geom)
1147
1817
{
1148
1818
/* UnionCascaded (single Collection of polygons) */
1155
1825
    gaiaPointPtr pt;
1156
1826
    gaiaLinestringPtr ln;
1157
1827
    gaiaPolygonPtr pg;
 
1828
    gaiaResetGeosMsg ();
1158
1829
    if (!geom)
1159
1830
        return NULL;
1160
1831
    if (gaiaIsToxic (geom))
1205
1876
}
1206
1877
 
1207
1878
GAIAGEO_DECLARE gaiaGeomCollPtr
 
1879
gaiaUnionCascaded_r (const void *p_cache, gaiaGeomCollPtr geom)
 
1880
{
 
1881
/* UnionCascaded (single Collection of polygons) */
 
1882
    GEOSGeometry *g1;
 
1883
    GEOSGeometry *g2;
 
1884
    gaiaGeomCollPtr result;
 
1885
    int pts = 0;
 
1886
    int lns = 0;
 
1887
    int pgs = 0;
 
1888
    gaiaPointPtr pt;
 
1889
    gaiaLinestringPtr ln;
 
1890
    gaiaPolygonPtr pg;
 
1891
    struct splite_internal_cache *cache =
 
1892
        (struct splite_internal_cache *) p_cache;
 
1893
    GEOSContextHandle_t handle = NULL;
 
1894
    if (cache == NULL)
 
1895
        return NULL;
 
1896
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
1897
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
1898
        return NULL;
 
1899
    handle = cache->GEOS_handle;
 
1900
    if (handle == NULL)
 
1901
        return NULL;
 
1902
    gaiaResetGeosMsg_r (cache);
 
1903
    if (!geom)
 
1904
        return NULL;
 
1905
    if (gaiaIsToxic_r (cache, geom))
 
1906
        return NULL;
 
1907
 
 
1908
/* testing if geom only contains Polygons */
 
1909
    pt = geom->FirstPoint;
 
1910
    while (pt)
 
1911
      {
 
1912
          pts++;
 
1913
          pt = pt->Next;
 
1914
      }
 
1915
    ln = geom->FirstLinestring;
 
1916
    while (ln)
 
1917
      {
 
1918
          lns++;
 
1919
          ln = ln->Next;
 
1920
      }
 
1921
    pg = geom->FirstPolygon;
 
1922
    while (pg)
 
1923
      {
 
1924
          pgs++;
 
1925
          pg = pg->Next;
 
1926
      }
 
1927
    if (pts || lns)
 
1928
        return NULL;
 
1929
    if (!pgs)
 
1930
        return NULL;
 
1931
 
 
1932
    g1 = gaiaToGeos_r (cache, geom);
 
1933
    g2 = GEOSUnionCascaded_r (handle, g1);
 
1934
    GEOSGeom_destroy_r (handle, g1);
 
1935
    if (!g2)
 
1936
        return NULL;
 
1937
    if (geom->DimensionModel == GAIA_XY_Z)
 
1938
        result = gaiaFromGeos_XYZ_r (cache, g2);
 
1939
    else if (geom->DimensionModel == GAIA_XY_M)
 
1940
        result = gaiaFromGeos_XYM_r (cache, g2);
 
1941
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
1942
        result = gaiaFromGeos_XYZM_r (cache, g2);
 
1943
    else
 
1944
        result = gaiaFromGeos_XY_r (cache, g2);
 
1945
    GEOSGeom_destroy_r (handle, g2);
 
1946
    if (result == NULL)
 
1947
        return NULL;
 
1948
    result->Srid = geom->Srid;
 
1949
    return result;
 
1950
}
 
1951
 
 
1952
GAIAGEO_DECLARE gaiaGeomCollPtr
1208
1953
gaiaGeometryDifference (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
1209
1954
{
1210
1955
/* builds a new geometry representing the "spatial difference" of GEOM-1 and GEOM-2 */
1212
1957
    GEOSGeometry *g1;
1213
1958
    GEOSGeometry *g2;
1214
1959
    GEOSGeometry *g3;
 
1960
    gaiaResetGeosMsg ();
1215
1961
    if (!geom1 || !geom2)
1216
1962
        return NULL;
1217
1963
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
1239
1985
}
1240
1986
 
1241
1987
GAIAGEO_DECLARE gaiaGeomCollPtr
 
1988
gaiaGeometryDifference_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
1989
                          gaiaGeomCollPtr geom2)
 
1990
{
 
1991
/* builds a new geometry representing the "spatial difference" of GEOM-1 and GEOM-2 */
 
1992
    gaiaGeomCollPtr geo;
 
1993
    GEOSGeometry *g1;
 
1994
    GEOSGeometry *g2;
 
1995
    GEOSGeometry *g3;
 
1996
    struct splite_internal_cache *cache =
 
1997
        (struct splite_internal_cache *) p_cache;
 
1998
    GEOSContextHandle_t handle = NULL;
 
1999
    if (cache == NULL)
 
2000
        return NULL;
 
2001
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2002
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2003
        return NULL;
 
2004
    handle = cache->GEOS_handle;
 
2005
    if (handle == NULL)
 
2006
        return NULL;
 
2007
    gaiaResetGeosMsg_r (cache);
 
2008
    if (!geom1 || !geom2)
 
2009
        return NULL;
 
2010
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
2011
        return NULL;
 
2012
    g1 = gaiaToGeos_r (cache, geom1);
 
2013
    g2 = gaiaToGeos_r (cache, geom2);
 
2014
    g3 = GEOSDifference_r (handle, g1, g2);
 
2015
    GEOSGeom_destroy_r (handle, g1);
 
2016
    GEOSGeom_destroy_r (handle, g2);
 
2017
    if (!g3)
 
2018
        return NULL;
 
2019
    if (geom1->DimensionModel == GAIA_XY_Z)
 
2020
        geo = gaiaFromGeos_XYZ_r (cache, g3);
 
2021
    else if (geom1->DimensionModel == GAIA_XY_M)
 
2022
        geo = gaiaFromGeos_XYM_r (cache, g3);
 
2023
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
 
2024
        geo = gaiaFromGeos_XYZM_r (cache, g3);
 
2025
    else
 
2026
        geo = gaiaFromGeos_XY_r (cache, g3);
 
2027
    GEOSGeom_destroy_r (handle, g3);
 
2028
    if (geo == NULL)
 
2029
        return NULL;
 
2030
    geo->Srid = geom1->Srid;
 
2031
    return geo;
 
2032
}
 
2033
 
 
2034
GAIAGEO_DECLARE gaiaGeomCollPtr
1242
2035
gaiaGeometrySymDifference (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
1243
2036
{
1244
2037
/* builds a new geometry representing the "spatial symmetric difference" of GEOM-1 and GEOM-2 */
1246
2039
    GEOSGeometry *g1;
1247
2040
    GEOSGeometry *g2;
1248
2041
    GEOSGeometry *g3;
 
2042
    gaiaResetGeosMsg ();
1249
2043
    if (!geom1 || !geom2)
1250
2044
        return NULL;
1251
2045
    if (gaiaIsToxic (geom1) || gaiaIsToxic (geom2))
1273
2067
}
1274
2068
 
1275
2069
GAIAGEO_DECLARE gaiaGeomCollPtr
 
2070
gaiaGeometrySymDifference_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
2071
                             gaiaGeomCollPtr geom2)
 
2072
{
 
2073
/* builds a new geometry representing the "spatial symmetric difference" of GEOM-1 and GEOM-2 */
 
2074
    gaiaGeomCollPtr geo;
 
2075
    GEOSGeometry *g1;
 
2076
    GEOSGeometry *g2;
 
2077
    GEOSGeometry *g3;
 
2078
    struct splite_internal_cache *cache =
 
2079
        (struct splite_internal_cache *) p_cache;
 
2080
    GEOSContextHandle_t handle = NULL;
 
2081
    if (cache == NULL)
 
2082
        return NULL;
 
2083
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2084
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2085
        return NULL;
 
2086
    handle = cache->GEOS_handle;
 
2087
    if (handle == NULL)
 
2088
        return NULL;
 
2089
    gaiaResetGeosMsg_r (cache);
 
2090
    if (!geom1 || !geom2)
 
2091
        return NULL;
 
2092
    if (gaiaIsToxic_r (cache, geom1) || gaiaIsToxic_r (cache, geom2))
 
2093
        return NULL;
 
2094
    g1 = gaiaToGeos_r (cache, geom1);
 
2095
    g2 = gaiaToGeos_r (cache, geom2);
 
2096
    g3 = GEOSSymDifference_r (handle, g1, g2);
 
2097
    GEOSGeom_destroy_r (handle, g1);
 
2098
    GEOSGeom_destroy_r (handle, g2);
 
2099
    if (!g3)
 
2100
        return NULL;
 
2101
    if (geom1->DimensionModel == GAIA_XY_Z)
 
2102
        geo = gaiaFromGeos_XYZ_r (cache, g3);
 
2103
    else if (geom1->DimensionModel == GAIA_XY_M)
 
2104
        geo = gaiaFromGeos_XYM_r (cache, g3);
 
2105
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
 
2106
        geo = gaiaFromGeos_XYZM_r (cache, g3);
 
2107
    else
 
2108
        geo = gaiaFromGeos_XY_r (cache, g3);
 
2109
    GEOSGeom_destroy_r (handle, g3);
 
2110
    if (geo == NULL)
 
2111
        return NULL;
 
2112
    geo->Srid = geom1->Srid;
 
2113
    return geo;
 
2114
}
 
2115
 
 
2116
GAIAGEO_DECLARE gaiaGeomCollPtr
1276
2117
gaiaBoundary (gaiaGeomCollPtr geom)
1277
2118
{
1278
2119
/* builds a new geometry representing the combinatorial boundary of GEOM */
1279
2120
    gaiaGeomCollPtr geo;
1280
2121
    GEOSGeometry *g1;
1281
2122
    GEOSGeometry *g2;
 
2123
    gaiaResetGeosMsg ();
1282
2124
    if (!geom)
1283
2125
        return NULL;
1284
2126
    if (gaiaIsToxic (geom))
1303
2145
    return geo;
1304
2146
}
1305
2147
 
 
2148
GAIAGEO_DECLARE gaiaGeomCollPtr
 
2149
gaiaBoundary_r (const void *p_cache, gaiaGeomCollPtr geom)
 
2150
{
 
2151
/* builds a new geometry representing the combinatorial boundary of GEOM */
 
2152
    gaiaGeomCollPtr geo;
 
2153
    GEOSGeometry *g1;
 
2154
    GEOSGeometry *g2;
 
2155
    struct splite_internal_cache *cache =
 
2156
        (struct splite_internal_cache *) p_cache;
 
2157
    GEOSContextHandle_t handle = NULL;
 
2158
    if (cache == NULL)
 
2159
        return NULL;
 
2160
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2161
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2162
        return NULL;
 
2163
    handle = cache->GEOS_handle;
 
2164
    if (handle == NULL)
 
2165
        return NULL;
 
2166
    gaiaResetGeosMsg_r (cache);
 
2167
    if (!geom)
 
2168
        return NULL;
 
2169
    if (gaiaIsToxic_r (cache, geom))
 
2170
        return NULL;
 
2171
    g1 = gaiaToGeos_r (cache, geom);
 
2172
    g2 = GEOSBoundary_r (handle, g1);
 
2173
    GEOSGeom_destroy_r (handle, g1);
 
2174
    if (!g2)
 
2175
        return NULL;
 
2176
    if (geom->DimensionModel == GAIA_XY_Z)
 
2177
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
2178
    else if (geom->DimensionModel == GAIA_XY_M)
 
2179
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
2180
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
2181
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
2182
    else
 
2183
        geo = gaiaFromGeos_XY_r (cache, g2);
 
2184
    GEOSGeom_destroy_r (handle, g2);
 
2185
    if (geo == NULL)
 
2186
        return NULL;
 
2187
    geo->Srid = geom->Srid;
 
2188
    return geo;
 
2189
}
 
2190
 
1306
2191
GAIAGEO_DECLARE int
1307
2192
gaiaGeomCollCentroid (gaiaGeomCollPtr geom, double *x, double *y)
1308
2193
{
1310
2195
    gaiaGeomCollPtr geo;
1311
2196
    GEOSGeometry *g1;
1312
2197
    GEOSGeometry *g2;
 
2198
    gaiaResetGeosMsg ();
1313
2199
    if (!geom)
1314
2200
        return 0;
1315
2201
    if (gaiaIsToxic (geom))
1344
2230
}
1345
2231
 
1346
2232
GAIAGEO_DECLARE int
 
2233
gaiaGeomCollCentroid_r (const void *p_cache, gaiaGeomCollPtr geom, double *x,
 
2234
                        double *y)
 
2235
{
 
2236
/* returns a Point representing the centroid for this Geometry */
 
2237
    gaiaGeomCollPtr geo;
 
2238
    GEOSGeometry *g1;
 
2239
    GEOSGeometry *g2;
 
2240
    struct splite_internal_cache *cache =
 
2241
        (struct splite_internal_cache *) p_cache;
 
2242
    GEOSContextHandle_t handle = NULL;
 
2243
    if (cache == NULL)
 
2244
        return 0;
 
2245
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2246
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2247
        return 0;
 
2248
    handle = cache->GEOS_handle;
 
2249
    if (handle == NULL)
 
2250
        return 0;
 
2251
    gaiaResetGeosMsg_r (cache);
 
2252
    if (!geom)
 
2253
        return 0;
 
2254
    if (gaiaIsToxic_r (cache, geom))
 
2255
      {
 
2256
          return 0;
 
2257
      }
 
2258
    g1 = gaiaToGeos_r (cache, geom);
 
2259
    g2 = GEOSGetCentroid_r (handle, g1);
 
2260
    GEOSGeom_destroy_r (handle, g1);
 
2261
    if (!g2)
 
2262
        return 0;
 
2263
    if (geom->DimensionModel == GAIA_XY_Z)
 
2264
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
2265
    else if (geom->DimensionModel == GAIA_XY_M)
 
2266
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
2267
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
2268
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
2269
    else
 
2270
        geo = gaiaFromGeos_XY_r (cache, g2);
 
2271
    GEOSGeom_destroy_r (handle, g2);
 
2272
    if (geo == NULL)
 
2273
        return 0;
 
2274
    if (geo->FirstPoint)
 
2275
      {
 
2276
          *x = geo->FirstPoint->X;
 
2277
          *y = geo->FirstPoint->Y;
 
2278
          gaiaFreeGeomColl (geo);
 
2279
          return 1;
 
2280
      }
 
2281
    gaiaFreeGeomColl (geo);
 
2282
    return 0;
 
2283
}
 
2284
 
 
2285
GAIAGEO_DECLARE int
1347
2286
gaiaGetPointOnSurface (gaiaGeomCollPtr geom, double *x, double *y)
1348
2287
{
1349
2288
/* returns a Point guaranteed to lie on the Surface */
1350
2289
    gaiaGeomCollPtr geo;
1351
2290
    GEOSGeometry *g1;
1352
2291
    GEOSGeometry *g2;
 
2292
    gaiaResetGeosMsg ();
1353
2293
    if (!geom)
1354
2294
        return 0;
1355
2295
    if (gaiaIsToxic (geom))
1384
2324
}
1385
2325
 
1386
2326
GAIAGEO_DECLARE int
 
2327
gaiaGetPointOnSurface_r (const void *p_cache, gaiaGeomCollPtr geom, double *x,
 
2328
                         double *y)
 
2329
{
 
2330
/* returns a Point guaranteed to lie on the Surface */
 
2331
    gaiaGeomCollPtr geo;
 
2332
    GEOSGeometry *g1;
 
2333
    GEOSGeometry *g2;
 
2334
    struct splite_internal_cache *cache =
 
2335
        (struct splite_internal_cache *) p_cache;
 
2336
    GEOSContextHandle_t handle = NULL;
 
2337
    if (cache == NULL)
 
2338
        return 0;
 
2339
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2340
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2341
        return 0;
 
2342
    handle = cache->GEOS_handle;
 
2343
    if (handle == NULL)
 
2344
        return 0;
 
2345
    gaiaResetGeosMsg_r (cache);
 
2346
    if (!geom)
 
2347
        return 0;
 
2348
    if (gaiaIsToxic_r (cache, geom))
 
2349
      {
 
2350
          return 0;
 
2351
      }
 
2352
    g1 = gaiaToGeos_r (cache, geom);
 
2353
    g2 = GEOSPointOnSurface_r (handle, g1);
 
2354
    GEOSGeom_destroy_r (handle, g1);
 
2355
    if (!g2)
 
2356
        return 0;
 
2357
    if (geom->DimensionModel == GAIA_XY_Z)
 
2358
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
2359
    else if (geom->DimensionModel == GAIA_XY_M)
 
2360
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
2361
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
2362
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
2363
    else
 
2364
        geo = gaiaFromGeos_XY_r (cache, g2);
 
2365
    GEOSGeom_destroy_r (handle, g2);
 
2366
    if (geo == NULL)
 
2367
        return 0;
 
2368
    if (geo->FirstPoint)
 
2369
      {
 
2370
          *x = geo->FirstPoint->X;
 
2371
          *y = geo->FirstPoint->Y;
 
2372
          gaiaFreeGeomColl (geo);
 
2373
          return 1;
 
2374
      }
 
2375
    gaiaFreeGeomColl (geo);
 
2376
    return 0;
 
2377
}
 
2378
 
 
2379
GAIAGEO_DECLARE int
1387
2380
gaiaIsSimple (gaiaGeomCollPtr geom)
1388
2381
{
1389
2382
/* checks if this GEOMETRYCOLLECTION is a simple one */
1390
2383
    int ret;
1391
2384
    GEOSGeometry *g;
 
2385
    gaiaResetGeosMsg ();
1392
2386
    if (!geom)
1393
2387
        return -1;
1394
2388
    if (gaiaIsToxic (geom))
1402
2396
}
1403
2397
 
1404
2398
GAIAGEO_DECLARE int
 
2399
gaiaIsSimple_r (const void *p_cache, gaiaGeomCollPtr geom)
 
2400
{
 
2401
/* checks if this GEOMETRYCOLLECTION is a simple one */
 
2402
    int ret;
 
2403
    GEOSGeometry *g;
 
2404
    struct splite_internal_cache *cache =
 
2405
        (struct splite_internal_cache *) p_cache;
 
2406
    GEOSContextHandle_t handle = NULL;
 
2407
    if (cache == NULL)
 
2408
        return -1;
 
2409
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2410
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2411
        return -1;
 
2412
    handle = cache->GEOS_handle;
 
2413
    if (handle == NULL)
 
2414
        return -1;
 
2415
    gaiaResetGeosMsg_r (cache);
 
2416
    if (!geom)
 
2417
        return -1;
 
2418
    if (gaiaIsToxic_r (cache, geom))
 
2419
        return -1;
 
2420
    g = gaiaToGeos_r (cache, geom);
 
2421
    ret = GEOSisSimple_r (handle, g);
 
2422
    GEOSGeom_destroy_r (handle, g);
 
2423
    if (ret == 2)
 
2424
        return -1;
 
2425
    return ret;
 
2426
}
 
2427
 
 
2428
GAIAGEO_DECLARE int
1405
2429
gaiaIsRing (gaiaLinestringPtr line)
1406
2430
{
1407
2431
/* checks if this LINESTRING can be a valid RING */
1414
2438
    double z;
1415
2439
    double m;
1416
2440
    GEOSGeometry *g;
 
2441
    gaiaResetGeosMsg ();
1417
2442
    if (!line)
1418
2443
        return -1;
1419
2444
    if (line->DimensionModel == GAIA_XY_Z)
1477
2502
}
1478
2503
 
1479
2504
GAIAGEO_DECLARE int
 
2505
gaiaIsRing_r (const void *p_cache, gaiaLinestringPtr line)
 
2506
{
 
2507
/* checks if this LINESTRING can be a valid RING */
 
2508
    gaiaGeomCollPtr geo;
 
2509
    gaiaLinestringPtr line2;
 
2510
    int ret;
 
2511
    int iv;
 
2512
    double x;
 
2513
    double y;
 
2514
    double z;
 
2515
    double m;
 
2516
    GEOSGeometry *g;
 
2517
    struct splite_internal_cache *cache =
 
2518
        (struct splite_internal_cache *) p_cache;
 
2519
    GEOSContextHandle_t handle = NULL;
 
2520
    if (cache == NULL)
 
2521
        return -1;
 
2522
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2523
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2524
        return -1;
 
2525
    handle = cache->GEOS_handle;
 
2526
    if (handle == NULL)
 
2527
        return -1;
 
2528
    gaiaResetGeosMsg_r (cache);
 
2529
    if (!line)
 
2530
        return -1;
 
2531
    if (line->DimensionModel == GAIA_XY_Z)
 
2532
        geo = gaiaAllocGeomCollXYZ ();
 
2533
    else if (line->DimensionModel == GAIA_XY_M)
 
2534
        geo = gaiaAllocGeomCollXYM ();
 
2535
    else if (line->DimensionModel == GAIA_XY_Z_M)
 
2536
        geo = gaiaAllocGeomCollXYZM ();
 
2537
    else
 
2538
        geo = gaiaAllocGeomColl ();
 
2539
    line2 = gaiaAddLinestringToGeomColl (geo, line->Points);
 
2540
    for (iv = 0; iv < line2->Points; iv++)
 
2541
      {
 
2542
          z = 0.0;
 
2543
          m = 0.0;
 
2544
          if (line->DimensionModel == GAIA_XY_Z)
 
2545
            {
 
2546
                gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
 
2547
            }
 
2548
          else if (line->DimensionModel == GAIA_XY_M)
 
2549
            {
 
2550
                gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
 
2551
            }
 
2552
          else if (line->DimensionModel == GAIA_XY_Z_M)
 
2553
            {
 
2554
                gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
 
2555
            }
 
2556
          else
 
2557
            {
 
2558
                gaiaGetPoint (line->Coords, iv, &x, &y);
 
2559
            }
 
2560
          if (line2->DimensionModel == GAIA_XY_Z)
 
2561
            {
 
2562
                gaiaSetPointXYZ (line2->Coords, iv, x, y, z);
 
2563
            }
 
2564
          else if (line2->DimensionModel == GAIA_XY_M)
 
2565
            {
 
2566
                gaiaSetPointXYM (line2->Coords, iv, x, y, m);
 
2567
            }
 
2568
          else if (line2->DimensionModel == GAIA_XY_Z_M)
 
2569
            {
 
2570
                gaiaSetPointXYZM (line2->Coords, iv, x, y, z, m);
 
2571
            }
 
2572
          else
 
2573
            {
 
2574
                gaiaSetPoint (line2->Coords, iv, x, y);
 
2575
            }
 
2576
      }
 
2577
    if (gaiaIsToxic_r (cache, geo))
 
2578
      {
 
2579
          gaiaFreeGeomColl (geo);
 
2580
          return -1;
 
2581
      }
 
2582
    g = gaiaToGeos_r (cache, geo);
 
2583
    gaiaFreeGeomColl (geo);
 
2584
    ret = GEOSisRing_r (handle, g);
 
2585
    GEOSGeom_destroy_r (handle, g);
 
2586
    if (ret == 2)
 
2587
        return -1;
 
2588
    return ret;
 
2589
}
 
2590
 
 
2591
GAIAGEO_DECLARE int
1480
2592
gaiaIsValid (gaiaGeomCollPtr geom)
1481
2593
{
1482
2594
/* checks if this GEOMETRYCOLLECTION is a valid one */
1498
2610
}
1499
2611
 
1500
2612
GAIAGEO_DECLARE int
1501
 
gaiaIsClosedGeom (gaiaGeomCollPtr geom)
 
2613
gaiaIsValid_r (const void *p_cache, gaiaGeomCollPtr geom)
 
2614
{
 
2615
/* checks if this GEOMETRYCOLLECTION is a valid one */
 
2616
    int ret;
 
2617
    GEOSGeometry *g;
 
2618
    struct splite_internal_cache *cache =
 
2619
        (struct splite_internal_cache *) p_cache;
 
2620
    GEOSContextHandle_t handle = NULL;
 
2621
    if (cache == NULL)
 
2622
        return -1;
 
2623
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2624
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2625
        return -1;
 
2626
    handle = cache->GEOS_handle;
 
2627
    if (handle == NULL)
 
2628
        return -1;
 
2629
    gaiaResetGeosMsg_r (cache);
 
2630
    if (!geom)
 
2631
        return -1;
 
2632
    if (gaiaIsToxic_r (cache, geom))
 
2633
        return 0;
 
2634
    if (gaiaIsNotClosedGeomColl_r (cache, geom))
 
2635
        return 0;
 
2636
    g = gaiaToGeos_r (cache, geom);
 
2637
    ret = GEOSisValid_r (handle, g);
 
2638
    GEOSGeom_destroy_r (handle, g);
 
2639
    if (ret == 2)
 
2640
        return -1;
 
2641
    return ret;
 
2642
}
 
2643
 
 
2644
GAIAGEO_DECLARE char *
 
2645
gaiaIsValidReason (gaiaGeomCollPtr geom)
 
2646
{
 
2647
/* return a TEXT string stating if a Geometry is valid
 
2648
/ and if not valid, a reason why */
 
2649
    char *text;
 
2650
    int len;
 
2651
    const char *str;
 
2652
    char *gstr;
 
2653
    GEOSGeometry *g;
 
2654
    gaiaResetGeosMsg ();
 
2655
    if (!geom)
 
2656
      {
 
2657
          str = "Invalid: NULL Geometry";
 
2658
          len = strlen (str);
 
2659
          text = malloc (len + 1);
 
2660
          strcpy (text, str);
 
2661
          return text;
 
2662
      }
 
2663
    if (gaiaIsToxic (geom))
 
2664
      {
 
2665
          str = "Invalid: Toxic Geometry ... too few points";
 
2666
          len = strlen (str);
 
2667
          text = malloc (len + 1);
 
2668
          strcpy (text, str);
 
2669
          return text;
 
2670
      }
 
2671
    if (gaiaIsNotClosedGeomColl (geom))
 
2672
      {
 
2673
          str = "Invalid: Unclosed Rings were detected";
 
2674
          len = strlen (str);
 
2675
          text = malloc (len + 1);
 
2676
          strcpy (text, str);
 
2677
          return text;
 
2678
      }
 
2679
    g = gaiaToGeos (geom);
 
2680
    gstr = GEOSisValidReason (g);
 
2681
    GEOSGeom_destroy (g);
 
2682
    if (gstr == NULL)
 
2683
        return NULL;
 
2684
    len = strlen (gstr);
 
2685
    text = malloc (len + 1);
 
2686
    strcpy (text, gstr);
 
2687
    GEOSFree (gstr);
 
2688
    return text;
 
2689
}
 
2690
 
 
2691
GAIAGEO_DECLARE char *
 
2692
gaiaIsValidReason_r (const void *p_cache, gaiaGeomCollPtr geom)
 
2693
{
 
2694
/* return a TEXT string stating if a Geometry is valid
 
2695
/ and if not valid, a reason why */
 
2696
    char *text;
 
2697
    int len;
 
2698
    const char *str;
 
2699
    char *gstr;
 
2700
    GEOSGeometry *g;
 
2701
    struct splite_internal_cache *cache =
 
2702
        (struct splite_internal_cache *) p_cache;
 
2703
    GEOSContextHandle_t handle = NULL;
 
2704
    if (cache == NULL)
 
2705
        return NULL;
 
2706
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2707
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2708
        return NULL;
 
2709
    handle = cache->GEOS_handle;
 
2710
    if (handle == NULL)
 
2711
        return NULL;
 
2712
    gaiaResetGeosMsg_r (cache);
 
2713
    if (!geom)
 
2714
      {
 
2715
          str = "Invalid: NULL Geometry";
 
2716
          len = strlen (str);
 
2717
          text = malloc (len + 1);
 
2718
          strcpy (text, str);
 
2719
          return text;
 
2720
      }
 
2721
    if (gaiaIsToxic (geom))
 
2722
      {
 
2723
          str = "Invalid: Toxic Geometry ... too few points";
 
2724
          len = strlen (str);
 
2725
          text = malloc (len + 1);
 
2726
          strcpy (text, str);
 
2727
          return text;
 
2728
      }
 
2729
    if (gaiaIsNotClosedGeomColl (geom))
 
2730
      {
 
2731
          str = "Invalid: Unclosed Rings were detected";
 
2732
          len = strlen (str);
 
2733
          text = malloc (len + 1);
 
2734
          strcpy (text, str);
 
2735
          return text;
 
2736
      }
 
2737
    g = gaiaToGeos_r (cache, geom);
 
2738
    gstr = GEOSisValidReason_r (handle, g);
 
2739
    GEOSGeom_destroy_r (handle, g);
 
2740
    if (gstr == NULL)
 
2741
        return NULL;
 
2742
    len = strlen (gstr);
 
2743
    text = malloc (len + 1);
 
2744
    strcpy (text, gstr);
 
2745
    GEOSFree_r (handle, gstr);
 
2746
    return text;
 
2747
}
 
2748
 
 
2749
GAIAGEO_DECLARE gaiaGeomCollPtr
 
2750
gaiaIsValidDetail (gaiaGeomCollPtr geom)
 
2751
{
 
2752
/* return a Geometry detail causing a Geometry to be invalid */
 
2753
    char *reason = NULL;
 
2754
    GEOSGeometry *g;
 
2755
    GEOSGeometry *d = NULL;
 
2756
    gaiaGeomCollPtr detail;
 
2757
    gaiaResetGeosMsg ();
 
2758
    if (!geom)
 
2759
        return NULL;
 
2760
    if (gaiaIsToxic (geom))
 
2761
        return NULL;
 
2762
    if (gaiaIsNotClosedGeomColl (geom))
 
2763
        return NULL;
 
2764
    g = gaiaToGeos (geom);
 
2765
    GEOSisValidDetail (g, 0, &reason, &d);
 
2766
    GEOSGeom_destroy (g);
 
2767
    if (reason != NULL)
 
2768
        GEOSFree (reason);
 
2769
    if (d == NULL)
 
2770
        return NULL;
 
2771
    detail = gaiaFromGeos_XY (d);
 
2772
    GEOSGeom_destroy (d);
 
2773
    return detail;
 
2774
}
 
2775
 
 
2776
GAIAGEO_DECLARE gaiaGeomCollPtr
 
2777
gaiaIsValidDetail_r (const void *p_cache, gaiaGeomCollPtr geom)
 
2778
{
 
2779
/* return a Geometry detail causing a Geometry to be invalid */
 
2780
    char *reason = NULL;
 
2781
    GEOSGeometry *g;
 
2782
    GEOSGeometry *d = NULL;
 
2783
    gaiaGeomCollPtr detail;
 
2784
    struct splite_internal_cache *cache =
 
2785
        (struct splite_internal_cache *) p_cache;
 
2786
    GEOSContextHandle_t handle = NULL;
 
2787
    if (cache == NULL)
 
2788
        return NULL;
 
2789
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2790
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2791
        return NULL;
 
2792
    handle = cache->GEOS_handle;
 
2793
    if (handle == NULL)
 
2794
        return NULL;
 
2795
    gaiaResetGeosMsg_r (cache);
 
2796
    if (!geom)
 
2797
        return NULL;
 
2798
    if (gaiaIsToxic (geom))
 
2799
        return NULL;
 
2800
    if (gaiaIsNotClosedGeomColl (geom))
 
2801
        return NULL;
 
2802
    g = gaiaToGeos_r (cache, geom);
 
2803
    GEOSisValidDetail_r (handle, g, 0, &reason, &d);
 
2804
    GEOSGeom_destroy_r (handle, g);
 
2805
    if (reason != NULL)
 
2806
        GEOSFree_r (handle, reason);
 
2807
    if (d == NULL)
 
2808
        return NULL;
 
2809
    detail = gaiaFromGeos_XY_r (cache, d);
 
2810
    GEOSGeom_destroy_r (handle, d);
 
2811
    return detail;
 
2812
}
 
2813
 
 
2814
GAIAGEO_DECLARE int
 
2815
gaiaIsClosedGeom_r (const void *cache, gaiaGeomCollPtr geom)
1502
2816
{
1503
2817
/* checks if this geometry is a closed linestring (or multilinestring) */
1504
2818
    int ret = 0;
1505
2819
    gaiaLinestringPtr ln;
 
2820
    if (cache != NULL)
 
2821
        gaiaResetGeosMsg_r (cache);
1506
2822
    if (!geom)
1507
2823
        return -1;
1508
 
    if (gaiaIsToxic (geom))
 
2824
    if (cache != NULL)
 
2825
        ret = gaiaIsToxic_r (cache, geom);
 
2826
    else
 
2827
        ret = gaiaIsToxic (geom);
 
2828
    if (ret)
1509
2829
        return 0;
1510
2830
    ln = geom->FirstLinestring;
1511
2831
    while (ln)
1549
2869
    return ret;
1550
2870
}
1551
2871
 
 
2872
GAIAGEO_DECLARE int
 
2873
gaiaIsClosedGeom (gaiaGeomCollPtr geom)
 
2874
{
 
2875
    gaiaResetGeosMsg ();
 
2876
    return gaiaIsClosedGeom_r (NULL, geom);
 
2877
}
 
2878
 
1552
2879
GAIAGEO_DECLARE gaiaGeomCollPtr
1553
2880
gaiaGeomCollSimplify (gaiaGeomCollPtr geom, double tolerance)
1554
2881
{
1556
2883
    gaiaGeomCollPtr geo;
1557
2884
    GEOSGeometry *g1;
1558
2885
    GEOSGeometry *g2;
 
2886
    gaiaResetGeosMsg ();
1559
2887
    if (!geom)
1560
2888
        return NULL;
1561
2889
    if (gaiaIsToxic (geom))
1581
2909
}
1582
2910
 
1583
2911
GAIAGEO_DECLARE gaiaGeomCollPtr
 
2912
gaiaGeomCollSimplify_r (const void *p_cache, gaiaGeomCollPtr geom,
 
2913
                        double tolerance)
 
2914
{
 
2915
/* builds a simplified geometry using the Douglas-Peuker algorihtm */
 
2916
    gaiaGeomCollPtr geo;
 
2917
    GEOSGeometry *g1;
 
2918
    GEOSGeometry *g2;
 
2919
    struct splite_internal_cache *cache =
 
2920
        (struct splite_internal_cache *) p_cache;
 
2921
    GEOSContextHandle_t handle = NULL;
 
2922
    if (cache == NULL)
 
2923
        return NULL;
 
2924
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
2925
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
2926
        return NULL;
 
2927
    handle = cache->GEOS_handle;
 
2928
    if (handle == NULL)
 
2929
        return NULL;
 
2930
    gaiaResetGeosMsg_r (cache);
 
2931
    if (!geom)
 
2932
        return NULL;
 
2933
    if (gaiaIsToxic_r (cache, geom))
 
2934
        return NULL;
 
2935
    g1 = gaiaToGeos_r (cache, geom);
 
2936
    g2 = GEOSSimplify_r (handle, g1, tolerance);
 
2937
    GEOSGeom_destroy_r (handle, g1);
 
2938
    if (!g2)
 
2939
        return NULL;
 
2940
    if (geom->DimensionModel == GAIA_XY_Z)
 
2941
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
2942
    else if (geom->DimensionModel == GAIA_XY_M)
 
2943
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
2944
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
2945
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
2946
    else
 
2947
        geo = gaiaFromGeos_XY_r (cache, g2);
 
2948
    GEOSGeom_destroy_r (handle, g2);
 
2949
    if (geo == NULL)
 
2950
        return NULL;
 
2951
    geo->Srid = geom->Srid;
 
2952
    return geo;
 
2953
}
 
2954
 
 
2955
GAIAGEO_DECLARE gaiaGeomCollPtr
1584
2956
gaiaGeomCollSimplifyPreserveTopology (gaiaGeomCollPtr geom, double tolerance)
1585
2957
{
1586
2958
/* builds a simplified geometry using the Douglas-Peuker algorihtm [preserving topology] */
1587
2959
    gaiaGeomCollPtr geo;
1588
2960
    GEOSGeometry *g1;
1589
2961
    GEOSGeometry *g2;
 
2962
    gaiaResetGeosMsg ();
1590
2963
    if (!geom)
1591
2964
        return NULL;
1592
2965
    if (gaiaIsToxic (geom))
1612
2985
}
1613
2986
 
1614
2987
GAIAGEO_DECLARE gaiaGeomCollPtr
 
2988
gaiaGeomCollSimplifyPreserveTopology_r (const void *p_cache,
 
2989
                                        gaiaGeomCollPtr geom, double tolerance)
 
2990
{
 
2991
/* builds a simplified geometry using the Douglas-Peuker algorihtm [preserving topology] */
 
2992
    gaiaGeomCollPtr geo;
 
2993
    GEOSGeometry *g1;
 
2994
    GEOSGeometry *g2;
 
2995
    struct splite_internal_cache *cache =
 
2996
        (struct splite_internal_cache *) p_cache;
 
2997
    GEOSContextHandle_t handle = NULL;
 
2998
    if (cache == NULL)
 
2999
        return NULL;
 
3000
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3001
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3002
        return NULL;
 
3003
    handle = cache->GEOS_handle;
 
3004
    if (handle == NULL)
 
3005
        return NULL;
 
3006
    gaiaResetGeosMsg_r (cache);
 
3007
    if (!geom)
 
3008
        return NULL;
 
3009
    if (gaiaIsToxic_r (cache, geom))
 
3010
        return NULL;
 
3011
    g1 = gaiaToGeos_r (cache, geom);
 
3012
    g2 = GEOSTopologyPreserveSimplify_r (handle, g1, tolerance);
 
3013
    GEOSGeom_destroy_r (handle, g1);
 
3014
    if (!g2)
 
3015
        return NULL;
 
3016
    if (geom->DimensionModel == GAIA_XY_Z)
 
3017
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
3018
    else if (geom->DimensionModel == GAIA_XY_M)
 
3019
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
3020
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
3021
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
3022
    else
 
3023
        geo = gaiaFromGeos_XY_r (cache, g2);
 
3024
    GEOSGeom_destroy_r (handle, g2);
 
3025
    if (geo == NULL)
 
3026
        return NULL;
 
3027
    geo->Srid = geom->Srid;
 
3028
    return geo;
 
3029
}
 
3030
 
 
3031
GAIAGEO_DECLARE gaiaGeomCollPtr
1615
3032
gaiaConvexHull (gaiaGeomCollPtr geom)
1616
3033
{
1617
3034
/* builds a geometry that is the convex hull of GEOM */
1618
3035
    gaiaGeomCollPtr geo;
1619
3036
    GEOSGeometry *g1;
1620
3037
    GEOSGeometry *g2;
 
3038
    gaiaResetGeosMsg ();
1621
3039
    if (!geom)
1622
3040
        return NULL;
1623
3041
    if (gaiaIsToxic (geom))
1643
3061
}
1644
3062
 
1645
3063
GAIAGEO_DECLARE gaiaGeomCollPtr
 
3064
gaiaConvexHull_r (const void *p_cache, gaiaGeomCollPtr geom)
 
3065
{
 
3066
/* builds a geometry that is the convex hull of GEOM */
 
3067
    gaiaGeomCollPtr geo;
 
3068
    GEOSGeometry *g1;
 
3069
    GEOSGeometry *g2;
 
3070
    struct splite_internal_cache *cache =
 
3071
        (struct splite_internal_cache *) p_cache;
 
3072
    GEOSContextHandle_t handle = NULL;
 
3073
    if (cache == NULL)
 
3074
        return NULL;
 
3075
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3076
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3077
        return NULL;
 
3078
    handle = cache->GEOS_handle;
 
3079
    if (handle == NULL)
 
3080
        return NULL;
 
3081
    gaiaResetGeosMsg_r (cache);
 
3082
    if (!geom)
 
3083
        return NULL;
 
3084
    if (gaiaIsToxic_r (cache, geom))
 
3085
        return NULL;
 
3086
    g1 = gaiaToGeos_r (cache, geom);
 
3087
    g2 = GEOSConvexHull_r (handle, g1);
 
3088
    GEOSGeom_destroy_r (handle, g1);
 
3089
    if (!g2)
 
3090
        return NULL;
 
3091
    if (geom->DimensionModel == GAIA_XY_Z)
 
3092
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
3093
    else if (geom->DimensionModel == GAIA_XY_M)
 
3094
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
3095
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
3096
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
3097
    else
 
3098
        geo = gaiaFromGeos_XY_r (cache, g2);
 
3099
    GEOSGeom_destroy_r (handle, g2);
 
3100
    if (geo == NULL)
 
3101
        return NULL;
 
3102
    geo->Srid = geom->Srid;
 
3103
    return geo;
 
3104
}
 
3105
 
 
3106
GAIAGEO_DECLARE gaiaGeomCollPtr
1646
3107
gaiaGeomCollBuffer (gaiaGeomCollPtr geom, double radius, int points)
1647
3108
{
1648
3109
/* builds a geometry that is the GIS buffer of GEOM */
1649
3110
    gaiaGeomCollPtr geo;
1650
3111
    GEOSGeometry *g1;
1651
3112
    GEOSGeometry *g2;
 
3113
    gaiaResetGeosMsg ();
1652
3114
    if (!geom)
1653
3115
        return NULL;
1654
3116
    if (gaiaIsToxic (geom))
1673
3135
    return geo;
1674
3136
}
1675
3137
 
 
3138
GAIAGEO_DECLARE gaiaGeomCollPtr
 
3139
gaiaGeomCollBuffer_r (const void *p_cache, gaiaGeomCollPtr geom, double radius,
 
3140
                      int points)
 
3141
{
 
3142
/* builds a geometry that is the GIS buffer of GEOM */
 
3143
    gaiaGeomCollPtr geo;
 
3144
    GEOSGeometry *g1;
 
3145
    GEOSGeometry *g2;
 
3146
    struct splite_internal_cache *cache =
 
3147
        (struct splite_internal_cache *) p_cache;
 
3148
    GEOSContextHandle_t handle = NULL;
 
3149
    if (cache == NULL)
 
3150
        return NULL;
 
3151
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3152
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3153
        return NULL;
 
3154
    handle = cache->GEOS_handle;
 
3155
    if (handle == NULL)
 
3156
        return NULL;
 
3157
    gaiaResetGeosMsg_r (cache);
 
3158
    if (!geom)
 
3159
        return NULL;
 
3160
    if (gaiaIsToxic_r (cache, geom))
 
3161
        return NULL;
 
3162
    g1 = gaiaToGeos_r (cache, geom);
 
3163
    g2 = GEOSBuffer_r (handle, g1, radius, points);
 
3164
    GEOSGeom_destroy_r (handle, g1);
 
3165
    if (!g2)
 
3166
        return NULL;
 
3167
    if (geom->DimensionModel == GAIA_XY_Z)
 
3168
        geo = gaiaFromGeos_XYZ_r (cache, g2);
 
3169
    else if (geom->DimensionModel == GAIA_XY_M)
 
3170
        geo = gaiaFromGeos_XYM_r (cache, g2);
 
3171
    else if (geom->DimensionModel == GAIA_XY_Z_M)
 
3172
        geo = gaiaFromGeos_XYZM_r (cache, g2);
 
3173
    else
 
3174
        geo = gaiaFromGeos_XY_r (cache, g2);
 
3175
    GEOSGeom_destroy_r (handle, g2);
 
3176
    if (geo == NULL)
 
3177
        return NULL;
 
3178
    geo->Srid = geom->Srid;
 
3179
    return geo;
 
3180
}
 
3181
 
1676
3182
static void
1677
 
auxFromGeosPolygon (const GEOSGeometry * geos, gaiaGeomCollPtr result)
 
3183
auxFromGeosPolygon (GEOSContextHandle_t handle, const GEOSGeometry * geos,
 
3184
                    gaiaGeomCollPtr result)
1678
3185
{
1679
3186
/* converting a Polygon from GEOS to SpatiaLite */
1680
3187
    const GEOSGeometry *geos_ring;
1690
3197
    gaiaPolygonPtr pg;
1691
3198
    gaiaRingPtr rng;
1692
3199
 
1693
 
    geos_ring = GEOSGetExteriorRing (geos);
1694
 
    interiors = GEOSGetNumInteriorRings (geos);
1695
 
    coords = GEOSGeom_getCoordSeq (geos_ring);
1696
 
    GEOSCoordSeq_getDimensions (coords, &geos_dims);
1697
 
    GEOSCoordSeq_getSize (coords, &pts);
 
3200
    if (handle != NULL)
 
3201
      {
 
3202
          geos_ring = GEOSGetExteriorRing_r (handle, geos);
 
3203
          interiors = GEOSGetNumInteriorRings_r (handle, geos);
 
3204
          coords = GEOSGeom_getCoordSeq_r (handle, geos_ring);
 
3205
          GEOSCoordSeq_getDimensions_r (handle, coords, &geos_dims);
 
3206
          GEOSCoordSeq_getSize_r (handle, coords, &pts);
 
3207
      }
 
3208
    else
 
3209
      {
 
3210
          geos_ring = GEOSGetExteriorRing (geos);
 
3211
          interiors = GEOSGetNumInteriorRings (geos);
 
3212
          coords = GEOSGeom_getCoordSeq (geos_ring);
 
3213
          GEOSCoordSeq_getDimensions (coords, &geos_dims);
 
3214
          GEOSCoordSeq_getSize (coords, &pts);
 
3215
      }
1698
3216
 
1699
3217
    pg = gaiaAddPolygonToGeomColl (result, pts, interiors);
1700
3218
/* setting up the Exterior ring */
1703
3221
      {
1704
3222
          if (geos_dims == 3)
1705
3223
            {
1706
 
                GEOSCoordSeq_getX (coords, iv, &x);
1707
 
                GEOSCoordSeq_getY (coords, iv, &y);
1708
 
                GEOSCoordSeq_getZ (coords, iv, &z);
 
3224
                if (handle != NULL)
 
3225
                  {
 
3226
                      GEOSCoordSeq_getX_r (handle, coords, iv, &x);
 
3227
                      GEOSCoordSeq_getY_r (handle, coords, iv, &y);
 
3228
                      GEOSCoordSeq_getZ_r (handle, coords, iv, &z);
 
3229
                  }
 
3230
                else
 
3231
                  {
 
3232
                      GEOSCoordSeq_getX (coords, iv, &x);
 
3233
                      GEOSCoordSeq_getY (coords, iv, &y);
 
3234
                      GEOSCoordSeq_getZ (coords, iv, &z);
 
3235
                  }
1709
3236
            }
1710
3237
          else
1711
3238
            {
1712
 
                GEOSCoordSeq_getX (coords, iv, &x);
1713
 
                GEOSCoordSeq_getY (coords, iv, &y);
 
3239
                if (handle != NULL)
 
3240
                  {
 
3241
                      GEOSCoordSeq_getX_r (handle, coords, iv, &x);
 
3242
                      GEOSCoordSeq_getY_r (handle, coords, iv, &y);
 
3243
                  }
 
3244
                else
 
3245
                  {
 
3246
                      GEOSCoordSeq_getX (coords, iv, &x);
 
3247
                      GEOSCoordSeq_getY (coords, iv, &y);
 
3248
                  }
1714
3249
                z = 0.0;
1715
3250
            }
1716
3251
          if (rng->DimensionModel == GAIA_XY_Z)
1734
3269
    for (ib = 0; ib < interiors; ib++)
1735
3270
      {
1736
3271
          /* setting up any interior ring */
1737
 
          geos_ring = GEOSGetInteriorRingN (geos, ib);
1738
 
          coords = GEOSGeom_getCoordSeq (geos_ring);
1739
 
          GEOSCoordSeq_getDimensions (coords, &geos_dims);
1740
 
          GEOSCoordSeq_getSize (coords, &pts);
 
3272
          if (handle != NULL)
 
3273
            {
 
3274
                geos_ring = GEOSGetInteriorRingN_r (handle, geos, ib);
 
3275
                coords = GEOSGeom_getCoordSeq_r (handle, geos_ring);
 
3276
                GEOSCoordSeq_getDimensions_r (handle, coords, &geos_dims);
 
3277
                GEOSCoordSeq_getSize_r (handle, coords, &pts);
 
3278
            }
 
3279
          else
 
3280
            {
 
3281
                geos_ring = GEOSGetInteriorRingN (geos, ib);
 
3282
                coords = GEOSGeom_getCoordSeq (geos_ring);
 
3283
                GEOSCoordSeq_getDimensions (coords, &geos_dims);
 
3284
                GEOSCoordSeq_getSize (coords, &pts);
 
3285
            }
1741
3286
          rng = gaiaAddInteriorRing (pg, ib, pts);
1742
3287
          for (iv = 0; iv < (int) pts; iv++)
1743
3288
            {
1744
3289
                if (geos_dims == 3)
1745
3290
                  {
1746
 
                      GEOSCoordSeq_getX (coords, iv, &x);
1747
 
                      GEOSCoordSeq_getY (coords, iv, &y);
1748
 
                      GEOSCoordSeq_getZ (coords, iv, &z);
 
3291
                      if (handle != NULL)
 
3292
                        {
 
3293
                            GEOSCoordSeq_getX_r (handle, coords, iv, &x);
 
3294
                            GEOSCoordSeq_getY_r (handle, coords, iv, &y);
 
3295
                            GEOSCoordSeq_getZ_r (handle, coords, iv, &z);
 
3296
                        }
 
3297
                      else
 
3298
                        {
 
3299
                            GEOSCoordSeq_getX (coords, iv, &x);
 
3300
                            GEOSCoordSeq_getY (coords, iv, &y);
 
3301
                            GEOSCoordSeq_getZ (coords, iv, &z);
 
3302
                        }
1749
3303
                  }
1750
3304
                else
1751
3305
                  {
1752
 
                      GEOSCoordSeq_getX (coords, iv, &x);
1753
 
                      GEOSCoordSeq_getY (coords, iv, &y);
 
3306
                      if (handle != NULL)
 
3307
                        {
 
3308
                            GEOSCoordSeq_getX_r (handle, coords, iv, &x);
 
3309
                            GEOSCoordSeq_getY_r (handle, coords, iv, &y);
 
3310
                        }
 
3311
                      else
 
3312
                        {
 
3313
                            GEOSCoordSeq_getX (coords, iv, &x);
 
3314
                            GEOSCoordSeq_getY (coords, iv, &y);
 
3315
                        }
1754
3316
                      z = 0.0;
1755
3317
                  }
1756
3318
                if (rng->DimensionModel == GAIA_XY_Z)
1774
3336
}
1775
3337
 
1776
3338
static void
1777
 
auxGeosMbr (const GEOSCoordSequence * cs, unsigned int pts, double *min_x,
1778
 
            double *min_y, double *max_x, double *max_y)
 
3339
auxGeosMbr (GEOSContextHandle_t handle, const GEOSCoordSequence * cs,
 
3340
            unsigned int pts, double *min_x, double *min_y, double *max_x,
 
3341
            double *max_y)
1779
3342
{
1780
3343
/* computing the MBR */
1781
3344
    int iv;
1787
3350
    *max_y = 0 - DBL_MAX;
1788
3351
    for (iv = 0; iv < (int) pts; iv++)
1789
3352
      {
1790
 
          GEOSCoordSeq_getX (cs, iv, &x);
1791
 
          GEOSCoordSeq_getY (cs, iv, &y);
 
3353
          if (handle != NULL)
 
3354
            {
 
3355
                GEOSCoordSeq_getX_r (handle, cs, iv, &x);
 
3356
                GEOSCoordSeq_getY_r (handle, cs, iv, &y);
 
3357
            }
 
3358
          else
 
3359
            {
 
3360
                GEOSCoordSeq_getX (cs, iv, &x);
 
3361
                GEOSCoordSeq_getY (cs, iv, &y);
 
3362
            }
1792
3363
          if (x < *min_x)
1793
3364
              *min_x = x;
1794
3365
          if (x > *max_x)
1800
3371
      }
1801
3372
}
1802
3373
 
1803
 
GAIAGEO_DECLARE gaiaGeomCollPtr
1804
 
gaiaPolygonize (gaiaGeomCollPtr geom, int force_multi)
 
3374
static gaiaGeomCollPtr
 
3375
gaiaPolygonizeCommon (const void *cache, GEOSContextHandle_t handle,
 
3376
                      gaiaGeomCollPtr geom, int force_multi)
1805
3377
{
1806
3378
/* attempts to rearrange a generic Geometry into a (multi)polygon */
1807
3379
    int ig;
1841
3413
    double max_x2;
1842
3414
    double min_y2;
1843
3415
    double max_y2;
 
3416
    int ret;
1844
3417
 
1845
3418
    if (!geom)
1846
3419
        return NULL;
1847
 
    if (gaiaIsToxic (geom))
 
3420
    if (cache != NULL)
 
3421
        ret = gaiaIsToxic_r (cache, geom);
 
3422
    else
 
3423
        ret = gaiaIsToxic (geom);
 
3424
    if (ret)
1848
3425
        return NULL;
1849
3426
    pt = geom->FirstPoint;
1850
3427
    while (pt)
1886
3463
    ln = geom->FirstLinestring;
1887
3464
    while (ln)
1888
3465
      {
1889
 
          cs = GEOSCoordSeq_create (ln->Points, geos_dims);
 
3466
          if (handle != NULL)
 
3467
              cs = GEOSCoordSeq_create_r (handle, ln->Points, geos_dims);
 
3468
          else
 
3469
              cs = GEOSCoordSeq_create (ln->Points, geos_dims);
1890
3470
          for (iv = 0; iv < ln->Points; iv++)
1891
3471
            {
1892
3472
                /* exterior ring segments */
1909
3489
                  }
1910
3490
                if (geos_dims == 3)
1911
3491
                  {
1912
 
                      GEOSCoordSeq_setX (cs, iv, x);
1913
 
                      GEOSCoordSeq_setY (cs, iv, y);
1914
 
                      GEOSCoordSeq_setZ (cs, iv, z);
 
3492
                      if (handle != NULL)
 
3493
                        {
 
3494
                            GEOSCoordSeq_setX_r (handle, cs, iv, x);
 
3495
                            GEOSCoordSeq_setY_r (handle, cs, iv, y);
 
3496
                            GEOSCoordSeq_setZ_r (handle, cs, iv, z);
 
3497
                        }
 
3498
                      else
 
3499
                        {
 
3500
                            GEOSCoordSeq_setX (cs, iv, x);
 
3501
                            GEOSCoordSeq_setY (cs, iv, y);
 
3502
                            GEOSCoordSeq_setZ (cs, iv, z);
 
3503
                        }
1915
3504
                  }
1916
3505
                else
1917
3506
                  {
1918
 
                      GEOSCoordSeq_setX (cs, iv, x);
1919
 
                      GEOSCoordSeq_setY (cs, iv, y);
 
3507
                      if (handle != NULL)
 
3508
                        {
 
3509
                            GEOSCoordSeq_setX_r (handle, cs, iv, x);
 
3510
                            GEOSCoordSeq_setY_r (handle, cs, iv, y);
 
3511
                        }
 
3512
                      else
 
3513
                        {
 
3514
                            GEOSCoordSeq_setX (cs, iv, x);
 
3515
                            GEOSCoordSeq_setY (cs, iv, y);
 
3516
                        }
1920
3517
                  }
1921
3518
            }
1922
 
          *p_item++ = GEOSGeom_createLineString (cs);
 
3519
          if (handle != NULL)
 
3520
              *p_item++ = GEOSGeom_createLineString_r (handle, cs);
 
3521
          else
 
3522
              *p_item++ = GEOSGeom_createLineString (cs);
1923
3523
          ln = ln->Next;
1924
3524
      }
1925
3525
 
1926
3526
/* calling GEOSPolygonize */
1927
 
    geos = GEOSPolygonize (geos_list, lns);
 
3527
    if (handle != NULL)
 
3528
        geos = GEOSPolygonize_r (handle, geos_list, lns);
 
3529
    else
 
3530
        geos = GEOSPolygonize (geos_list, lns);
1928
3531
    if (geos == NULL)
1929
3532
        goto cleanup;
1930
3533
 
1939
3542
/
1940
3543
*/
1941
3544
    error = 0;
1942
 
    items = GEOSGetNumGeometries (geos);
 
3545
    if (handle != NULL)
 
3546
        items = GEOSGetNumGeometries_r (handle, geos);
 
3547
    else
 
3548
        items = GEOSGetNumGeometries (geos);
1943
3549
    for (ig = 0; ig < items; ig++)
1944
3550
      {
1945
3551
          /* looping on elementaty GEOS geometries */
1946
 
          geos_item = GEOSGetGeometryN (geos, ig);
1947
 
          if (GEOSGeomTypeId (geos_item) != GEOS_POLYGON)
1948
 
            {
1949
 
                /* not a Polygon ... ouch ... */
1950
 
                error = 1;
1951
 
                goto cleanup;
 
3552
          if (handle != NULL)
 
3553
            {
 
3554
                geos_item = GEOSGetGeometryN_r (handle, geos, ig);
 
3555
                if (GEOSGeomTypeId_r (handle, geos_item) != GEOS_POLYGON)
 
3556
                  {
 
3557
                      /* not a Polygon ... ouch ... */
 
3558
                      error = 1;
 
3559
                      goto cleanup;
 
3560
                  }
 
3561
            }
 
3562
          else
 
3563
            {
 
3564
                geos_item = GEOSGetGeometryN (geos, ig);
 
3565
                if (GEOSGeomTypeId (geos_item) != GEOS_POLYGON)
 
3566
                  {
 
3567
                      /* not a Polygon ... ouch ... */
 
3568
                      error = 1;
 
3569
                      goto cleanup;
 
3570
                  }
1952
3571
            }
1953
3572
      }
1954
3573
 
1959
3578
    for (ig = 0; ig < items; ig++)
1960
3579
      {
1961
3580
          /* looping on elementaty GEOS Polygons */
1962
 
          geos_item = GEOSGetGeometryN (geos, ig);
1963
 
          interiors = GEOSGetNumInteriorRings (geos_item);
 
3581
          if (handle != NULL)
 
3582
            {
 
3583
                geos_item = GEOSGetGeometryN_r (handle, geos, ig);
 
3584
                interiors = GEOSGetNumInteriorRings_r (handle, geos_item);
 
3585
            }
 
3586
          else
 
3587
            {
 
3588
                geos_item = GEOSGetGeometryN (geos, ig);
 
3589
                interiors = GEOSGetNumInteriorRings (geos_item);
 
3590
            }
1964
3591
          for (ib = 0; ib < interiors; ib++)
1965
3592
            {
1966
3593
                /* looping on any interior ring */
1967
 
                geos_ring = GEOSGetInteriorRingN (geos_item, ib);
1968
 
                coords = GEOSGeom_getCoordSeq (geos_ring);
1969
 
                GEOSCoordSeq_getSize (coords, &pts1);
1970
 
                auxGeosMbr (coords, pts1, &min_x1, &min_y1, &max_x1, &max_y1);
 
3594
                if (handle != NULL)
 
3595
                  {
 
3596
                      geos_ring =
 
3597
                          GEOSGetInteriorRingN_r (handle, geos_item, ib);
 
3598
                      coords = GEOSGeom_getCoordSeq_r (handle, geos_ring);
 
3599
                      GEOSCoordSeq_getSize_r (handle, coords, &pts1);
 
3600
                  }
 
3601
                else
 
3602
                  {
 
3603
                      geos_ring = GEOSGetInteriorRingN (geos_item, ib);
 
3604
                      coords = GEOSGeom_getCoordSeq (geos_ring);
 
3605
                      GEOSCoordSeq_getSize (coords, &pts1);
 
3606
                  }
 
3607
                auxGeosMbr (handle, coords, pts1, &min_x1, &min_y1, &max_x1,
 
3608
                            &max_y1);
1971
3609
                for (iv = 0; iv < items; iv++)
1972
3610
                  {
1973
3611
                      if (iv == ig)
1980
3618
                            /* skipping any already invalid Polygon */
1981
3619
                            continue;
1982
3620
                        }
1983
 
                      geos_item2 = GEOSGetGeometryN (geos, iv);
1984
 
                      if (GEOSGetNumInteriorRings (geos_item2) > 0)
1985
 
                        {
1986
 
                            /* this Polygon contains holes [surely valid] */
1987
 
                            continue;
1988
 
                        }
1989
 
                      geos_ring = GEOSGetExteriorRing (geos_item2);
1990
 
                      coords = GEOSGeom_getCoordSeq (geos_ring);
1991
 
                      GEOSCoordSeq_getSize (coords, &pts2);
 
3621
                      if (handle != NULL)
 
3622
                        {
 
3623
                            geos_item2 = GEOSGetGeometryN_r (handle, geos, iv);
 
3624
                            if (GEOSGetNumInteriorRings_r (handle, geos_item2) >
 
3625
                                0)
 
3626
                              {
 
3627
                                  /* this Polygon contains holes [surely valid] */
 
3628
                                  continue;
 
3629
                              }
 
3630
                            geos_ring =
 
3631
                                GEOSGetExteriorRing_r (handle, geos_item2);
 
3632
                            coords = GEOSGeom_getCoordSeq_r (handle, geos_ring);
 
3633
                            GEOSCoordSeq_getSize_r (handle, coords, &pts2);
 
3634
                        }
 
3635
                      else
 
3636
                        {
 
3637
                            geos_item2 = GEOSGetGeometryN (geos, iv);
 
3638
                            if (GEOSGetNumInteriorRings (geos_item2) > 0)
 
3639
                              {
 
3640
                                  /* this Polygon contains holes [surely valid] */
 
3641
                                  continue;
 
3642
                              }
 
3643
                            geos_ring = GEOSGetExteriorRing (geos_item2);
 
3644
                            coords = GEOSGeom_getCoordSeq (geos_ring);
 
3645
                            GEOSCoordSeq_getSize (coords, &pts2);
 
3646
                        }
1992
3647
                      if (pts1 == pts2)
1993
3648
                        {
1994
 
                            auxGeosMbr (coords, pts2, &min_x2, &min_y2, &max_x2,
1995
 
                                        &max_y2);
 
3649
                            auxGeosMbr (handle, coords, pts2, &min_x2, &min_y2,
 
3650
                                        &max_x2, &max_y2);
1996
3651
                            if (min_x1 == min_x2 && min_y1 == min_y2
1997
3652
                                && max_x1 == max_x2 && max_y1 == max_y2)
1998
3653
                              {
2022
3677
    for (ig = 0; ig < items; ig++)
2023
3678
      {
2024
3679
          /* looping on GEOS Polygons */
2025
 
          geos_item = GEOSGetGeometryN (geos, ig);
 
3680
          if (handle != NULL)
 
3681
              geos_item = GEOSGetGeometryN_r (handle, geos, ig);
 
3682
          else
 
3683
              geos_item = GEOSGetGeometryN (geos, ig);
2026
3684
          if (valid_polygons[ig] == 'Y')
2027
 
              auxFromGeosPolygon (geos_item, result);
 
3685
              auxFromGeosPolygon (handle, geos_item, result);
2028
3686
      }
2029
3687
 
2030
3688
  cleanup:
2037
3695
          for (iv = 0; iv < lns; iv++)
2038
3696
            {
2039
3697
                if (*p_item != NULL)
2040
 
                    GEOSGeom_destroy (*p_item);
 
3698
                  {
 
3699
                      if (handle != NULL)
 
3700
                          GEOSGeom_destroy_r (handle, *p_item);
 
3701
                      else
 
3702
                          GEOSGeom_destroy (*p_item);
 
3703
                  }
2041
3704
                p_item++;
2042
3705
            }
2043
3706
          p_item = (GEOSGeometry **) geos_list;
2044
3707
          free (p_item);
2045
3708
      }
2046
3709
    if (geos != NULL)
2047
 
        GEOSGeom_destroy (geos);
 
3710
      {
 
3711
          if (handle != NULL)
 
3712
              GEOSGeom_destroy_r (handle, geos);
 
3713
          else
 
3714
              GEOSGeom_destroy (geos);
 
3715
      }
2048
3716
    if (error || result->FirstPolygon == NULL)
2049
3717
      {
2050
3718
          gaiaFreeGeomColl (result);
2053
3721
    return result;
2054
3722
}
2055
3723
 
2056
 
#ifdef GEOS_ADVANCED            /* GEOS advanced features */
2057
 
 
2058
 
GAIAGEO_DECLARE gaiaGeomCollPtr
2059
 
gaiaOffsetCurve (gaiaGeomCollPtr geom, double radius, int points,
2060
 
                 int left_right)
2061
 
{
2062
 
/*
2063
 
// builds a geometry that is the OffsetCurve of GEOM 
2064
 
// (which is expected to be of the LINESTRING type)
2065
 
//
2066
 
*/
2067
 
    gaiaGeomCollPtr geo;
2068
 
    GEOSGeometry *g1;
2069
 
    GEOSGeometry *g2;
2070
 
    gaiaPointPtr pt;
2071
 
    gaiaLinestringPtr ln;
2072
 
    gaiaPolygonPtr pg;
2073
 
    int pts = 0;
2074
 
    int lns = 0;
2075
 
    int pgs = 0;
2076
 
    int closed = 0;
2077
 
    if (!geom)
2078
 
        return NULL;
2079
 
 
2080
 
/* checking the input geometry for validity */
2081
 
    pt = geom->FirstPoint;
2082
 
    while (pt)
2083
 
      {
2084
 
          /* counting how many POINTs are there */
2085
 
          pts++;
2086
 
          pt = pt->Next;
2087
 
      }
2088
 
    ln = geom->FirstLinestring;
2089
 
    while (ln)
2090
 
      {
2091
 
          /* counting how many LINESTRINGs are there */
2092
 
          lns++;
2093
 
          if (gaiaIsClosed (ln))
2094
 
              closed++;
2095
 
          ln = ln->Next;
2096
 
      }
2097
 
    pg = geom->FirstPolygon;
2098
 
    while (pg)
2099
 
      {
2100
 
          /* counting how many POLYGON are there */
2101
 
          pgs++;
2102
 
          pg = pg->Next;
2103
 
      }
2104
 
    if (pts > 0 || pgs > 0 || lns > 1 || closed > 0)
2105
 
        return NULL;
2106
 
 
2107
 
/* all right: this one simply is a LINESTRING */
2108
 
    geom->DeclaredType = GAIA_LINESTRING;
2109
 
 
2110
 
    g1 = gaiaToGeos (geom);
2111
 
    g2 = GEOSSingleSidedBuffer (g1, radius, points, GEOSBUF_JOIN_ROUND, 5.0,
2112
 
                                left_right);
2113
 
    GEOSGeom_destroy (g1);
2114
 
    if (!g2)
2115
 
        return NULL;
2116
 
    if (geom->DimensionModel == GAIA_XY_Z)
2117
 
        geo = gaiaFromGeos_XYZ (g2);
2118
 
    else if (geom->DimensionModel == GAIA_XY_M)
2119
 
        geo = gaiaFromGeos_XYM (g2);
2120
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
2121
 
        geo = gaiaFromGeos_XYZM (g2);
2122
 
    else
2123
 
        geo = gaiaFromGeos_XY (g2);
2124
 
    GEOSGeom_destroy (g2);
2125
 
    if (geo == NULL)
2126
 
        return NULL;
2127
 
    geo->Srid = geom->Srid;
2128
 
    return geo;
2129
 
}
2130
 
 
2131
 
GAIAGEO_DECLARE gaiaGeomCollPtr
2132
 
gaiaSingleSidedBuffer (gaiaGeomCollPtr geom, double radius, int points,
2133
 
                       int left_right)
2134
 
{
2135
 
/*
2136
 
// builds a geometry that is the SingleSided BUFFER of GEOM 
2137
 
// (which is expected to be of the LINESTRING type)
2138
 
//
2139
 
*/
2140
 
    gaiaGeomCollPtr geo;
2141
 
    GEOSGeometry *g1;
2142
 
    GEOSGeometry *g2;
2143
 
    GEOSBufferParams *params = NULL;
2144
 
    gaiaPointPtr pt;
2145
 
    gaiaLinestringPtr ln;
2146
 
    gaiaPolygonPtr pg;
2147
 
    int pts = 0;
2148
 
    int lns = 0;
2149
 
    int pgs = 0;
2150
 
    int closed = 0;
2151
 
    if (!geom)
2152
 
        return NULL;
2153
 
 
2154
 
/* checking the input geometry for validity */
2155
 
    pt = geom->FirstPoint;
2156
 
    while (pt)
2157
 
      {
2158
 
          /* counting how many POINTs are there */
2159
 
          pts++;
2160
 
          pt = pt->Next;
2161
 
      }
2162
 
    ln = geom->FirstLinestring;
2163
 
    while (ln)
2164
 
      {
2165
 
          /* counting how many LINESTRINGs are there */
2166
 
          lns++;
2167
 
          if (gaiaIsClosed (ln))
2168
 
              closed++;
2169
 
          ln = ln->Next;
2170
 
      }
2171
 
    pg = geom->FirstPolygon;
2172
 
    while (pg)
2173
 
      {
2174
 
          /* counting how many POLYGON are there */
2175
 
          pgs++;
2176
 
          pg = pg->Next;
2177
 
      }
2178
 
    if (pts > 0 || pgs > 0 || lns > 1 || closed > 0)
2179
 
        return NULL;
2180
 
 
2181
 
/* all right: this one simply is a LINESTRING */
2182
 
    geom->DeclaredType = GAIA_LINESTRING;
2183
 
 
2184
 
    g1 = gaiaToGeos (geom);
2185
 
/* setting up Buffer params */
2186
 
    params = GEOSBufferParams_create ();
2187
 
    GEOSBufferParams_setJoinStyle (params, GEOSBUF_JOIN_ROUND);
2188
 
    GEOSBufferParams_setMitreLimit (params, 5.0);
2189
 
    GEOSBufferParams_setQuadrantSegments (params, points);
2190
 
    GEOSBufferParams_setSingleSided (params, 1);
2191
 
 
2192
 
/* creating the SingleSided Buffer */
2193
 
    if (left_right == 0)
2194
 
      {
2195
 
          /* right-sided requires NEGATIVE radius */
2196
 
          radius *= -1.0;
2197
 
      }
2198
 
    g2 = GEOSBufferWithParams (g1, params, radius);
2199
 
    GEOSGeom_destroy (g1);
2200
 
    GEOSBufferParams_destroy (params);
2201
 
    if (!g2)
2202
 
        return NULL;
2203
 
    if (geom->DimensionModel == GAIA_XY_Z)
2204
 
        geo = gaiaFromGeos_XYZ (g2);
2205
 
    else if (geom->DimensionModel == GAIA_XY_M)
2206
 
        geo = gaiaFromGeos_XYM (g2);
2207
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
2208
 
        geo = gaiaFromGeos_XYZM (g2);
2209
 
    else
2210
 
        geo = gaiaFromGeos_XY (g2);
2211
 
    GEOSGeom_destroy (g2);
2212
 
    if (geo == NULL)
2213
 
        return NULL;
2214
 
    geo->Srid = geom->Srid;
2215
 
    return geo;
2216
 
}
2217
 
 
2218
 
GAIAGEO_DECLARE int
2219
 
gaiaHausdorffDistance (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2,
2220
 
                       double *xdist)
2221
 
{
2222
 
/* 
2223
 
/ computes the (discrete) Hausdorff distance intercurring 
2224
 
/ between GEOM-1 and GEOM-2 
2225
 
*/
2226
 
    double dist;
2227
 
    int ret;
2228
 
    GEOSGeometry *g1;
2229
 
    GEOSGeometry *g2;
2230
 
    if (!geom1 || !geom2)
2231
 
        return 0;
2232
 
    g1 = gaiaToGeos (geom1);
2233
 
    g2 = gaiaToGeos (geom2);
2234
 
    ret = GEOSHausdorffDistance (g1, g2, &dist);
2235
 
    GEOSGeom_destroy (g1);
2236
 
    GEOSGeom_destroy (g2);
2237
 
    if (ret)
2238
 
        *xdist = dist;
2239
 
    return ret;
2240
 
}
2241
 
 
2242
 
static gaiaGeomCollPtr
2243
 
geom_as_lines (gaiaGeomCollPtr geom)
2244
 
{
2245
 
/* transforms a Geometry into a LINESTRING/MULTILINESTRING (if possible) */
2246
 
    gaiaGeomCollPtr result;
2247
 
    gaiaLinestringPtr ln;
2248
 
    gaiaLinestringPtr new_ln;
2249
 
    gaiaPolygonPtr pg;
2250
 
    gaiaRingPtr rng;
2251
 
    int iv;
2252
 
    int ib;
2253
 
    double x;
2254
 
    double y;
2255
 
    double z;
2256
 
    double m;
2257
 
 
2258
 
    if (!geom)
2259
 
        return NULL;
2260
 
    if (geom->FirstPoint != NULL)
2261
 
      {
2262
 
          /* invalid: GEOM contains at least one POINT */
2263
 
          return NULL;
2264
 
      }
2265
 
 
2266
 
    switch (geom->DimensionModel)
2267
 
      {
2268
 
      case GAIA_XY_Z_M:
2269
 
          result = gaiaAllocGeomCollXYZM ();
2270
 
          break;
2271
 
      case GAIA_XY_Z:
2272
 
          result = gaiaAllocGeomCollXYZ ();
2273
 
          break;
2274
 
      case GAIA_XY_M:
2275
 
          result = gaiaAllocGeomCollXYM ();
2276
 
          break;
2277
 
      default:
2278
 
          result = gaiaAllocGeomColl ();
2279
 
          break;
2280
 
      };
2281
 
    result->Srid = geom->Srid;
2282
 
    ln = geom->FirstLinestring;
2283
 
    while (ln)
2284
 
      {
2285
 
          /* copying any Linestring */
2286
 
          new_ln = gaiaAddLinestringToGeomColl (result, ln->Points);
2287
 
          for (iv = 0; iv < ln->Points; iv++)
2288
 
            {
2289
 
                if (ln->DimensionModel == GAIA_XY_Z)
2290
 
                  {
2291
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2292
 
                      gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
2293
 
                  }
2294
 
                else if (ln->DimensionModel == GAIA_XY_M)
2295
 
                  {
2296
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2297
 
                      gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
2298
 
                  }
2299
 
                else if (ln->DimensionModel == GAIA_XY_Z_M)
2300
 
                  {
2301
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2302
 
                      gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
2303
 
                  }
2304
 
                else
2305
 
                  {
2306
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
2307
 
                      gaiaSetPoint (new_ln->Coords, iv, x, y);
2308
 
                  }
2309
 
            }
2310
 
          ln = ln->Next;
2311
 
      }
2312
 
    pg = geom->FirstPolygon;
2313
 
    while (pg)
2314
 
      {
2315
 
          /* copying any Polygon Ring (as Linestring) */
2316
 
          rng = pg->Exterior;
2317
 
          new_ln = gaiaAddLinestringToGeomColl (result, rng->Points);
2318
 
          for (iv = 0; iv < rng->Points; iv++)
2319
 
            {
2320
 
                /* exterior Ring */
2321
 
                if (rng->DimensionModel == GAIA_XY_Z)
2322
 
                  {
2323
 
                      gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2324
 
                      gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
2325
 
                  }
2326
 
                else if (rng->DimensionModel == GAIA_XY_M)
2327
 
                  {
2328
 
                      gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2329
 
                      gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
2330
 
                  }
2331
 
                else if (rng->DimensionModel == GAIA_XY_Z_M)
2332
 
                  {
2333
 
                      gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2334
 
                      gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
2335
 
                  }
2336
 
                else
2337
 
                  {
2338
 
                      gaiaGetPoint (rng->Coords, iv, &x, &y);
2339
 
                      gaiaSetPoint (new_ln->Coords, iv, x, y);
2340
 
                  }
2341
 
            }
2342
 
          for (ib = 0; ib < pg->NumInteriors; ib++)
2343
 
            {
2344
 
                rng = pg->Interiors + ib;
2345
 
                new_ln = gaiaAddLinestringToGeomColl (result, rng->Points);
2346
 
                for (iv = 0; iv < rng->Points; iv++)
2347
 
                  {
2348
 
                      /* any interior Ring */
2349
 
                      if (rng->DimensionModel == GAIA_XY_Z)
2350
 
                        {
2351
 
                            gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
2352
 
                            gaiaSetPointXYZ (new_ln->Coords, iv, x, y, z);
2353
 
                        }
2354
 
                      else if (rng->DimensionModel == GAIA_XY_M)
2355
 
                        {
2356
 
                            gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
2357
 
                            gaiaSetPointXYM (new_ln->Coords, iv, x, y, m);
2358
 
                        }
2359
 
                      else if (rng->DimensionModel == GAIA_XY_Z_M)
2360
 
                        {
2361
 
                            gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
2362
 
                            gaiaSetPointXYZM (new_ln->Coords, iv, x, y, z, m);
2363
 
                        }
2364
 
                      else
2365
 
                        {
2366
 
                            gaiaGetPoint (rng->Coords, iv, &x, &y);
2367
 
                            gaiaSetPoint (new_ln->Coords, iv, x, y);
2368
 
                        }
2369
 
                  }
2370
 
            }
2371
 
          pg = pg->Next;
2372
 
      }
2373
 
    return result;
2374
 
}
2375
 
 
2376
 
static void
2377
 
add_shared_linestring (gaiaGeomCollPtr geom, gaiaDynamicLinePtr dyn)
2378
 
{
2379
 
/* adding a LINESTRING from Dynamic Line */
2380
 
    int count = 0;
2381
 
    gaiaLinestringPtr ln;
2382
 
    gaiaPointPtr pt;
2383
 
    int iv;
2384
 
 
2385
 
    if (!geom)
2386
 
        return;
2387
 
    if (!dyn)
2388
 
        return;
2389
 
    pt = dyn->First;
2390
 
    while (pt)
2391
 
      {
2392
 
          /* counting how many Points are there */
2393
 
          count++;
2394
 
          pt = pt->Next;
2395
 
      }
2396
 
    if (count == 0)
2397
 
        return;
2398
 
    ln = gaiaAddLinestringToGeomColl (geom, count);
2399
 
    iv = 0;
2400
 
    pt = dyn->First;
2401
 
    while (pt)
2402
 
      {
2403
 
          /* copying points into the LINESTRING */
2404
 
          if (ln->DimensionModel == GAIA_XY_Z)
2405
 
            {
2406
 
                gaiaSetPointXYZ (ln->Coords, iv, pt->X, pt->Y, pt->Z);
2407
 
            }
2408
 
          else if (ln->DimensionModel == GAIA_XY_M)
2409
 
            {
2410
 
                gaiaSetPointXYM (ln->Coords, iv, pt->X, pt->Y, pt->M);
2411
 
            }
2412
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
2413
 
            {
2414
 
                gaiaSetPointXYZM (ln->Coords, iv, pt->X, pt->Y, pt->Z, pt->M);
2415
 
            }
2416
 
          else
2417
 
            {
2418
 
                gaiaSetPoint (ln->Coords, iv, pt->X, pt->Y);
2419
 
            }
2420
 
          iv++;
2421
 
          pt = pt->Next;
2422
 
      }
2423
 
}
2424
 
 
2425
 
static void
2426
 
append_shared_path (gaiaDynamicLinePtr dyn, gaiaLinestringPtr ln, int order)
2427
 
{
2428
 
/* appends a Shared Path item to Dynamic Line */
2429
 
    int iv;
2430
 
    double x;
2431
 
    double y;
2432
 
    double z;
2433
 
    double m;
2434
 
    if (order)
2435
 
      {
2436
 
          /* reversed order */
2437
 
          for (iv = ln->Points - 1; iv >= 0; iv--)
2438
 
            {
2439
 
                if (ln->DimensionModel == GAIA_XY_Z)
2440
 
                  {
2441
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2442
 
                      if (x == dyn->Last->X && y == dyn->Last->Y
2443
 
                          && z == dyn->Last->Z)
2444
 
                          ;
2445
 
                      else
2446
 
                          gaiaAppendPointZToDynamicLine (dyn, x, y, z);
2447
 
                  }
2448
 
                else if (ln->DimensionModel == GAIA_XY_M)
2449
 
                  {
2450
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2451
 
                      if (x == dyn->Last->X && y == dyn->Last->Y
2452
 
                          && m == dyn->Last->M)
2453
 
                          ;
2454
 
                      else
2455
 
                          gaiaAppendPointMToDynamicLine (dyn, x, y, m);
2456
 
                  }
2457
 
                else if (ln->DimensionModel == GAIA_XY_Z_M)
2458
 
                  {
2459
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2460
 
                      if (x == dyn->Last->X && y == dyn->Last->Y
2461
 
                          && z == dyn->Last->Z && m == dyn->Last->M)
2462
 
                          ;
2463
 
                      else
2464
 
                          gaiaAppendPointZMToDynamicLine (dyn, x, y, z, m);
2465
 
                  }
2466
 
                else
2467
 
                  {
2468
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
2469
 
                      if (x == dyn->Last->X && y == dyn->Last->Y)
2470
 
                          ;
2471
 
                      else
2472
 
                          gaiaAppendPointToDynamicLine (dyn, x, y);
2473
 
                  }
2474
 
            }
2475
 
      }
2476
 
    else
2477
 
      {
2478
 
          /* conformant order */
2479
 
          for (iv = 0; iv < ln->Points; iv++)
2480
 
            {
2481
 
                if (ln->DimensionModel == GAIA_XY_Z)
2482
 
                  {
2483
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2484
 
                      if (x == dyn->Last->X && y == dyn->Last->Y
2485
 
                          && z == dyn->Last->Z)
2486
 
                          ;
2487
 
                      else
2488
 
                          gaiaAppendPointZToDynamicLine (dyn, x, y, z);
2489
 
                  }
2490
 
                else if (ln->DimensionModel == GAIA_XY_M)
2491
 
                  {
2492
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2493
 
                      if (x == dyn->Last->X && y == dyn->Last->Y
2494
 
                          && m == dyn->Last->M)
2495
 
                          ;
2496
 
                      else
2497
 
                          gaiaAppendPointMToDynamicLine (dyn, x, y, m);
2498
 
                  }
2499
 
                else if (ln->DimensionModel == GAIA_XY_Z_M)
2500
 
                  {
2501
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2502
 
                      if (x == dyn->Last->X && y == dyn->Last->Y
2503
 
                          && z == dyn->Last->Z && m == dyn->Last->M)
2504
 
                          ;
2505
 
                      else
2506
 
                          gaiaAppendPointZMToDynamicLine (dyn, x, y, z, m);
2507
 
                  }
2508
 
                else
2509
 
                  {
2510
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
2511
 
                      if (x == dyn->Last->X && y == dyn->Last->Y)
2512
 
                          ;
2513
 
                      else
2514
 
                          gaiaAppendPointToDynamicLine (dyn, x, y);
2515
 
                  }
2516
 
            }
2517
 
      }
2518
 
}
2519
 
 
2520
 
static void
2521
 
prepend_shared_path (gaiaDynamicLinePtr dyn, gaiaLinestringPtr ln, int order)
2522
 
{
2523
 
/* prepends a Shared Path item to Dynamic Line */
2524
 
    int iv;
2525
 
    double x;
2526
 
    double y;
2527
 
    double z;
2528
 
    double m;
2529
 
    if (order)
2530
 
      {
2531
 
          /* reversed order */
2532
 
          for (iv = 0; iv < ln->Points; iv++)
2533
 
            {
2534
 
                if (ln->DimensionModel == GAIA_XY_Z)
2535
 
                  {
2536
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2537
 
                      if (x == dyn->First->X && y == dyn->First->Y
2538
 
                          && z == dyn->First->Z)
2539
 
                          ;
2540
 
                      else
2541
 
                          gaiaPrependPointZToDynamicLine (dyn, x, y, z);
2542
 
                  }
2543
 
                else if (ln->DimensionModel == GAIA_XY_M)
2544
 
                  {
2545
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2546
 
                      if (x == dyn->First->X && y == dyn->First->Y
2547
 
                          && m == dyn->First->M)
2548
 
                          ;
2549
 
                      else
2550
 
                          gaiaPrependPointMToDynamicLine (dyn, x, y, m);
2551
 
                  }
2552
 
                else if (ln->DimensionModel == GAIA_XY_Z_M)
2553
 
                  {
2554
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2555
 
                      if (x == dyn->First->X && y == dyn->First->Y
2556
 
                          && z == dyn->First->Z && m == dyn->First->M)
2557
 
                          ;
2558
 
                      else
2559
 
                          gaiaPrependPointZMToDynamicLine (dyn, x, y, z, m);
2560
 
                  }
2561
 
                else
2562
 
                  {
2563
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
2564
 
                      if (x == dyn->First->X && y == dyn->First->Y)
2565
 
                          ;
2566
 
                      else
2567
 
                          gaiaPrependPointToDynamicLine (dyn, x, y);
2568
 
                  }
2569
 
            }
2570
 
      }
2571
 
    else
2572
 
      {
2573
 
          /* conformant order */
2574
 
          for (iv = ln->Points - 1; iv >= 0; iv--)
2575
 
            {
2576
 
                if (ln->DimensionModel == GAIA_XY_Z)
2577
 
                  {
2578
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2579
 
                      if (x == dyn->First->X && y == dyn->First->Y
2580
 
                          && z == dyn->First->Z)
2581
 
                          ;
2582
 
                      else
2583
 
                          gaiaPrependPointZToDynamicLine (dyn, x, y, z);
2584
 
                  }
2585
 
                else if (ln->DimensionModel == GAIA_XY_M)
2586
 
                  {
2587
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2588
 
                      if (x == dyn->First->X && y == dyn->First->Y
2589
 
                          && m == dyn->First->M)
2590
 
                          ;
2591
 
                      else
2592
 
                          gaiaPrependPointMToDynamicLine (dyn, x, y, m);
2593
 
                  }
2594
 
                else if (ln->DimensionModel == GAIA_XY_Z_M)
2595
 
                  {
2596
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
2597
 
                      if (x == dyn->First->X && y == dyn->First->Y
2598
 
                          && z == dyn->First->Z && m == dyn->First->M)
2599
 
                          ;
2600
 
                      else
2601
 
                          gaiaPrependPointZMToDynamicLine (dyn, x, y, z, m);
2602
 
                  }
2603
 
                else
2604
 
                  {
2605
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
2606
 
                      if (x == dyn->First->X && y == dyn->First->Y)
2607
 
                          ;
2608
 
                      else
2609
 
                          gaiaPrependPointToDynamicLine (dyn, x, y);
2610
 
                  }
2611
 
            }
2612
 
      }
2613
 
}
2614
 
 
2615
 
static gaiaGeomCollPtr
2616
 
arrange_shared_paths (gaiaGeomCollPtr geom)
2617
 
{
2618
 
/* final aggregation step for shared paths */
2619
 
    gaiaLinestringPtr ln;
2620
 
    gaiaLinestringPtr *ln_array;
2621
 
    gaiaGeomCollPtr result;
2622
 
    gaiaDynamicLinePtr dyn;
2623
 
    int count;
2624
 
    int i;
2625
 
    int i2;
2626
 
    int iv;
2627
 
    double x;
2628
 
    double y;
2629
 
    double z;
2630
 
    double m;
2631
 
    int ok;
2632
 
    int ok2;
2633
 
 
2634
 
    if (!geom)
2635
 
        return NULL;
2636
 
    count = 0;
2637
 
    ln = geom->FirstLinestring;
2638
 
    while (ln)
2639
 
      {
2640
 
          /* counting how many Linestrings are there */
2641
 
          count++;
2642
 
          ln = ln->Next;
2643
 
      }
2644
 
    if (count == 0)
2645
 
        return NULL;
2646
 
 
2647
 
    ln_array = malloc (sizeof (gaiaLinestringPtr) * count);
2648
 
    i = 0;
2649
 
    ln = geom->FirstLinestring;
2650
 
    while (ln)
2651
 
      {
2652
 
          /* populating the Linestring references array */
2653
 
          ln_array[i++] = ln;
2654
 
          ln = ln->Next;
2655
 
      }
2656
 
 
2657
 
/* allocating a new Geometry [MULTILINESTRING] */
2658
 
    switch (geom->DimensionModel)
2659
 
      {
2660
 
      case GAIA_XY_Z_M:
2661
 
          result = gaiaAllocGeomCollXYZM ();
2662
 
          break;
2663
 
      case GAIA_XY_Z:
2664
 
          result = gaiaAllocGeomCollXYZ ();
2665
 
          break;
2666
 
      case GAIA_XY_M:
2667
 
          result = gaiaAllocGeomCollXYM ();
2668
 
          break;
2669
 
      default:
2670
 
          result = gaiaAllocGeomColl ();
2671
 
          break;
2672
 
      };
2673
 
    result->Srid = geom->Srid;
2674
 
    result->DeclaredType = GAIA_MULTILINESTRING;
2675
 
 
2676
 
    ok = 1;
2677
 
    while (ok)
2678
 
      {
2679
 
          /* looping until we have processed any input item */
2680
 
          ok = 0;
2681
 
          for (i = 0; i < count; i++)
2682
 
            {
2683
 
                if (ln_array[i] != NULL)
2684
 
                  {
2685
 
                      /* starting a new LINESTRING */
2686
 
                      dyn = gaiaAllocDynamicLine ();
2687
 
                      ln = ln_array[i];
2688
 
                      ln_array[i] = NULL;
2689
 
                      for (iv = 0; iv < ln->Points; iv++)
2690
 
                        {
2691
 
                            /* inserting the 'seed' path */
2692
 
                            if (ln->DimensionModel == GAIA_XY_Z)
2693
 
                              {
2694
 
                                  gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
2695
 
                                  gaiaAppendPointZToDynamicLine (dyn, x, y, z);
2696
 
                              }
2697
 
                            else if (ln->DimensionModel == GAIA_XY_M)
2698
 
                              {
2699
 
                                  gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
2700
 
                                  gaiaAppendPointMToDynamicLine (dyn, x, y, m);
2701
 
                              }
2702
 
                            else if (ln->DimensionModel == GAIA_XY_Z_M)
2703
 
                              {
2704
 
                                  gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z,
2705
 
                                                    &m);
2706
 
                                  gaiaAppendPointZMToDynamicLine (dyn, x, y, z,
2707
 
                                                                  m);
2708
 
                              }
2709
 
                            else
2710
 
                              {
2711
 
                                  gaiaGetPoint (ln->Coords, iv, &x, &y);
2712
 
                                  gaiaAppendPointToDynamicLine (dyn, x, y);
2713
 
                              }
2714
 
                        }
2715
 
                      ok2 = 1;
2716
 
                      while (ok2)
2717
 
                        {
2718
 
                            /* looping until we have checked any other item */
2719
 
                            ok2 = 0;
2720
 
                            for (i2 = 0; i2 < count; i2++)
2721
 
                              {
2722
 
                                  /* expanding the 'seed' path */
2723
 
                                  if (ln_array[i2] == NULL)
2724
 
                                      continue;
2725
 
                                  ln = ln_array[i2];
2726
 
                                  /* checking the first vertex */
2727
 
                                  iv = 0;
2728
 
                                  if (ln->DimensionModel == GAIA_XY_Z)
2729
 
                                    {
2730
 
                                        gaiaGetPointXYZ (ln->Coords, iv, &x, &y,
2731
 
                                                         &z);
2732
 
                                    }
2733
 
                                  else if (ln->DimensionModel == GAIA_XY_M)
2734
 
                                    {
2735
 
                                        gaiaGetPointXYM (ln->Coords, iv, &x, &y,
2736
 
                                                         &m);
2737
 
                                    }
2738
 
                                  else if (ln->DimensionModel == GAIA_XY_Z_M)
2739
 
                                    {
2740
 
                                        gaiaGetPointXYZM (ln->Coords, iv, &x,
2741
 
                                                          &y, &z, &m);
2742
 
                                    }
2743
 
                                  else
2744
 
                                    {
2745
 
                                        gaiaGetPoint (ln->Coords, iv, &x, &y);
2746
 
                                    }
2747
 
                                  if (x == dyn->Last->X && y == dyn->Last->Y)
2748
 
                                    {
2749
 
                                        /* appending this item to the 'seed' (conformant order) */
2750
 
                                        append_shared_path (dyn, ln, 0);
2751
 
                                        ln_array[i2] = NULL;
2752
 
                                        ok2 = 1;
2753
 
                                        continue;
2754
 
                                    }
2755
 
                                  if (x == dyn->First->X && y == dyn->First->Y)
2756
 
                                    {
2757
 
                                        /* prepending this item to the 'seed' (reversed order) */
2758
 
                                        prepend_shared_path (dyn, ln, 1);
2759
 
                                        ln_array[i2] = NULL;
2760
 
                                        ok2 = 1;
2761
 
                                        continue;
2762
 
                                    }
2763
 
                                  /* checking the last vertex */
2764
 
                                  iv = ln->Points - 1;
2765
 
                                  if (ln->DimensionModel == GAIA_XY_Z)
2766
 
                                    {
2767
 
                                        gaiaGetPointXYZ (ln->Coords, iv, &x, &y,
2768
 
                                                         &z);
2769
 
                                    }
2770
 
                                  else if (ln->DimensionModel == GAIA_XY_M)
2771
 
                                    {
2772
 
                                        gaiaGetPointXYM (ln->Coords, iv, &x, &y,
2773
 
                                                         &m);
2774
 
                                    }
2775
 
                                  else if (ln->DimensionModel == GAIA_XY_Z_M)
2776
 
                                    {
2777
 
                                        gaiaGetPointXYZM (ln->Coords, iv, &x,
2778
 
                                                          &y, &z, &m);
2779
 
                                    }
2780
 
                                  else
2781
 
                                    {
2782
 
                                        gaiaGetPoint (ln->Coords, iv, &x, &y);
2783
 
                                    }
2784
 
                                  if (x == dyn->Last->X && y == dyn->Last->Y)
2785
 
                                    {
2786
 
                                        /* appending this item to the 'seed' (reversed order) */
2787
 
                                        append_shared_path (dyn, ln, 1);
2788
 
                                        ln_array[i2] = NULL;
2789
 
                                        ok2 = 1;
2790
 
                                        continue;
2791
 
                                    }
2792
 
                                  if (x == dyn->First->X && y == dyn->First->Y)
2793
 
                                    {
2794
 
                                        /* prepending this item to the 'seed' (conformant order) */
2795
 
                                        prepend_shared_path (dyn, ln, 0);
2796
 
                                        ln_array[i2] = NULL;
2797
 
                                        ok2 = 1;
2798
 
                                        continue;
2799
 
                                    }
2800
 
                              }
2801
 
                        }
2802
 
                      add_shared_linestring (result, dyn);
2803
 
                      gaiaFreeDynamicLine (dyn);
2804
 
                      ok = 1;
2805
 
                      break;
2806
 
                  }
2807
 
            }
2808
 
      }
2809
 
    free (ln_array);
2810
 
    return result;
2811
 
}
2812
 
 
2813
 
GAIAGEO_DECLARE gaiaGeomCollPtr
2814
 
gaiaSharedPaths (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
2815
 
{
2816
 
/*
2817
 
// builds a geometry containing Shared Paths commons to GEOM1 & GEOM2 
2818
 
// (which are expected to be of the LINESTRING/MULTILINESTRING type)
2819
 
//
2820
 
*/
2821
 
    gaiaGeomCollPtr geo;
2822
 
    gaiaGeomCollPtr result;
2823
 
    gaiaGeomCollPtr line1;
2824
 
    gaiaGeomCollPtr line2;
2825
 
    GEOSGeometry *g1;
2826
 
    GEOSGeometry *g2;
2827
 
    GEOSGeometry *g3;
2828
 
    if (!geom1)
2829
 
        return NULL;
2830
 
    if (!geom2)
2831
 
        return NULL;
2832
 
/* transforming input geoms as Lines */
2833
 
    line1 = geom_as_lines (geom1);
2834
 
    line2 = geom_as_lines (geom2);
2835
 
    if (line1 == NULL || line2 == NULL)
2836
 
      {
2837
 
          if (line1)
2838
 
              gaiaFreeGeomColl (line1);
2839
 
          if (line2)
2840
 
              gaiaFreeGeomColl (line2);
2841
 
          return NULL;
2842
 
      }
2843
 
 
2844
 
    g1 = gaiaToGeos (line1);
2845
 
    g2 = gaiaToGeos (line2);
2846
 
    gaiaFreeGeomColl (line1);
2847
 
    gaiaFreeGeomColl (line2);
2848
 
    g3 = GEOSSharedPaths (g1, g2);
2849
 
    GEOSGeom_destroy (g1);
2850
 
    GEOSGeom_destroy (g2);
2851
 
    if (!g3)
2852
 
        return NULL;
2853
 
    if (geom1->DimensionModel == GAIA_XY_Z)
2854
 
        geo = gaiaFromGeos_XYZ (g3);
2855
 
    else if (geom1->DimensionModel == GAIA_XY_M)
2856
 
        geo = gaiaFromGeos_XYM (g3);
2857
 
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
2858
 
        geo = gaiaFromGeos_XYZM (g3);
2859
 
    else
2860
 
        geo = gaiaFromGeos_XY (g3);
2861
 
    GEOSGeom_destroy (g3);
2862
 
    if (geo == NULL)
2863
 
        return NULL;
2864
 
    geo->Srid = geom1->Srid;
2865
 
    result = arrange_shared_paths (geo);
2866
 
    gaiaFreeGeomColl (geo);
2867
 
    return result;
 
3724
GAIAGEO_DECLARE gaiaGeomCollPtr
 
3725
gaiaPolygonize (gaiaGeomCollPtr geom, int force_multi)
 
3726
{
 
3727
/* attempts to rearrange a generic Geometry into a (multi)polygon */
 
3728
    gaiaResetGeosMsg ();
 
3729
    return gaiaPolygonizeCommon (NULL, NULL, geom, force_multi);
 
3730
}
 
3731
 
 
3732
GAIAGEO_DECLARE gaiaGeomCollPtr
 
3733
gaiaPolygonize_r (const void *p_cache, gaiaGeomCollPtr geom, int force_multi)
 
3734
{
 
3735
/* attempts to rearrange a generic Geometry into a (multi)polygon */
 
3736
    struct splite_internal_cache *cache =
 
3737
        (struct splite_internal_cache *) p_cache;
 
3738
    GEOSContextHandle_t handle = NULL;
 
3739
    if (cache == NULL)
 
3740
        return NULL;
 
3741
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3742
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3743
        return NULL;
 
3744
    handle = cache->GEOS_handle;
 
3745
    if (handle == NULL)
 
3746
        return NULL;
 
3747
    gaiaResetGeosMsg_r (cache);
 
3748
    return gaiaPolygonizeCommon (cache, handle, geom, force_multi);
2868
3749
}
2869
3750
 
2870
3751
GAIAGEO_DECLARE int
2874
3755
    int ret;
2875
3756
    GEOSGeometry *g1;
2876
3757
    GEOSGeometry *g2;
 
3758
    gaiaResetGeosMsg ();
2877
3759
    if (!geom1 || !geom2)
2878
3760
        return -1;
2879
3761
 
2892
3774
}
2893
3775
 
2894
3776
GAIAGEO_DECLARE int
2895
 
gaiaGeomCollPreparedCovers (void *p_cache, gaiaGeomCollPtr geom1,
 
3777
gaiaGeomCollCovers_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
3778
                      gaiaGeomCollPtr geom2)
 
3779
{
 
3780
/* checks if geom1 "spatially covers" geom2 */
 
3781
    int ret;
 
3782
    GEOSGeometry *g1;
 
3783
    GEOSGeometry *g2;
 
3784
    struct splite_internal_cache *cache =
 
3785
        (struct splite_internal_cache *) p_cache;
 
3786
    GEOSContextHandle_t handle = NULL;
 
3787
    if (cache == NULL)
 
3788
        return -1;
 
3789
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3790
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3791
        return -1;
 
3792
    handle = cache->GEOS_handle;
 
3793
    if (handle == NULL)
 
3794
        return -1;
 
3795
    gaiaResetGeosMsg_r (cache);
 
3796
    if (!geom1 || !geom2)
 
3797
        return -1;
 
3798
 
 
3799
/* quick check based on MBRs comparison */
 
3800
    if (!splite_mbr_contains (geom1, geom2))
 
3801
        return 0;
 
3802
 
 
3803
    g1 = gaiaToGeos_r (cache, geom1);
 
3804
    g2 = gaiaToGeos_r (cache, geom2);
 
3805
    ret = GEOSCovers_r (handle, g1, g2);
 
3806
    GEOSGeom_destroy_r (handle, g1);
 
3807
    GEOSGeom_destroy_r (handle, g2);
 
3808
    if (ret == 2)
 
3809
        return -1;
 
3810
    return ret;
 
3811
}
 
3812
 
 
3813
GAIAGEO_DECLARE int
 
3814
gaiaGeomCollPreparedCovers (const void *p_cache, gaiaGeomCollPtr geom1,
2896
3815
                            unsigned char *blob1, int size1,
2897
3816
                            gaiaGeomCollPtr geom2, unsigned char *blob2,
2898
3817
                            int size2)
2905
3824
    GEOSGeometry *g1;
2906
3825
    GEOSGeometry *g2;
2907
3826
    gaiaGeomCollPtr geom;
 
3827
    GEOSContextHandle_t handle = NULL;
 
3828
    if (cache == NULL)
 
3829
        return -1;
 
3830
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3831
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3832
        return -1;
 
3833
    handle = cache->GEOS_handle;
 
3834
    if (handle == NULL)
 
3835
        return -1;
 
3836
    gaiaResetGeosMsg_r (cache);
2908
3837
    if (!geom1 || !geom2)
2909
3838
        return -1;
2910
3839
 
2916
3845
    if (evalGeosCache
2917
3846
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
2918
3847
      {
2919
 
          g2 = gaiaToGeos (geom);
 
3848
          g2 = gaiaToGeos_r (cache, geom);
2920
3849
          if (geom == geom2)
2921
 
              ret = GEOSPreparedCovers (gPrep, g2);
 
3850
              ret = GEOSPreparedCovers_r (handle, gPrep, g2);
2922
3851
          else
2923
 
              ret = GEOSPreparedCoveredBy (gPrep, g2);
2924
 
          GEOSGeom_destroy (g2);
 
3852
              ret = GEOSPreparedCoveredBy_r (handle, gPrep, g2);
 
3853
          GEOSGeom_destroy_r (handle, g2);
2925
3854
          if (ret == 2)
2926
3855
              return -1;
2927
3856
          return ret;
2928
3857
      }
2929
3858
 
2930
 
    g1 = gaiaToGeos (geom1);
2931
 
    g2 = gaiaToGeos (geom2);
2932
 
    ret = GEOSCovers (g1, g2);
2933
 
    GEOSGeom_destroy (g1);
2934
 
    GEOSGeom_destroy (g2);
 
3859
    g1 = gaiaToGeos_r (cache, geom1);
 
3860
    g2 = gaiaToGeos_r (cache, geom2);
 
3861
    ret = GEOSCovers_r (handle, g1, g2);
 
3862
    GEOSGeom_destroy_r (handle, g1);
 
3863
    GEOSGeom_destroy_r (handle, g2);
2935
3864
    if (ret == 2)
2936
3865
        return -1;
2937
3866
    return ret;
2944
3873
    int ret;
2945
3874
    GEOSGeometry *g1;
2946
3875
    GEOSGeometry *g2;
 
3876
    gaiaResetGeosMsg ();
2947
3877
    if (!geom1 || !geom2)
2948
3878
        return -1;
2949
3879
 
2962
3892
}
2963
3893
 
2964
3894
GAIAGEO_DECLARE int
2965
 
gaiaGeomCollPreparedCoveredBy (void *p_cache, gaiaGeomCollPtr geom1,
 
3895
gaiaGeomCollCoveredBy_r (const void *p_cache, gaiaGeomCollPtr geom1,
 
3896
                         gaiaGeomCollPtr geom2)
 
3897
{
 
3898
/* checks if geom1 is "spatially covered by" geom2 */
 
3899
    int ret;
 
3900
    GEOSGeometry *g1;
 
3901
    GEOSGeometry *g2;
 
3902
    struct splite_internal_cache *cache =
 
3903
        (struct splite_internal_cache *) p_cache;
 
3904
    GEOSContextHandle_t handle = NULL;
 
3905
    if (cache == NULL)
 
3906
        return -1;
 
3907
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3908
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3909
        return -1;
 
3910
    handle = cache->GEOS_handle;
 
3911
    if (handle == NULL)
 
3912
        return -1;
 
3913
    gaiaResetGeosMsg_r (cache);
 
3914
    if (!geom1 || !geom2)
 
3915
        return -1;
 
3916
 
 
3917
/* quick check based on MBRs comparison */
 
3918
    if (!splite_mbr_within (geom1, geom2))
 
3919
        return 0;
 
3920
 
 
3921
    g1 = gaiaToGeos_r (cache, geom1);
 
3922
    g2 = gaiaToGeos_r (cache, geom2);
 
3923
    ret = GEOSCoveredBy_r (handle, g1, g2);
 
3924
    GEOSGeom_destroy_r (handle, g1);
 
3925
    GEOSGeom_destroy_r (handle, g2);
 
3926
    if (ret == 2)
 
3927
        return -1;
 
3928
    return ret;
 
3929
}
 
3930
 
 
3931
GAIAGEO_DECLARE int
 
3932
gaiaGeomCollPreparedCoveredBy (const void *p_cache, gaiaGeomCollPtr geom1,
2966
3933
                               unsigned char *blob1, int size1,
2967
3934
                               gaiaGeomCollPtr geom2, unsigned char *blob2,
2968
3935
                               int size2)
2974
3941
    GEOSPreparedGeometry *gPrep;
2975
3942
    GEOSGeometry *g1;
2976
3943
    GEOSGeometry *g2;
 
3944
    GEOSContextHandle_t handle = NULL;
2977
3945
    gaiaGeomCollPtr geom;
 
3946
    gaiaResetGeosMsg ();
 
3947
    if (cache == NULL)
 
3948
        return -1;
 
3949
    if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
 
3950
        || cache->magic2 != SPATIALITE_CACHE_MAGIC2)
 
3951
        return -1;
 
3952
    handle = cache->GEOS_handle;
 
3953
    if (handle == NULL)
 
3954
        return -1;
 
3955
    gaiaResetGeosMsg_r (cache);
2978
3956
    if (!geom1 || !geom2)
2979
3957
        return -1;
2980
3958
 
2986
3964
    if (evalGeosCache
2987
3965
        (cache, geom1, blob1, size1, geom2, blob2, size2, &gPrep, &geom))
2988
3966
      {
2989
 
          g2 = gaiaToGeos (geom);
 
3967
          g2 = gaiaToGeos_r (cache, geom);
2990
3968
          if (geom == geom2)
2991
 
              ret = GEOSPreparedCoveredBy (gPrep, g2);
 
3969
              ret = GEOSPreparedCoveredBy_r (handle, gPrep, g2);
2992
3970
          else
2993
 
              ret = GEOSPreparedCovers (gPrep, g2);
2994
 
          GEOSGeom_destroy (g2);
 
3971
              ret = GEOSPreparedCovers_r (handle, gPrep, g2);
 
3972
          GEOSGeom_destroy_r (handle, g2);
2995
3973
          if (ret == 2)
2996
3974
              return -1;
2997
3975
          return ret;
2998
3976
      }
2999
3977
 
3000
 
    g1 = gaiaToGeos (geom1);
3001
 
    g2 = gaiaToGeos (geom2);
3002
 
    ret = GEOSCoveredBy (g1, g2);
3003
 
    GEOSGeom_destroy (g1);
3004
 
    GEOSGeom_destroy (g2);
 
3978
    g1 = gaiaToGeos_r (cache, geom1);
 
3979
    g2 = gaiaToGeos_r (cache, geom2);
 
3980
    ret = GEOSCoveredBy_r (handle, g1, g2);
 
3981
    GEOSGeom_destroy_r (handle, g1);
 
3982
    GEOSGeom_destroy_r (handle, g2);
3005
3983
    if (ret == 2)
3006
3984
        return -1;
3007
3985
    return ret;
3008
3986
}
3009
3987
 
3010
 
GAIAGEO_DECLARE gaiaGeomCollPtr
3011
 
gaiaLineInterpolatePoint (gaiaGeomCollPtr geom, double fraction)
3012
 
{
3013
 
/*
3014
 
 * attempts to intepolate a point on line at dist "fraction" 
3015
 
 *
3016
 
 * the fraction is expressed into the range from 0.0 to 1.0
3017
 
 */
3018
 
    int pts = 0;
3019
 
    int lns = 0;
3020
 
    int pgs = 0;
3021
 
    gaiaGeomCollPtr result;
3022
 
    gaiaPointPtr pt;
3023
 
    gaiaLinestringPtr ln;
3024
 
    gaiaPolygonPtr pg;
3025
 
    GEOSGeometry *g;
3026
 
    GEOSGeometry *g_pt;
3027
 
    double length;
3028
 
    double projection;
3029
 
    if (!geom)
3030
 
        return NULL;
3031
 
 
3032
 
/* checking if a single Linestring has been passed */
3033
 
    pt = geom->FirstPoint;
3034
 
    while (pt)
3035
 
      {
3036
 
          pts++;
3037
 
          pt = pt->Next;
3038
 
      }
3039
 
    ln = geom->FirstLinestring;
3040
 
    while (ln)
3041
 
      {
3042
 
          lns++;
3043
 
          ln = ln->Next;
3044
 
      }
3045
 
    pg = geom->FirstPolygon;
3046
 
    while (pg)
3047
 
      {
3048
 
          pgs++;
3049
 
          pg = pg->Next;
3050
 
      }
3051
 
    if (pts == 0 && lns == 1 && pgs == 0)
3052
 
        ;
3053
 
    else
3054
 
        return NULL;
3055
 
 
3056
 
    g = gaiaToGeos (geom);
3057
 
    if (GEOSLength (g, &length))
3058
 
      {
3059
 
          /* transforming fraction to length */
3060
 
          if (fraction < 0.0)
3061
 
              fraction = 0.0;
3062
 
          if (fraction > 1.0)
3063
 
              fraction = 1.0;
3064
 
          projection = length * fraction;
3065
 
      }
3066
 
    else
3067
 
      {
3068
 
          GEOSGeom_destroy (g);
3069
 
          return NULL;
3070
 
      }
3071
 
    g_pt = GEOSInterpolate (g, projection);
3072
 
    GEOSGeom_destroy (g);
3073
 
    if (!g_pt)
3074
 
        return NULL;
3075
 
    if (geom->DimensionModel == GAIA_XY_Z)
3076
 
        result = gaiaFromGeos_XYZ (g_pt);
3077
 
    else if (geom->DimensionModel == GAIA_XY_M)
3078
 
        result = gaiaFromGeos_XYM (g_pt);
3079
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
3080
 
        result = gaiaFromGeos_XYZM (g_pt);
3081
 
    else
3082
 
        result = gaiaFromGeos_XY (g_pt);
3083
 
    GEOSGeom_destroy (g_pt);
3084
 
    if (result == NULL)
3085
 
        return NULL;
3086
 
    result->Srid = geom->Srid;
3087
 
    return result;
3088
 
}
3089
 
 
3090
 
GAIAGEO_DECLARE gaiaGeomCollPtr
3091
 
gaiaLineInterpolateEquidistantPoints (gaiaGeomCollPtr geom, double distance)
3092
 
{
3093
 
/*
3094
 
 * attempts to intepolate a set of points on line at regular distances 
3095
 
 */
3096
 
    int pts = 0;
3097
 
    int lns = 0;
3098
 
    int pgs = 0;
3099
 
    gaiaGeomCollPtr result;
3100
 
    gaiaGeomCollPtr xpt;
3101
 
    gaiaPointPtr pt;
3102
 
    gaiaLinestringPtr ln;
3103
 
    gaiaPolygonPtr pg;
3104
 
    GEOSGeometry *g;
3105
 
    GEOSGeometry *g_pt;
3106
 
    double length;
3107
 
    double current_length = 0.0;
3108
 
    if (!geom)
3109
 
        return NULL;
3110
 
    if (distance <= 0.0)
3111
 
        return NULL;
3112
 
 
3113
 
/* checking if a single Linestring has been passed */
3114
 
    pt = geom->FirstPoint;
3115
 
    while (pt)
3116
 
      {
3117
 
          pts++;
3118
 
          pt = pt->Next;
3119
 
      }
3120
 
    ln = geom->FirstLinestring;
3121
 
    while (ln)
3122
 
      {
3123
 
          lns++;
3124
 
          ln = ln->Next;
3125
 
      }
3126
 
    pg = geom->FirstPolygon;
3127
 
    while (pg)
3128
 
      {
3129
 
          pgs++;
3130
 
          pg = pg->Next;
3131
 
      }
3132
 
    if (pts == 0 && lns == 1 && pgs == 0)
3133
 
        ;
3134
 
    else
3135
 
        return NULL;
3136
 
 
3137
 
    g = gaiaToGeos (geom);
3138
 
    if (GEOSLength (g, &length))
3139
 
      {
3140
 
          if (length <= distance)
3141
 
            {
3142
 
                /* the line is too short to apply interpolation */
3143
 
                GEOSGeom_destroy (g);
3144
 
                return NULL;
3145
 
            }
3146
 
      }
3147
 
    else
3148
 
      {
3149
 
          GEOSGeom_destroy (g);
3150
 
          return NULL;
3151
 
      }
3152
 
 
3153
 
/* creating the MultiPoint [always supporting M] */
3154
 
    if (geom->DimensionModel == GAIA_XY_Z
3155
 
        || geom->DimensionModel == GAIA_XY_Z_M)
3156
 
        result = gaiaAllocGeomCollXYZM ();
3157
 
    else
3158
 
        result = gaiaAllocGeomCollXYM ();
3159
 
    if (result == NULL)
3160
 
      {
3161
 
          GEOSGeom_destroy (g);
3162
 
          return NULL;
3163
 
      }
3164
 
 
3165
 
    while (1)
3166
 
      {
3167
 
          /* increasing the current distance */
3168
 
          current_length += distance;
3169
 
          if (current_length >= length)
3170
 
              break;
3171
 
          /* interpolating a point */
3172
 
          g_pt = GEOSInterpolate (g, current_length);
3173
 
          if (!g_pt)
3174
 
              goto error;
3175
 
          if (geom->DimensionModel == GAIA_XY_Z)
3176
 
            {
3177
 
                xpt = gaiaFromGeos_XYZ (g_pt);
3178
 
                if (!xpt)
3179
 
                    goto error;
3180
 
                pt = xpt->FirstPoint;
3181
 
                if (!pt)
3182
 
                    goto error;
3183
 
                gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z,
3184
 
                                            current_length);
3185
 
            }
3186
 
          else if (geom->DimensionModel == GAIA_XY_M)
3187
 
            {
3188
 
                xpt = gaiaFromGeos_XYM (g_pt);
3189
 
                if (!xpt)
3190
 
                    goto error;
3191
 
                pt = xpt->FirstPoint;
3192
 
                if (!pt)
3193
 
                    goto error;
3194
 
                gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y,
3195
 
                                           current_length);
3196
 
            }
3197
 
          else if (geom->DimensionModel == GAIA_XY_Z_M)
3198
 
            {
3199
 
                xpt = gaiaFromGeos_XYZM (g_pt);
3200
 
                if (!xpt)
3201
 
                    goto error;
3202
 
                pt = xpt->FirstPoint;
3203
 
                if (!pt)
3204
 
                    goto error;
3205
 
                gaiaAddPointToGeomCollXYZM (result, pt->X, pt->Y, pt->Z,
3206
 
                                            current_length);
3207
 
            }
3208
 
          else
3209
 
            {
3210
 
                xpt = gaiaFromGeos_XY (g_pt);
3211
 
                if (!xpt)
3212
 
                    goto error;
3213
 
                pt = xpt->FirstPoint;
3214
 
                if (!pt)
3215
 
                    goto error;
3216
 
                gaiaAddPointToGeomCollXYM (result, pt->X, pt->Y,
3217
 
                                           current_length);
3218
 
            }
3219
 
          GEOSGeom_destroy (g_pt);
3220
 
          gaiaFreeGeomColl (xpt);
3221
 
      }
3222
 
    GEOSGeom_destroy (g);
3223
 
    result->Srid = geom->Srid;
3224
 
    result->DeclaredType = GAIA_MULTIPOINT;
3225
 
    return result;
3226
 
 
3227
 
  error:
3228
 
    if (g_pt)
3229
 
        GEOSGeom_destroy (g_pt);
3230
 
    GEOSGeom_destroy (g);
3231
 
    gaiaFreeGeomColl (result);
3232
 
    return NULL;
3233
 
}
3234
 
 
3235
 
GAIAGEO_DECLARE double
3236
 
gaiaLineLocatePoint (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
3237
 
{
3238
 
/* 
3239
 
 * attempts to compute the location of the closest point on LineString 
3240
 
 * to the given Point, as a fraction of total 2d line length 
3241
 
 *
3242
 
 * the fraction is expressed into the range from 0.0 to 1.0
3243
 
 */
3244
 
    int pts1 = 0;
3245
 
    int lns1 = 0;
3246
 
    int pgs1 = 0;
3247
 
    int pts2 = 0;
3248
 
    int lns2 = 0;
3249
 
    int pgs2 = 0;
3250
 
    double length;
3251
 
    double projection;
3252
 
    double result;
3253
 
    gaiaPointPtr pt;
3254
 
    gaiaLinestringPtr ln;
3255
 
    gaiaPolygonPtr pg;
3256
 
    GEOSGeometry *g1;
3257
 
    GEOSGeometry *g2;
3258
 
    if (!geom1 || !geom2)
3259
 
        return -1.0;
3260
 
 
3261
 
/* checking if a single Linestring has been passed */
3262
 
    pt = geom1->FirstPoint;
3263
 
    while (pt)
3264
 
      {
3265
 
          pts1++;
3266
 
          pt = pt->Next;
3267
 
      }
3268
 
    ln = geom1->FirstLinestring;
3269
 
    while (ln)
3270
 
      {
3271
 
          lns1++;
3272
 
          ln = ln->Next;
3273
 
      }
3274
 
    pg = geom1->FirstPolygon;
3275
 
    while (pg)
3276
 
      {
3277
 
          pgs1++;
3278
 
          pg = pg->Next;
3279
 
      }
3280
 
    if (pts1 == 0 && lns1 >= 1 && pgs1 == 0)
3281
 
        ;
3282
 
    else
3283
 
        return -1.0;
3284
 
 
3285
 
/* checking if a single Point has been passed */
3286
 
    pt = geom2->FirstPoint;
3287
 
    while (pt)
3288
 
      {
3289
 
          pts2++;
3290
 
          pt = pt->Next;
3291
 
      }
3292
 
    ln = geom2->FirstLinestring;
3293
 
    while (ln)
3294
 
      {
3295
 
          lns2++;
3296
 
          ln = ln->Next;
3297
 
      }
3298
 
    pg = geom2->FirstPolygon;
3299
 
    while (pg)
3300
 
      {
3301
 
          pgs2++;
3302
 
          pg = pg->Next;
3303
 
      }
3304
 
    if (pts2 == 1 && lns2 == 0 && pgs2 == 0)
3305
 
        ;
3306
 
    else
3307
 
        return -1.0;
3308
 
 
3309
 
    g1 = gaiaToGeos (geom1);
3310
 
    g2 = gaiaToGeos (geom2);
3311
 
    projection = GEOSProject (g1, g2);
3312
 
    if (GEOSLength (g1, &length))
3313
 
      {
3314
 
          /* normalizing as a fraction between 0.0 and 1.0 */
3315
 
          result = projection / length;
3316
 
      }
3317
 
    else
3318
 
        result = -1.0;
3319
 
    GEOSGeom_destroy (g1);
3320
 
    GEOSGeom_destroy (g2);
3321
 
    return result;
3322
 
}
3323
 
 
3324
 
GAIAGEO_DECLARE gaiaGeomCollPtr
3325
 
gaiaLineSubstring (gaiaGeomCollPtr geom, double start_fraction,
3326
 
                   double end_fraction)
3327
 
{
3328
 
/* 
3329
 
 * attempts to build a new Linestring being a substring of the input one starting 
3330
 
 * and ending at the given fractions of total 2d length 
3331
 
 */
3332
 
    int pts = 0;
3333
 
    int lns = 0;
3334
 
    int pgs = 0;
3335
 
    gaiaGeomCollPtr result;
3336
 
    gaiaPointPtr pt;
3337
 
    gaiaLinestringPtr ln;
3338
 
    gaiaLinestringPtr out;
3339
 
    gaiaPolygonPtr pg;
3340
 
    GEOSGeometry *g;
3341
 
    GEOSGeometry *g_start;
3342
 
    GEOSGeometry *g_end;
3343
 
    GEOSCoordSequence *cs;
3344
 
    const GEOSCoordSequence *in_cs;
3345
 
    GEOSGeometry *segm;
3346
 
    double length;
3347
 
    double total = 0.0;
3348
 
    double start;
3349
 
    double end;
3350
 
    int iv;
3351
 
    int i_start = -1;
3352
 
    int i_end = -1;
3353
 
    int points;
3354
 
    double x;
3355
 
    double y;
3356
 
    double z;
3357
 
    double m;
3358
 
    unsigned int dims;
3359
 
    if (!geom)
3360
 
        return NULL;
3361
 
 
3362
 
/* checking if a single Linestring has been passed */
3363
 
    pt = geom->FirstPoint;
3364
 
    while (pt)
3365
 
      {
3366
 
          pts++;
3367
 
          pt = pt->Next;
3368
 
      }
3369
 
    ln = geom->FirstLinestring;
3370
 
    while (ln)
3371
 
      {
3372
 
          lns++;
3373
 
          ln = ln->Next;
3374
 
      }
3375
 
    pg = geom->FirstPolygon;
3376
 
    while (pg)
3377
 
      {
3378
 
          pgs++;
3379
 
          pg = pg->Next;
3380
 
      }
3381
 
    if (pts == 0 && lns == 1 && pgs == 0)
3382
 
        ;
3383
 
    else
3384
 
        return NULL;
3385
 
 
3386
 
    if (start_fraction < 0.0)
3387
 
        start_fraction = 0.0;
3388
 
    if (start_fraction > 1.0)
3389
 
        start_fraction = 1.0;
3390
 
    if (end_fraction < 0.0)
3391
 
        end_fraction = 0.0;
3392
 
    if (end_fraction > 1.0)
3393
 
        end_fraction = 1.0;
3394
 
    if (start_fraction >= end_fraction)
3395
 
        return NULL;
3396
 
    g = gaiaToGeos (geom);
3397
 
    if (GEOSLength (g, &length))
3398
 
      {
3399
 
          start = length * start_fraction;
3400
 
          end = length * end_fraction;
3401
 
      }
3402
 
    else
3403
 
      {
3404
 
          GEOSGeom_destroy (g);
3405
 
          return NULL;
3406
 
      }
3407
 
    g_start = GEOSInterpolate (g, start);
3408
 
    g_end = GEOSInterpolate (g, end);
3409
 
    GEOSGeom_destroy (g);
3410
 
    if (!g_start || !g_end)
3411
 
        return NULL;
3412
 
 
3413
 
/* identifying first and last valid vertex */
3414
 
    ln = geom->FirstLinestring;
3415
 
    for (iv = 0; iv < ln->Points; iv++)
3416
 
      {
3417
 
 
3418
 
          double x0;
3419
 
          double y0;
3420
 
          switch (ln->DimensionModel)
3421
 
            {
3422
 
            case GAIA_XY_Z:
3423
 
                gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3424
 
                break;
3425
 
            case GAIA_XY_M:
3426
 
                gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3427
 
                break;
3428
 
            case GAIA_XY_Z_M:
3429
 
                gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3430
 
                break;
3431
 
            default:
3432
 
                gaiaGetPoint (ln->Coords, iv, &x, &y);
3433
 
                break;
3434
 
            };
3435
 
 
3436
 
          if (iv > 0)
3437
 
            {
3438
 
                cs = GEOSCoordSeq_create (2, 2);
3439
 
                GEOSCoordSeq_setX (cs, 0, x0);
3440
 
                GEOSCoordSeq_setY (cs, 0, y0);
3441
 
                GEOSCoordSeq_setX (cs, 1, x);
3442
 
                GEOSCoordSeq_setY (cs, 1, y);
3443
 
                segm = GEOSGeom_createLineString (cs);
3444
 
                GEOSLength (segm, &length);
3445
 
                total += length;
3446
 
                GEOSGeom_destroy (segm);
3447
 
                if (total > start && i_start < 0)
3448
 
                    i_start = iv;
3449
 
                if (total < end)
3450
 
                    i_end = iv;
3451
 
            }
3452
 
          x0 = x;
3453
 
          y0 = y;
3454
 
      }
3455
 
    if (i_start < 0 || i_end < 0)
3456
 
      {
3457
 
          i_start = -1;
3458
 
          i_end = -1;
3459
 
          points = 2;
3460
 
      }
3461
 
    else
3462
 
        points = i_end - i_start + 3;
3463
 
 
3464
 
/* creating the output geometry */
3465
 
    switch (ln->DimensionModel)
3466
 
      {
3467
 
      case GAIA_XY_Z:
3468
 
          result = gaiaAllocGeomCollXYZ ();
3469
 
          break;
3470
 
      case GAIA_XY_M:
3471
 
          result = gaiaAllocGeomCollXYM ();
3472
 
          break;
3473
 
      case GAIA_XY_Z_M:
3474
 
          result = gaiaAllocGeomCollXYZM ();
3475
 
          break;
3476
 
      default:
3477
 
          result = gaiaAllocGeomColl ();
3478
 
          break;
3479
 
      };
3480
 
    result->Srid = geom->Srid;
3481
 
    out = gaiaAddLinestringToGeomColl (result, points);
3482
 
 
3483
 
/* start vertex */
3484
 
    points = 0;
3485
 
    in_cs = GEOSGeom_getCoordSeq (g_start);
3486
 
    GEOSCoordSeq_getDimensions (in_cs, &dims);
3487
 
    if (dims == 3)
3488
 
      {
3489
 
          GEOSCoordSeq_getX (in_cs, 0, &x);
3490
 
          GEOSCoordSeq_getY (in_cs, 0, &y);
3491
 
          GEOSCoordSeq_getZ (in_cs, 0, &z);
3492
 
          m = 0.0;
3493
 
      }
3494
 
    else
3495
 
      {
3496
 
          GEOSCoordSeq_getX (in_cs, 0, &x);
3497
 
          GEOSCoordSeq_getY (in_cs, 0, &y);
3498
 
          z = 0.0;
3499
 
          m = 0.0;
3500
 
      }
3501
 
    GEOSGeom_destroy (g_start);
3502
 
    switch (out->DimensionModel)
3503
 
      {
3504
 
      case GAIA_XY_Z:
3505
 
          gaiaSetPointXYZ (out->Coords, points, x, y, z);
3506
 
          break;
3507
 
      case GAIA_XY_M:
3508
 
          gaiaSetPointXYM (out->Coords, points, x, y, 0.0);
3509
 
          break;
3510
 
      case GAIA_XY_Z_M:
3511
 
          gaiaSetPointXYZM (out->Coords, points, x, y, z, 0.0);
3512
 
          break;
3513
 
      default:
3514
 
          gaiaSetPoint (out->Coords, points, x, y);
3515
 
          break;
3516
 
      };
3517
 
    points++;
3518
 
 
3519
 
    if (i_start < 0 || i_end < 0)
3520
 
        ;
3521
 
    else
3522
 
      {
3523
 
          for (iv = i_start; iv <= i_end; iv++)
3524
 
            {
3525
 
                z = 0.0;
3526
 
                m = 0.0;
3527
 
                switch (ln->DimensionModel)
3528
 
                  {
3529
 
                  case GAIA_XY_Z:
3530
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3531
 
                      break;
3532
 
                  case GAIA_XY_M:
3533
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3534
 
                      break;
3535
 
                  case GAIA_XY_Z_M:
3536
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3537
 
                      break;
3538
 
                  default:
3539
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
3540
 
                      break;
3541
 
                  };
3542
 
                switch (out->DimensionModel)
3543
 
                  {
3544
 
                  case GAIA_XY_Z:
3545
 
                      gaiaSetPointXYZ (out->Coords, points, x, y, z);
3546
 
                      break;
3547
 
                  case GAIA_XY_M:
3548
 
                      gaiaSetPointXYM (out->Coords, points, x, y, 0.0);
3549
 
                      break;
3550
 
                  case GAIA_XY_Z_M:
3551
 
                      gaiaSetPointXYZM (out->Coords, points, x, y, z, 0.0);
3552
 
                      break;
3553
 
                  default:
3554
 
                      gaiaSetPoint (out->Coords, points, x, y);
3555
 
                      break;
3556
 
                  };
3557
 
                points++;
3558
 
            }
3559
 
      }
3560
 
 
3561
 
/* end vertex */
3562
 
    in_cs = GEOSGeom_getCoordSeq (g_end);
3563
 
    GEOSCoordSeq_getDimensions (in_cs, &dims);
3564
 
    if (dims == 3)
3565
 
      {
3566
 
          GEOSCoordSeq_getX (in_cs, 0, &x);
3567
 
          GEOSCoordSeq_getY (in_cs, 0, &y);
3568
 
          GEOSCoordSeq_getZ (in_cs, 0, &z);
3569
 
          m = 0.0;
3570
 
      }
3571
 
    else
3572
 
      {
3573
 
          GEOSCoordSeq_getX (in_cs, 0, &x);
3574
 
          GEOSCoordSeq_getY (in_cs, 0, &y);
3575
 
          z = 0.0;
3576
 
          m = 0.0;
3577
 
      }
3578
 
    GEOSGeom_destroy (g_end);
3579
 
    switch (out->DimensionModel)
3580
 
      {
3581
 
      case GAIA_XY_Z:
3582
 
          gaiaSetPointXYZ (out->Coords, points, x, y, z);
3583
 
          break;
3584
 
      case GAIA_XY_M:
3585
 
          gaiaSetPointXYM (out->Coords, points, x, y, 0.0);
3586
 
          break;
3587
 
      case GAIA_XY_Z_M:
3588
 
          gaiaSetPointXYZM (out->Coords, points, x, y, z, 0.0);
3589
 
          break;
3590
 
      default:
3591
 
          gaiaSetPoint (out->Coords, points, x, y);
3592
 
          break;
3593
 
      };
3594
 
    return result;
3595
 
}
3596
 
 
3597
 
static GEOSGeometry *
3598
 
buildGeosPoints (const gaiaGeomCollPtr gaia)
3599
 
{
3600
 
/* converting a GAIA Geometry into a GEOS Geometry of POINTS */
3601
 
    int pts = 0;
3602
 
    unsigned int dims;
3603
 
    int iv;
3604
 
    int ib;
3605
 
    int nItem;
3606
 
    double x;
3607
 
    double y;
3608
 
    double z;
3609
 
    double m;
3610
 
    gaiaPointPtr pt;
3611
 
    gaiaLinestringPtr ln;
3612
 
    gaiaPolygonPtr pg;
3613
 
    gaiaRingPtr rng;
3614
 
    GEOSGeometry *geos;
3615
 
    GEOSGeometry *geos_item;
3616
 
    GEOSGeometry **geos_coll;
3617
 
    GEOSCoordSequence *cs;
3618
 
    if (!gaia)
3619
 
        return NULL;
3620
 
    pt = gaia->FirstPoint;
3621
 
    while (pt)
3622
 
      {
3623
 
          /* counting how many POINTs are there */
3624
 
          pts++;
3625
 
          pt = pt->Next;
3626
 
      }
3627
 
    ln = gaia->FirstLinestring;
3628
 
    while (ln)
3629
 
      {
3630
 
          /* counting how many POINTs are there */
3631
 
          pts += ln->Points;
3632
 
          ln = ln->Next;
3633
 
      }
3634
 
    pg = gaia->FirstPolygon;
3635
 
    while (pg)
3636
 
      {
3637
 
          /* counting how many POINTs are there */
3638
 
          rng = pg->Exterior;
3639
 
          pts += rng->Points - 1;       /* exterior ring */
3640
 
          for (ib = 0; ib < pg->NumInteriors; ib++)
3641
 
            {
3642
 
                /* interior ring */
3643
 
                rng = pg->Interiors + ib;
3644
 
                pts += rng->Points - 1;
3645
 
            }
3646
 
          pg = pg->Next;
3647
 
      }
3648
 
    if (pts == 0)
3649
 
        return NULL;
3650
 
    switch (gaia->DimensionModel)
3651
 
      {
3652
 
      case GAIA_XY_Z:
3653
 
      case GAIA_XY_Z_M:
3654
 
          dims = 3;
3655
 
          break;
3656
 
      default:
3657
 
          dims = 2;
3658
 
          break;
3659
 
      };
3660
 
    nItem = 0;
3661
 
    geos_coll = malloc (sizeof (GEOSGeometry *) * (pts));
3662
 
    pt = gaia->FirstPoint;
3663
 
    while (pt)
3664
 
      {
3665
 
          cs = GEOSCoordSeq_create (1, dims);
3666
 
          switch (pt->DimensionModel)
3667
 
            {
3668
 
            case GAIA_XY_Z:
3669
 
            case GAIA_XY_Z_M:
3670
 
                GEOSCoordSeq_setX (cs, 0, pt->X);
3671
 
                GEOSCoordSeq_setY (cs, 0, pt->Y);
3672
 
                GEOSCoordSeq_setZ (cs, 0, pt->Z);
3673
 
                break;
3674
 
            default:
3675
 
                GEOSCoordSeq_setX (cs, 0, pt->X);
3676
 
                GEOSCoordSeq_setY (cs, 0, pt->Y);
3677
 
                break;
3678
 
            };
3679
 
          geos_item = GEOSGeom_createPoint (cs);
3680
 
          *(geos_coll + nItem++) = geos_item;
3681
 
          pt = pt->Next;
3682
 
      }
3683
 
    ln = gaia->FirstLinestring;
3684
 
    while (ln)
3685
 
      {
3686
 
          for (iv = 0; iv < ln->Points; iv++)
3687
 
            {
3688
 
                z = 0.0;
3689
 
                m = 0.0;
3690
 
                switch (ln->DimensionModel)
3691
 
                  {
3692
 
                  case GAIA_XY_Z:
3693
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3694
 
                      break;
3695
 
                  case GAIA_XY_M:
3696
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3697
 
                      break;
3698
 
                  case GAIA_XY_Z_M:
3699
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3700
 
                      break;
3701
 
                  default:
3702
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
3703
 
                      break;
3704
 
                  };
3705
 
                cs = GEOSCoordSeq_create (1, dims);
3706
 
                if (dims == 3)
3707
 
                  {
3708
 
                      GEOSCoordSeq_setX (cs, 0, x);
3709
 
                      GEOSCoordSeq_setY (cs, 0, y);
3710
 
                      GEOSCoordSeq_setZ (cs, 0, z);
3711
 
                  }
3712
 
                else
3713
 
                  {
3714
 
                      GEOSCoordSeq_setX (cs, 0, x);
3715
 
                      GEOSCoordSeq_setY (cs, 0, y);
3716
 
                  }
3717
 
                geos_item = GEOSGeom_createPoint (cs);
3718
 
                *(geos_coll + nItem++) = geos_item;
3719
 
            }
3720
 
          ln = ln->Next;
3721
 
      }
3722
 
    pg = gaia->FirstPolygon;
3723
 
    while (pg)
3724
 
      {
3725
 
          rng = pg->Exterior;
3726
 
          for (iv = 1; iv < rng->Points; iv++)
3727
 
            {
3728
 
                /* exterior ring */
3729
 
                z = 0.0;
3730
 
                m = 0.0;
3731
 
                switch (rng->DimensionModel)
3732
 
                  {
3733
 
                  case GAIA_XY_Z:
3734
 
                      gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3735
 
                      break;
3736
 
                  case GAIA_XY_M:
3737
 
                      gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3738
 
                      break;
3739
 
                  case GAIA_XY_Z_M:
3740
 
                      gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3741
 
                      break;
3742
 
                  default:
3743
 
                      gaiaGetPoint (rng->Coords, iv, &x, &y);
3744
 
                      break;
3745
 
                  };
3746
 
                cs = GEOSCoordSeq_create (1, dims);
3747
 
                if (dims == 3)
3748
 
                  {
3749
 
                      GEOSCoordSeq_setX (cs, 0, x);
3750
 
                      GEOSCoordSeq_setY (cs, 0, y);
3751
 
                      GEOSCoordSeq_setZ (cs, 0, z);
3752
 
                  }
3753
 
                else
3754
 
                  {
3755
 
                      GEOSCoordSeq_setX (cs, 0, x);
3756
 
                      GEOSCoordSeq_setY (cs, 0, y);
3757
 
                  }
3758
 
                geos_item = GEOSGeom_createPoint (cs);
3759
 
                *(geos_coll + nItem++) = geos_item;
3760
 
            }
3761
 
          for (ib = 0; ib < pg->NumInteriors; ib++)
3762
 
            {
3763
 
                /* interior ring */
3764
 
                rng = pg->Interiors + ib;
3765
 
                for (iv = 1; iv < rng->Points; iv++)
3766
 
                  {
3767
 
                      /* exterior ring */
3768
 
                      z = 0.0;
3769
 
                      m = 0.0;
3770
 
                      switch (rng->DimensionModel)
3771
 
                        {
3772
 
                        case GAIA_XY_Z:
3773
 
                            gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3774
 
                            break;
3775
 
                        case GAIA_XY_M:
3776
 
                            gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3777
 
                            break;
3778
 
                        case GAIA_XY_Z_M:
3779
 
                            gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3780
 
                            break;
3781
 
                        default:
3782
 
                            gaiaGetPoint (rng->Coords, iv, &x, &y);
3783
 
                            break;
3784
 
                        };
3785
 
                      cs = GEOSCoordSeq_create (1, dims);
3786
 
                      if (dims == 3)
3787
 
                        {
3788
 
                            GEOSCoordSeq_setX (cs, 0, x);
3789
 
                            GEOSCoordSeq_setY (cs, 0, y);
3790
 
                            GEOSCoordSeq_setZ (cs, 0, z);
3791
 
                        }
3792
 
                      else
3793
 
                        {
3794
 
                            GEOSCoordSeq_setX (cs, 0, x);
3795
 
                            GEOSCoordSeq_setY (cs, 0, y);
3796
 
                        }
3797
 
                      geos_item = GEOSGeom_createPoint (cs);
3798
 
                      *(geos_coll + nItem++) = geos_item;
3799
 
                  }
3800
 
            }
3801
 
          pg = pg->Next;
3802
 
      }
3803
 
    geos = GEOSGeom_createCollection (GEOS_MULTIPOINT, geos_coll, pts);
3804
 
    free (geos_coll);
3805
 
    GEOSSetSRID (geos, gaia->Srid);
3806
 
    return geos;
3807
 
}
3808
 
 
3809
 
static GEOSGeometry *
3810
 
buildGeosSegments (const gaiaGeomCollPtr gaia)
3811
 
{
3812
 
/* converting a GAIA Geometry into a GEOS Geometry of SEGMENTS */
3813
 
    int segms = 0;
3814
 
    unsigned int dims;
3815
 
    int iv;
3816
 
    int ib;
3817
 
    int nItem;
3818
 
    double x;
3819
 
    double y;
3820
 
    double z;
3821
 
    double m;
3822
 
    double x0;
3823
 
    double y0;
3824
 
    double z0;
3825
 
    gaiaLinestringPtr ln;
3826
 
    gaiaPolygonPtr pg;
3827
 
    gaiaRingPtr rng;
3828
 
    GEOSGeometry *geos;
3829
 
    GEOSGeometry *geos_item;
3830
 
    GEOSGeometry **geos_coll;
3831
 
    GEOSCoordSequence *cs;
3832
 
    if (!gaia)
3833
 
        return NULL;
3834
 
    ln = gaia->FirstLinestring;
3835
 
    while (ln)
3836
 
      {
3837
 
          /* counting how many SEGMENTs are there */
3838
 
          segms += ln->Points - 1;
3839
 
          ln = ln->Next;
3840
 
      }
3841
 
    pg = gaia->FirstPolygon;
3842
 
    while (pg)
3843
 
      {
3844
 
          /* counting how many SEGMENTs are there */
3845
 
          rng = pg->Exterior;
3846
 
          segms += rng->Points - 1;     /* exterior ring */
3847
 
          for (ib = 0; ib < pg->NumInteriors; ib++)
3848
 
            {
3849
 
                /* interior ring */
3850
 
                rng = pg->Interiors + ib;
3851
 
                segms += rng->Points - 1;
3852
 
            }
3853
 
          pg = pg->Next;
3854
 
      }
3855
 
    if (segms == 0)
3856
 
        return NULL;
3857
 
    switch (gaia->DimensionModel)
3858
 
      {
3859
 
      case GAIA_XY_Z:
3860
 
      case GAIA_XY_Z_M:
3861
 
          dims = 3;
3862
 
          break;
3863
 
      default:
3864
 
          dims = 2;
3865
 
          break;
3866
 
      };
3867
 
    nItem = 0;
3868
 
    geos_coll = malloc (sizeof (GEOSGeometry *) * (segms));
3869
 
    ln = gaia->FirstLinestring;
3870
 
    while (ln)
3871
 
      {
3872
 
          for (iv = 0; iv < ln->Points; iv++)
3873
 
            {
3874
 
                z = 0.0;
3875
 
                m = 0.0;
3876
 
                switch (ln->DimensionModel)
3877
 
                  {
3878
 
                  case GAIA_XY_Z:
3879
 
                      gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
3880
 
                      break;
3881
 
                  case GAIA_XY_M:
3882
 
                      gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
3883
 
                      break;
3884
 
                  case GAIA_XY_Z_M:
3885
 
                      gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
3886
 
                      break;
3887
 
                  default:
3888
 
                      gaiaGetPoint (ln->Coords, iv, &x, &y);
3889
 
                      break;
3890
 
                  };
3891
 
                if (iv > 0)
3892
 
                  {
3893
 
                      cs = GEOSCoordSeq_create (2, dims);
3894
 
                      if (dims == 3)
3895
 
                        {
3896
 
                            GEOSCoordSeq_setX (cs, 0, x0);
3897
 
                            GEOSCoordSeq_setY (cs, 0, y0);
3898
 
                            GEOSCoordSeq_setZ (cs, 0, z0);
3899
 
                            GEOSCoordSeq_setX (cs, 1, x);
3900
 
                            GEOSCoordSeq_setY (cs, 1, y);
3901
 
                            GEOSCoordSeq_setZ (cs, 1, z);
3902
 
                        }
3903
 
                      else
3904
 
                        {
3905
 
                            GEOSCoordSeq_setX (cs, 0, x0);
3906
 
                            GEOSCoordSeq_setY (cs, 0, y0);
3907
 
                            GEOSCoordSeq_setX (cs, 1, x);
3908
 
                            GEOSCoordSeq_setY (cs, 1, y);
3909
 
                        }
3910
 
                      geos_item = GEOSGeom_createLineString (cs);
3911
 
                      *(geos_coll + nItem++) = geos_item;
3912
 
                  }
3913
 
                x0 = x;
3914
 
                y0 = y;
3915
 
                z0 = z;
3916
 
            }
3917
 
          ln = ln->Next;
3918
 
      }
3919
 
    pg = gaia->FirstPolygon;
3920
 
    while (pg)
3921
 
      {
3922
 
          rng = pg->Exterior;
3923
 
          for (iv = 0; iv < rng->Points; iv++)
3924
 
            {
3925
 
                /* exterior ring */
3926
 
                z = 0.0;
3927
 
                m = 0.0;
3928
 
                switch (rng->DimensionModel)
3929
 
                  {
3930
 
                  case GAIA_XY_Z:
3931
 
                      gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3932
 
                      break;
3933
 
                  case GAIA_XY_M:
3934
 
                      gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3935
 
                      break;
3936
 
                  case GAIA_XY_Z_M:
3937
 
                      gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3938
 
                      break;
3939
 
                  default:
3940
 
                      gaiaGetPoint (rng->Coords, iv, &x, &y);
3941
 
                      break;
3942
 
                  };
3943
 
                if (iv > 0)
3944
 
                  {
3945
 
                      cs = GEOSCoordSeq_create (2, dims);
3946
 
                      if (dims == 3)
3947
 
                        {
3948
 
                            GEOSCoordSeq_setX (cs, 0, x0);
3949
 
                            GEOSCoordSeq_setY (cs, 0, y0);
3950
 
                            GEOSCoordSeq_setZ (cs, 0, z0);
3951
 
                            GEOSCoordSeq_setX (cs, 1, x);
3952
 
                            GEOSCoordSeq_setY (cs, 1, y);
3953
 
                            GEOSCoordSeq_setZ (cs, 1, z);
3954
 
                        }
3955
 
                      else
3956
 
                        {
3957
 
                            GEOSCoordSeq_setX (cs, 0, x0);
3958
 
                            GEOSCoordSeq_setY (cs, 0, y0);
3959
 
                            GEOSCoordSeq_setX (cs, 1, x);
3960
 
                            GEOSCoordSeq_setY (cs, 1, y);
3961
 
                        }
3962
 
                      geos_item = GEOSGeom_createLineString (cs);
3963
 
                      *(geos_coll + nItem++) = geos_item;
3964
 
                  }
3965
 
                x0 = x;
3966
 
                y0 = y;
3967
 
                z0 = z;
3968
 
            }
3969
 
          for (ib = 0; ib < pg->NumInteriors; ib++)
3970
 
            {
3971
 
                /* interior ring */
3972
 
                rng = pg->Interiors + ib;
3973
 
                for (iv = 0; iv < rng->Points; iv++)
3974
 
                  {
3975
 
                      /* exterior ring */
3976
 
                      z = 0.0;
3977
 
                      m = 0.0;
3978
 
                      switch (rng->DimensionModel)
3979
 
                        {
3980
 
                        case GAIA_XY_Z:
3981
 
                            gaiaGetPointXYZ (rng->Coords, iv, &x, &y, &z);
3982
 
                            break;
3983
 
                        case GAIA_XY_M:
3984
 
                            gaiaGetPointXYM (rng->Coords, iv, &x, &y, &m);
3985
 
                            break;
3986
 
                        case GAIA_XY_Z_M:
3987
 
                            gaiaGetPointXYZM (rng->Coords, iv, &x, &y, &z, &m);
3988
 
                            break;
3989
 
                        default:
3990
 
                            gaiaGetPoint (rng->Coords, iv, &x, &y);
3991
 
                            break;
3992
 
                        };
3993
 
                      if (iv > 0)
3994
 
                        {
3995
 
                            cs = GEOSCoordSeq_create (2, dims);
3996
 
                            if (dims == 3)
3997
 
                              {
3998
 
                                  GEOSCoordSeq_setX (cs, 0, x0);
3999
 
                                  GEOSCoordSeq_setY (cs, 0, y0);
4000
 
                                  GEOSCoordSeq_setZ (cs, 0, z0);
4001
 
                                  GEOSCoordSeq_setX (cs, 1, x);
4002
 
                                  GEOSCoordSeq_setY (cs, 1, y);
4003
 
                                  GEOSCoordSeq_setZ (cs, 1, z);
4004
 
                              }
4005
 
                            else
4006
 
                              {
4007
 
                                  GEOSCoordSeq_setX (cs, 0, x0);
4008
 
                                  GEOSCoordSeq_setY (cs, 0, y0);
4009
 
                                  GEOSCoordSeq_setX (cs, 1, x);
4010
 
                                  GEOSCoordSeq_setY (cs, 1, y);
4011
 
                              }
4012
 
                            geos_item = GEOSGeom_createLineString (cs);
4013
 
                            *(geos_coll + nItem++) = geos_item;
4014
 
                        }
4015
 
                      x0 = x;
4016
 
                      y0 = y;
4017
 
                      z0 = z;
4018
 
                  }
4019
 
            }
4020
 
          pg = pg->Next;
4021
 
      }
4022
 
    geos = GEOSGeom_createCollection (GEOS_MULTILINESTRING, geos_coll, segms);
4023
 
    free (geos_coll);
4024
 
    GEOSSetSRID (geos, gaia->Srid);
4025
 
    return geos;
4026
 
}
4027
 
 
4028
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4029
 
gaiaShortestLine (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
4030
 
{
4031
 
/* attempts to compute the the shortest line between two geometries */
4032
 
    GEOSGeometry *g1_points;
4033
 
    GEOSGeometry *g1_segments;
4034
 
    const GEOSGeometry *g1_item;
4035
 
    GEOSGeometry *g2_points;
4036
 
    GEOSGeometry *g2_segments;
4037
 
    const GEOSGeometry *g2_item;
4038
 
    const GEOSCoordSequence *cs;
4039
 
    GEOSGeometry *g_pt;
4040
 
    gaiaGeomCollPtr result;
4041
 
    gaiaLinestringPtr ln;
4042
 
    int nItems1;
4043
 
    int nItems2;
4044
 
    int it1;
4045
 
    int it2;
4046
 
    unsigned int dims;
4047
 
    double x_ini;
4048
 
    double y_ini;
4049
 
    double z_ini;
4050
 
    double x_fin;
4051
 
    double y_fin;
4052
 
    double z_fin;
4053
 
    double dist;
4054
 
    double min_dist = DBL_MAX;
4055
 
    double projection;
4056
 
    if (!geom1 || !geom2)
4057
 
        return NULL;
4058
 
 
4059
 
    g1_points = buildGeosPoints (geom1);
4060
 
    g1_segments = buildGeosSegments (geom1);
4061
 
    g2_points = buildGeosPoints (geom2);
4062
 
    g2_segments = buildGeosSegments (geom2);
4063
 
 
4064
 
    if (g1_points && g2_points)
4065
 
      {
4066
 
          /* computing distances between POINTs */
4067
 
          nItems1 = GEOSGetNumGeometries (g1_points);
4068
 
          nItems2 = GEOSGetNumGeometries (g2_points);
4069
 
          for (it1 = 0; it1 < nItems1; it1++)
4070
 
            {
4071
 
                g1_item = GEOSGetGeometryN (g1_points, it1);
4072
 
                for (it2 = 0; it2 < nItems2; it2++)
4073
 
                  {
4074
 
                      g2_item = GEOSGetGeometryN (g2_points, it2);
4075
 
                      if (GEOSDistance (g1_item, g2_item, &dist))
4076
 
                        {
4077
 
                            if (dist < min_dist)
4078
 
                              {
4079
 
                                  /* saving min-dist points */
4080
 
                                  min_dist = dist;
4081
 
                                  cs = GEOSGeom_getCoordSeq (g1_item);
4082
 
                                  GEOSCoordSeq_getDimensions (cs, &dims);
4083
 
                                  if (dims == 3)
4084
 
                                    {
4085
 
                                        GEOSCoordSeq_getX (cs, 0, &x_ini);
4086
 
                                        GEOSCoordSeq_getY (cs, 0, &y_ini);
4087
 
                                        GEOSCoordSeq_getZ (cs, 0, &z_ini);
4088
 
                                    }
4089
 
                                  else
4090
 
                                    {
4091
 
                                        GEOSCoordSeq_getX (cs, 0, &x_ini);
4092
 
                                        GEOSCoordSeq_getY (cs, 0, &y_ini);
4093
 
                                        z_ini = 0.0;
4094
 
                                    }
4095
 
                                  cs = GEOSGeom_getCoordSeq (g2_item);
4096
 
                                  GEOSCoordSeq_getDimensions (cs, &dims);
4097
 
                                  if (dims == 3)
4098
 
                                    {
4099
 
                                        GEOSCoordSeq_getX (cs, 0, &x_fin);
4100
 
                                        GEOSCoordSeq_getY (cs, 0, &y_fin);
4101
 
                                        GEOSCoordSeq_getZ (cs, 0, &z_fin);
4102
 
                                    }
4103
 
                                  else
4104
 
                                    {
4105
 
                                        GEOSCoordSeq_getX (cs, 0, &x_fin);
4106
 
                                        GEOSCoordSeq_getY (cs, 0, &y_fin);
4107
 
                                        z_fin = 0.0;
4108
 
                                    }
4109
 
                              }
4110
 
                        }
4111
 
                  }
4112
 
            }
4113
 
      }
4114
 
 
4115
 
    if (g1_points && g2_segments)
4116
 
      {
4117
 
          /* computing distances between POINTs (g1) and SEGMENTs (g2) */
4118
 
          nItems1 = GEOSGetNumGeometries (g1_points);
4119
 
          nItems2 = GEOSGetNumGeometries (g2_segments);
4120
 
          for (it1 = 0; it1 < nItems1; it1++)
4121
 
            {
4122
 
                g1_item = GEOSGetGeometryN (g1_points, it1);
4123
 
                for (it2 = 0; it2 < nItems2; it2++)
4124
 
                  {
4125
 
                      g2_item = GEOSGetGeometryN (g2_segments, it2);
4126
 
                      if (GEOSDistance (g1_item, g2_item, &dist))
4127
 
                        {
4128
 
                            if (dist < min_dist)
4129
 
                              {
4130
 
                                  /* saving min-dist points */
4131
 
                                  projection = GEOSProject (g2_item, g1_item);
4132
 
                                  g_pt = GEOSInterpolate (g2_item, projection);
4133
 
                                  if (g_pt)
4134
 
                                    {
4135
 
                                        min_dist = dist;
4136
 
                                        cs = GEOSGeom_getCoordSeq (g1_item);
4137
 
                                        GEOSCoordSeq_getDimensions (cs, &dims);
4138
 
                                        if (dims == 3)
4139
 
                                          {
4140
 
                                              GEOSCoordSeq_getX (cs, 0, &x_ini);
4141
 
                                              GEOSCoordSeq_getY (cs, 0, &y_ini);
4142
 
                                              GEOSCoordSeq_getZ (cs, 0, &z_ini);
4143
 
                                          }
4144
 
                                        else
4145
 
                                          {
4146
 
                                              GEOSCoordSeq_getX (cs, 0, &x_ini);
4147
 
                                              GEOSCoordSeq_getY (cs, 0, &y_ini);
4148
 
                                              z_ini = 0.0;
4149
 
                                          }
4150
 
                                        cs = GEOSGeom_getCoordSeq (g_pt);
4151
 
                                        GEOSCoordSeq_getDimensions (cs, &dims);
4152
 
                                        if (dims == 3)
4153
 
                                          {
4154
 
                                              GEOSCoordSeq_getX (cs, 0, &x_fin);
4155
 
                                              GEOSCoordSeq_getY (cs, 0, &y_fin);
4156
 
                                              GEOSCoordSeq_getZ (cs, 0, &z_fin);
4157
 
                                          }
4158
 
                                        else
4159
 
                                          {
4160
 
                                              GEOSCoordSeq_getX (cs, 0, &x_fin);
4161
 
                                              GEOSCoordSeq_getY (cs, 0, &y_fin);
4162
 
                                              z_fin = 0.0;
4163
 
                                          }
4164
 
                                        GEOSGeom_destroy (g_pt);
4165
 
                                    }
4166
 
                              }
4167
 
                        }
4168
 
                  }
4169
 
            }
4170
 
      }
4171
 
 
4172
 
    if (g1_segments && g2_points)
4173
 
      {
4174
 
          /* computing distances between SEGMENTs (g1) and POINTs (g2) */
4175
 
          nItems1 = GEOSGetNumGeometries (g1_segments);
4176
 
          nItems2 = GEOSGetNumGeometries (g2_points);
4177
 
          for (it1 = 0; it1 < nItems1; it1++)
4178
 
            {
4179
 
                g1_item = GEOSGetGeometryN (g1_segments, it1);
4180
 
                for (it2 = 0; it2 < nItems2; it2++)
4181
 
                  {
4182
 
                      g2_item = GEOSGetGeometryN (g2_points, it2);
4183
 
                      if (GEOSDistance (g1_item, g2_item, &dist))
4184
 
                        {
4185
 
                            if (dist < min_dist)
4186
 
                              {
4187
 
                                  /* saving min-dist points */
4188
 
                                  projection = GEOSProject (g1_item, g2_item);
4189
 
                                  g_pt = GEOSInterpolate (g1_item, projection);
4190
 
                                  if (g_pt)
4191
 
                                    {
4192
 
                                        min_dist = dist;
4193
 
                                        cs = GEOSGeom_getCoordSeq (g_pt);
4194
 
                                        GEOSCoordSeq_getDimensions (cs, &dims);
4195
 
                                        if (dims == 3)
4196
 
                                          {
4197
 
                                              GEOSCoordSeq_getX (cs, 0, &x_ini);
4198
 
                                              GEOSCoordSeq_getY (cs, 0, &y_ini);
4199
 
                                              GEOSCoordSeq_getZ (cs, 0, &z_ini);
4200
 
                                          }
4201
 
                                        else
4202
 
                                          {
4203
 
                                              GEOSCoordSeq_getX (cs, 0, &x_ini);
4204
 
                                              GEOSCoordSeq_getY (cs, 0, &y_ini);
4205
 
                                              z_ini = 0.0;
4206
 
                                          }
4207
 
                                        cs = GEOSGeom_getCoordSeq (g2_item);
4208
 
                                        GEOSCoordSeq_getDimensions (cs, &dims);
4209
 
                                        if (dims == 3)
4210
 
                                          {
4211
 
                                              GEOSCoordSeq_getX (cs, 0, &x_fin);
4212
 
                                              GEOSCoordSeq_getY (cs, 0, &y_fin);
4213
 
                                              GEOSCoordSeq_getZ (cs, 0, &z_fin);
4214
 
                                          }
4215
 
                                        else
4216
 
                                          {
4217
 
                                              GEOSCoordSeq_getX (cs, 0, &x_fin);
4218
 
                                              GEOSCoordSeq_getY (cs, 0, &y_fin);
4219
 
                                              z_fin = 0.0;
4220
 
                                          }
4221
 
                                        GEOSGeom_destroy (g_pt);
4222
 
                                    }
4223
 
                              }
4224
 
                        }
4225
 
                  }
4226
 
            }
4227
 
      }
4228
 
    if (g1_points)
4229
 
        GEOSGeom_destroy (g1_points);
4230
 
    if (g1_segments)
4231
 
        GEOSGeom_destroy (g1_segments);
4232
 
    if (g2_points)
4233
 
        GEOSGeom_destroy (g2_points);
4234
 
    if (g2_segments)
4235
 
        GEOSGeom_destroy (g2_segments);
4236
 
    if (min_dist == DBL_MAX || min_dist <= 0.0)
4237
 
        return NULL;
4238
 
 
4239
 
/* building the shortest line */
4240
 
    switch (geom1->DimensionModel)
4241
 
      {
4242
 
      case GAIA_XY_Z:
4243
 
          result = gaiaAllocGeomCollXYZ ();
4244
 
          break;
4245
 
      case GAIA_XY_M:
4246
 
          result = gaiaAllocGeomCollXYM ();
4247
 
          break;
4248
 
      case GAIA_XY_Z_M:
4249
 
          result = gaiaAllocGeomCollXYZM ();
4250
 
          break;
4251
 
      default:
4252
 
          result = gaiaAllocGeomColl ();
4253
 
          break;
4254
 
      };
4255
 
    result->Srid = geom1->Srid;
4256
 
    ln = gaiaAddLinestringToGeomColl (result, 2);
4257
 
    switch (ln->DimensionModel)
4258
 
      {
4259
 
      case GAIA_XY_Z:
4260
 
          gaiaSetPointXYZ (ln->Coords, 0, x_ini, y_ini, z_ini);
4261
 
          gaiaSetPointXYZ (ln->Coords, 1, x_fin, y_fin, z_fin);
4262
 
          break;
4263
 
      case GAIA_XY_M:
4264
 
          gaiaSetPointXYM (ln->Coords, 0, x_ini, y_ini, 0.0);
4265
 
          gaiaSetPointXYM (ln->Coords, 1, x_fin, y_fin, 0.0);
4266
 
          break;
4267
 
      case GAIA_XY_Z_M:
4268
 
          gaiaSetPointXYZM (ln->Coords, 0, x_ini, y_ini, z_ini, 0.0);
4269
 
          gaiaSetPointXYZM (ln->Coords, 1, x_fin, y_fin, z_fin, 0.0);
4270
 
          break;
4271
 
      default:
4272
 
          gaiaSetPoint (ln->Coords, 0, x_ini, y_ini);
4273
 
          gaiaSetPoint (ln->Coords, 1, x_fin, y_fin);
4274
 
          break;
4275
 
      };
4276
 
    return result;
4277
 
}
4278
 
 
4279
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4280
 
gaiaSnap (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2, double tolerance)
4281
 
{
4282
 
/* attempts to "snap" geom1 on geom2 using the given tolerance */
4283
 
    GEOSGeometry *g1;
4284
 
    GEOSGeometry *g2;
4285
 
    GEOSGeometry *g3;
4286
 
    gaiaGeomCollPtr result;
4287
 
    if (!geom1 || !geom2)
4288
 
        return NULL;
4289
 
 
4290
 
    g1 = gaiaToGeos (geom1);
4291
 
    g2 = gaiaToGeos (geom2);
4292
 
    g3 = GEOSSnap (g1, g2, tolerance);
4293
 
    GEOSGeom_destroy (g1);
4294
 
    GEOSGeom_destroy (g2);
4295
 
    if (!g3)
4296
 
        return NULL;
4297
 
    if (geom1->DimensionModel == GAIA_XY_Z)
4298
 
        result = gaiaFromGeos_XYZ (g3);
4299
 
    else if (geom1->DimensionModel == GAIA_XY_M)
4300
 
        result = gaiaFromGeos_XYM (g3);
4301
 
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
4302
 
        result = gaiaFromGeos_XYZM (g3);
4303
 
    else
4304
 
        result = gaiaFromGeos_XY (g3);
4305
 
    GEOSGeom_destroy (g3);
4306
 
    if (result == NULL)
4307
 
        return NULL;
4308
 
    result->Srid = geom1->Srid;
4309
 
    return result;
4310
 
}
4311
 
 
4312
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4313
 
gaiaLineMerge (gaiaGeomCollPtr geom)
4314
 
{
4315
 
/* attempts to reassemble lines from a collection of sparse fragments */
4316
 
    GEOSGeometry *g1;
4317
 
    GEOSGeometry *g2;
4318
 
    gaiaGeomCollPtr result;
4319
 
    if (!geom)
4320
 
        return NULL;
4321
 
    if (gaiaIsToxic (geom))
4322
 
        return NULL;
4323
 
 
4324
 
    g1 = gaiaToGeos (geom);
4325
 
    g2 = GEOSLineMerge (g1);
4326
 
    GEOSGeom_destroy (g1);
4327
 
    if (!g2)
4328
 
        return NULL;
4329
 
    if (geom->DimensionModel == GAIA_XY_Z)
4330
 
        result = gaiaFromGeos_XYZ (g2);
4331
 
    else if (geom->DimensionModel == GAIA_XY_M)
4332
 
        result = gaiaFromGeos_XYM (g2);
4333
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
4334
 
        result = gaiaFromGeos_XYZM (g2);
4335
 
    else
4336
 
        result = gaiaFromGeos_XY (g2);
4337
 
    GEOSGeom_destroy (g2);
4338
 
    if (result == NULL)
4339
 
        return NULL;
4340
 
    result->Srid = geom->Srid;
4341
 
    return result;
4342
 
}
4343
 
 
4344
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4345
 
gaiaUnaryUnion (gaiaGeomCollPtr geom)
4346
 
{
4347
 
/* Unary Union (single Collection) */
4348
 
    GEOSGeometry *g1;
4349
 
    GEOSGeometry *g2;
4350
 
    gaiaGeomCollPtr result;
4351
 
    if (!geom)
4352
 
        return NULL;
4353
 
    if (gaiaIsToxic (geom))
4354
 
        return NULL;
4355
 
    g1 = gaiaToGeos (geom);
4356
 
    g2 = GEOSUnaryUnion (g1);
4357
 
    GEOSGeom_destroy (g1);
4358
 
    if (!g2)
4359
 
        return NULL;
4360
 
    if (geom->DimensionModel == GAIA_XY_Z)
4361
 
        result = gaiaFromGeos_XYZ (g2);
4362
 
    else if (geom->DimensionModel == GAIA_XY_M)
4363
 
        result = gaiaFromGeos_XYM (g2);
4364
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
4365
 
        result = gaiaFromGeos_XYZM (g2);
4366
 
    else
4367
 
        result = gaiaFromGeos_XY (g2);
4368
 
    GEOSGeom_destroy (g2);
4369
 
    if (result == NULL)
4370
 
        return NULL;
4371
 
    result->Srid = geom->Srid;
4372
 
    return result;
4373
 
}
4374
 
 
4375
 
static void
4376
 
rotateRingBeforeCut (gaiaLinestringPtr ln, gaiaPointPtr node)
4377
 
{
4378
 
/* rotating a Ring, so to ensure that Start/End points match the node */
4379
 
    int io = 0;
4380
 
    int iv;
4381
 
    int copy = 0;
4382
 
    int base_idx = -1;
4383
 
    double x;
4384
 
    double y;
4385
 
    double z;
4386
 
    double m;
4387
 
    gaiaLinestringPtr new_ln = NULL;
4388
 
 
4389
 
    if (ln->DimensionModel == GAIA_XY_Z)
4390
 
        new_ln = gaiaAllocLinestringXYZ (ln->Points);
4391
 
    else if (ln->DimensionModel == GAIA_XY_M)
4392
 
        new_ln = gaiaAllocLinestringXYM (ln->Points);
4393
 
    else if (ln->DimensionModel == GAIA_XY_Z_M)
4394
 
        new_ln = gaiaAllocLinestringXYZM (ln->Points);
4395
 
    else
4396
 
        new_ln = gaiaAllocLinestring (ln->Points);
4397
 
 
4398
 
/* first pass */
4399
 
    for (iv = 0; iv < ln->Points; iv++)
4400
 
      {
4401
 
          if (ln->DimensionModel == GAIA_XY_Z)
4402
 
            {
4403
 
                gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4404
 
            }
4405
 
          else if (ln->DimensionModel == GAIA_XY_M)
4406
 
            {
4407
 
                gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4408
 
            }
4409
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4410
 
            {
4411
 
                gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4412
 
            }
4413
 
          else
4414
 
            {
4415
 
                gaiaGetPoint (ln->Coords, iv, &x, &y);
4416
 
            }
4417
 
          if (!copy)            /* CAZZO */
4418
 
            {
4419
 
                if (ln->DimensionModel == GAIA_XY_Z
4420
 
                    || ln->DimensionModel == GAIA_XY_Z_M)
4421
 
                  {
4422
 
                      if (node->X == x && node->Y == y && node->Z == z)
4423
 
                        {
4424
 
                            base_idx = iv;
4425
 
                            copy = 1;
4426
 
                        }
4427
 
                  }
4428
 
                else if (node->X == x && node->Y == y)
4429
 
                  {
4430
 
                      base_idx = iv;
4431
 
                      copy = 1;
4432
 
                  }
4433
 
            }
4434
 
          if (copy)
4435
 
            {
4436
 
                /* copying points */
4437
 
                if (ln->DimensionModel == GAIA_XY_Z)
4438
 
                  {
4439
 
                      gaiaSetPointXYZ (new_ln->Coords, io, x, y, z);
4440
 
                  }
4441
 
                else if (ln->DimensionModel == GAIA_XY_M)
4442
 
                  {
4443
 
                      gaiaSetPointXYM (new_ln->Coords, io, x, y, m);
4444
 
                  }
4445
 
                else if (ln->DimensionModel == GAIA_XY_Z_M)
4446
 
                  {
4447
 
                      gaiaSetPointXYZM (new_ln->Coords, io, x, y, z, m);
4448
 
                  }
4449
 
                else
4450
 
                  {
4451
 
                      gaiaSetPoint (new_ln->Coords, io, x, y);
4452
 
                  }
4453
 
                io++;
4454
 
            }
4455
 
      }
4456
 
    if (base_idx <= 0)
4457
 
      {
4458
 
          gaiaFreeLinestring (new_ln);
4459
 
          return;
4460
 
      }
4461
 
 
4462
 
/* second pass */
4463
 
    for (iv = 1; iv <= base_idx; iv++)
4464
 
      {
4465
 
          if (ln->DimensionModel == GAIA_XY_Z)
4466
 
            {
4467
 
                gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4468
 
            }
4469
 
          else if (ln->DimensionModel == GAIA_XY_M)
4470
 
            {
4471
 
                gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4472
 
            }
4473
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4474
 
            {
4475
 
                gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4476
 
            }
4477
 
          else
4478
 
            {
4479
 
                gaiaGetPoint (ln->Coords, iv, &x, &y);
4480
 
            }
4481
 
          if (ln->DimensionModel == GAIA_XY_Z)
4482
 
            {
4483
 
                gaiaSetPointXYZ (new_ln->Coords, io, x, y, z);
4484
 
            }
4485
 
          else if (ln->DimensionModel == GAIA_XY_M)
4486
 
            {
4487
 
                gaiaSetPointXYM (new_ln->Coords, io, x, y, m);
4488
 
            }
4489
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4490
 
            {
4491
 
                gaiaSetPointXYZM (new_ln->Coords, io, x, y, z, m);
4492
 
            }
4493
 
          else
4494
 
            {
4495
 
                gaiaSetPoint (new_ln->Coords, io, x, y);
4496
 
            }
4497
 
          io++;
4498
 
      }
4499
 
 
4500
 
/* copying back */
4501
 
    for (iv = 0; iv < new_ln->Points; iv++)
4502
 
      {
4503
 
          if (ln->DimensionModel == GAIA_XY_Z)
4504
 
            {
4505
 
                gaiaGetPointXYZ (new_ln->Coords, iv, &x, &y, &z);
4506
 
            }
4507
 
          else if (ln->DimensionModel == GAIA_XY_M)
4508
 
            {
4509
 
                gaiaGetPointXYM (new_ln->Coords, iv, &x, &y, &m);
4510
 
            }
4511
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4512
 
            {
4513
 
                gaiaGetPointXYZM (new_ln->Coords, iv, &x, &y, &z, &m);
4514
 
            }
4515
 
          else
4516
 
            {
4517
 
                gaiaGetPoint (new_ln->Coords, iv, &x, &y);
4518
 
            }
4519
 
          if (ln->DimensionModel == GAIA_XY_Z)
4520
 
            {
4521
 
                gaiaSetPointXYZ (ln->Coords, iv, x, y, z);
4522
 
            }
4523
 
          else if (ln->DimensionModel == GAIA_XY_M)
4524
 
            {
4525
 
                gaiaSetPointXYM (ln->Coords, iv, x, y, m);
4526
 
            }
4527
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4528
 
            {
4529
 
                gaiaSetPointXYZM (ln->Coords, iv, x, y, z, m);
4530
 
            }
4531
 
          else
4532
 
            {
4533
 
                gaiaSetPoint (ln->Coords, iv, x, y);
4534
 
            }
4535
 
      }
4536
 
    gaiaFreeLinestring (new_ln);
4537
 
}
4538
 
 
4539
 
static void
4540
 
extractSubLine (gaiaGeomCollPtr result, gaiaLinestringPtr ln, int i_start,
4541
 
                int i_end)
4542
 
{
4543
 
/* extracting s SubLine */
4544
 
    int iv;
4545
 
    int io = 0;
4546
 
    int pts = i_end - i_start + 1;
4547
 
    gaiaLinestringPtr new_ln = NULL;
4548
 
    double x;
4549
 
    double y;
4550
 
    double z;
4551
 
    double m;
4552
 
 
4553
 
    new_ln = gaiaAddLinestringToGeomColl (result, pts);
4554
 
 
4555
 
    for (iv = i_start; iv <= i_end; iv++)
4556
 
      {
4557
 
          if (ln->DimensionModel == GAIA_XY_Z)
4558
 
            {
4559
 
                gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4560
 
            }
4561
 
          else if (ln->DimensionModel == GAIA_XY_M)
4562
 
            {
4563
 
                gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4564
 
            }
4565
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4566
 
            {
4567
 
                gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4568
 
            }
4569
 
          else
4570
 
            {
4571
 
                gaiaGetPoint (ln->Coords, iv, &x, &y);
4572
 
            }
4573
 
          if (ln->DimensionModel == GAIA_XY_Z)
4574
 
            {
4575
 
                gaiaSetPointXYZ (new_ln->Coords, io, x, y, z);
4576
 
            }
4577
 
          else if (ln->DimensionModel == GAIA_XY_M)
4578
 
            {
4579
 
                gaiaSetPointXYM (new_ln->Coords, io, x, y, m);
4580
 
            }
4581
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4582
 
            {
4583
 
                gaiaSetPointXYZM (new_ln->Coords, io, x, y, z, m);
4584
 
            }
4585
 
          else
4586
 
            {
4587
 
                gaiaSetPoint (new_ln->Coords, io, x, y);
4588
 
            }
4589
 
          io++;
4590
 
      }
4591
 
}
4592
 
 
4593
 
static void
4594
 
cutLineAtNodes (gaiaLinestringPtr ln, gaiaPointPtr pt_base,
4595
 
                gaiaGeomCollPtr result)
4596
 
{
4597
 
/* attempts to cut a single Line accordingly to given nodes */
4598
 
    int closed = 0;
4599
 
    int match = 0;
4600
 
    int iv;
4601
 
    int i_start;
4602
 
    double x;
4603
 
    double y;
4604
 
    double z;
4605
 
    double m;
4606
 
    gaiaPointPtr pt;
4607
 
    gaiaPointPtr node = NULL;
4608
 
 
4609
 
    if (gaiaIsClosed (ln))
4610
 
        closed = 1;
4611
 
/* pre-check */
4612
 
    for (iv = 0; iv < ln->Points; iv++)
4613
 
      {
4614
 
          if (ln->DimensionModel == GAIA_XY_Z)
4615
 
            {
4616
 
                gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4617
 
            }
4618
 
          else if (ln->DimensionModel == GAIA_XY_M)
4619
 
            {
4620
 
                gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4621
 
            }
4622
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4623
 
            {
4624
 
                gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4625
 
            }
4626
 
          else
4627
 
            {
4628
 
                gaiaGetPoint (ln->Coords, iv, &x, &y);
4629
 
            }
4630
 
          pt = pt_base;
4631
 
          while (pt)
4632
 
            {
4633
 
                if (ln->DimensionModel == GAIA_XY_Z
4634
 
                    || ln->DimensionModel == GAIA_XY_Z_M)
4635
 
                  {
4636
 
                      if (pt->X == x && pt->Y == y && pt->Z == z)
4637
 
                        {
4638
 
                            node = pt;
4639
 
                            match++;
4640
 
                        }
4641
 
                  }
4642
 
                else if (pt->X == x && pt->Y == y)
4643
 
                  {
4644
 
                      node = pt;
4645
 
                      match++;
4646
 
                  }
4647
 
                pt = pt->Next;
4648
 
            }
4649
 
      }
4650
 
 
4651
 
    if (closed && node)
4652
 
        rotateRingBeforeCut (ln, node);
4653
 
 
4654
 
    i_start = 0;
4655
 
    for (iv = 1; iv < ln->Points - 1; iv++)
4656
 
      {
4657
 
          /* identifying sub-linestrings */
4658
 
          if (ln->DimensionModel == GAIA_XY_Z)
4659
 
            {
4660
 
                gaiaGetPointXYZ (ln->Coords, iv, &x, &y, &z);
4661
 
            }
4662
 
          else if (ln->DimensionModel == GAIA_XY_M)
4663
 
            {
4664
 
                gaiaGetPointXYM (ln->Coords, iv, &x, &y, &m);
4665
 
            }
4666
 
          else if (ln->DimensionModel == GAIA_XY_Z_M)
4667
 
            {
4668
 
                gaiaGetPointXYZM (ln->Coords, iv, &x, &y, &z, &m);
4669
 
            }
4670
 
          else
4671
 
            {
4672
 
                gaiaGetPoint (ln->Coords, iv, &x, &y);
4673
 
            }
4674
 
          match = 0;
4675
 
          pt = pt_base;
4676
 
          while (pt)
4677
 
            {
4678
 
                if (ln->DimensionModel == GAIA_XY_Z
4679
 
                    || ln->DimensionModel == GAIA_XY_Z_M)
4680
 
                  {
4681
 
                      if (pt->X == x && pt->Y == y && pt->Z == z)
4682
 
                        {
4683
 
                            match = 1;
4684
 
                            break;
4685
 
                        }
4686
 
                  }
4687
 
                else if (pt->X == x && pt->Y == y)
4688
 
                  {
4689
 
                      match = 1;
4690
 
                      break;
4691
 
                  }
4692
 
                pt = pt->Next;
4693
 
            }
4694
 
          if (match)
4695
 
            {
4696
 
                /* cutting the line */
4697
 
                extractSubLine (result, ln, i_start, iv);
4698
 
                i_start = iv;
4699
 
            }
4700
 
      }
4701
 
    if (i_start != 0 && i_start != ln->Points - 1)
4702
 
      {
4703
 
          /* extracting the last SubLine */
4704
 
          extractSubLine (result, ln, i_start, ln->Points - 1);
4705
 
      }
4706
 
    else
4707
 
      {
4708
 
          /* cloning the untouched Line */
4709
 
          extractSubLine (result, ln, 0, ln->Points - 1);
4710
 
      }
4711
 
}
4712
 
 
4713
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4714
 
gaiaLinesCutAtNodes (gaiaGeomCollPtr geom1, gaiaGeomCollPtr geom2)
4715
 
{
4716
 
/* attempts to cut lines accordingly to nodes */
4717
 
    int pts1 = 0;
4718
 
    int lns1 = 0;
4719
 
    int pgs1 = 0;
4720
 
    int pts2 = 0;
4721
 
    int lns2 = 0;
4722
 
    int pgs2 = 0;
4723
 
    gaiaPointPtr pt;
4724
 
    gaiaLinestringPtr ln;
4725
 
    gaiaPolygonPtr pg;
4726
 
    gaiaGeomCollPtr result = NULL;
4727
 
 
4728
 
    if (!geom1)
4729
 
        return NULL;
4730
 
    if (!geom2)
4731
 
        return NULL;
4732
 
 
4733
 
/* both Geometryes should have identical Dimensions */
4734
 
    if (geom1->DimensionModel != geom2->DimensionModel)
4735
 
        return NULL;
4736
 
 
4737
 
    pt = geom1->FirstPoint;
4738
 
    while (pt)
4739
 
      {
4740
 
          pts1++;
4741
 
          pt = pt->Next;
4742
 
      }
4743
 
    ln = geom1->FirstLinestring;
4744
 
    while (ln)
4745
 
      {
4746
 
          lns1++;
4747
 
          ln = ln->Next;
4748
 
      }
4749
 
    pg = geom1->FirstPolygon;
4750
 
    while (pg)
4751
 
      {
4752
 
          pgs1++;
4753
 
          pg = pg->Next;
4754
 
      }
4755
 
    pt = geom2->FirstPoint;
4756
 
    while (pt)
4757
 
      {
4758
 
          pts2++;
4759
 
          pt = pt->Next;
4760
 
      }
4761
 
    ln = geom2->FirstLinestring;
4762
 
    while (ln)
4763
 
      {
4764
 
          lns2++;
4765
 
          ln = ln->Next;
4766
 
      }
4767
 
    pg = geom2->FirstPolygon;
4768
 
    while (pg)
4769
 
      {
4770
 
          pgs2++;
4771
 
          pg = pg->Next;
4772
 
      }
4773
 
 
4774
 
/* the first Geometry is expected to contain one or more Linestring(s) */
4775
 
    if (pts1 == 0 && lns1 > 0 && pgs1 == 0)
4776
 
        ;
4777
 
    else
4778
 
        return NULL;
4779
 
/* the second Geometry is expected to contain one or more Point(s) */
4780
 
    if (pts2 > 0 && lns2 == 0 && pgs2 == 0)
4781
 
        ;
4782
 
    else
4783
 
        return NULL;
4784
 
 
4785
 
/* attempting to cut Lines accordingly to Nodes */
4786
 
    if (geom1->DimensionModel == GAIA_XY_Z)
4787
 
        result = gaiaAllocGeomCollXYZ ();
4788
 
    else if (geom1->DimensionModel == GAIA_XY_M)
4789
 
        result = gaiaAllocGeomCollXYM ();
4790
 
    else if (geom1->DimensionModel == GAIA_XY_Z_M)
4791
 
        result = gaiaAllocGeomCollXYZM ();
4792
 
    else
4793
 
        result = gaiaAllocGeomColl ();
4794
 
    ln = geom1->FirstLinestring;
4795
 
    while (ln)
4796
 
      {
4797
 
          cutLineAtNodes (ln, geom2->FirstPoint, result);
4798
 
          ln = ln->Next;
4799
 
      }
4800
 
    if (result->FirstLinestring == NULL)
4801
 
      {
4802
 
          gaiaFreeGeomColl (result);
4803
 
          return NULL;
4804
 
      }
4805
 
    result->Srid = geom1->Srid;
4806
 
    return result;
4807
 
}
4808
 
 
4809
 
#endif /* end GEOS advanced features */
4810
 
 
4811
 
#ifdef GEOS_TRUNK               /* GEOS experimental features */
4812
 
 
4813
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4814
 
gaiaDelaunayTriangulation (gaiaGeomCollPtr geom, double tolerance,
4815
 
                           int only_edges)
4816
 
{
4817
 
/* Delaunay Triangulation */
4818
 
    GEOSGeometry *g1;
4819
 
    GEOSGeometry *g2;
4820
 
    gaiaGeomCollPtr result;
4821
 
    if (!geom)
4822
 
        return NULL;
4823
 
    g1 = gaiaToGeos (geom);
4824
 
    g2 = GEOSDelaunayTriangulation (g1, tolerance, only_edges);
4825
 
    GEOSGeom_destroy (g1);
4826
 
    if (!g2)
4827
 
        return NULL;
4828
 
    if (geom->DimensionModel == GAIA_XY_Z)
4829
 
        result = gaiaFromGeos_XYZ (g2);
4830
 
    else if (geom->DimensionModel == GAIA_XY_M)
4831
 
        result = gaiaFromGeos_XYM (g2);
4832
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
4833
 
        result = gaiaFromGeos_XYZM (g2);
4834
 
    else
4835
 
        result = gaiaFromGeos_XY (g2);
4836
 
    GEOSGeom_destroy (g2);
4837
 
    if (result == NULL)
4838
 
        return NULL;
4839
 
    result->Srid = geom->Srid;
4840
 
    if (only_edges)
4841
 
        result->DeclaredType = GAIA_MULTILINESTRING;
4842
 
    else
4843
 
        result->DeclaredType = GAIA_MULTIPOLYGON;
4844
 
    return result;
4845
 
}
4846
 
 
4847
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4848
 
gaiaVoronojDiagram (gaiaGeomCollPtr geom, double extra_frame_size,
4849
 
                    double tolerance, int only_edges)
4850
 
{
4851
 
/* Voronoj Diagram */
4852
 
    GEOSGeometry *g1;
4853
 
    GEOSGeometry *g2;
4854
 
    gaiaGeomCollPtr result;
4855
 
    gaiaPolygonPtr pg;
4856
 
    int pgs = 0;
4857
 
    int errs = 0;
4858
 
    void *voronoj;
4859
 
    if (!geom)
4860
 
        return NULL;
4861
 
    g1 = gaiaToGeos (geom);
4862
 
    g2 = GEOSDelaunayTriangulation (g1, tolerance, 0);
4863
 
    GEOSGeom_destroy (g1);
4864
 
    if (!g2)
4865
 
        return NULL;
4866
 
    if (geom->DimensionModel == GAIA_XY_Z)
4867
 
        result = gaiaFromGeos_XYZ (g2);
4868
 
    else if (geom->DimensionModel == GAIA_XY_M)
4869
 
        result = gaiaFromGeos_XYM (g2);
4870
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
4871
 
        result = gaiaFromGeos_XYZM (g2);
4872
 
    else
4873
 
        result = gaiaFromGeos_XY (g2);
4874
 
    GEOSGeom_destroy (g2);
4875
 
    if (result == NULL)
4876
 
        return NULL;
4877
 
    pg = result->FirstPolygon;
4878
 
    while (pg)
4879
 
      {
4880
 
          /* counting how many triangles are in Delaunay */
4881
 
          if (delaunay_triangle_check (pg))
4882
 
              pgs++;
4883
 
          else
4884
 
              errs++;
4885
 
          pg = pg->Next;
4886
 
      }
4887
 
    if (pgs == 0 || errs)
4888
 
      {
4889
 
          gaiaFreeGeomColl (result);
4890
 
          return NULL;
4891
 
      }
4892
 
 
4893
 
/* building the Voronoj Diagram from Delaunay */
4894
 
    voronoj = voronoj_build (pgs, result->FirstPolygon, extra_frame_size);
4895
 
    gaiaFreeGeomColl (result);
4896
 
 
4897
 
/* creating the Geometry representing Voronoj */
4898
 
    if (geom->DimensionModel == GAIA_XY_Z)
4899
 
        result = gaiaAllocGeomCollXYZ ();
4900
 
    else if (geom->DimensionModel == GAIA_XY_M)
4901
 
        result = gaiaAllocGeomCollXYM ();
4902
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
4903
 
        result = gaiaAllocGeomCollXYZM ();
4904
 
    else
4905
 
        result = gaiaAllocGeomColl ();
4906
 
    result = voronoj_export (voronoj, result, only_edges);
4907
 
    voronoj_free (voronoj);
4908
 
 
4909
 
    result->Srid = geom->Srid;
4910
 
    if (only_edges)
4911
 
        result->DeclaredType = GAIA_MULTILINESTRING;
4912
 
    else
4913
 
        result->DeclaredType = GAIA_MULTIPOLYGON;
4914
 
    return result;
4915
 
}
4916
 
 
4917
 
GAIAGEO_DECLARE gaiaGeomCollPtr
4918
 
gaiaConcaveHull (gaiaGeomCollPtr geom, double factor, double tolerance,
4919
 
                 int allow_holes)
4920
 
{
4921
 
/* Concave Hull */
4922
 
    GEOSGeometry *g1;
4923
 
    GEOSGeometry *g2;
4924
 
    gaiaGeomCollPtr result;
4925
 
    gaiaGeomCollPtr concave_hull;
4926
 
    gaiaPolygonPtr pg;
4927
 
    int pgs = 0;
4928
 
    int errs = 0;
4929
 
    if (!geom)
4930
 
        return NULL;
4931
 
    g1 = gaiaToGeos (geom);
4932
 
    g2 = GEOSDelaunayTriangulation (g1, tolerance, 0);
4933
 
    GEOSGeom_destroy (g1);
4934
 
    if (!g2)
4935
 
        return NULL;
4936
 
    if (geom->DimensionModel == GAIA_XY_Z)
4937
 
        result = gaiaFromGeos_XYZ (g2);
4938
 
    else if (geom->DimensionModel == GAIA_XY_M)
4939
 
        result = gaiaFromGeos_XYM (g2);
4940
 
    else if (geom->DimensionModel == GAIA_XY_Z_M)
4941
 
        result = gaiaFromGeos_XYZM (g2);
4942
 
    else
4943
 
        result = gaiaFromGeos_XY (g2);
4944
 
    GEOSGeom_destroy (g2);
4945
 
    if (result == NULL)
4946
 
        return NULL;
4947
 
    pg = result->FirstPolygon;
4948
 
    while (pg)
4949
 
      {
4950
 
          /* counting how many triangles are in Delaunay */
4951
 
          if (delaunay_triangle_check (pg))
4952
 
              pgs++;
4953
 
          else
4954
 
              errs++;
4955
 
          pg = pg->Next;
4956
 
      }
4957
 
    if (pgs == 0 || errs)
4958
 
      {
4959
 
          gaiaFreeGeomColl (result);
4960
 
          return NULL;
4961
 
      }
4962
 
 
4963
 
/* building the Concave Hull from Delaunay */
4964
 
    concave_hull =
4965
 
        concave_hull_build (result->FirstPolygon, geom->DimensionModel, factor,
4966
 
                            allow_holes);
4967
 
    gaiaFreeGeomColl (result);
4968
 
    if (!concave_hull)
4969
 
        return NULL;
4970
 
    result = concave_hull;
4971
 
 
4972
 
    result->Srid = geom->Srid;
4973
 
    return result;
4974
 
}
4975
 
 
4976
 
#endif /* end GEOS experimental features */
4977
 
 
4978
3988
#endif /* end including GEOS */