Complete Note#1 in the http://wiki.osgeo.org/wiki/GEOS_Provenance_Review to get out...

[geos.git] / doc / example.cpp

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2001-2002 Vivid Solutions Inc.
 * Copyright (C) 2005 Refractions Research Inc.
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation. 
 * See the COPYING file for more information.
 *
 *********************************************************************
 *
 * This file should document by example usage of the GEOS library.
 * It could actually be a live discuss-by-example board for
 * architectural design choices.
 *
 *                      --strk;
 * 
 * DEBUGGING TIPS:
 *  use -D__USE_MALLOC at compile time for gcc 2.91, 2.95, 3.0 and 3.1
 *  and GLIBCXX_FORCE_NEW or GLIBCPP_FORCE_NEW at run time with gcc 3.2.2+
 *  to force libstdc++ avoid caching memory. This should remove some
 *  obscure reports from memory checkers like valgrind.
 *
 **********************************************************************/

#include <geos/geom/PrecisionModel.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/Point.h>
#include <geos/geom/LinearRing.h>
#include <geos/geom/LineString.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/CoordinateArraySequence.h>
#include <geos/geom/IntersectionMatrix.h>
#include <geos/io/WKBReader.h>
#include <geos/io/WKBWriter.h>
#include <geos/io/WKTWriter.h>
#include <geos/util/GeometricShapeFactory.h>
#include <geos/geom/util/SineStarFactory.h>
#include <geos/util/GEOSException.h>
#include <geos/util/IllegalArgumentException.h>
#include <geos/opLinemerge.h>
#include <geos/opPolygonize.h>
#include <vector>
#include <sstream>
#include <iomanip>
#include <cstdlib> // exit()

#ifndef M_PI
#define M_PI       3.14159265358979323846
#endif

// Set to 0 to skip section
#define GEOMETRIC_SHAPES 1
#define RELATIONAL_OPERATORS 1
#define COMBINATIONS 1
#define UNARY_OPERATIONS 1
#define LINEMERGE 1
#define POLYGONIZE 1


using namespace std;
using namespace geos;
using namespace geos::geom;
using namespace geos::operation::polygonize;
using namespace geos::operation::linemerge;
using geos::util::GEOSException;
using geos::util::IllegalArgumentException;


// Prototypes
void wkt_print_geoms(vector<Geometry *> *geoms);


// This object will be used to construct our geometries.
// It might be bypassed by directly call geometry constructors,
// but that would be boring because you'd need to specify
// a PrecisionModel and a SRID everytime: those infos are
// cached inside a GeometryFactory object.
GeometryFactory *global_factory;

//#define DEBUG_STREAM_STATE 1


