File Coverage

blib/lib/Geo/OGC/Geometry.pm

Criterion	Covered	Total	%
statement	716	1471	48.6
branch	220	604	36.4
condition	36	129	27.9
subroutine	124	269	46.1
pod	34	42	80.9
total	1130	2515	44.9

=pod

=head1 NAME

Geo::OGC::Geometry - Perl extension for OGC simple feature geometries

=head1 SYNOPSIS

    use Geo::OGC::Geometry;

    my $point = Geo::OGC::Point->new(Text => 'POINT 1 1');

=head1 DESCRIPTION

OGC simple feature geometries is a class hierarchy for geographic
information. This module conforms to the document OGC 06-103r4, which
is the current document describing the standard. The document is
available from L.

This module is currently mostly an interface and storage classes and
does not contain implementations of most of the methods for testing
spatial relations etc.

=head1 CLASSES

=head2 Geo::OGC::Geometry

Geo::OGC::Geometry is the root class and it represents an arbitrary
geospatial geometry.

=cut

package Geo::OGC::Geometry;

use strict;
use POSIX;
use Carp;

BEGIN {
    use Exporter 'import';
    use vars qw /$SNAP_DISTANCE_SQR/;
    POSIX::setlocale( &POSIX::LC_NUMERIC, "C" ); # Force to "C" locale
    our %EXPORT_TAGS = ( 'all' => [ qw( &ccw &intersect &distance_point_line_sqr
					$SNAP_DISTANCE_SQR ) ] );
    our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );
    our $VERSION = '0.06';
    $SNAP_DISTANCE_SQR = 1E-6;
}

=pod

=over

=item new(%params)

Create a new Geometry object or a new object belonging to a concrete
subclass of Geometry if initialization data is given.

%param is a hash containing named parameters. The keywords are
(subclasses of Geometry may add new keywords):

=over

=item Text

A well-known text to initialize an object of a concrete class.

=item SRID

Sets the Spatial Reference System ID of the object, default is -1.

=item Precision

If specified, has the effect that numeric comparisons in the Equals
method is is preceded with a rounding operation (using sprintf "%.pe",
where p is the Precision-1, i.e the number of meaningful numbers is
Precision). Affects also AsText.

=back

=cut

sub new {
    my($package, %params) = @_;
    return parse_wkt($params{Text}) if exists $params{Text};
    my $self = {};
    bless $self => (ref($package) or $package);
    $self->init(%params);
    return $self;
}

# methods beginning with lower case letter are private
# and thus not in pod

# Override in new classes but call $self->SUPER::init(%params);
sub init {
    my($self, %params) = @_;
    $self->{SRID} = exists $params{SRID} ? $params{SRID} : -1;
    $self->{Precision} = $params{Precision}-1 if exists $params{Precision};
}

# Copies the attributes of this object to the other object,
# which is a newly created object of same class.
sub copy {
    my($self, $clone) = @_;
    $clone->{SRID} = $self->{SRID};
    $clone->{Precision} = $self->{Precision} if exists $self->{Precision};
}

# parse well known text and construct respective geometry
sub parse_wkt {
    my $text = shift;
    my $self;
    my $type;
    my $prefix;
    for my $t (qw/POINT MULTIPOINT LINESTRING MULTILINESTRING
    LINEARRING POLYGON POLYHEDRALSURFACE MULTIPOLYGON
    GEOMETRYCOLLECTION/) {
        next unless $text =~ /^\s*$t/i;
        $type = $t;
        $text =~ s/^\s*$t\s*//i;
        if ($text =~ /^EMPTY/i) {
            $text = undef;
        } else {
            $text =~ s/^([ZM]*)\s*\(\s*//i;
            $prefix = lc($1);
            $text =~ s/\s*\)\s*$//;
        }
        last;
    }
    if ($type eq 'POINT') {
        return Geo::OGC::Point->new unless $text;
        my @point = split /\s+/, $text;
        $self = Geo::OGC::Point->new("point$prefix"=>\@point);
    } elsif ($type eq 'MULTIPOINT') {
        return Geo::OGC::MultiPoint->new unless $text;
        $text =~ s/^\(\s*//;
        $text =~ s/\)\s*$//;
        my @points = split /\s*,\s*/, $text;
        $self = Geo::OGC::MultiPoint->new();
        for my $p (@points) {
            $self->AddGeometry(parse_wkt("POINT $prefix ($p)"));
        }
    } elsif ($type eq 'LINESTRING') {
        return Geo::OGC::LineString->new unless $text;
        $text =~ s/^\(\s*//;
        $text =~ s/\)\s*$//;
        my @points = split /\s*,\s*/, $text;
        $self = Geo::OGC::LineString->new();
        for my $p (@points) {
            $self->AddPoint(parse_wkt("POINT $prefix ($p)"));
        }
    } elsif ($type eq 'MULTILINESTRING') {
        return Geo::OGC::MultiLineString->new unless $text;
        $text =~ s/^\(\s*//;
        $text =~ s/\)\s*$//;
        my @strings = split /\)\s*,\s*\(/, $text;
        $self = Geo::OGC::MultiLineString->new();
        for my $s (@strings) {
            $self->AddGeometry(parse_wkt("LINESTRING $prefix ($s)"));
        }
    } elsif ($type eq 'LINEARRING') {
        return Geo::OGC::LinearRing->new unless $text;
        $text =~ s/^\(\s*//;
        $text =~ s/\)\s*$//;
        my @points = split /\s*,\s*/, $text;
        $self = Geo::OGC::LinearRing->new();
        for my $p (@points) {
            $self->AddPoint(parse_wkt("POINT $prefix ($p)"));
        }
    } elsif ($type eq 'POLYGON') {
        return Geo::OGC::Polygon->new unless $text;
        $text =~ s/^\(\s*//;
        $text =~ s/\)\s*$//;
        my @rings = split /\)\s*,\s*\(/, $text;
        $self = Geo::OGC::Polygon->new();
        $self->ExteriorRing(parse_wkt("LINEARRING $prefix (".shift(@rings).")"));
        for my $r (@rings) {
            $self->AddInteriorRing(parse_wkt("LINEARRING $prefix ($r)"));
        }
    } elsif ($type eq 'POLYHEDRALSURFACE') {
        return Geo::OGC::PolyhedralSurface->new unless $text;
        $text =~ s/^\(\s*//;
        $text =~ s/\)\s*$//;
        my @patches = split /\)\s*,\s*\(/, $text;
        $self = Geo::OGC::PolyhedralSurface->new();
        for my $p (@patches) {
            $self->AddPatch(parse_wkt("POLYGON $prefix ($p)"));
        }
    } elsif ($type eq 'MULTIPOLYGON') {
        return Geo::OGC::MultiPolygon->new unless $text;
        $text =~ s/^\(\s*\(\s*//;
        $text =~ s/\)\s*\)\s*$//;
        my @polygons = split /\)\s*\)\s*,\s*\(\s*\(/, $text;
        $self = Geo::OGC::MultiPolygon->new();
        for my $p (@polygons) {
            $self->AddGeometry(parse_wkt("POLYGON $prefix (($p))"));
        }
    } elsif ($type eq 'GEOMETRYCOLLECTION') {
        return Geo::OGC::GeometryCollection->new unless $text;
        my @b = $text =~ /,([A-Z])/g;
        unshift @b,'';
        my @geometries = split /,[A-Z]/, $text;
        $self = Geo::OGC::GeometryCollection->new();
        for my $i (0..$#geometries) {
            $self->AddGeometry(parse_wkt($b[$i].$geometries[$i]));
        }
    } else {
        my $b = substr $text, 0, 20;
        croak "can't parse text beginning as: $b";
    }
    return $self;
}

# counterclockwise from Sedgewick: Algorithms in C
sub ccw {
    my($x0, $y0, $x1, $y1, $x2, $y2) = @_;
    my $dx1 = $x1 - $x0;
    my $dy1 = $y1 - $y0;
    my $dx2 = $x2 - $x0;
    my $dy2 = $y2 - $y0;
    return +1 if $dx1*$dy2 > $dy1*$dx2;
    return -1 if $dx1*$dy2 < $dy1*$dx2;
    return -1 if ($dx1*$dx2 < 0) or ($dy1*$dy2 < 0);
    return +1 if ($dx1*$dx1+$dy1*$dy1) < ($dx2*$dx2+$dy2*$dy2);
    return 0;
}

# Test intersection of two lines from Sedgewick: Algorithms in C
sub intersect {
    my($x11, $y11, $x12, $y12, $x21, $y21, $x22, $y22) = @_;
    return ((ccw($x11, $y11, $x12, $y12, $x21, $y21)
	     *ccw($x11, $y11, $x12, $y12, $x22, $y22)) <= 0)
	&& ((ccw($x21, $y21, $x22, $y22, $x11, $y11)
	     *ccw($x21, $y21, $x22, $y22, $x12, $y12)) <= 0);
}

sub intersection_point {
    my($x11, $y11, $x12, $y12, $x21, $y21, $x22, $y22) = @_;
    my $dy1 = $y12 - $y11;
    my $dx1 = $x12 - $x11;
    my $dy2 = $y22 - $y21;
    my $dx2 = $x22 - $x21;
    # (dy1*dx2 - dy2*dx1)*x = dx1*dx2*(y21-y11) - dy2*dx1*x21 + dy1*dx2*x11
    # (dy1*dx1 - dy2*dx2)*y = dy1*dy2*(x21-x11) - dy1*dx2*y21 + dy2*dx1*y11
    my $x = ($dx1*$dx2*($y21-$y11) - $dy2*$dx1*$x21 + $dy1*$dx2*$x11)/($dy1*$dx2 - $dy2*$dx1);
    my $y = ($dy1*$dy2*($x21-$x11) - $dy1*$dx2*$y21 + $dy2*$dx1*$y11)/($dy1*$dx1 - $dy2*$dx2);
    return ($x, $y);
}

