Wednesday November 05 2008 at 01:00
(WARNING: I have completely rewritten my centroid code since this article when I discovered that the algorithm I had been supplied by a Tino failed in one important case. Instead of fixing the existing algorithm I completely rewrote it and have also added support for polygons and mutilpoint geometries. I will edit all centroid related articles some time soon.)
This blog follows on from the first “Centroid Shoot Out” and the comment left by Andy. It covers two issues. Firstly it publishes the original centroid code written by Tino Delbourgo when he worked for a company based in Hobart, Tasmania called Geometry Pty Ltd ; secondly, it shows how the PL/SQL version of the Java code handles holes in polygons, multipart polygons and the generation of multipoint centroids.
Original Java Code.
The original Java code is:
private static void calcParaCentroid(Geometry geometry, Double centroidX, Double centroidY)
{
Enumeration segs = ((Polygon)geometry).getAllSegments();
int segCount = 0;
Envelope envp = geometry.getEnvelope();
// Go halfway along bottom edge of bounding box to get our candidate point (x,y)
double x = (envp.getMinX()+envp.getMaxX())/2.0;
double y = envp.getMinY();
// Go through each line segment in turn, working out crossing points from the
// halfway along the bottom edge of the bounding box in a line due north.
// We only need to keep the lowest and second lowest of these points.
double lowestCrossing = 0.0;
double secondLowestCrossing = 0.0;
while (segs.hasMoreElements()) {
segCount++;
if (debugPinP)
System.out.println("Considering polygon segment #"+segCount);
LinearSegment seg = (LinearSegment)segs.nextElement();
double coords[] = seg.getCoordArray();
double x1 = coords[0];
double y1 = coords[1];
double x2 = 0.0;
double y2 = 0.0;
for (int i=2; i<coords.length; i+=2) {
double ycrossing = 0.0;
x2 = coords[i];
y2 = coords[i+1];
if (debugPinP)
System.out.println(" Considering line segment: ("+x1+","+y1+") to ("+x2+","+y2+")");
if ((x1<x && x2<x)  (x1>x && x2>x)  (y1<y && y2<y)) {
// do nothing  segment is either wholly to left, to right or below our point
} else if (y1>=y && y2>=y) {
// the line segment is above our point with one end to the left and the other to the right
//  this is a definite crossing
if ((x1<=x && x<=x2)  (x2<=x && x<=x1))
ycrossing = y1+((y2y1)/(x2x1))*(xx1);
} else if (x1==x && x2==x) {
// do nothing  we are comparing with a vertical line (no crossing possible)
} else if (y1==y && y2==y) {
// we are dealing with a horizontal line going through our point  a definite crossing
ycrossing = y;
} else {
// we are dealing with a line above our point with our point being within the line's bounding box
//  we need to project up to the line to see if it crosses
double ytest = y1+((y2y1)/(x2x1))*(xx1);
if (ytest>y)
ycrossing = ytest;
}
// If it crosses, then see if it is the crossing with the lowest or second lowest yvalue
if (ycrossing!=0.0) {
if ( ycrossing lowestCrossing ) {
lowestCrossing = ycrossing;
} else if ( ycrossing secondLowestCrossing ) {
secondLowestCrossing = ycrossing;
} else if (lowestCrossing==0.0 
ycrossing < lowestCrossing) {
secondLowestCrossing = lowestCrossing;
lowestCrossing = ycrossing;
} else if (secondLowestCrossing==0.0 
ycrossing < secondLowestCrossing) {
secondLowestCrossing = ycrossing;
}
}
// Go round the loop again with next segment
x1 = x2;
y1 = y2;
}
}
if (secondLowestCrossing==0.0  lowestCrossing==0.0)
throw new RuntimeException("Paracentroid calculation failed, couldn't find two crossing points up from centre bottom edge of bounding box");
centroidX = new Double(x);
centroidY = new Double((lowestCrossing+secondLowestCrossing)/2.0);
}
You will note that the code requires some functions in the old sdoapi.jar from 9i (the api was changed with the 10g release). The function it needs from this jar file is one that “vectorizes” an sdo_geometry polygon into simple start/end 2 point linestrings (or vectors): getAllSegments().
I have looked at deploying this code inside the Oracle database JVM with PL/SQL wrappers. It is possible but because the sdoapi.jar file is deprecated I decided the only long term solution was to use the Java Topology Suite and to build a vectorisation function using it. I have not time to do this. There are also other issues such as the conversion of Sdo_Geometry/JGeometry to JTS via GeoTools (I think that the current converter does not support circular arcs).
When I converted the Java to PL/SQL, I had to build my own vectoriser (see my article on spatial pipelining). I also added in support for rectangles, circles and compound linestrings/polygons with three point circular arcs. I added the ability to choose the largest part of a multipart geometry object (SdoGType x007). All these could be added to the Java code but have not been done. In creating the PL/SQL version I found a few centroid use cases that the original Java code didn’t implement. I have added these to the PL/SQL version (compiled and tested) and the Java version (never compiled or tested).
MultiPart Polygons: Single Centroid
As indicated in the previous paragraph the PL/SQL implementation supports generating a centroid for multipart polygons (SDO_GTYPE = x007). The method I use is simple: firstly, choose the largest part; secondly, generate the centroid as normal. The algorithm used gets the minimum bounding rectangle of each outer shell and selects the largest. Of course, this may still not be correct as in certain situations area would be better. Perhaps one day I will add this in.
MultiPart Polygons: Multiple Centroids
The PL/SQL version has support for generating a multipoint sdo_geometry (SDO_GTYPE = x005) holding the centroid of each and every outer shell (EType 1003) of a multipart polygon (Sdo_Gtype x007). The algorithm simple extracts each outer shell as a single (2003) polygon and passes it to the standard centroid function: each centroid created is appended to a multipoint sdo_geometry object. Once all parts have been processed the resultant multipoint geometry is returned.
Holes
There is little that I can say other than the code will never put a centroid into any of the holes (SDO_GTYPE = 2003) of any part of a polygon of any type.
Data and Diagram
Putting it all together. I have created a single, complex, multipart polygon (SDO_GTYPE = 2007) object (see below). For this polygon I have generated a centroid from the standard Oracle function sdo_geom.sdo_centroid, my own geom.sdo_centroid function and finally I have generated a single multipoint geometry via my geom.sdo_multi_centroid function.
codesys@XE> CREATE TABLE MultiCentroidShootOut AS
SELECT
SDO_GEOMETRY(2007, 28355, NULL,
SDO_ELEM_INFO_ARRAY(
1, 1003, 1,
25, 2003, 1,
35, 2003, 1,
45, 2003, 1,
63, 1003, 1,
87, 1003, 1,
99, 1003, 1),
SDO_ORDINATE_ARRAY(
17.0158371, 492.10068,
60.6568627, 602.4868,
186.445701, 612.75528,
378.979638, 712.87293,
740.943439, 712.87293,
669.064103, 515.20475,
835.926848, 284.16403,
879.567873, 140.40535,
645.96003, 16.188914,
299.398944, 13.6217949,
76.0595777, 171.21078,
17.0158371, 492.10068,
386.680995, 576.81561,
396.949472, 366.31184,
602.319005, 489.53356,
566.379336, 651.26207,
386.680995, 576.81561,
91.4622926, 474.13084,
168.475867, 399.68439,
307.100302, 481.8322,
181.311463, 558.84578,
91.4622926, 474.13084,
176.177225, 209.71757,
299.398944, 86.495852,
561.245098, 30.019231,
692.168175, 94.19721,
789.718703, 171.21078,
633.124434, 340.64065,
376.412519, 255.92572,
214.684012, 361.1776,
176.177225, 209.71757,
69, 9.5,
206, 86.5,
397, 185.5,
553, 189.5,
698, 143.5,
920, 7.5,
853, 105.5,
704, 259.5,
537, 307.5,
403, 271.5,
183, 134.5,
69, 9.5,
412.352187, 158.37519,
468.828808, 225.12029,
592.050528, 237.95588,
661.362745, 202.01621,
568.946456, 89.062971,
412.352187, 158.37519,
886, 424.5,
999, 357.5,
1153, 208.5,
1201, 41.5,
1165, 92.5,
1028, 312.5,
903, 426.5,
980, 289.5,
1079, 98.5,
1083, 57.5,
1037, 202.5,
886, 424.5)) as geom
from dual;
codesys@XE> select sdo_geom.sdo_centroid(geom,0.05) from MultiCentroidShootout;
SDO_GEOM.SDO_CENTROID(GEOM,0.05)(SDO_GTYPE, SDO_SRID, SDO_POINT(X, Y, Z), SDO_ELEM_INFO, SDO_ORDINATES)

SDO_GEOMETRY(2001, 28355, SDO_POINT_TYPE(545.679984, 227.04292, NULL), NULL, NULL)
codesys@XE> select geom.sdo_centroid(geom,0.05) from MultiCentroidShootout;
GEOM.SDO_CENTROID(GEOM,0.05)(SDO_GTYPE, SDO_SRID, SDO_POINT(X, Y, Z), SDO_ELEM_INFO, SDO_ORDINATES)

SDO_GEOMETRY(2001, 28355, SDO_POINT_TYPE(448.3, 657.6, NULL), NULL, NULL)
codesys@XE> select geom.sdo_multi_centroid(geom,0.05) from MultiCentroidShootout;
GEOM.SDO_MULTI_CENTROID(GEOM,0.05)(SDO_GTYPE, SDO_SRID, SDO_POINT(X, Y, Z), SDO_ELEM_INFO, SDO_ORDINATES)

SDO_GEOMETRY(2005, 28355, NULL, SDO_ELEM_INFO_ARRAY(1, 1, 4), SDO_ORDINATE_ARRAY(448.3, 657.6, 494.5, 242.1, 536.9, 167.8, 1043.5, 248.2))
Hope all this helps in understanding how the centroid functions in my PL/SQL packages work.

Comment