//
// This function tests writing and reading WKB
// TODO:
//      - compare input and output geometries for equality
//      - remove debugging lines (on stream state)
//
void WKBtest(vector<Geometry*>*geoms)
{
        stringstream s(ios_base::binary|ios_base::in|ios_base::out);
        io::WKBReader wkbReader(*global_factory);
        io::WKBWriter wkbWriter;
        Geometry *gout;

#if DEBUG_STREAM_STATE
        cout<<"WKBtest: machine byte order: "<<BYTE_ORDER<<endl;
#endif


        unsigned int ngeoms=geoms->size();
        for (unsigned int i=0; i<ngeoms; ++i)
        {
                Geometry *gin = (*geoms)[i];

#if DEBUG_STREAM_STATE
                cout<<"State of stream before WRITE: ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<
                        " good:"<<s.good()<<
                        " eof:"<<s.eof()<<
                        " bad:"<<s.bad()<<
                        " fail:"<<s.fail()<<endl; 
#endif

#if DEBUG_STREAM_STATE
                cout<<"State of stream after SEEKP: ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<
                        " good:"<<s.good()<<
                        " eof:"<<s.eof()<<
                        " bad:"<<s.bad()<<
                        " fail:"<<s.fail()<<endl; 
#endif

                wkbWriter.write(*gin, s);
#if DEBUG_STREAM_STATE
                cout<<"wkbWriter wrote and reached ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<endl; 

                cout<<"State of stream before DUMP: ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<
                        " good:"<<s.good()<<
                        " eof:"<<s.eof()<<
                        " bad:"<<s.bad()<<
                        " fail:"<<s.fail()<<endl; 
#endif

#if DEBUG_STREAM_STATE
                cout<<"State of stream after DUMP: ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<
                        " good:"<<s.good()<<
                        " eof:"<<s.eof()<<
                        " bad:"<<s.bad()<<
                        " fail:"<<s.fail()<<endl; 
#endif

                s.seekg(0, ios::beg); // rewind reader pointer

#if DEBUG_STREAM_STATE
                cout<<"State of stream before READ: ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<
                        " good:"<<s.good()<<
                        " eof:"<<s.eof()<<
                        " bad:"<<s.bad()<<
                        " fail:"<<s.fail()<<endl; 
#endif

                gout = wkbReader.read(s);

#if DEBUG_STREAM_STATE
                cout<<"State of stream after READ: ";
                cout<<"p:"<<s.tellp()<<" g:"<<s.tellg()<<
                        " good:"<<s.good()<<
                        " eof:"<<s.eof()<<
                        " bad:"<<s.bad()<<
                        " fail:"<<s.fail()<<endl; 
#endif

                gin->normalize();
                gout->normalize();
                int failed = gin->compareTo(gout);
                if ( failed ) cout<<"{"<<i<<"} (WKB) ";
                else cout<<"["<<i<<"] (WKB) ";

                io::WKBReader::printHEX(s, cout);
                cout<<endl;

                if ( failed ) {
                        io::WKTWriter wkt;
                        cout<<"  IN: "<<wkt.write(gin)<<endl;
                        cout<<" OUT: "<<wkt.write(gout)<<endl;
                }

                s.seekp(0, ios::beg); // rewind writer pointer

                delete gout;
        }

}


// This function will print given geometries in WKT
// format to stdout. As a side-effect, will test WKB
// output and input, using the WKBtest function.
void
wkt_print_geoms(vector<Geometry *> *geoms)
{
        WKBtest(geoms); // test WKB parser

        // WKT-print given geometries
        io::WKTWriter *wkt = new io::WKTWriter();
        for (unsigned int i=0; i<geoms->size(); i++) {
                const Geometry *g = (*geoms)[i];
                string tmp=wkt->write(g);
                cout<<"["<<i<<"] (WKT) "<<tmp<<endl;
        }
        delete wkt;
}

// This is the simpler geometry you can get: a point.
Point *
create_point(double x, double y)
{
        Coordinate c(x, y);
        Point *p = global_factory->createPoint(c);
        return p;
}

// This function will create a LinearString
// geometry with the shape of the letter U
// having top-left corner at given coordinates
// and 'side' height and width 
LineString *
create_ushaped_linestring(double xoffset, double yoffset, double side)
{
        // We will use a coordinate list to build the linestring
        CoordinateSequence *cl = new CoordinateArraySequence();

        cl->add(Coordinate(xoffset, yoffset));
        cl->add(Coordinate(xoffset, yoffset+side));
        cl->add(Coordinate(xoffset+side, yoffset+side));
        cl->add(Coordinate(xoffset+side, yoffset));

        // Now that we have a CoordinateSequence we can create 
        // the linestring.
        // The newly created LineString will take ownership
        // of the CoordinateSequence.
        LineString *ls = global_factory->createLineString(cl);

        // This is what you do if you want the new LineString
        // to make a copy of your CoordinateSequence:
        // LineString *ls = global_factory->createLineString(*cl);
        
        return ls; // our LineString
}

// This function will create a LinearRing
// geometry rapresenting a square with the given origin 
// and side 
LinearRing *
create_square_linearring(double xoffset, double yoffset, double side)
{
        // We will use a coordinate list to build the linearring
        CoordinateSequence *cl = new CoordinateArraySequence();

        cl->add(Coordinate(xoffset, yoffset));
        cl->add(Coordinate(xoffset, yoffset+side));
        cl->add(Coordinate(xoffset+side, yoffset+side));
        cl->add(Coordinate(xoffset+side, yoffset));
        cl->add(Coordinate(xoffset, yoffset));

        // Now that we have a CoordinateSequence we can create 
        // the linearring.
        // The newly created LinearRing will take ownership
        // of the CoordinateSequence.
        LinearRing *lr = global_factory->createLinearRing(cl);

        // This is what you do if you want the new LinearRing
        // to make a copy of your CoordinateSequence:
        // LinearRing *lr = global_factory->createLinearRing(*cl);
        
        return lr; // our LinearRing
}

