[geos-devel] intersect/intersection() problem under Win

[Christian Pschierer]

Hi all,

we are currently writing a little application which checks geometries for
true intersections and returns the area of the intersection.

...
if( ( g1->intersects( g2 ) == TRUE ) &&
    ( g1->touches( g2 ) == FALSE ) ){
   intersectionArea = g1->intersection( g2 )->getArea();
...

On Linux this works pretty well. On Windows XP (Visual Studio 6.0) and
Irix (MIPSpro Compilers: Version 7.3.1.3m) we found 2 polygons which cause
an Exception in intersection().

We traced the error back to the file
 geos/source/graph/DirectedEdgeStar.cpp:178

if (state==LINKING_TO_OUTGOING) {
   if (firstOut==NULL)
      throw new TopologyException("no outgoing dirEdge
      found",&(getCoordinate()));
...

Any ideas what causes this problem and how to solve it? Our guess is, that
some kind of numerical precision problem leads to this behavior, as Linux
computes an area of 0 for in this case.

Below you can find a modified SimpleWKTTester-program which can reproduce
this error on Windows and Irix.


Greetings
Christian Pschierer


P.S.:
 The TopologyException-Class is missing the toString()-Method. 
Can someone with write-Access to the CVS-Repository please copy this
method from source/util/GEOSException.cpp to
source/geom/TopologyException.cpp? Thank you.


WKTIn:
POLYGON(( -79.63495899324135500000 43.68317101967772900000,
 -79.63535258791171100000 43.68308676161363700000,
 -79.63709417370184000000 43.68469600959869600000,
 -79.63582188599052400000 43.68383905676894600000,
 -79.63495899324135500000 43.68317101967772900000))
POLYGON(( -79.63544261367904700000 43.68354543050776300000,
 -79.63538454138728200000 43.68359271151414400000,
 -79.63528343774432000000 43.68352100382917100000,
 -79.63533917745682800000 43.68346535191968600000,
 -79.63544261367904700000 43.68354543050776300000))



Main.cpp:

#include <iostream>
#include <fstream>

#include "headers/io.h"

using namespace std;
using namespace geos;

int main(int argc, char** argv)
{
	 
	try {
		ifstream in("WKTIn");
		string instr;
		WKTReader *r = new WKTReader(new GeometryFactory(new
PrecisionModel(),10));
		Geometry *g1;
		Geometry *g2;
		Geometry *intersectg1g2;
		bool intersect;
		bool touch;

		getline(in,instr);
		g1=r->read(instr);
		getline(in,instr);
		g2=r->read(instr);

		intersect=g1->intersects(g2);
		if(intersect)
			cout << "g1 and g2 intersect!" << endl;
		else
			cout << "g1 and g2 do not intersect!" << endl;

		touch=g1->touches(g2);
		if(touch)
			cout << "g1 and g2 touch!" << endl;
		else
			cout << "g1 and g2 do not touch!" << endl;

		///////////////////////////////////////////
		// This leads to an exception under WinXP and Irix

		cout << "Before intersection" << endl;
		intersectg1g2=g1->intersection( g2 );
		cout << "After intersection" << endl;
		
		///////////////////////////////////////////

		if( intersectg1g2 ) {
			cout << "The intersection has an area of " <<
intersectg1g2->getArea() << endl;
		}
	}
	catch (GEOSException *ge) {
		cout << "Caught GEOSException!" << endl;
		cout << ge->toString() << endl;
	}
	catch (TopologyException *te) {
		cout << "Caught TopologyException!" << endl;
	}

	return 0;
}