[geos-devel] symDifference, malfunction or misunderstanding?

[Frederik Ramm]

Dear all,

    apologies for using up about 95% of bandwidth on this list recently. 
I hope that my numerous queries, and your answers, will at least help 
others in my situation in the future through Google et al.

I have a problem with symDifference (in geos-3.0.0 as shipped with 
Ubuntu jaunty, but using the 3.2.0rc2 doesn't change things). Perhaps it 
is a misunderstanding on my part.

I have two symmetric linear rings formed like adjoining angular 
brackets, like this:

+------+------+
|      |(a)   |
|   +--+--+   |
|   |     |   |
|   +--+--+   |
|      |(b)   |
+------+------+

They touch in the middle. Here's the WKT:

LINEARRING (1 4, 3 4, 3 3, 2 3, 2 2, 3 2, 3 1, 1 1, 1 4)
LINEARRING (3 4, 5 4, 5 1, 3 1, 3 2, 4 2, 4 3, 3 3, 3 4)

If I compute the intersection between both, using x->intersection(y), 
then I get what one would expect, namely the two segments labelled "a" 
and "b" in the above diagram:

MULTILINESTRING ((3.00 4.00, 3.00 3.00), (3.00 2.00, 3.00 1.00))

When I do a "symDifference" between the two, I would expect to get the 
original shape minus the "a" and "b" segments, something like

+------+------+
|             |
|   +--+--+   |
|   |     |   |
|   +--+--+   |
|             |
+------+------+

Instead, I get a four-part Multilinestring that contains the original 
shape minus the segment "b", but including the segment "a":

MULTILINESTRING ((1.00 4.00, 3.00 4.00), (3.00 3.00, 2.00 3.00, 2.00 
2.00, 3.00 2.00), (3.00 1.00, 1.00 1.00, 1.00 4.00), (3.00 4.00, 5.00 
4.00, 5.00 1.00, 3.00 1.00), (3.00 2.00, 4.00 2.00, 4.00 3.00, 3.00 3.00))

Just to be on the safe side, I then combined both linear rings into a 
multilinestring (using Union), and then computed the difference between 
this union and the intersection - again expecting to get the shape minus 
"a" and "b", and again getting the shape minus "a" but including "b".

Is this a bug in GEOS, or a bug in my approach?

I'm attaching the full C++ test file.

Bye
Frederik

#include <stdio.h>
#include <geos/geom/PrecisionModel.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Geometry.h>
#include <geos/io/WKTReader.h>

using namespace geos::geom;
using namespace geos;

GeometryFactory *global_factory;

int main(int argc, char **argv)
{
     PrecisionModel *pm = new PrecisionModel(2.0, 0, 0);
     global_factory = new GeometryFactory(pm, -1);
     delete pm;

     io::WKTReader *reader = new io::WKTReader(global_factory);
     Geometry *a = reader->read(
         "LINEARRING (1 4, 3 4, 3 3, 2 3, 2 2, 3 2, 3 1, 1 1, 1 4)"
         );
     Geometry *b = reader->read(
         "LINEARRING (3 4, 5 4, 5 1, 3 1, 3 2, 4 2, 4 3, 3 3, 3 4)"
         );

     Geometry *inter = a->intersection(b);
     Geometry *diff = a->symDifference(b);

     printf("intersection: %s\n", inter->toString().c_str());
     printf("sym difference: %s\n", diff->toString().c_str());
}

-- 
Frederik Ramm  ##  eMail frederik at remote.org  ##  N49°00'09" E008°23'33"