# Compute the distance of a point to a line.
sub distance_point_line_sqr {
    my($x, $y, $x1, $y1, $x2, $y2) = @_;
    my $dx = $x2-$x1;
    my $dy = $y2-$y1;
    my $l2 = $dx*$dx + $dy*$dy;
    my $u = (($x - $x1) * $dx + ($y - $y1) * $dy) / $l2;
    if ($u < 0) { # distance to point 1
	return (($x-$x1)*($x-$x1) + ($y-$y1)*($y-$y1));
    } elsif ($u > 1) { # distance to point 2
	return (($x-$x2)*($x-$x2) + ($y-$y2)*($y-$y2));
    } else {
	my $ix = $x1 + $u * $dx;
	my $iy = $y1 + $u * $dy;
	return (($x-$ix)*($x-$ix) + ($y-$iy)*($y-$iy));
    }
}

=pod

=item Clone()

Clones this object, i.e., creates a spatially exact copy.

=cut

sub Clone {
    my($self) = @_;
    my $clone = Geo::OGC::Geometry::new($self);
    $self->copy($clone);
    return $clone;
}

=pod

=item Precision($precision)

Sets or gets the precision.

Not in the specification.

=cut

sub Precision {
    my($self, $Precision) = @_;
    defined $Precision ?
	$self->{Precision} = $Precision-1 : $self->{Precision}+1;
}

=pod

=item Dimension()

The dimension (2 or 3) of this geometric object. In non-homogeneous
collections, this will return the largest topological dimension of the
contained objects.

=cut

sub Dimension {
    my($self) = @_;
    croak "Dimension method for class ".ref($self)." is not implemented yet";
}

=pod

=item GeometryType()

Returns the geometry type of this object.

=cut

sub GeometryType {
    my($self) = @_;
    croak "GeometryType method for class ".ref($self)." is not implemented yet";
}

=pod

=item SRID($srid)

Get or set Spatial Reference System ID for this object.

=cut

sub SRID {
    my($self, $SRID) = @_;
    defined $SRID ?
	$self->{SRID} = $SRID : $self->{SRID};
}

=pod

=item Envelope()

Compute the minimum bounding box for this geometry as a ring.

=cut

sub Envelope {
    my($self) = @_;
    croak "Envelope method for class ".ref($self)." is not implemented yet";
}

# A helper method used by AsText
sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    croak "as_text method for class ".ref($self)." is not implemented yet";
}

=pod

=item AsText()

This object in Well-known Text representation.

=cut

sub AsText {
    my($self) = @_;
    return uc($self->GeometryType).' '.'EMPTY' if $self->IsEmpty;
    return $self->as_text(1, 1);
}

=pod

=item AsBinary()

This object in Well-known Binary representation.

=cut

sub AsBinary {
    my($self) = @_;
    croak "AsBinary method for class ".ref($self)." is not implemented yet";
}

=pod

=item IsEmpty()

Returns true if this object is empty, i.e. contains no points.

=cut

sub IsEmpty {
    my($self) = @_;
    croak "IsEmpty method for class ".ref($self)." is not implemented yet";
}

=pod

=item IsSimple()

Returns true if this geometric object has no anomalous geometric
points, such as self intersection or self tangency.

=cut

sub IsSimple {
    my($self) = @_;
    croak "IsSimple method for class ".ref($self)." is not implemented yet";
}

=pod

=item Is3D()

Returns true if this geometric object has z coordinate values.

=cut

sub Is3D {
    my($self) = @_;
    croak "Is3D method for class ".ref($self)." is not implemented yet";
}

=pod

=item IsMeasured()

Returns true if this object has m coordinate values.

=cut

sub IsMeasured {
    my($self) = @_;
    croak "IsMeasured method for class ".ref($self)." is not implemented yet";
}

=pod

=item Boundary()

Returns the closure of the combinatorial boundary of this object.

=cut

sub Boundary {
    my($self) = @_;
    croak "Boundary method for class ".ref($self)." is not implemented yet";
}

=pod

=item Equals($geometry)

Returns true if this object is "spatially equal" to the other
geometry.

=cut

sub Equals {
    my($self, $geom) = @_;
    croak "Equals method for class ".ref($self)." is not implemented yet";
}

=pod

=item Disjoint($geometry)

Returns true if this object is "spatially disjoint" from the other geometry.

=cut

sub Disjoint {
    my($self, $geom) = @_;
    croak "Disjoint method for class ".ref($self)." is not implemented yet";
}

=pod

=item Intersects($geometry)

Returns true if this object "spatially intersects" the other geometry.

=cut

sub Intersects {
    my($self, $geom) = @_;
    croak "Intersects method for class ".ref($self)." is not implemented yet";
}

=pod

=item Touches($geometry)

Returns true if this object "spatially touches" the other geometry.

=cut

sub Touches {
    my($self, $geom) = @_;
    croak "Touches method for class ".ref($self)." is not implemented yet";
}

=pod

=item Crosses($geometry)

Returns true if this object "spatially crosses" the other geometry.

=cut

sub Crosses {
    my($self, $geom) = @_;
    croak "Crosses method for class ".ref($self)." is not implemented yet";
}

=pod

=item Within($geometry)

Returns true if this object is "spatially within" the other geometry.

=cut

sub Within {
    my($self, $geom) = @_;
    croak "Within method for class ".ref($self)." is not implemented yet";
}

=pod

=item Contains($geometry)

Returns true if this object "spatially contains" the other geometry.

=cut

sub Contains {
    my($self, $geom) = @_;
    croak "Contains method for class ".ref($self)." is not implemented yet";
}

=pod

=item Overlaps($geometry)

Returns true if this object "spatially overlaps" the other geometry.

=cut

sub Overlaps {
    my($self, $geom) = @_;
    croak "Overlaps method for class ".ref($self)." is not implemented yet";
}

=pod

=item Relate($geometry, $intersectionPatternMatrix)

Returns true if this object is spatially related to the other geometry
by testing for intersections between the interior, boundary and
exterior of the two geometric objects as specified by the values in
the intersectionPatternMatrix.  This returns FALSE if all the tested
intersections are empty except exterior (this) intersect exterior
(another).

=cut

sub Relate {
    my($self, $geom, $int_pat) = @_;
    croak "Relate method for class ".ref($self)." is not implemented yet";
}

=pod

=item LocateAlong($mValue)

Returns a derived geometry collection value that matches the
specified m coordinate value.

=cut

sub LocateAlong {
    my($self, $mValue) = @_;
    croak "LocateAlong method for class ".ref($self)." is not implemented yet";
}

=pod

=item LocateBetween($mStart, $mEnd)

Returns a derived geometry collection value that matches the specified
range of m coordinate values inclusively.

=cut

sub LocateBetween {
    my($self, $mStart, $mEnd) = @_;
    croak "LocateBetween method for class ".ref($self)." is not implemented yet";
}

=pod

=item Distance($geometry)

Returns the shortest distance between any two points between this object and the other geometry.

=cut

sub Distance {
    my($self, $geom) = @_;
    croak "Distance method for class ".ref($self)." is not implemented yet";
}

=pod

=item Buffer($distance)

Returns a surface whose points are closer than or at the given distance from this object.

=cut

sub Buffer {
    my($self, $distance) = @_;
    croak "Buffer method for class ".ref($self)." is not implemented yet";
}

=pod

=item ConvexHull()

Returns a convex hull of this object.

=cut

sub ConvexHull {
    my($self) = @_;
    croak "ConvexHull method for class ".ref($self)." is not implemented yet";
}

=pod

=item Intersection($geometry)

Returns a point set intersection of this object with the other geometry.

=cut

sub Intersection {
    my($self, $geom) = @_;
    croak "Intersection method for class ".ref($self)." is not implemented yet";
}

=pod

=item Union($geometry)

Returns a point set union of this object with the other geometry.

=cut

sub Union {
    my($self, $geom) = @_;
    croak "Union method for class ".ref($self)." is not implemented yet";
}

=pod

=item Difference($geometry)

Returns a point set difference of this object with the other geometry.

=cut

sub Difference {
    my($self, $geom) = @_;
    croak "Difference method for class ".ref($self)." is not implemented yet";
}

=pod

=item SymDifference($geometry)

Returns a point set symmetric difference of this object with the other geometry.

=cut

sub SymDifference {
    my($self, $geom) = @_;
    croak "SymDifference method for class ".ref($self)." is not implemented yet";
}

=pod

=item MakeCollection()

Creates a collection which contains this geometry.

=cut

sub MakeCollection {
    my($self) = @_;
    croak "MakeCollection method for class ".ref($self)." is not implemented";
}

=pod

=item ApplyTransformation($transformation)

transf = A point transformation method which will be applied for all
the points in the geometry as:

    ($new_x, $new_y, $new_z) = $transformation->($x, $y, $z)

Not in the specification.

=cut

sub ApplyTransformation {
    my($self, $transf) = @_;
    croak "ApplyTransformation method for class ".ref($self)." is not implemented";
}

=pod

=back

=head2 Geo::OGC::SpatialReferenceSystem

=cut

package Geo::OGC::SpatialReferenceSystem;

use strict;
use Carp;

sub new {
    my($package, %params) = @_;
    my $self = {};
    bless $self => (ref($package) or $package);
}

=pod

=head2 Geo::OGC::Point

A 0-dimensional geometric object.

=cut

package Geo::OGC::Point;

use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;

our @ISA = qw( Geo::OGC::Geometry );

=pod

=over

=item new()

A new point object can be created (besides using WKT) with

    $point = Geo::OGC::Point->new($x, $y);
    $point = Geo::OGC::Point->new($x, $y, $z);
    $point = Geo::OGC::Point->new(point => [$x, $y]);
    $point = Geo::OGC::Point->new(point => [$x, $y, $z]);
    $point = Geo::OGC::Point->new(point => [$x, $y, $z, $m]);
    $point = Geo::OGC::Point->new(pointz => [$x, $y, $z]);
    $point = Geo::OGC::Point->new(pointz => [$x, $y, $z, $m]);
    $point = Geo::OGC::Point->new(pointm => [$x, $y, $m]);
    $point = Geo::OGC::Point->new(pointm => [$x, $y, $z, $m]);
    $point = Geo::OGC::Point->new(pointzm => [$x, $y, $z, $m]);
    $point = Geo::OGC::Point->new(X => $x, Y => $y);
    $point = Geo::OGC::Point->new(X => $x, Y => $y, Z => $z);
    $point = Geo::OGC::Point->new(X => $x, Y => $y, Z => $z, M => $m);