// This function will create a Polygon
// geometry rapresenting a square with the given origin 
// and side and with a central hole 1/3 sided.
Polygon *
create_square_polygon(double xoffset, double yoffset, double side)
{
        // We need a LinearRing for the polygon shell 
        LinearRing *outer = create_square_linearring(xoffset,yoffset,side);

        // And another for the hole 
        LinearRing *inner = create_square_linearring(xoffset+(side/3),
                        yoffset+(side/3),(side/3));
        
        // If we need to specify any hole, we do it using
        // a vector of Geometry pointers (I don't know why
        // not LinearRings)
        vector<Geometry *> *holes = new vector<Geometry *>;

        // We add the newly created geometry to the vector
        // of holes.
        holes->push_back(inner);

        // And finally we call the polygon constructor.
        // Both the outer LinearRing and the vector of holes
        // will be referenced by the resulting Polygon object,
        // thus we CANNOT delete them, neither the holes, nor
        // the vector containing their pointers, nor the outer
        // LinearRing. Everything will be deleted at Polygon
        // deletion time (this is inconsistent with LinearRing
        // behaviour... what should we do?).
        Polygon *poly = global_factory->createPolygon(outer, holes);

        return poly;
}

//
// This function will create a GeometryCollection
// containing copies of all Geometries in given vector.
//
GeometryCollection *
create_simple_collection(vector<Geometry *> *geoms)
{
        return global_factory->createGeometryCollection(*geoms);
        // if you wanted to transfer ownership of vector end
        // its elements you should have call:
        // return global_factory->createGeometryCollection(geoms);
}

//
// This function uses GeometricShapeFactory to render
// a circle having given center and radius
//
Polygon *
create_circle(double centerX, double centerY, double radius)
{
        geos::util::GeometricShapeFactory shapefactory(global_factory);
        shapefactory.setCentre(Coordinate(centerX, centerY));
        shapefactory.setSize(radius);
        // same as:
        //      shapefactory.setHeight(radius);
        //      shapefactory.setWidth(radius);
        return shapefactory.createCircle();
}

//
// This function uses GeometricShapeFactory to render
// an ellipse having given center and axis size
//
Polygon *
create_ellipse(double centerX, double centerY, double width, double height)
{
        geos::util::GeometricShapeFactory shapefactory(global_factory);
        shapefactory.setCentre(Coordinate(centerX, centerY));
        shapefactory.setHeight(height);
        shapefactory.setWidth(width);
        return shapefactory.createCircle();
}

//
// This function uses GeometricShapeFactory to render
// a rectangle having lower-left corner at given coordinates
// and given sizes.
//
Polygon *
create_rectangle(double llX, double llY, double width, double height)
{
        geos::util::GeometricShapeFactory shapefactory(global_factory);
        shapefactory.setBase(Coordinate(llX, llY));
        shapefactory.setHeight(height);
        shapefactory.setWidth(width);
        shapefactory.setNumPoints(4); // we don't need more then 4 points for a rectangle...
        // can use setSize for a square
        return shapefactory.createRectangle();
}

//
// This function uses GeometricShapeFactory to render
// an arc having lower-left corner at given coordinates,
// given sizes and given angles. 
//
LineString *
create_arc(double llX, double llY, double width, double height, double startang, double endang)
{
        geos::util::GeometricShapeFactory shapefactory(global_factory);
        shapefactory.setBase(Coordinate(llX, llY));
        shapefactory.setHeight(height);
        shapefactory.setWidth(width);
        // shapefactory.setNumPoints(100); // the default (100 pts)
        // can use setSize for a square
        return shapefactory.createArc(startang, endang);
}

auto_ptr<Polygon>
create_sinestar(double cx, double cy, double size, int nArms, double armLenRat)
{
        geos::geom::util::SineStarFactory fact(global_factory);
        fact.setCentre(Coordinate(cx, cy));
        fact.setSize(size);
        fact.setNumPoints(nArms*5);
  fact.setArmLengthRatio(armLenRat);
        fact.setNumArms(nArms);
        return fact.createSineStar();
}

