Multi-Centroid Shootout

(WARNING: I have completely re-written 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 re-wrote it and have also added support for polygons and mutil-point 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, multi-part polygons and the generation of multi-point 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 half-way 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
    // half-way 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+((y2-y1)/(x2-x1))*(x-x1);
          } 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+((y2-y1)/(x2-x1))*(x-x1);
             if (ytest>y)
                ycrossing = ytest;
          }

          // If it crosses, then see if it is the crossing with the lowest or second lowest y-value
          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("Para-centroid 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 multi-part 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).

Multi-Part Polygons: Single Centroid

As indicated in the previous paragraph the PL/SQL implementation supports generating a centroid for multi-part 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.

Multi-Part Polygons: Multiple Centroids

The PL/SQL version has support for generating a multi-point sdo_geometry (SDO_GTYPE = x005) holding the centroid of each and every outer shell (EType 1003) of a multi-part 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 multi-point sdo_geometry object. Once all parts have been processed the resultant multi-point 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, multi-part 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 multi-point 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))

!/images/17.gif!

Hope all this helps in understanding how the centroid functions in my PL/SQL packages work.