=cut

sub new {
    my $package = shift;
    my %params;
    if (@_ == 2 and !($_[0] =~ /^[XYZMpP]/)) { # allow syntax Point->new($x, $y);
	$params{X} = shift;
	$params{Y} = shift;
    } elsif (@_ == 3) { # allow syntax Point->new($x, $y, $z);
	$params{X} = shift;
	$params{Y} = shift;
	$params{Z} = shift;
    } else {
	%params = @_;
    }
    # support comma as decimal point, and space in numbers
    for my $k (keys %params) {
	if (ref($params{$k})) {
	    for my $p (@{$params{$k}}) {
		$p =~ s/,/./g;
		$p =~ s/ //g;
		#print STDERR "point: $_\n";
	    }
	} else {
	    $params{$k} =~ s/,/./g;
	    $params{$k} =~ s/ //g;
	    #print STDERR "point: $_ => $params{$_}\n";
	}
    }
    my $self = Geo::OGC::Geometry::new($package, %params);
    return $self;
}

sub init {
    my($self, %params) = @_;
    $self->SUPER::init(%params);
    # +0 catches non-numeric error, if warnings are on
    if ($params{point}) {
	$self->{X} = $params{point}[0]+0;
	$self->{Y} = $params{point}[1]+0;
	$self->{Z} = $params{point}[2]+0 if @{$params{point}} > 2;
	$self->{M} = $params{point}[3]+0 if @{$params{point}} > 3;
    } elsif ($params{pointz}) {
	$self->{X} = $params{pointz}[0]+0;
	$self->{Y} = $params{pointz}[1]+0;
	$self->{Z} = $params{pointz}[2]+0;
	if (@{$params{pointz}} == 4) {
	    $self->{M} = $params{pointz}[3]+0;
	}
    } elsif ($params{pointm}) {
	$self->{X} = $params{pointm}[0]+0;
	$self->{Y} = $params{pointm}[1]+0;
	if (@{$params{pointm}} == 3) {
	    $self->{M} = $params{pointm}[2]+0;
	} elsif (@{$params{pointm}} == 4) {
	    $self->{Z} = $params{pointm}[2]+0;
	    $self->{M} = $params{pointm}[3]+0;
	}
    } elsif ($params{pointzm}) {
	$self->{X} = $params{pointm}[0]+0;
	$self->{Y} = $params{pointm}[1]+0;
	$self->{Z} = $params{pointm}[2]+0;
	$self->{M} = $params{pointm}[3]+0;
    } else {
	$self->{X} = $params{X}+0 if exists $params{X};
	$self->{Y} = $params{Y}+0 if exists $params{Y};
	$self->{Z} = $params{Z}+0 if exists $params{Z};
	$self->{M} = $params{M}+0 if exists $params{M};
    }
}

sub copy {
    my($self, $clone) = @_;
    $self->SUPER::copy($clone);
    for my $a (qw/X Y Z M/) {
	$clone->{$a} = $self->{$a} if exists $self->{$a};
    }
}

# Return a reference to an anonymous array that contains the point data.
# Note that there is no difference between [x,y,z] and [x,y,m]
sub point {
    my($self) = @_;
    my @point = ($self->{X}, $self->{Y});
    push @point, $self->{Z} if exists $self->{Z};
    push @point, $self->{M} if exists $self->{M};
    return [@point];
}

sub GeometryType {
    return 'Point';
}

sub Dimension {
    return 0;
}

sub Clone {
    my($self) = @_;
    my $m = exists $self->{M} ? 'm' : '';
    return Geo::OGC::Point::new($self, "point$m" => $self->point);
}

sub IsEmpty {
    my($self) = @_;
    return !(exists $self->{X});
}

# A point is always simple.
sub IsSimple {
    my($self) = @_;
    return 1;
}

sub Is3D {
    my($self) = @_;
    return exists $self->{Z};
}

sub IsMeasured {
    my($self) = @_;
    return exists $self->{M};
}

sub Boundary {
    my($self) = @_;
    return $self->Clone;
}

sub X {
    my($self, $X) = @_;
    defined $X ?
	$self->{X} = $X : $self->{X};
}

sub Y {
    my($self, $Y) = @_;
    defined $Y ?
	$self->{Y} = $Y : $self->{Y};
}

=pod

=item Z()

Get or set the z coordinate.

Setting is not in the specification.

=cut

sub Z {
    my($self, $Z) = @_;
    defined $Z ?
	$self->{Z} = $Z : (exists $self->{Z} ? $self->{Z} : undef);
}

=pod

=item M()

Get or set the measure.

Setting is not in the specification.

=cut

sub M {
    my($self, $M) = @_;
    defined $M ?
	$self->{M} = $M : (exists $self->{M} ? $self->{M} : undef);
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my @coords;
    my $ZM = exists $self->{Z} ? 'Z' : '';
    if (exists $self->{Precision}) {
	for my $a (qw/X Y Z/) {
	    last unless exists $self->{$a};
	    my $s = sprintf('%.'.$self->{Precision}.'e', $self->{$a});
	    push @coords, $s;
	}
    } else {
	for my $a (qw/X Y Z/) {
	    push @coords, $self->{$a} if exists $self->{$a};
	}
    }
    if (exists $self->{M}) {
	push @coords, $self->{M};
	$ZM .= 'M';
    }
    my $text = join(' ', @coords);
    $text =~ s/,/./g; # setting POSIX numeric locale does not seem to work??
    $text = '('.$text.')' if $force_parens;
    $text = "POINT$ZM ".$text if $include_tag;
    return $text;
}

# what should we do with z?
sub Equals {
    my($self, $geom) = @_;
    return 0 unless $geom->isa('Geo::OGC::Point');
    if (exists $self->{Precision}) {
	for my $a (qw/X Y Z/) {
	    last unless exists $self->{$a} and exists $geom->{$a};
	    my $s = sprintf('%.'.$self->{Precision}.'e', $self->{$a});
	    my $g = sprintf('%.'.$self->{Precision}.'e', $geom->{$a});
	    return 0 if $s != $g;
	}
	return 1;
    }
    return (($self->{X}-$geom->{X})**2+($self->{Y}-$geom->{Y})**2) < $SNAP_DISTANCE_SQR;
}

sub DistanceToLineStringSqr {
    my($self, $linestring) = @_;
    my($x, $y) = ($self->{X}, $self->{Y});
    my $p1 = $linestring->{Points}[0];
    return unless $p1;
    my $p2 = $linestring->{Points}[1];
    return (($x-$p1->{X})**2+($y-$p1->{Y})**2) unless $p2;
    my $distance = distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y});
    for my $i (2..$#{$linestring->{Points}}) {
	$p1 = $p2;
	$p2 = $linestring->{Points}[$i];
	my $d = distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y});
	$distance = $d if $d < $distance;
    }
    return $distance;
}

sub Distance {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return sqrt(($self->{X}-$geom->{X})**2 + ($self->{Y}-$geom->{Y})**2);
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	return sqrt($self->DistanceToLineStringSqr($geom));
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	if ($geom->{ExteriorRing}->IsPointIn($self)) {
	    for my $ring (@{$geom->{InteriorRings}}) {
		return sqrt($self->DistanceToLineStringSqr($ring)) if $ring->IsPointIn($self);
	    }
	    return 0;
	} else {
	    return sqrt($self->DistanceToLineStringSqr($geom->{ExteriorRing}));
	}
    } elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
	my $dist = $self->Distance($geom->{Geometries}[0]);
	for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
	    my $d = $self->Distance($g);
	    $dist = $d if $d < $dist;
	}
	return $dist;
    } else {
	croak "can't compute distance between a ".ref($geom)." and a point";
    }
}

# what should this be?
sub Envelope {
    my($self) = @_;
    my $r = Geo::OGC::LinearRing->new;
    $r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
    $r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
    $r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
    $r->AddPoint(Geo::OGC::Point->new($self->{X}, $self->{Y}));
    $r->Close;
    return $r;
}

sub Area {
    return 0;
}

sub Intersection {
    my($self, $geom) = @_;
    return $self->Clone if $self->Within($geom);
    return undef;
}

sub Within {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return $self->Equals($geom) ? 1 : 0;
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	return $self->DistanceToLineStringSqr($geom) < $SNAP_DISTANCE_SQR ? 1 : 0;
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	if (!($geom->{ExteriorRing}->IsPointStricktlyOut($self))) {
	    for my $ring (@{$geom->{InteriorRing}}) {
		return 0 if $ring->IsPointStricktlyIn($self);
	    }
	    return 1;
	}
	return 0;
    } elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
	for my $g (@{$geom->{Geometries}}) {
	    return 1 if $self->Within($g);
	}
	return 0;
    } else {
	croak "point within ".ref($geom)." is not implemented yet";
    }
}

sub MakeCollection {
    my($self) = @_;
    my $coll = Geo::OGC::MultiPoint->new;
    $coll->AddGeometry($self);
    return $coll;
}

sub ApplyTransformation {
    my($self, $transf) = @_;
    if (@_ > 2) {
	($self->{X}, $self->{Y}, $self->{Z}) = $transf->($self->{X}, $self->{Y}, $self->{Z});
    } else {
	($self->{X}, $self->{Y}) = $transf->($self->{X}, $self->{Y});
    }
}

sub ClosestVertex {
    my($self, $x, $y) = @_;
    return (($self->{X}-$x)**2 + ($self->{Y}-$y)**2);
}

sub VertexAt {
    my $self = shift;
    return ($self);
}

sub ClosestPoint {
    my($self, $x, $y) = @_;
    return (($self->{X}-$x)**2 + ($self->{Y}-$y)**2);
}

sub AddVertex {
}

sub DeleteVertex {
}

=pod

=back

=head2 Geo::OGC::Curve

A 1-dimensional geometric object.

=cut

package Geo::OGC::Curve;

use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;

our @ISA = qw( Geo::OGC::Geometry );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Geometry::new($package, %params);
    return $self;
}