// Start reading here
void do_all()
{
        vector<Geometry *> *geoms = new vector<Geometry *>;
        vector<Geometry *> *newgeoms;
        
        // Define a precision model using 0,0 as the reference origin
        // and 2.0 as coordinates scale.
        PrecisionModel *pm = new PrecisionModel(2.0, 0, 0);

        // Initialize global factory with defined PrecisionModel
        // and a SRID of -1 (undefined).
        global_factory = new GeometryFactory(pm, -1);

        // We do not need PrecisionMode object anymore, it has
        // been copied to global_factory private storage
        delete pm;

////////////////////////////////////////////////////////////////////////
// GEOMETRY CREATION
////////////////////////////////////////////////////////////////////////

        // Read function bodies to see the magic behind them
        geoms->push_back(create_point(150, 350));
        geoms->push_back(create_square_linearring(0,0,100));
        geoms->push_back(create_ushaped_linestring(60,60,100));
        geoms->push_back(create_square_linearring(0,0,100));
        geoms->push_back(create_square_polygon(0,200,300));
        geoms->push_back(create_square_polygon(0,250,300));
        geoms->push_back(create_simple_collection(geoms));

#if GEOMETRIC_SHAPES
        // These ones use a GeometricShapeFactory
        geoms->push_back(create_circle(0, 0, 10));
        geoms->push_back(create_ellipse(0, 0, 8, 12));
        geoms->push_back(create_rectangle(-5, -5, 10, 10)); // a square
        geoms->push_back(create_rectangle(-5, -5, 10, 20)); // a rectangle
        // The upper-right quarter of a vertical ellipse
        geoms->push_back(create_arc(0, 0, 10, 20, 0, M_PI/2));
        geoms->push_back(create_sinestar(10, 10, 100, 5, 2).release()); // a sine star
#endif

        // Print all geoms.
        cout<<"--------HERE ARE THE BASE GEOMS ----------"<<endl;
        wkt_print_geoms(geoms);


#if UNARY_OPERATIONS

////////////////////////////////////////////////////////////////////////
// UNARY OPERATIONS
////////////////////////////////////////////////////////////////////////

        /////////////////////////////////////////////
        // CENTROID
        /////////////////////////////////////////////
        
        // Find centroid of each base geometry
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g = (*geoms)[i];
                newgeoms->push_back( g->getCentroid() );
        }

        // Print all convex hulls
        cout<<endl<<"------- AND HERE ARE THEIR CENTROIDS -----"<<endl;
        wkt_print_geoms(newgeoms);

        // Delete the centroids
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;
        
        /////////////////////////////////////////////
        // BUFFER
        /////////////////////////////////////////////
        
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g = (*geoms)[i];
                try {
                        Geometry *g2 = g->buffer(10);
                        newgeoms->push_back(g2);
                }
                catch (const GEOSException& exc) {
                        cerr <<"GEOS Exception: geometry "<<i<<"->buffer(10): "<<exc.what()<<"\n";
                }
        }

        cout<<endl<<"--------HERE COMES THE BUFFERED GEOMS ----------"<<endl;
        wkt_print_geoms(newgeoms);

        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;
        
        /////////////////////////////////////////////
        // CONVEX HULL
        /////////////////////////////////////////////
        
        // Make convex hulls of geometries
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g = (*geoms)[i];
                newgeoms->push_back( g->convexHull() );
        }

        // Print all convex hulls
        cout<<endl<<"--------HERE COMES THE HULLS----------"<<endl;
        wkt_print_geoms(newgeoms);

        // Delete the hulls
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;

#endif // UNARY_OPERATIONS

#if RELATIONAL_OPERATORS

////////////////////////////////////////////////////////////////////////
// RELATIONAL OPERATORS
////////////////////////////////////////////////////////////////////////

