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