sub init {
    my($self, %params) = @_;
    $self->SUPER::init(%params);
    $self->{Points} = [];
    if ($params{points}) {
	for my $p (@{$params{points}}) {
	    $self->AddPoint(Geo::OGC::Point->new(point=>$p));
	}
    } elsif ($params{pointsm}) {
	for my $p (@{$params{pointsm}}) {
	    $self->AddPoint(Geo::OGC::Point->new(pointm=>$p));
	}
    }
}

sub copy {
    my($self, $clone) = @_;
    $self->SUPER::copy($clone);
    for my $p (@{$self->{Points}}) {
	$clone->AddPoint($p->Clone);
    }
}

sub IsEmpty {
    my($self) = @_;
    return 1 if !$self->{Points} or @{$self->{Points}} == 0;
    return 0;
}

sub GeometryType {
    return 'Curve';
}

sub Dimension {
    return 1;
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text = join(',', map {$_->as_text} @{$self->{Points}});
    $text = '('.$text.')';
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' '.$text if $include_tag;
    return $text;
}

=pod

=over

=item AddPoint($point, $i)

Adds a point to the end (N+1) of the curve by default.

point = A Point object
i [optional] = The location in the sequence (1..N+1) where to add the Point.

Not in the specification.

=cut

sub AddPoint {
    my($self, $point, $i) = @_;
    croak 'usage: Curve->AddPoint($point) '
	unless $point and $point->isa('Geo::OGC::Point');
    my $points = $self->{Points};
    if (defined $i) {
	my $temp = $points->[$i-1];
	splice @$points, $i-1, 1, ($point, $temp);
    } else {
	push @$points, $point;
    }
}

=pod

=item DeletePoint($i)

Delete a point from the curve.

i = The location in the sequence (1..N) from where to delete the Point.

Not in the specification.

=cut

sub DeletePoint {
    my($self, $i) = @_;
    my $points = $self->{Points};
    splice @$points, $i-1, 1;
}

sub StartPoint {
    my($self) = @_;
    my $points = $self->{Points};
    return $points->[0] if @$points;
}

sub EndPoint {
    my($self) = @_;
    my $points = $self->{Points};
    return $points->[$#$points] if @$points;
}

=pod

=item NumPoints()

Return the number of points in the sequence.

=cut

sub NumPoints {
    my($self) = @_;
    scalar(@{$self->{Points}});
}

=pod

=item PointN($N, $point)

Get or set the point of the sequence.

N = A location in the sequence.
point [optional] = A Point object, if defined sets the point to index N.

The first point has the index 1 as OGC SF SQL conformance test uses 1-based indexing.

=cut

sub PointN {
    my($self, $N, $point) = @_;
    my $points = $self->{Points};
    $points->[$N-1] = $point if defined $point;
    return $points->[$N-1];
}

sub Is3D {
    my($self) = @_;
    for my $p (@{$self->{Points}}) {
	return 1 if $p->Is3D;
    }
    return 0;
}

sub IsMeasured {
    my($self) = @_;
    for my $p (@{$self->{Points}}) {
	return 1 if $p->IsMeasured;
    }
    return 0;
}

sub IsClosed {
    my($self) = @_;
    $self->StartPoint()->Equals($self->EndPoint());
}

=pod

=item Close()

Close the curve by adding the first point also as the last point.

Not in the specification.

=cut

sub Close {
    my($self) = @_;
    push @{$self->{Points}}, $self->{Points}[0];
}

=pod

=item IsRing($upgrade)

Tests whether this curve is a ring, i.e., closed and simple.

upgrade [optional, not in the specification] = Upgrades this curve into a Ring if this really could be a ring.

=cut

sub IsRing {
    my($self, $upgrade) = @_;
    my $ret = ($self->IsClosed and $self->IsSimple);
    bless($self, 'Geo::OGC::LinearRing') if $ret and $upgrade;
    return $ret;
}

# should use Precision if one exists
sub Equals {
    my($self, $geom) = @_;
    return 0 unless $geom->isa('Geo::OGC::Curve');
    return 0 unless $#{$self->{Points}} == $#{$geom->{Points}};
    for my $i (0..$#{$self->{Points}}) {
	return 0 unless $self->{Points}[$i]->Equals($geom->{Points}[$i]);
    }
    return 1;
}

sub Area {
    return 0;
}

sub MakeCollection {
    my($self) = @_;
    my $coll = Geo::OGC::MultiCurve->new;
    $coll->AddGeometry($self);
    return $coll;
}

sub ApplyTransformation {
    my($self, $transf) = @_;
    for my $p (@{$self->{Points}}) {
	$p->ApplyTransformation($transf);
    }
}

=pod

=item Reverse()

Reverse the order of the points in the sequence.

Not in the specification.

=cut

sub Reverse {
    my($self) = @_;
    @{$self->{Points}} = reverse @{$self->{Points}};
}

sub ClosestVertex {
    my($self, $x, $y) = @_;
    return unless @{$self->{Points}};
    my($dmin) = $self->{Points}[0]->ClosestVertex($x, $y);
    my $i = 0;
    for my $j (1..$#{$self->{Points}}) {
	my($d) = $self->{Points}[$j]->ClosestVertex($x, $y);
	($i, $dmin) = ($j, $d) if $d < $dmin;
    }
    return ($i, $dmin);
}

sub VertexAt {
    my($self, $i) = @_;
    return ($self->{Points}[0], $self->{Points}[$#{$self->{Points}}])
	if (($i == 0 or $i == $#{$self->{Points}}) and $self->isa('Geo::OGC::LinearRing'));
    return ($self->{Points}[$i]);
}

sub _closest_point {
    my($x0, $y0, $x1, $y1, $x, $y) = @_;
    my $ab2 = ($x1-$x0)*($x1-$x0) + ($y1-$y0)*($y1-$y0);
    my $ap_ab = ($x-$x0)*($x1-$x0) + ($y-$y0)*($y1-$y0);
    my $t = $ap_ab/$ab2;
    if ($t < 0) {$t = 0} elsif ($t > 1) {$t = 1}
    my $xp = $x0+$t*($x1-$x0);
    my $yp = $y0+$t*($y1-$y0);
    return ($xp, $yp, ($x-$xp)*($x-$xp)+($y-$yp)*($y-$yp));
}

sub ClosestPoint {
    my($self, $x, $y) = @_;
    return unless @{$self->{Points}};
    my($i, $pmin, $dmin);
    for my $j (1..$#{$self->{Points}}) {
	my($xp, $yp, $d) =
	    _closest_point($self->{Points}[$j-1]{X}, $self->{Points}[$j-1]{Y},
			   $self->{Points}[$j]{X}, $self->{Points}[$j]{Y}, $x, $y);
	($i, $pmin, $dmin) = ($j, Geo::OGC::Point->new($xp, $yp), $d)
	    if (!defined($dmin) or $d  < $dmin);
    }
    return ($i, $pmin, $dmin)
}

sub AddVertex {
    my($self, $i, $p) = @_;
    splice @{$self->{Points}}, $i, 0, $p;
}

sub DeleteVertex {
    my($self, $i) = @_;
    splice @{$self->{Points}}, $i, 1;
}

=pod

=back

=head2 Geo::OGC::LineString

=cut

package Geo::OGC::LineString;

use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;

our @ISA = qw( Geo::OGC::Curve );

# de Berg et al p. 25
sub FindIntersections {
    # @_ contains a list of linestrings that make up the S
    my @linestrings = @_;
    #my $precision =

    # event queue
    my @Q = (); # [ymax,index1,index2], ..., index1 is index to @_ index2 is index to Points
    for my $index1 (0 .. $#linestrings) {
	my $s = $linestrings[$index1]->{Points};
	for my $index2 (0 .. $#$s-1) {
	    my($y1, $y2) = ( $s->[$index2]{Y}, $s->[$index2+1]{Y} );
	    if ($y1 > $y2) {
		push @Q, [$y1, $index1, $index2];
		push @Q, [$y2];
	    } else {
		push @Q, [$y2, $index1, $index2];
		push @Q, [$y1];
	    }
	}
    }

    # process event points in descending ymax order
    @Q = sort {$b->[0] <=> $a->[0]} @Q;

    my $T = Tree::Binary::Search->new(); # the status structure
    $T->setComparisonFunction
	(sub {
	    my($a, $b) = @_; # two keys

	});

    my $i = 0;
    while ($i < @Q) {
	my $j = $i+1;
	while ($j < @Q and sqrt(($Q[$i][0]-$Q[$j][0])**2) < $SNAP_DISTANCE_SQR) {
	    $j++;
	}
	# $i .. $j-1 are the event points in @Q
	for my $k ($i .. $j-1) {
	    if (@{$Q[$k]} == 3) {
		my $i1 = $Q[$k][1];
		my $i2 = $Q[$k][2];
		my $s = $linestrings[$i1]->{Points};
		$T->insert( "$i1,$i2" => [$s->[$i2]{X}, $s->[$i2]{Y}, $s->[$i2+1]{X}, $s->[$i2+1]{Y}] );
	    }
	}
    }
}

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Curve::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'LineString';
}

sub IsSimple {
    my($self) = @_;
    my $edges = @{$self->{Points}} - 1;
    return 1 if $edges < 2;
    my $closed = $self->IsClosed;
    my $simple = 1;
    for my $i (0..$edges-1) {
	# check a case of self tangency
	return 0 if $i < $edges-1 and $self->{Points}[$i+2]->Equals($self->{Points}[$i]);
	for my $j ($i+2..$edges-1) {
	    next if $closed and $i == 0 and $j == $edges-1;
	    return 0 if intersect
		(
		 $self->{Points}[$i]{X}, $self->{Points}[$i]{Y},
		 $self->{Points}[$i+1]{X}, $self->{Points}[$i+1]{Y},
		 $self->{Points}[$j]{X}, $self->{Points}[$j]{Y},
		 $self->{Points}[$j+1]{X},$self->{Points}[$j+1]{Y}
		 );
	}
    }
    return 1;
}

sub Envelope {
    my($self) = @_;
    my($minx, $miny, $maxx, $maxy);
    for my $p (@{$self->{Points}}) {
	$minx = $p->{X} if !defined($minx) or $minx > $p->{X};
	$miny = $p->{Y} if !defined($miny) or $miny > $p->{Y};
	$maxx = $p->{X} if !defined($maxx) or $maxx > $p->{X};
	$maxy = $p->{Y} if !defined($maxy) or $maxy > $p->{Y};
    }
    my $r = Geo::OGC::LinearRing->new;
    $r->AddPoint(Geo::OGC::Point->new($minx, $miny));
    $r->AddPoint(Geo::OGC::Point->new($maxx, $miny));
    $r->AddPoint(Geo::OGC::Point->new($maxx, $maxy));
    $r->AddPoint(Geo::OGC::Point->new($minx, $maxy));
    $r->Close;
    return $r;
}

=pod

=over

=item Length()

The length of this LineString in its associated spatial reference.

Currently computed as a simple euclidean distance.

=cut

sub Length {
    my($self) = @_;
    my $l = 0;
    my($x0, $y0) = ($self->{Points}[0]{X}, $self->{Points}[0]{Y});
    for my $i (1..$#{$self->{Points}}) {
	my($x1, $y1) = ($self->{Points}[$i]{X}, $self->{Points}[$i]{Y});
	$l += sqrt(($x1 - $x0)**2+($y1 - $y0)**2);
	($x0, $y0) = ($x1, $y1);
    }
    return $l;
}

sub Distance {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return $geom->DistanceToLineStringSqr($self);
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	my $dist;
	for my $p (@{$self->{Points}}) {
	    my $d = $p->DistanceToLineStringSqr($geom);
	    $dist = $d if !(defined $dist) or $d < $dist;
	}
	return $dist;
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	my $dist;
	for my $p (@{$self->{Points}}) {
	    my $d = $p->Distance($geom);
	    $dist = $d if !(defined $dist) or $d < $dist;
	}
	return $dist;
    } elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
	my $dist = $self->Distance($geom->{Geometries}[0]);
	for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
	    my $d = $self->Distance($g);
	    $dist = $d if $d < $dist;
	}
	return $dist;
    } else {
	croak "can't compute distance between a ".ref($geom)." and a line string";
    }
}

sub LinesWhereWithin {
    my($self, $point) = @_;
    my($x, $y) = ($point->{X}, $point->{Y});
    my @ret;
    my $p1 = $self->{Points}[0];
    return @ret unless $p1;
    my $p2 = $self->{Points}[1];
    return @ret unless $p1;
    push @ret, 1 if
	distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y}) < $SNAP_DISTANCE_SQR;
    for my $i (2..$#{$self->{Points}}) {
	$p1 = $p2;
	$p2 = $self->{Points}[$i];
	push @ret, $i if
	    distance_point_line_sqr($x, $y, $p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y}) < $SNAP_DISTANCE_SQR;
    }
    return @ret;
}