cout<<"-------------------------------------------------------------------------------"<<endl;
cout<<"RELATIONAL OPERATORS"<<endl;
cout<<"-------------------------------------------------------------------------------"<<endl;

        /////////////////////////////////////////////
        // DISJOINT
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"   DISJOINT   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->disjoint(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }
        
        /////////////////////////////////////////////
        // TOUCHES
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"    TOUCHES   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->touches(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // INTERSECTS
        /////////////////////////////////////////////

        cout<<endl;
        cout<<" INTERSECTS   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->intersects(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // CROSSES
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"    CROSSES   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->crosses(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // WITHIN
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"     WITHIN   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->within(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // CONTAINS
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"   CONTAINS   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->contains(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // OVERLAPS
        /////////////////////////////////////////////
        
        cout<<endl;
        cout<<"   OVERLAPS   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->overlaps(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // RELATE
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"     RELATE   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        IntersectionMatrix *im=NULL;
                        try {
                                // second argument is intersectionPattern
                                string pattern = "212101212";
                                if ( g1->relate(g2, pattern) ) cout<<" 1\t";
                                else cout<<" 0\t";
                                
                                // get the intersectionMatrix itself
                                im=g1->relate(g2);
                                delete im; // delete afterwards
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // EQUALS
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"     EQUALS   ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                if ( g1->equals(g2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // EQUALS_EXACT
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"EQUALS_EXACT  ";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                // second argument is a tolerance
                                if ( g1->equalsExact(g2, 0.5) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

        /////////////////////////////////////////////
        // IS_WITHIN_DISTANCE
        /////////////////////////////////////////////

        cout<<endl;
        cout<<"IS_WITHIN_DIST";
        for (unsigned int i=0; i<geoms->size(); i++) {
                cout<<"\t["<<i<<"]";
        }
        cout<<endl;
        for (unsigned int i=0; i<geoms->size(); i++) {
                Geometry *g1 = (*geoms)[i];
                cout<<"      ["<<i<<"]\t";
                for (unsigned int j=0; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                // second argument is the distance
                                if ( g1->isWithinDistance(g2,2) ) cout<<" 1\t";
                                else cout<<" 0\t";
                        }
                        // Geometry Collection is not a valid argument
                        catch (const IllegalArgumentException& exc) {
                                cout<<" X\t";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
                cout<<endl;
        }

#endif // RELATIONAL_OPERATORS

#if COMBINATIONS

////////////////////////////////////////////////////////////////////////
// COMBINATIONS
////////////////////////////////////////////////////////////////////////

cout<<endl;
cout<<"-------------------------------------------------------------------------------"<<endl;
cout<<"COMBINATIONS"<<endl;
cout<<"-------------------------------------------------------------------------------"<<endl;

        /////////////////////////////////////////////
        // UNION
        /////////////////////////////////////////////

        // Make unions of all geoms
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size()-1; i++) {
                Geometry *g1 = (*geoms)[i];
                for (unsigned int j=i+1; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                Geometry *g3 = g1->Union(g2);
                                newgeoms->push_back(g3);
                        }
                        // It's illegal to union a collection ...
                        catch (const IllegalArgumentException& ill) {
                                //cerr <<ill.toString()<<"\n";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
        }

        // Print all unions
        cout<<endl<<"----- AND HERE ARE SOME UNION COMBINATIONS ------"<<endl;
        wkt_print_geoms(newgeoms);

        // Delete the resulting geoms
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;
        

        /////////////////////////////////////////////
        // INTERSECTION
        /////////////////////////////////////////////
        
        // Compute intersection of adhiacent geometries
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size()-1; i++) {
                Geometry *g1 = (*geoms)[i];
                for (unsigned int j=i+1; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                Geometry *g3 = g1->intersection(g2);
                                newgeoms->push_back(g3);
                        }
                        // Collection are illegal as intersection argument
                        catch (const IllegalArgumentException& ill) {
                                //cerr <<ill.toString()<<"\n";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
        }

        cout<<endl<<"----- HERE ARE SOME INTERSECTIONS COMBINATIONS ------"<<endl;
        wkt_print_geoms(newgeoms);

        // Delete the resulting geoms
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;

        /////////////////////////////////////////////
        // DIFFERENCE
        /////////////////////////////////////////////
        
        // Compute difference of adhiacent geometries
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size()-1; i++) {
                Geometry *g1 = (*geoms)[i];
                for (unsigned int j=i+1; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                Geometry *g3 = g1->difference(g2);
                                newgeoms->push_back(g3);
                        }
                        // Collection are illegal as difference argument
                        catch (const IllegalArgumentException& ill) {
                                //cerr <<ill.toString()<<"\n";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
        }

        cout<<endl<<"----- HERE ARE SOME DIFFERENCE COMBINATIONS ------"<<endl;
        wkt_print_geoms(newgeoms);

        // Delete the resulting geoms
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;
        
        /////////////////////////////////////////////
        // SYMMETRIC DIFFERENCE
        /////////////////////////////////////////////
        
        // Compute symmetric difference of adhiacent geometries
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<geoms->size()-1; i++) {
                Geometry *g1 = (*geoms)[i];
                for (unsigned int j=i+1; j<geoms->size(); j++) {
                        Geometry *g2 = (*geoms)[j];
                        try {
                                Geometry *g3 = g1->symDifference(g2);
                                newgeoms->push_back(g3);
                        }
                        // Collection are illegal as symdifference argument
                        catch (const IllegalArgumentException& ill) {
                                //cerr <<ill.toString()<<"\n";
                        }
                        catch (const std::exception& exc) {
                                cerr<<exc.what()<<endl;
                        }
                }
        }

        cout<<endl<<"----- HERE ARE SYMMETRIC DIFFERENCES ------"<<endl;
        wkt_print_geoms(newgeoms);

        // Delete the resulting geoms
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;

#endif // COMBINATIONS

#if LINEMERGE
        
        /////////////////////////////////////////////
        // LINEMERGE
        /////////////////////////////////////////////
        LineMerger lm;
        lm.add(geoms);
        vector<LineString *> *mls = lm.getMergedLineStrings();
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<mls->size(); i++)
                newgeoms->push_back((*mls)[i]);
        delete mls;

        cout<<endl<<"----- HERE IS THE LINEMERGE OUTPUT ------"<<endl;
        wkt_print_geoms(newgeoms);
        
        // Delete the resulting geoms
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;

#endif // LINEMERGE

#if POLYGONIZE

        /////////////////////////////////////////////
        // POLYGONIZE
        /////////////////////////////////////////////
        Polygonizer plgnzr;
        plgnzr.add(geoms);
        vector<Polygon *> *polys = plgnzr.getPolygons();
        newgeoms = new vector<Geometry *>;
        for (unsigned int i=0; i<polys->size(); i++)
                newgeoms->push_back((*polys)[i]);
        delete polys;

        cout<<endl<<"----- HERE IS POLYGONIZE OUTPUT ------"<<endl;
        wkt_print_geoms(newgeoms);
        
        // Delete the resulting geoms
        for (unsigned int i=0; i<newgeoms->size(); i++) {
                delete (*newgeoms)[i];
        }
        delete newgeoms;

#endif // POLYGONIZE

        /////////////////////////////////////////////
        // CLEANUP
        /////////////////////////////////////////////

        // Delete base geometries 
        for (unsigned int i=0; i<geoms->size(); i++) {
                delete (*geoms)[i];
        }
        delete geoms;

        delete global_factory;
}

