31
31
import org .elasticsearch .geometry .Point ;
32
32
import org .elasticsearch .geometry .Polygon ;
33
33
import org .elasticsearch .geometry .Rectangle ;
34
+ import org .elasticsearch .search .aggregations .metrics .CompensatedSum ;
35
+
36
+ import static org .elasticsearch .common .geo .DimensionalShapeType .LINE ;
37
+ import static org .elasticsearch .common .geo .DimensionalShapeType .POINT ;
38
+ import static org .elasticsearch .common .geo .DimensionalShapeType .POLYGON ;
34
39
35
40
/**
36
41
* This class keeps a running Kahan-sum of coordinates
37
42
* that are to be averaged in {@link TriangleTreeWriter} for use
38
43
* as the centroid of a shape.
39
44
*/
40
45
public class CentroidCalculator {
41
- private double compX ;
42
- private double compY ;
43
- private double sumX ;
44
- private double sumY ;
45
- private double sumWeight ;
46
+ CompensatedSum compSumX ;
47
+ CompensatedSum compSumY ;
48
+ CompensatedSum compSumWeight ;
49
+
50
+ private CentroidCalculatorVisitor visitor ;
46
51
private DimensionalShapeType dimensionalShapeType ;
47
52
48
53
public CentroidCalculator (Geometry geometry ) {
49
- this .sumX = 0.0 ;
50
- this .compX = 0.0 ;
51
- this .sumY = 0.0 ;
52
- this .compY = 0.0 ;
53
- this .sumWeight = 0.0 ;
54
- CentroidCalculatorVisitor visitor = new CentroidCalculatorVisitor (this );
54
+ this .compSumX = new CompensatedSum (0 , 0 );
55
+ this .compSumY = new CompensatedSum (0 , 0 );
56
+ this .compSumWeight = new CompensatedSum (0 , 0 );
57
+ this .dimensionalShapeType = null ;
58
+ this .visitor = new CentroidCalculatorVisitor (this );
55
59
geometry .visit (visitor );
56
- this .dimensionalShapeType = DimensionalShapeType .forGeometry (geometry );
60
+ this .dimensionalShapeType = visitor .calculator .dimensionalShapeType ;
61
+ }
62
+
63
+ private CompensatedSum getCompSumX () {
64
+ return compSumX ;
65
+ }
66
+
67
+ private CompensatedSum getCompSumY () {
68
+ return compSumY ;
69
+ }
70
+
71
+ private CompensatedSum getCompSumWeight () {
72
+ return compSumWeight ;
57
73
}
58
74
59
75
/**
@@ -63,18 +79,19 @@ public CentroidCalculator(Geometry geometry) {
63
79
* @param y the y-coordinate of the point
64
80
* @param weight the associated weight of the coordinate
65
81
*/
66
- private void addCoordinate (double x , double y , double weight ) {
67
- double correctedX = weight * x - compX ;
68
- double newSumX = sumX + correctedX ;
69
- compX = (newSumX - sumX ) - correctedX ;
70
- sumX = newSumX ;
71
-
72
- double correctedY = weight * y - compY ;
73
- double newSumY = sumY + correctedY ;
74
- compY = (newSumY - sumY ) - correctedY ;
75
- sumY = newSumY ;
76
-
77
- sumWeight += weight ;
82
+ private void addCoordinate (double x , double y , double weight , DimensionalShapeType dimensionalShapeType ) {
83
+ if (this .dimensionalShapeType == null || this .dimensionalShapeType == dimensionalShapeType ) {
84
+ compSumX .add (x * weight );
85
+ compSumY .add (y * weight );
86
+ compSumWeight .add (weight );
87
+ this .dimensionalShapeType = dimensionalShapeType ;
88
+ } else if (dimensionalShapeType .compareTo (this .dimensionalShapeType ) > 0 ) {
89
+ // reset counters
90
+ compSumX .reset (x * weight , 0 );
91
+ compSumY .reset (y * weight , 0 );
92
+ compSumWeight .reset (weight , 0 );
93
+ this .dimensionalShapeType = dimensionalShapeType ;
94
+ }
78
95
}
79
96
80
97
/**
@@ -86,16 +103,17 @@ private void addCoordinate(double x, double y, double weight) {
86
103
* @param otherCalculator the other centroid calculator to add from
87
104
*/
88
105
public void addFrom (CentroidCalculator otherCalculator ) {
89
- int compared = DimensionalShapeType . COMPARATOR . compare ( dimensionalShapeType , otherCalculator .dimensionalShapeType );
106
+ int compared = dimensionalShapeType . compareTo ( otherCalculator .dimensionalShapeType );
90
107
if (compared < 0 ) {
91
- sumWeight = otherCalculator .sumWeight ;
92
108
dimensionalShapeType = otherCalculator .dimensionalShapeType ;
93
- sumX = otherCalculator .sumX ;
94
- sumY = otherCalculator .sumY ;
95
- compX = otherCalculator .compX ;
96
- compY = otherCalculator . compY ;
109
+ this . compSumX = otherCalculator .compSumX ;
110
+ this . compSumY = otherCalculator .compSumY ;
111
+ this . compSumWeight = otherCalculator .compSumWeight ;
112
+
97
113
} else if (compared == 0 ) {
98
- addCoordinate (otherCalculator .sumX , otherCalculator .sumY , otherCalculator .sumWeight );
114
+ this .compSumX .add (otherCalculator .getCompSumX ().value ());
115
+ this .compSumY .add (otherCalculator .getCompSumY ().value ());
116
+ this .compSumWeight .add (otherCalculator .getCompSumWeight ().value ());
99
117
} // else (compared > 0) do not modify centroid calculation since otherCalculator is of lower dimension than this calculator
100
118
}
101
119
@@ -104,22 +122,22 @@ public void addFrom(CentroidCalculator otherCalculator) {
104
122
*/
105
123
public double getX () {
106
124
// normalization required due to floating point precision errors
107
- return GeoUtils .normalizeLon (sumX / sumWeight );
125
+ return GeoUtils .normalizeLon (getCompSumX (). value () / getCompSumWeight (). value () );
108
126
}
109
127
110
128
/**
111
129
* @return the y-coordinate centroid
112
130
*/
113
131
public double getY () {
114
132
// normalization required due to floating point precision errors
115
- return GeoUtils .normalizeLat (sumY / sumWeight );
133
+ return GeoUtils .normalizeLat (getCompSumY (). value () / getCompSumWeight (). value () );
116
134
}
117
135
118
136
/**
119
137
* @return the sum of all the weighted coordinates summed in the calculator
120
138
*/
121
139
public double sumWeight () {
122
- return sumWeight ;
140
+ return getCompSumWeight (). value () ;
123
141
}
124
142
125
143
/**
@@ -152,62 +170,34 @@ public Void visit(GeometryCollection<?> collection) {
152
170
153
171
@ Override
154
172
public Void visit (Line line ) {
155
- // a line's centroid is calculated by summing the center of each
156
- // line segment weighted by the line segment's length in degrees
157
- for (int i = 0 ; i < line .length () - 1 ; i ++) {
158
- double diffX = line .getX (i ) - line .getX (i + 1 );
159
- double diffY = line .getY (i ) - line .getY (i + 1 );
160
- double x = (line .getX (i ) + line .getX (i + 1 )) / 2 ;
161
- double y = (line .getY (i ) + line .getY (i + 1 )) / 2 ;
162
- calculator .addCoordinate (x , y , Math .sqrt (diffX * diffX + diffY * diffY ));
173
+ if (calculator .dimensionalShapeType != POLYGON ) {
174
+ visitLine (line .length (), line ::getX , line ::getY );
163
175
}
164
176
return null ;
165
177
}
178
+
166
179
@ Override
167
180
public Void visit (LinearRing ring ) {
168
181
throw new IllegalArgumentException ("invalid shape type found [LinearRing] while calculating centroid" );
169
182
}
170
183
171
- private Void visit (LinearRing ring , boolean isHole ) {
172
- // implementation of calculation defined in
173
- // https://www.seas.upenn.edu/~sys502/extra_materials/Polygon%20Area%20and%20Centroid.pdf
174
- //
175
- // centroid of a ring is a weighted coordinate based on the ring's area.
176
- // the sign of the area is positive for the outer-shell of a polygon and negative for the holes
177
-
178
- int sign = isHole ? -1 : 1 ;
179
- double totalRingArea = 0.0 ;
180
- for (int i = 0 ; i < ring .length () - 1 ; i ++) {
181
- totalRingArea += (ring .getX (i ) * ring .getY (i + 1 )) - (ring .getX (i + 1 ) * ring .getY (i ));
182
- }
183
- totalRingArea = totalRingArea / 2 ;
184
-
185
- double sumX = 0.0 ;
186
- double sumY = 0.0 ;
187
- for (int i = 0 ; i < ring .length () - 1 ; i ++) {
188
- double twiceArea = (ring .getX (i ) * ring .getY (i + 1 )) - (ring .getX (i + 1 ) * ring .getY (i ));
189
- sumX += twiceArea * (ring .getX (i ) + ring .getX (i + 1 ));
190
- sumY += twiceArea * (ring .getY (i ) + ring .getY (i + 1 ));
191
- }
192
- double cX = sumX / (6 * totalRingArea );
193
- double cY = sumY / (6 * totalRingArea );
194
- calculator .addCoordinate (cX , cY , sign * Math .abs (totalRingArea ));
195
-
196
- return null ;
197
- }
198
184
199
185
@ Override
200
186
public Void visit (MultiLine multiLine ) {
201
- for (Line line : multiLine ) {
202
- visit (line );
187
+ if (calculator .getDimensionalShapeType () != POLYGON ) {
188
+ for (Line line : multiLine ) {
189
+ visit (line );
190
+ }
203
191
}
204
192
return null ;
205
193
}
206
194
207
195
@ Override
208
196
public Void visit (MultiPoint multiPoint ) {
209
- for (Point point : multiPoint ) {
210
- visit (point );
197
+ if (calculator .getDimensionalShapeType () == null || calculator .getDimensionalShapeType () == POINT ) {
198
+ for (Point point : multiPoint ) {
199
+ visit (point );
200
+ }
211
201
}
212
202
return null ;
213
203
}
@@ -222,16 +212,39 @@ public Void visit(MultiPolygon multiPolygon) {
222
212
223
213
@ Override
224
214
public Void visit (Point point ) {
225
- calculator .addCoordinate (point .getX (), point .getY (), 1.0 );
215
+ if (calculator .getDimensionalShapeType () == null || calculator .getDimensionalShapeType () == POINT ) {
216
+ visitPoint (point .getX (), point .getY ());
217
+ }
226
218
return null ;
227
219
}
228
220
229
221
@ Override
230
222
public Void visit (Polygon polygon ) {
231
- visit (polygon .getPolygon (), false );
223
+ // check area of polygon
224
+
225
+ double [] centroidX = new double [1 + polygon .getNumberOfHoles ()];
226
+ double [] centroidY = new double [1 + polygon .getNumberOfHoles ()];
227
+ double [] weight = new double [1 + polygon .getNumberOfHoles ()];
228
+ visitLinearRing (polygon .getPolygon ().length (), polygon .getPolygon ()::getX , polygon .getPolygon ()::getY , false ,
229
+ centroidX , centroidY , weight , 0 );
232
230
for (int i = 0 ; i < polygon .getNumberOfHoles (); i ++) {
233
- visit (polygon .getHole (i ), true );
231
+ visitLinearRing (polygon .getHole (i ).length (), polygon .getHole (i )::getX , polygon .getHole (i )::getY , true ,
232
+ centroidX , centroidY , weight , i + 1 );
234
233
}
234
+
235
+ double sumWeight = 0 ;
236
+ for (double w : weight ) {
237
+ sumWeight += w ;
238
+ }
239
+
240
+ if (sumWeight == 0 && calculator .dimensionalShapeType != POLYGON ) {
241
+ visitLine (polygon .getPolygon ().length (), polygon .getPolygon ()::getX , polygon .getPolygon ()::getY );
242
+ } else {
243
+ for (int i = 0 ; i < 1 + polygon .getNumberOfHoles (); i ++) {
244
+ calculator .addCoordinate (centroidX [i ], centroidY [i ], weight [i ], POLYGON );
245
+ }
246
+ }
247
+
235
248
return null ;
236
249
}
237
250
@@ -241,9 +254,73 @@ public Void visit(Rectangle rectangle) {
241
254
double sumY = rectangle .getMaxY () + rectangle .getMinY ();
242
255
double diffX = rectangle .getMaxX () - rectangle .getMinX ();
243
256
double diffY = rectangle .getMaxY () - rectangle .getMinY ();
244
- calculator .addCoordinate (sumX / 2 , sumY / 2 , Math .abs (diffX * diffY ));
257
+ if (diffX != 0 && diffY != 0 ) {
258
+ calculator .addCoordinate (sumX / 2 , sumY / 2 , Math .abs (diffX * diffY ), POLYGON );
259
+ } else if (diffX != 0 ) {
260
+ calculator .addCoordinate (sumX / 2 , rectangle .getMinY (), diffX , LINE );
261
+ } else if (diffY != 0 ) {
262
+ calculator .addCoordinate (rectangle .getMinX (), sumY / 2 , diffY , LINE );
263
+ } else {
264
+ visitPoint (rectangle .getMinX (), rectangle .getMinY ());
265
+ }
245
266
return null ;
246
267
}
268
+
269
+
270
+ private void visitPoint (double x , double y ) {
271
+ calculator .addCoordinate (x , y , 1.0 , POINT );
272
+ }
273
+
274
+ private void visitLine (int length , CoordinateSupplier x , CoordinateSupplier y ) {
275
+ // check line has length
276
+ double originDiffX = x .get (0 ) - x .get (1 );
277
+ double originDiffY = y .get (0 ) - y .get (1 );
278
+ if (originDiffX != 0 || originDiffY != 0 ) {
279
+ // a line's centroid is calculated by summing the center of each
280
+ // line segment weighted by the line segment's length in degrees
281
+ for (int i = 0 ; i < length - 1 ; i ++) {
282
+ double diffX = x .get (i ) - x .get (i + 1 );
283
+ double diffY = y .get (i ) - y .get (i + 1 );
284
+ double xAvg = (x .get (i ) + x .get (i + 1 )) / 2 ;
285
+ double yAvg = (y .get (i ) + y .get (i + 1 )) / 2 ;
286
+ double weight = Math .sqrt (diffX * diffX + diffY * diffY );
287
+ calculator .addCoordinate (xAvg , yAvg , weight , LINE );
288
+ }
289
+ } else {
290
+ visitPoint (x .get (0 ), y .get (0 ));
291
+ }
292
+ }
293
+
294
+ private void visitLinearRing (int length , CoordinateSupplier x , CoordinateSupplier y , boolean isHole ,
295
+ double [] centroidX , double [] centroidY , double [] weight , int idx ) {
296
+ // implementation of calculation defined in
297
+ // https://www.seas.upenn.edu/~sys502/extra_materials/Polygon%20Area%20and%20Centroid.pdf
298
+ //
299
+ // centroid of a ring is a weighted coordinate based on the ring's area.
300
+ // the sign of the area is positive for the outer-shell of a polygon and negative for the holes
301
+
302
+ int sign = isHole ? -1 : 1 ;
303
+ double totalRingArea = 0.0 ;
304
+ for (int i = 0 ; i < length - 1 ; i ++) {
305
+ totalRingArea += (x .get (i ) * y .get (i + 1 )) - (x .get (i + 1 ) * y .get (i ));
306
+ }
307
+ totalRingArea = totalRingArea / 2 ;
308
+
309
+ double sumX = 0.0 ;
310
+ double sumY = 0.0 ;
311
+ for (int i = 0 ; i < length - 1 ; i ++) {
312
+ double twiceArea = (x .get (i ) * y .get (i + 1 )) - (x .get (i + 1 ) * y .get (i ));
313
+ sumX += twiceArea * (x .get (i ) + x .get (i + 1 ));
314
+ sumY += twiceArea * (y .get (i ) + y .get (i + 1 ));
315
+ }
316
+ centroidX [idx ] = sumX / (6 * totalRingArea );
317
+ centroidY [idx ] = sumY / (6 * totalRingArea );
318
+ weight [idx ] = sign * Math .abs (totalRingArea );
319
+ }
247
320
}
248
321
322
+ @ FunctionalInterface
323
+ private interface CoordinateSupplier {
324
+ double get (int idx );
325
+ }
249
326
}
0 commit comments