sub Within {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	for my $p (@{$self->{Points}}) {
	    return 0 unless $p->Equals($geom);
	}
	return 1;
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	my @w1 = ();
	for my $p (@{$self->{Points}}) {
	    my @w2 = $geom->LinesWhereWithin($p);
	    return 0 unless @w2;
	    next unless @w1;
	    my $overlap = 0;
	    for my $w1 (@w1) {
		for my $w2 (@w2) {
		    $overlap = 1 if $w1 == $w2;
		}
	    }
	    return 0 unless $overlap;
	    @w1 = @w2;
	}
	return 1;
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	for my $p (@{$self->{Points}}) {
	    return 0 if $geom->{ExteriorRing}->IsPointStricktlyOut($p);
	    for my $ring (@{$geom->{InteriorRings}}) {
		return 0 if $ring->IsPointStricktlyIn($p);
	    }
	}
	my $i = $self->Intersection($geom->{ExteriorRing});
	for my $g (@{$i->{Geometries}}) {
	    next unless $g->isa('Geo::OGC::Line');
	    # does the line go out of the polygon?
	    # yes, if its start and end points are on different lines
	    my @s = $geom->{ExteriorRing}->LinesWhereWithin($g->StartPoint);
	    my @e = $geom->{ExteriorRing}->LinesWhereWithin($g->EndPoint);
	    my $overlap = 0;
	    for my $s (@s) {
		for my $e (@e) {
		    $overlap = 1 if $s == $e;
		}
	    }
	    return 0 unless $overlap;
	}
	for my $ring (@{$geom->{InteriorRings}}) {
	    my $i = $self->Intersection($ring);
	    for my $g (@{$i->{Geometries}}) {
		next unless $g->isa('Geo::OGC::Line');
		# does the line go into the interior ring?
		# yes, if its start and end points are on different lines
		my @s = $ring->LinesWhereWithin($g->StartPoint);
		my @e = $ring->LinesWhereWithin($g->EndPoint);
		my $overlap = 0;
		for my $s (@s) {
		    for my $e (@e) {
			$overlap = 1 if $s == $e;
		    }
		}
		return 0 unless $overlap;
	    }
	}
	return 1;
    } else {
	croak "linestring within ".ref($geom)." is not yet implemented";
    }
}

# assuming simple lines
# z and m coords!
sub Intersection {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return $geom->Clone if $geom->DistanceToLineStringSqr($self) < $SNAP_DISTANCE_SQR;
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	#my $i = Geo::OGC::GeometryCollection->new;
	my %i;
	my $index = 1;
	my $p1;
	for my $p2 (@{$self->{Points}}) {
	    $p1 = $p2, next unless $p1;
	    my $q1;
	    for my $q2 (@{$geom->{Points}}) {
		$q1 = $q2, next unless $q1;
		my @p = ($p1->{X}, $p1->{Y}, $p2->{X}, $p2->{Y});
		my @q = ($q1->{X}, $q1->{Y}, $q2->{X}, $q2->{Y});
		#print STDERR "lines: @p; @q\n";
		if (intersect( @p, @q )) {
		    my $p1q = distance_point_line_sqr(@p[0..1], @q);
		    my $p2q = distance_point_line_sqr(@p[2..3], @q);
		    my $q1p = distance_point_line_sqr(@q[0..1], @p);
		    my $q2p = distance_point_line_sqr(@q[2..3], @p);
		    #print STDERR "p1=q $p1q p2=q $p2q q1=p $q1p q2=p $q2p\n";
		    if ($p1q < $SNAP_DISTANCE_SQR) {
			if ($q1p < $SNAP_DISTANCE_SQR) {
			    if ($p1->Equals($q1)) {
				$i{$index++} = $p1->Clone;
			    } else {
				$i{$index++} = Geo::OGC::Line->new( points=>[[@p[0..1]],[@q[0..1]]] );
			    }
			}
			if ($q2p < $SNAP_DISTANCE_SQR) {
			    if ($p1->Equals($q2)) {
				$i{$index++} = $p1->Clone;
			    } else {
				$i{$index++} = Geo::OGC::Line->new( points=>[[@p[0..1]],[@q[2..3]]] );
			    }
			}
			if ($p2q < $SNAP_DISTANCE_SQR) {
			    $i{$index++} = Geo::OGC::Line->new( points=>[[@p[0..1]],[@p[2..3]]] );
			} else {
			    $i{$index++} = Geo::OGC::Point->new( @p[0..1] );
			}
		    } elsif ($p2q < $SNAP_DISTANCE_SQR) {
			if ($q1p < $SNAP_DISTANCE_SQR) {
			    if ($p2->Equals($q1)) {
				$i{$index++} = $p2->Clone;
			    } else {
				$i{$index++} = Geo::OGC::Line->new( points=>[[@p[2..3]],[@q[0..1]]] );
			    }
			}
			if ($q2p < $SNAP_DISTANCE_SQR) {
			    if ($p2->Equals($q2)) {
				$i{$index++} = $p2->Clone;
			    } else {
				$i{$index++} = Geo::OGC::Line->new( points=>[[@p[2..3]],[@q[2..3]]] );
			    }
			}
			$i{$index++} = Geo::OGC::Point->new( @p[2..3] );
		    } elsif ($q1p < $SNAP_DISTANCE_SQR) {
			if ($q2p < $SNAP_DISTANCE_SQR) {
			    $i{$index++} = Geo::OGC::Line->new( points=>[[@q[0..1]],[@q[2..3]]] );
			} else {
			    $i{$index++} = Geo::OGC::Point->new( @q[0..1] );
			}
		    } elsif ($q2p < $SNAP_DISTANCE_SQR) {
			$i{$index++} = Geo::OGC::Point->new( @q[2..3] );
		    } else {
			$i{$index++} = Geo::OGC::Point->new( Geo::OGC::Geometry::intersection_point(@p, @q) );
		    }
		}
		$q1 = $q2;
	    }
	    $p1 = $p2;
	}
	#$i->Simplify;

	# delete unnecessary points and lines
	# comparisons are done unnecessarily twice??
	for my $g1 (keys %i) {
	    for my $g2 (keys %i) {
		next if $g1 == $g2;
		if ($i{$g1}->Within($i{$g2})) {
		    #print STDERR "delete ",$i{$g1}->AsText," because it is in ",$i{$g2}->AsText,"\n";
		    delete $i{$g1};
		    last;
		}
	    }
	}

	my $i = Geo::OGC::GeometryCollection->new;
	for my $g1 (keys %i) {
	    $i->AddGeometry($i{$g1});
	}
	return $i;
    } else {
	croak "intersection between a ".ref($geom)." and a line string is not yet implemented";
    }
}

sub MakeCollection {
    my($self) = @_;
    my $coll = Geo::OGC::MultiLineString->new;
    $coll->AddGeometry($self);
    return $coll;
}

# from http://everything2.com/index.pl?node_id=859282
sub pt_to_seg_dist {
    # distance of p to segment p1-p2, v12 is the vector p1p2
    my ($p1, $v12, $p) = @_;

    my $m12 = $v12->[0] * $v12->[0] + $v12->[1] * $v12->[1];
    my $v1p = [];
    $v1p->[0] = $p->{X} - $p1->{X};
    $v1p->[1] = $p->{Y} - $p1->{Y};
    my $dot = $v1p->[0] * $v12->[0] + $v1p->[1] * $v12->[1];
    if ($dot <= 0.0)
    {
	return sqrt ($v1p->[0] * $v1p->[0] + $v1p->[1] * $v1p->[1]);
    }
    else
    {
	if ($dot >= $m12)
	{
	    $v1p->[0] = $v1p->[0] + $v12->[0];
	    $v1p->[1] = $v1p->[1] + $v12->[1];
	    return sqrt ($v1p->[0] * $v1p->[0] + $v1p->[1] * $v1p->[1]);
	}
	else
	{
	    my $slash = $v1p->[0] * $v12->[1] - $v1p->[1] * $v12->[0];
	    return abs ($slash / sqrt ($m12));
	}
    }
}