int
main()
{
        cout<<"GEOS "<<geosversion()<<" ported from JTS "<<jtsport()<<endl;
        try
        {
                do_all();
        }
        // All exception thrown by GEOS are subclasses of this
        // one, so this is a catch-all 
        catch (const GEOSException& exc)
        {
                cerr <<"GEOS Exception: "<<exc.what()<<"\n";
                exit(1);
        }
        catch (const exception &e)
        {
                cerr <<"Standard exception thrown: "<<e.what()<<endl;
                exit(1);
        }
        // and this is a catch-all non standard ;)
        catch (...)
        {
                cerr <<"unknown exception trown!\n";
                exit(1);
        }

        // Unload is no more necessary
        //io::Unload::Release();

        exit(0);
}

/**********************************************************************
 * $Log$
 * Revision 1.46  2006/04/09 11:07:54  mloskot
 * Small fixes in doc/example.cpp.
 *
 * Revision 1.45  2006/04/04 08:16:46  strk
 * Changed GEOSException hierarchy to be derived from std::runtime_exception.
 * Removed the GEOSException::toString redundant method (use ::what() instead)
 *
 * Revision 1.44  2006/03/28 15:19:22  strk
 * Added macros for sections skip (useful in debugging)
 *
 * Revision 1.43  2006/03/15 18:44:51  strk
 * Bug #60 - Missing <cmath> header in some files
 *
 * Revision 1.42  2006/03/06 21:27:39  strk
 * Cascading fixed after Unload definition moved to geos::io namespace
 *
 * Revision 1.41  2006/03/06 19:40:46  strk
 * geos::util namespace. New GeometryCollection::iterator interface, many cleanups.
 *
 * Revision 1.40  2006/03/06 15:23:14  strk
 * geos::io namespace
 *
 * Revision 1.39  2006/03/03 10:46:21  strk
 * Removed 'using namespace' from headers, added missing headers in .cpp files, removed useless includes in headers (bug#46)
 *
 * Revision 1.38  2006/02/19 19:46:49  strk
 * Packages <-> namespaces mapping for most GEOS internal code (uncomplete, but working). Dir-level libs for index/ subdirs.
 *
 * Revision 1.37  2006/02/09 15:52:47  strk
 * GEOSException derived from std::exception; always thrown and cought by const ref.
 *
 * Revision 1.36  2006/01/31 19:07:33  strk
 * - Renamed DefaultCoordinateSequence to CoordinateArraySequence.
 * - Moved GetNumGeometries() and GetGeometryN() interfaces
 *   from GeometryCollection to Geometry class.
 * - Added getAt(int pos, Coordinate &to) funtion to CoordinateSequence class.
 * - Reworked automake scripts to produce a static lib for each subdir and
 *   then link all subsystem's libs togheter
 * - Moved C-API in it's own top-level dir capi/
 * - Moved source/bigtest and source/test to tests/bigtest and test/xmltester
 * - Fixed PointLocator handling of LinearRings
 * - Changed CoordinateArrayFilter to reduce memory copies
 * - Changed UniqueCoordinateArrayFilter to reduce memory copies
 * - Added CGAlgorithms::isPointInRing() version working with
 *   Coordinate::ConstVect type (faster!)
 * - Ported JTS-1.7 version of ConvexHull with big attention to
 *   memory usage optimizations.
 * - Improved XMLTester output and user interface
 * - geos::geom::util namespace used for geom/util stuff
 * - Improved memory use in geos::geom::util::PolygonExtractor
 * - New ShortCircuitedGeometryVisitor class
 * - New operation/predicate package
 *
 * Revision 1.35  2005/11/30 11:27:07  strk
 * catch std::exception by ref
 *
 * Revision 1.34  2005/11/30 11:25:12  strk
 * includes cleanup
 *
 * Revision 1.33  2005/09/27 16:20:43  strk
 * Reverted previous change, fixed the bug by turning WKBReader.factory
 * into a reference rather then a real object. ABI still breaks, but API
 * at least is safe (we didn't release any WKB-aware package so breaking
 * ABI is not a big deal at this stage).
 *
 * Revision 1.32  2005/09/27 16:00:26  strk
 * Fixed bug in WKBReader destroying the GeometryFactory used in Geometry
 * construction. Changed it's definition to *require* a GeometryFactory
 * pointer parameter.
 *
 * Revision 1.31  2005/09/26 08:17:19  strk
 * Removed memory leak from WKB tester
 *
 * Revision 1.30  2005/09/03 21:26:42  strk
 * Reworked WKB I/O to avoid use of templates and make better use of STL
 *
 * Revision 1.29  2005/07/11 12:17:26  strk
 * Commented out useless include
 *
 * Revision 1.28  2005/04/29 16:36:28  strk
 * Made WKBReader use global_factory, for having WKB reads produce
 * same context of input geoms.
 *
 * Revision 1.27  2005/04/29 15:34:20  strk
 * Typedef'ed biostringstream, preferred parameter for
 * WKB parser templates.
 * Added << operator for biostringstream.
 * Typedef'ed WKBWriter and WKBReader to be parametrized by
 * biostringstream.
 * Added WKBtest in doc/example.cpp
 *
 * Revision 1.26  2004/12/08 13:54:43  strk
 * gcc warnings checked and fixed, general cleanups.
 *
 * Revision 1.25  2004/10/13 10:03:02  strk
 * Added missing linemerge and polygonize operation.
 * Bug fixes and leaks removal from the newly added modules and
 * planargraph (used by them).
 * Some comments and indentation changes.
 *
 * Revision 1.24  2004/07/22 16:58:01  strk
 * runtime version extractor functions split. geos::version() is now
 * geos::geosversion() and geos::jtsport()
 *
 * Revision 1.23  2004/07/17 09:19:32  strk
 * added GEOS version report
 *
 *********************************************************************/