@@ -45,7 +45,9 @@ namespace pal
45
45
46
46
47
47
PointSet::PointSet ()
48
- : holeOf( NULL )
48
+ : mGeos ( 0 )
49
+ , mOwnsGeom ( false )
50
+ , holeOf( NULL )
49
51
, parent( NULL )
50
52
, xmin( DBL_MAX )
51
53
, xmax( -DBL_MAX )
@@ -60,7 +62,9 @@ namespace pal
60
62
}
61
63
62
64
PointSet::PointSet ( int nbPoints, double *x, double *y )
63
- : cHullSize( 0 )
65
+ : mGeos ( 0 )
66
+ , mOwnsGeom ( false )
67
+ , cHullSize( 0 )
64
68
, holeOf( NULL )
65
69
, parent( NULL )
66
70
, xmin( DBL_MAX )
@@ -83,8 +87,12 @@ namespace pal
83
87
}
84
88
85
89
PointSet::PointSet ( double aX, double aY )
86
- : xmin( aX ), xmax( aY )
87
- , ymin( aX ), ymax( aY )
90
+ : mGeos ( 0 )
91
+ , mOwnsGeom ( false )
92
+ , xmin( aX )
93
+ , xmax( aY )
94
+ , ymin( aX )
95
+ , ymax( aY )
88
96
{
89
97
nbPoints = cHullSize = 1 ;
90
98
x = new double [1 ];
@@ -100,7 +108,9 @@ namespace pal
100
108
}
101
109
102
110
PointSet::PointSet ( PointSet &ps )
103
- : parent( 0 )
111
+ : mGeos ( 0 )
112
+ , mOwnsGeom ( false )
113
+ , parent( 0 )
104
114
, xmin( DBL_MAX )
105
115
, xmax( -DBL_MAX )
106
116
, ymin( DBL_MAX )
@@ -138,8 +148,27 @@ namespace pal
138
148
holeOf = ps.holeOf ;
139
149
}
140
150
151
+ void PointSet::createGeosGeom ()
152
+ {
153
+ GEOSContextHandle_t geosctxt = geosContext ();
154
+ GEOSCoordSequence *coord = GEOSCoordSeq_create_r ( geosctxt, nbPoints, 2 );
155
+ for ( int i = 0 ; i < nbPoints; ++i )
156
+ {
157
+ GEOSCoordSeq_setX_r ( geosctxt, coord, i, x[i] );
158
+ GEOSCoordSeq_setY_r ( geosctxt, coord, i, y[i] );
159
+ }
160
+ mGeos = GEOSGeom_createPolygon_r ( geosctxt, GEOSGeom_createLinearRing_r ( geosctxt, coord ), 0 , 0 );
161
+ mOwnsGeom = true ;
162
+ }
163
+
141
164
PointSet::~PointSet ()
142
165
{
166
+ if ( mOwnsGeom )
167
+ {
168
+ GEOSGeom_destroy_r ( geosContext (), mGeos );
169
+ mGeos = NULL ;
170
+ }
171
+
143
172
deleteCoords ();
144
173
145
174
if ( cHull )
@@ -156,9 +185,6 @@ namespace pal
156
185
y = NULL ;
157
186
}
158
187
159
-
160
-
161
-
162
188
PointSet* PointSet::extractShape ( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
163
189
{
164
190
@@ -842,65 +868,38 @@ namespace pal
842
868
}
843
869
844
870
845
-
846
871
void PointSet::getCentroid ( double &px, double &py, bool forceInside )
847
872
{
848
- // for explanation see this page:
849
- // http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/
873
+ if ( ! mGeos )
874
+ createGeosGeom ();
850
875
851
- int i, j;
852
- double cx = 0 , cy = 0 , A = 0 , tmp, sumx = 0 , sumy = 0 ;
853
- for ( i = 0 ; i < nbPoints; i++ )
854
- {
855
- j = i + 1 ; if ( j == nbPoints ) j = 0 ;
856
- tmp = (( x[i] - x[0 ] ) * ( y[j] - y[0 ] ) - ( x[j] - x[0 ] ) * ( y[i] - y[0 ] ) );
857
- cx += ( x[i] + x[j] - 2 * x[0 ] ) * tmp;
858
- cy += ( y[i] + y[j] - 2 * y[0 ] ) * tmp;
859
- A += tmp;
860
- sumx += x[i];
861
- sumy += y[i];
862
- }
876
+ if ( !mGeos )
877
+ return ;
863
878
864
- if ( A == 0 )
879
+ GEOSContextHandle_t geosctxt = geosContext ();
880
+ GEOSGeometry *centroidGeom = GEOSGetCentroid_r ( geosctxt, mGeos );
881
+ if ( centroidGeom )
865
882
{
866
- px = sumx / nbPoints;
867
- py = sumy / nbPoints;
868
- }
869
- else
870
- {
871
- px = cx / ( 3 * A ) + x[0 ];
872
- py = cy / ( 3 * A ) + y[0 ];
883
+ const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r ( geosctxt, centroidGeom );
884
+ GEOSCoordSeq_getX_r ( geosctxt, coordSeq, 0 , &px );
885
+ GEOSCoordSeq_getY_r ( geosctxt, coordSeq, 0 , &py );
873
886
}
874
887
875
888
// check if centroid inside in polygon
876
889
if ( forceInside && !isPointInPolygon ( nbPoints, x, y, px, py ) )
877
890
{
878
- GEOSContextHandle_t geosctxt = geosContext ();
879
- GEOSCoordSequence *coord = GEOSCoordSeq_create_r ( geosctxt, nbPoints, 2 );
880
-
881
- for ( int i = 0 ; i < nbPoints; ++i )
882
- {
883
- GEOSCoordSeq_setX_r ( geosctxt, coord, i, x[i] );
884
- GEOSCoordSeq_setY_r ( geosctxt, coord, i, y[i] );
885
- }
886
-
887
- GEOSGeometry *geom = GEOSGeom_createPolygon_r ( geosctxt, GEOSGeom_createLinearRing_r ( geosctxt, coord ), 0 , 0 );
891
+ GEOSGeometry *pointGeom = GEOSPointOnSurface_r ( geosctxt, mGeos );
888
892
889
- if ( geom )
893
+ if ( pointGeom )
890
894
{
891
- GEOSGeometry *pointGeom = GEOSPointOnSurface_r ( geosctxt, geom );
892
-
893
- if ( pointGeom )
894
- {
895
- const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r ( geosctxt, pointGeom );
896
- GEOSCoordSeq_getX_r ( geosctxt, coordSeq, 0 , &px );
897
- GEOSCoordSeq_getY_r ( geosctxt, coordSeq, 0 , &py );
898
-
899
- GEOSGeom_destroy_r ( geosctxt, pointGeom );
900
- }
901
- GEOSGeom_destroy_r ( geosctxt, geom );
895
+ const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r ( geosctxt, pointGeom );
896
+ GEOSCoordSeq_getX_r ( geosctxt, coordSeq, 0 , &px );
897
+ GEOSCoordSeq_getY_r ( geosctxt, coordSeq, 0 , &py );
898
+ GEOSGeom_destroy_r ( geosctxt, pointGeom );
902
899
}
903
900
}
901
+
902
+ GEOSGeom_destroy_r ( geosctxt, centroidGeom );
904
903
}
905
904
906
905
} // end namespace
0 commit comments