sub simplify_part {
    my($self, $first, $last, $simple, $tolerance) = @_;
    if ($last > $first + 1)
    {
	my $p1 = $self->{Points}[$first];
	my $vfl = [$self->{Points}[$last]{X} - $self->{Points}[$first]{X},
		   $self->{Points}[$last]{Y} - $self->{Points}[$first]{Y}];
	# find the intermediate point
	# furthest from the segment
	# connecting first and last
	my $b = $first+1;
	my $db = pt_to_seg_dist ($p1, $vfl, $self->{Points}[$b]);
	my $i = $b + 1;
	while ($i < $last)
	{
	    my $di = pt_to_seg_dist ($p1, $vfl, $self->{Points}[$i]);
	    if ($di > $db)
	    {
		$b = $i;
		$db = $di;
	    }
	    $i++;
	}
	# if the furthest distance beats the tolerance,
	# recursively simplify the rest of the array.
	if ($db >= $tolerance)
	{
	    simplify_part ($self, $first, $b, $simple, $tolerance);
	    $simple->AddPoint($self->{Points}[$b]);
	    simplify_part ($self, $b, $last, $simple, $tolerance);
	}
    }
}

# Simplifies the linestring using Douglas-Peucker
sub simplify {
    my($self, $tolerance) = @_;
    my $simple = Geo::OGC::LineString->new;
    $simple->AddPoint($self->StartPoint);
    simplify_part ($self, 0, $self->NumPoints-1, $simple, $tolerance);
    $simple->AddPoint($self->EndPoint);
    return $simple;
}

=pod

=back

=head2 Geo::OGC::Line

=cut

package Geo::OGC::Line;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::LineString );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::LineString::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'Line';
}

=pod

=head2 Geo::OGC::LinearRing

=cut

package Geo::OGC::LinearRing;

use strict;
use Carp;
use Geo::OGC::Geometry qw/:all/;

our @ISA = qw( Geo::OGC::LineString );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::LineString::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'LinearRing';
}

=pod

=over

=item IsPointIn($point)

Tests whether the given point is within the ring.

Uses the pnpoly algorithm from L.

Assumes a simple closed ring. May or may not return true if the point is on the border.

Not in the specification.

=cut

sub IsPointIn {
    my($self, $point) = @_;
    my($x, $y) = ($point->{X}, $point->{Y});
    my $c = 0;
    my $prev;
    for my $p (@{$self->{Points}}) {
	$prev = $p, next unless $prev;
	$c = !$c if (((( $p->{Y} <= $y ) && ( $y < $prev->{Y} )) ||
		      (( $prev->{Y} <= $y ) && ( $y < $p->{Y} ))) &&
		     ( $x < ( $prev->{X} - $p->{X} ) *
		       ( $y - $p->{Y} ) / ( $prev->{Y} - $p->{Y} ) + $p->{X} ));
	$prev = $p;
    }
    return $c;
}

# @note not on the border
sub IsPointStricktlyIn {
    my($self, $point) = @_;
    return 1 if ( $self->IsPointIn($point) and
		  !($point->DistanceToLineStringSqr($self) < $SNAP_DISTANCE_SQR) );
    return 0;
}

# @note not on the border
sub IsPointStricktlyOut {
    my($self, $point) = @_;
    return 1 unless ( $self->IsPointIn($point) or
		      $point->DistanceToLineStringSqr($self) < $SNAP_DISTANCE_SQR );
    return 0;
}

=pod

=item Area()

Compute the area of the ring. The area is computed as a simple euclidean area.

Not in the specification.

Assumes a simple closed ring

Returns the area as a negative number if the sense of the rotation of the ring is clockwise.

=cut

sub Area {
    my($self) = @_;
    my $area = 0;
    my $N = $#{$self->{Points}}-1; # skip the closing point
    my $j = 0;
    for my $i (0..$N) {
	$j++;
	$area += $self->{Points}[$i]{X} * $self->{Points}[$j]{Y};
	$area -= $self->{Points}[$i]{Y} * $self->{Points}[$j]{X};
    }
    return $area/2;
}

sub Centroid {
    my($self) = @_;
    my($area, $x, $y) = (0, 0, 0);
    my $N = $#{$self->{Points}}-1; # skip the closing point
    for my $i (0..$N) {
	my($xi, $yi, $xj, $yj) = ( $self->{Points}[$i]{X}, $self->{Points}[$i]{Y},
				   $self->{Points}[$i+1]{X}, $self->{Points}[$i+1]{Y} );
	my $b = $xi * $yj - $xj * $yi;
	$area += $b;
	$x += ($xi + $xj) * $b;
	$y += ($yi + $yj) * $b;
    }
    $x /= abs(3*$area); # 6 but $area is 2 * real area
    $y /= abs(3*$area);
    return Geo::OGC::Point->new($x, $y);
}

# Returns true if the points in this ring are arranged counterclockwise. Assumes a simple closed ring.
sub IsCCW {
    my($self) = @_;
    # find the northernmost point
    my $t = 0;
    my $N = $#{$self->{Points}}-1; # skip the closing point
    for my $i (1..$N) {
	$t = $i if $self->{Points}[$i]{Y} > $self->{Points}[$t]{Y};
    }
    # the previous point
    my $p = $t-1;
    $p = $N if $p < 0;
    # the next point
    my $n = $t+1;
    $n = 0 if $n > $N;
    return ccw($self->{Points}[$p]{X}, $self->{Points}[$p]{Y},
	       $self->{Points}[$t]{X}, $self->{Points}[$t]{Y},
	       $self->{Points}[$n]{X}, $self->{Points}[$n]{Y}) == 1;
}

# Makes clockwise from counterclockwise and vice versa.
sub Rotate {
    my($self) = @_;
    @{$self->{Points}} = reverse @{$self->{Points}};
}

=pod

=back

=head2 Geo::OGC::Surface

A 2-dimensional geometric object.

=cut

package Geo::OGC::Surface;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::Geometry );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Geometry::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'Surface';
}

sub Dimension {
    return 2;
}

sub Area {
    my($self) = @_;
    croak "Area method for class ".ref($self)." is not implemented yet";
}

sub Centroid {
    my($self) = @_;
    croak "Centroid method for class ".ref($self)." is not implemented yet";
}

sub PointOnSurface {
    my($self) = @_;
    croak "PointOnSurface method for class ".ref($self)." is not implemented yet";
}

sub MakeCollection {
    my($self) = @_;
    my $coll = Geo::OGC::MultiSurface->new;
    $coll->AddGeometry($self);
    return $coll;
}

=pod

=head2 Geo::OGC::Polygon

=cut

package Geo::OGC::Polygon;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::Surface );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Surface::new($package, %params);
    return $self;
}

sub init {
    my($self, %params) = @_;
    $self->SUPER::init(%params);
    $self->{ExteriorRing} = undef;
    $self->{InteriorRings} = [];
}

sub copy {
    my($self, $clone) = @_;
    $self->SUPER::copy($clone);
    $clone->ExteriorRing($self->{ExteriorRing}->Clone) if $self->{ExteriorRing};
    for my $r (@{$self->{InteriorRings}}) {
	$clone->AddInteriorRing($r->Clone);
    }
}

sub IsEmpty {
    my($self) = @_;
    return 1 unless $self->{ExteriorRing};
    return undef;
}

sub GeometryType {
    return 'Polygon';
}

# Test the rules that define valid polygons.
sub Assert {
    my($self) = @_;

    # a) Polygons are topologically closed;
    croak "not at least triangle" unless @{$self->{ExteriorRing}->{Points}} > 3;
    croak "exterior not closed" unless $self->{ExteriorRing}->IsClosed;

    # b) The boundary of a Polygon consists of a set of LinearRings
    # that make up its exterior and interior boundaries;
    for my $r (@{$self->{InteriorRings}}) {
	croak "an interior is not closed" unless $r->IsClosed;
    }

    # c) No two Rings in the boundary cross and the Rings in the
    # boundary of a Polygon may intersect at a Point but only as a
    # tangent
    for my $ring (@{$self->{InteriorRings}}) {
	for my $p (@{$ring->{Points}}) {
	    croak "point in interior not within exterior" unless $self->{ExteriorRing}->IsPointIn($p);
	    for my $r2 (@{$self->{InteriorRings}}) {
		next if $ring == $r2;
		croak "point in interior is within another interior" if $r2->IsPointIn($p);
	    }
	}
    }

    # d) A Polygon may not have cut lines, spikes or punctures
    croak "exterior is not simple" unless $self->{ExteriorRing}->IsSimple;
    for my $r (@{$self->{InteriorRings}}) {
	croak "an interior is not simple" unless $r->IsSimple;
    }

    # e) The interior of every Polygon is a connected point set

    # f) The exterior of a Polygon with 1 or more holes is not
    # connected. Each hole defines a connected component of the
    # exterior.

    for my $i (0..$#{$self->{InteriorRings}}) {
	my $r1 = $self->{InteriorRings}[$i];
	my $r2 = $r1->Intersection($self->{ExteriorRing});
	croak "an interior intersects too much with the exterior"
	    if @{$r2->{Geometries}} and (@{$r2->{Geometries}} > 1 or
					 !$r2->{Geometries}[0]->isa('Geo::OGC::Point'));
	for my $j ($i+1..$#{$self->{InteriorRings}}) {
	    my $r2 = $self->{InteriorRings}[$j];
	    $r2 = $r1->Intersection($r2);
	    croak "an interior intersects too much with another interior"
		if @{$r2->{Geometries}} and (@{$r2->{Geometries}} > 1 or
					     !$r2->{Geometries}[0]->isa('Geo::OGC::Point'));
	}
    }

    return 1;

}

sub Is3D {
    my($self) = @_;
    return 1 if $self->{ExteriorRing}->Is3D;
    for my $r (@{$self->{InteriorRings}}) {
	return 1 if $r->Is3D;
    }
    return 0;
}

sub IsMeasured {
    my($self) = @_;
    return 1 if $self->{ExteriorRing}->IsMeasured;
    for my $r (@{$self->{InteriorRings}}) {
	return 1 if $r->IsMeasured;
    }
    return 0;
}

sub AddInteriorRing {
    my($self, $ring, $i) = @_;
    croak 'usage: Polygon->AddInteriorRing($ring[, $i])'
	unless $ring and $ring->isa('Geo::OGC::LinearRing');
    my $rings = $self->{InteriorRings};
    $i = @$rings unless defined $i;
    if (@$rings) {
	my $temp = $rings->[$i-1];
	splice @$rings,$i-1,1,($temp, $ring);
    } else {
	push @$rings, $ring;
    }
}

sub ExteriorRing {
    my($self, $ring) = @_;
    if (defined $ring) {
	croak 'usage: Polygon->ExteriorRing($ring)'
	    unless $ring->isa('Geo::OGC::LinearRing');
	$self->{ExteriorRing} = $ring;
    } else {
	return $self->{ExteriorRing};
    }
}

sub Envelope {
    my($self) = @_;
    return $self->{ExteriorRing}->Envelope;
}

=pod

=over

=item NumInteriorRing()

Return the number of interior rings in the polygon.

=cut

sub NumInteriorRing {
    my($self) = @_;
    scalar(@{$self->{InteriorRings}});
}

sub InteriorRingN {
    my($self, $N, $ring) = @_;
    my $rings = $self->{InteriorRings};
    $rings->[$N-1] = $ring if defined $ring;
    return $rings->[$N-1] if @$rings;
}

# @note Assumes the order of the interior rings is the same.
sub Equals {
    my($self, $geom) = @_;
    return 0 unless $geom->isa('Geo::OGC::Polygon');
    return 0 unless @{$self->{InteriorRings}} == @{$geom->{InteriorRings}};
    return 0 unless $self->{ExteriorRing}->Equals($geom->{ExteriorRing});
    for my $i (0..$#{$self->{InteriorRings}}) {
	return 0 unless $self->{InteriorRings}[$i]->Equals($geom->{InteriorRings}[$i]);
    }
    return 1;
}

sub Area {
    my($self) = @_;
    my $a = $self->{ExteriorRing}->Area;
    for my $ring (@{$self->{InteriorRings}}) {
	$a -= $ring->Area;
    }
    return $a;
}

sub IsPointIn {
    my($self, $point) = @_;
    my $c = $self->{ExteriorRing}->IsPointIn($point);
    if ($c) {
	for my $ring (@{$self->{InteriorRings}}) {
	    $c = 0, last if $ring->IsPointIn($point);
	}
    }
    return $c;
}

sub Distance {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return $geom->Distance($self);
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	return $geom->Distance($self);
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	my $dist;
	for my $p (@{$self->{ExteriorRing}->{Points}}) {
	    if ($geom->{ExteriorRing}->IsPointIn($p)) {
		my $c = 1;
		for my $ring (@{$self->{InteriorRings}}) {
		    if ($ring->IsPointIn($p)) {
			my $d = $p->DistanceToLineStringSqr($ring);
			$dist = $d if !(defined $dist) or $d < $dist;
			$c = 0;
		    }
		}
		return 0 if $c;
	    } else {
		my $d = $p->DistanceToLineStringSqr($geom->{ExteriorRing});
		$dist = $d if !(defined $dist) or $d < $dist;
	    }
	}
	return $dist;
    } elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
	my $dist = $self->Distance($geom->{Geometries}[0]);
	for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
	    my $d = $self->Distance($g);
	    $dist = $d if $d < $dist;
	}
	return $dist;
    } else {
	croak "can't compute distance between a ".ref($geom)." and a polygon";
    }
}

sub Within {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	for my $p (@{$self->{ExteriorRing}->{Points}}) {
	    return 0 unless $p->Equals($geom);
	}
	return 1;
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	for my $p (@{$self->{ExteriorRing}->{Points}}) {
	    return 0 unless $p->Within($geom);
	}
	return 1;
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	croak "polygon within ".ref($geom)." is not yet implemented";
	# the exterior and interior rings must be completely within
	# and the other's interiors must not be within ...
	for my $p (@{$self->{ExteriorRing}->{Points}}) {
	    return 0 unless $p->Within($geom);
	}
	return 1;
    } else {
	croak "polygon within ".ref($geom)." is not yet implemented";
    }
}

sub Intersection {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return $geom->Intersection($self);
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	return $geom->Intersection($self);
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	croak "polygon within ".ref($geom)." is not yet implemented";

	# 1st intersection between outer rings $A (this) and $B (other)
	my $A = $self->{ExteriorRing};
	my $B = $geom->{ExteriorRing};
	my $p1 = $A->{Points}[0];
	my $c = $B->IsPointStricktlyIn($p1);

	my $i = 0;
	my $j;
	while (1) {
	    my $p = $A->{Points}[$i];
	    $j = $i, last if $B->IsPointStricktlyOut($p);
	    $i++;
	    last if $i == @{$A->{Points}};
	}
	if (defined $j) { # there is at least one point
	}
	# the list of
	my @A; # lines on A

    } else {
	croak "polygon within ".ref($geom)." is not yet implemented";
    }
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text .= $self->{ExteriorRing}->as_text if $self->{ExteriorRing};
    for my $ring (@{$self->{InteriorRings}}) {
	$text .= ', ';
	$text .= $ring->as_text;
    }
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
    return $text;
}

sub MakeCollection {
    my($self) = @_;
    my $coll = Geo::OGC::MultiPolygon->new;
    $coll->AddGeometry($self);
    return $coll;
}

sub ClosestVertex {
    my($self, $x, $y) = @_;
    return unless $self->{ExteriorRing};
    my($imin, $dmin) = $self->{ExteriorRing}->ClosestVertex($x, $y);
    my $iring = -1;
    my $r = 0;
    for my $ring (@{$self->{InteriorRings}}) {
	my($i, $d) = $ring->ClosestVertex($x, $y);
	($iring, $imin, $dmin) = ($r, $i, $d) if $d < $dmin;
	$r++;
    }
    return ($iring, $imin, $dmin);
}

sub VertexAt {
    my($self, $iring, $ivertex) = @_;
    return $self->{ExteriorRing}->VertexAt($ivertex) if $iring < 0;
    return $self->{InteriorRings}[$iring]->VertexAt($ivertex);
}

sub ClosestPoint {
    my($self, $x, $y) = @_;
    return unless $self->{ExteriorRing};
    my($imin, $pmin, $dmin) = $self->{ExteriorRing}->ClosestPoint($x, $y);
    my $iring = -1;
    my $r = 0;
    for my $ring (@{$self->{InteriorRings}}) {
	my($i, $p, $d) = $ring->ClosestPoint($x, $y);
	($iring, $imin, $pmin, $dmin) = ($r, $i, $p, $d) if $d < $dmin;
	$r++;
    }
    return ($iring, $imin, $pmin, $dmin);
}

sub AddVertex {
    my($self, $ring, $i, $p) = @_;
    $self->{ExteriorRing}->AddVertex($i, $p), return if $ring < 0;
    $self->{InteriorRings}[$ring]->AddVertex($i, $p);
}

sub DeleteVertex {
    my($self, $ring, $i) = @_;
    $self->{ExteriorRing}->DeleteVertex($i), return if $ring < 0;
    $self->{InteriorRings}[$ring]->DeleteVertex($i);
}

## @method LastPolygon()
# @brief Returns self
sub LastPolygon {
    my($self) = @_;
    return $self;
}

=pod

=back

=head2 Geo::OGC::Triangle

=cut

package Geo::OGC::Triangle;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::Polygon );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Polygon::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'Triangle';
}

=pod

=head2 Geo::OGC::PolyhedralSurface

=cut

package Geo::OGC::PolyhedralSurface;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::Surface );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Surface::new($package, %params);
    return $self;
}

sub init {
    my($self, %params) = @_;
    $self->SUPER::init(%params);
    $self->{Patches} = []; # polygon
    if ($params{patches}) {
	for my $p (@{$params{patches}}) {
	    $self->AddPatch(Geo::OGC::Polygon->new(points=>$p));
	}
    } elsif ($params{patchesm}) {
	for my $p (@{$params{patches}}) {
	    $self->AddPatch(Geo::OGC::Polygon->new(pointsm=>$p));
	}
    }
}

sub copy {
    my($self, $clone) = @_;
    $self->SUPER::copy($clone);
    for my $p (@{$self->{Patches}}){
	$clone->AddPatch($p->Clone);
    }
}

sub IsEmpty {
    my($self) = @_;
    return 1 if !$self->{Patches} or @{$self->{Patches}} == 0;
    return undef;
}

sub GeometryType {
    return 'PolyhedralSurface';
}

sub AddPatch {
    my($self, $patch, $i) = @_;
    croak 'usage: PolyhedralSurface->AddPatch($patch[, $i])'
	unless $patch and $patch->isa('Geo::OGC::Polygon');
    my $patches = $self->{Patches};
    $i = @$patches unless defined $i;
    if (@$patches) {
	my $temp = $patches->[$i-1];
	splice @$patches,$i-1,1,($temp, $patch);
    } else {
	push @$patches, $patch;
    }
}

sub NumPatches {
    my($self) = @_;
    @{$self->{Patches}};
}

sub PatchN {
    my($self, $N, $patch) = @_;
    my $patches = $self->{Patches};
    $patches->[$N-1] = $patch if defined $patch;
    return $patches->[$N-1] if @$patches;
}

sub BoundingPolygons {
    my($self, $p) = @_;
    croak "BoundingPolygons method for class ".ref($self)." is not implemented yet";
}

sub IsClosed {
    my($self) = @_;
    croak "IsClosed method for class ".ref($self)." is not implemented yet";
}

sub IsMeasured {
    my($self) = @_;
    for my $p (@{$self->{Patches}}) {
	return 1 if $p->IsMeasured;
    }
    return 0;
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text = '';
    for my $patch (@{$self->{Patches}}) {
	$text .= ', ';
	$text .= $patch->as_text;
    }
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
    return $text;
}

=pod

=head2 Geo::OGC::TIN

=cut

package Geo::OGC::TIN;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::PolyhedralSurface );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Surface::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'TIN';
}

=pod

=head2 Geo::OGC::GeometryCollection

=cut

package Geo::OGC::GeometryCollection;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::Geometry );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::Geometry::new($package, %params);
    return $self;
}

sub init {
    my($self, %params) = @_;
    $self->SUPER::init(%params);
    $self->{Geometries} = [];
}

sub copy {
    my($self, $clone) = @_;
    $self->SUPER::copy($clone);
    for my $g (@{$self->{Geometries}}) {
	$clone->AddGeometry($g->Clone);
    }
}

sub IsEmpty {
    my($self) = @_;
    return 1 if !$self->{Geometries} or @{$self->{Geometries}} == 0;
    return undef;
}

sub GeometryType {
    return 'GeometryCollection';
}

sub Dimension {
    my($self) = @_;
    my $dim;
    for my $g (@{$self->{Geometries}}) {
	my $d = $g->Dimension;
	$dim = $d if !(defined $dim) or $d > $dim;
    }
    return $dim;
}

sub Is3D {
    my($self) = @_;
    for my $g (@{$self->{Geometries}}) {
	return 1 if $g->Is3D;
    }
    return 0;
}

sub IsMeasured {
    my($self) = @_;
    for my $g (@{$self->{Geometries}}) {
	return 1 if $g->IsMeasured;
    }
    return 0;
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text = join(',', map {$_->as_text(1, 1)} @{$self->{Geometries}});
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
    return $text;
}

sub ElementType {
    return 'Geometry';
}

sub AddGeometry {
    my($self, $geometry, $i) = @_;
    croak 'usage: GeometryCollection->AddGeometry($geometry[, $i])'
	unless $geometry and $geometry->isa('Geo::OGC::Geometry');
    my $geometries = $self->{Geometries};
    $i = @$geometries unless defined $i;
    if (@$geometries) {
	my $temp = $geometries->[$i-1];
	splice @$geometries,$i-1,1,($temp, $geometry);
    } else {
	push @$geometries, $geometry;
    }
}

=pod

=over

=item NumGeometries()

Return the number of geometries in this collection.

=cut

sub NumGeometries {
    my($self) = @_;
    @{$self->{Geometries}};
}

=pod

=item GeometryN($n)

Return the Nth geometry from the collection (the index of the first geometry is 1).

=cut

sub GeometryN {
    my($self, $N, $geometry) = @_;
    my $geometries = $self->{Geometries};
    $geometries->[$N-1] = $geometry if defined $geometry;
    return $geometries->[$N-1] if @$geometries;
}

sub Envelope {
    my($self) = @_;
    my($minx, $miny, $maxx, $maxy);
    for my $g (@{$self->{Geometries}}) {
	my $e = $g->Envelope;
	my $min = $e->PointN(1);
	my $max = $e->PointN(3);
	$minx = $min->{X} if !defined($minx) or $minx > $min->{X};
	$miny = $min->{Y} if !defined($miny) or $miny > $min->{Y};
	$maxx = $max->{X} if !defined($maxx) or $maxx > $max->{X};
	$maxy = $max->{Y} if !defined($maxy) or $maxy > $max->{Y};
    }
    my $r = new Geo::OGC::LinearRing;
    $r->AddPoint(Geo::OGC::Point->new($minx, $miny));
    $r->AddPoint(Geo::OGC::Point->new($maxx, $miny));
    $r->AddPoint(Geo::OGC::Point->new($maxx, $maxy));
    $r->AddPoint(Geo::OGC::Point->new($minx, $maxy));
    $r->Close;
    return $r;
}

# Assumes the order is the same.
sub Equals {
    my($self, $geom) = @_;
    return 0 unless $geom->isa('Geo::OGC::GeometryCollection');
    return 0 unless @{$self->{Geometries}} == @{$geom->{Geometries}};
    for my $i (0..$#{$self->{Geometries}}) {
	return 0 unless $self->{Geometries}[$i]->Equals($geom->{Geometries}[$i]);
    }
    return 1;
}

sub Distance {
    my($self, $geom) = @_;
    if ($geom->isa('Geo::OGC::Point')) {
	return $geom->Distance($self);
    } elsif ($geom->isa('Geo::OGC::LineString')) {
	return $geom->Distance($self);
    } elsif ($geom->isa('Geo::OGC::Polygon')) {
	return $geom->Distance($self);
    } elsif ($geom->isa('Geo::OGC::GeometryCollection')) {
	my $dist = $self->Distance($geom->{Geometries}[0]);
	for my $g (@{$geom->{Geometries}}[1..$#{$geom->{Geometries}}]) {
	    my $d = $g->Distance($self);
	    $dist = $d if $d < $dist;
	}
	return $dist;
    } else {
	croak "can't compute distance between a ".ref($geom)." and a geometry collection";
    }
}

sub ClosestVertex {
    my($self, $x, $y) = @_;
    return unless @{$self->{Geometries}};
    my @rmin = $self->{Geometries}[0]->ClosestVertex($x, $y);
    my $imin = 0;
    my $i = 0;
    for my $g (@{$self->{Geometries}}) {
	$i++, next if $i == 0;
	my @r = $g->ClosestVertex($x, $y);
	if ($r[$#r] < $rmin[$#rmin]) {
	    @rmin = @r;
	    $imin = $i;
	}
	$i++;
    }
    return ($imin, @rmin);
}

sub VertexAt {
    my $self = shift;
    my $i = shift;
    return $self->{Geometries}[$i]->VertexAt(@_);
}

sub ClosestPoint {
    my($self, $x, $y) = @_;
    return unless @{$self->{Geometries}};
    my @rmin = $self->{Geometries}[0]->ClosestPoint($x, $y);
    my $imin = 0;
    my $i = 0;
    for my $g (@{$self->{Geometries}}) {
	$i++, next if $i == 0;
	my @r = $g->ClosestPoint($x, $y);
	if ($r[$#r] < $rmin[$#rmin]) {
	    @rmin = @r;
	    $imin = $i;
	}
	$i++;
    }
    return ($imin, @rmin);
}

sub AddVertex {
    my $self = shift;
    my $i = shift;
    $self->{Geometries}[$i]->AddVertex(@_);
}

sub DeleteVertex {
    my $self = shift;
    my $i = shift;
    $self->{Geometries}[$i]->DeleteVertex(@_);
}

sub LastPolygon {
    my($self) = @_;
    for (my $i = $#{$self->{Geometries}}; $i >= 0; $i--) {
	my $polygon = $self->{Geometries}[$i]->LastPolygon;
	return $polygon if $polygon;
    }
}

=pod

=back

=head2 Geo::OGC::MultiSurface

=cut

package Geo::OGC::MultiSurface;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::GeometryCollection );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::GeometryCollection::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'MultiSurface';
}

sub ElementType {
    return 'Surface';
}

sub Area {
    my($self) = @_;
    croak "Area method for class ".ref($self)." is not implemented yet";
}

sub Centroid {
    my($self) = @_;
    croak "Centroid method for class ".ref($self)." is not implemented yet";
}

sub PointOnSurface {
    my($self) = @_;
    croak "PointOnSurface method for class ".ref($self)." is not implemented yet";
}

=pod

=head2 Geo::OGC::MultiCurve

=cut

package Geo::OGC::MultiCurve;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::GeometryCollection );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::GeometryCollection::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'MultiCurve';
}

sub ElementType {
    return 'Curve';
}

=pod

=head2 Geo::OGC::MultiPoint

=cut

package Geo::OGC::MultiPoint;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::GeometryCollection );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::GeometryCollection::new($package, %params);
    return $self;
}

sub init {
    my($self, %params) = @_;
    $self->SUPER::init(%params);
    if ($params{points}) {
	for my $p (@{$params{points}}) {
	    $self->AddGeometry(Geo::OGC::Point->new(point=>$p));
	}
    } elsif ($params{pointsm}) {
	for my $p (@{$params{pointsm}}) {
	    $self->AddGeometry(Geo::OGC::Point->new(pointm=>$p));
	}
    }
}

sub GeometryType {
    return 'MultiPoint';
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text = join(',', map {$_->as_text(1)} @{$self->{Geometries}});
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
    return $text;
}

sub ElementType {
    return 'Point';
}

sub Boundary {
    my($self) = @_;
    return Geo::OGC::GeometryCollection->new();
}

=pod

=head2 Geo::OGC::MultiPolygon

=cut

package Geo::OGC::MultiPolygon;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::MultiSurface );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::MultiSurface::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'MultiPolygon';
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text = join(',', map {$_->as_text(1)} @{$self->{Geometries}});
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
    return $text;
}

sub ElementType {
    return 'Polygon';
}

=pod

=over

=item LastPolygon()

Returns last (latest added) polygon or undef.

=cut

sub LastPolygon {
    return undef;
}

=pod

=back

=head2 Geo::OGC::MultiLineString

=cut

package Geo::OGC::MultiLineString;

use strict;
use Carp;

our @ISA = qw( Geo::OGC::MultiCurve );

sub new {
    my($package, %params) = @_;
    my $self = Geo::OGC::MultiCurve::new($package, %params);
    return $self;
}

sub GeometryType {
    return 'MultiLineString';
}

sub as_text {
    my($self, $force_parens, $include_tag) = @_;
    my $text = join(',', map {$_->as_text(1)} @{$self->{Geometries}});
    my $ZM = $self->Is3D ? 'Z' : '';
    $ZM .= 'M' if $self->IsMeasured;
    $text = uc($self->GeometryType).$ZM.' ('.$text.')' if $include_tag;
    return $text;
}

sub ElementType {
    return 'LineString';
}

=pod

=head1 ACKNOWLEDGEMENTS

The OpenGIS (r) Implementation Standard for Geographic information -
Simple feature access - Part 1: Common architecture was heavily used
in preparation of this module.

=head1 REPOSITORY

L

=head1 AUTHOR

Ari Jolma L

=head1 COPYRIGHT AND LICENSE

Copyright (c) 2008- Ari Jolma

This program is free software; you can redistribute it and/or modify
it under the same terms as Perl 5 itself.

=cut

1;
__END__