[geos-devel] Computational Geometry Problem + Exception
Jo
doublebyte at gmail.com
Mon Jul 20 18:16:26 EDT 2009
Hi,
Sorry to bother you again, but I am getting sometimes this exception from
GEOS and I don't have a clue what it is:
test: CoordinateArraySequence.cpp:105: virtual const geos::geom::Coordinate&
geos::geom::CoordinateArraySequence::getAt(size_t) const: Assertion
`pos<vect->size()' failed.
Backtracing it, I discovered that it comes from the ISCCW algorithm:
#0 0xb7119c66 in raise () from /lib/libc.so.6
#1 0xb711b571 in abort () from /lib/libc.so.6
#2 0xb7112e60 in __assert_fail () from /lib/libc.so.6
#3 0xb7884ee6 in geos::geom::CoordinateArraySequence::getAt (this=0x0,
pos=7117)
at
/usr/local/lib/gcc/i686-pc-linux-gnu/3.4.6/../../../../include/c++/3.4.6/bits/stl_vector.h:375
#4 0xb7876b79 in geos::algorithm::CGAlgorithms::isCCW (ring=0x804f2b8) at
CGAlgorithms.cpp:175
Does anybody has any ideas in what can be causing this?
Strangely it breaks in different files in my two computers: both using GEOS,
vs 3.2.0 from svn.
Any help on this would be much appreciated!
cheers,
Jo
2009/7/19 <geos-devel-request at lists.osgeo.org>
> Send geos-devel mailing list submissions to
> geos-devel at lists.osgeo.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.osgeo.org/mailman/listinfo/geos-devel
> or, via email, send a message with subject or body 'help' to
> geos-devel-request at lists.osgeo.org
>
> You can reach the person managing the list at
> geos-devel-owner at lists.osgeo.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of geos-devel digest..."
>
>
> Today's Topics:
>
> 1. geos svn (Jo)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Sun, 19 Jul 2009 23:10:01 +0100
> From: Jo <doublebyte at gmail.com>
> Subject: [geos-devel] geos svn
> To: geos-devel at lists.osgeo.org
> Message-ID:
> <23ab5f0a0907191510s1d567375qdf13e491e1e06933 at mail.gmail.com>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Latest SVN version of geos seems to be incomplete: no libtool and
> lt-main.sh.
> Also, the nightly release appears to be a broken tar file (at least I
> cannot
> open it in my pc)...
> cheers,
> Jo
>
> 2009/7/19 <geos-devel-request at lists.osgeo.org>
>
> > Send geos-devel mailing list submissions to
> > geos-devel at lists.osgeo.org
> >
> > To subscribe or unsubscribe via the World Wide Web, visit
> > http://lists.osgeo.org/mailman/listinfo/geos-devel
> > or, via email, send a message with subject or body 'help' to
> > geos-devel-request at lists.osgeo.org
> >
> > You can reach the person managing the list at
> > geos-devel-owner at lists.osgeo.org
> >
> > When replying, please edit your Subject line so it is more specific
> > than "Re: Contents of geos-devel digest..."
> >
> >
> > Today's Topics:
> >
> > 1. Computational Geometry Problem (Jo)
> >
> >
> > ----------------------------------------------------------------------
> >
> > Message: 1
> > Date: Sun, 19 Jul 2009 22:56:53 +0100
> > From: Jo <doublebyte at gmail.com>
> > Subject: [geos-devel] Computational Geometry Problem
> > To: geos-devel at lists.osgeo.org
> > Message-ID:
> > <23ab5f0a0907191456l5f8d648i4dce9666abf3dcf8 at mail.gmail.com>
> > Content-Type: text/plain; charset="iso-8859-1"
> >
> > Hi Sanak,
> > Thanks a lot for the modified code: it worked on perfection!
> > Here is a screenshot of a run with my shapefile:
> >
> > http://ladybug.no-ip.org/files/inCircleFinal.png
> >
> > Have a good rest of the weekend!
> > Jo
> >
> > 2009/7/19 <geos-devel-request at lists.osgeo.org>
> >
> > > Send geos-devel mailing list submissions to
> > > geos-devel at lists.osgeo.org
> > >
> > > To subscribe or unsubscribe via the World Wide Web, visit
> > > http://lists.osgeo.org/mailman/listinfo/geos-devel
> > > or, via email, send a message with subject or body 'help' to
> > > geos-devel-request at lists.osgeo.org
> > >
> > > You can reach the person managing the list at
> > > geos-devel-owner at lists.osgeo.org
> > >
> > > When replying, please edit your Subject line so it is more specific
> > > than "Re: Contents of geos-devel digest..."
> > >
> > >
> > > Today's Topics:
> > >
> > > 1. Re: Re: Computational Geometry Problem (Sanak Goe)
> > >
> > >
> > > ----------------------------------------------------------------------
> > >
> > > Message: 1
> > > Date: Sun, 19 Jul 2009 21:24:34 +0900
> > > From: Sanak Goe <geosanak at gmail.com>
> > > Subject: Re: [geos-devel] Re: Computational Geometry Problem
> > > To: GEOS Development List <geos-devel at lists.osgeo.org>
> > > Message-ID:
> > > <5f9be0a0907190524n5a8de434y38731cf04e7770d0 at mail.gmail.com>
> > > Content-Type: text/plain; charset="iso-8859-1"
> > >
> > > Skipped content of type multipart/alternative-------------- next part
> > > --------------
> > > #include <locale.h>
> > > #include <iostream>
> > > #include <vector>
> > > #include <algorithm>
> > > #include <stdexcept>
> > > #include <map> // added by sanak 2009.07.19
> > > #ifdef _MSC_VER
> > > #include <io.h>
> > > #endif // _MSC_VER
> > >
> > > //Geos
> > > #include <geos_c.h>
> > >
> > > #include <geos/algorithm/CGAlgorithms.h>
> > > #include <geos/algorithm/HCoordinate.h>
> > > #include <geos/geom/Coordinate.h>
> > > #include <geos/geom/CoordinateArraySequence.h>
> > > #include <geos/geom/LineString.h>
> > > #include <geos/geom/MultiPoint.h>
> > > #include <geos/geom/Point.h>
> > > #include <geos/geom/Polygon.h>
> > > #include <geos/io/WKTReader.h>
> > > #include <geos/io/WKTWriter.h>
> > >
> > > // <-- modified by sanak 2009.07.19
> > > //const double epsilonArea = 100.0;
> > > const double epsilonRadius = 0.1;
> > > // --> modified by sanak 2009.07.19
> > >
> > > // <-- Boost code
> > > template < class BidirectionalIterator >
> > > bool next_combination ( BidirectionalIterator first1 ,
> > > BidirectionalIterator last1 ,
> > > BidirectionalIterator first2 ,
> > > BidirectionalIterator last2 )
> > > {
> > > if (( first1 == last1 ) || ( first2 == last2 )) {
> > > return false ;
> > > }
> > > BidirectionalIterator m1 = last1 ;
> > > BidirectionalIterator m2 = last2 ; --m2;
> > > while (--m1 != first1 && !(* m1 < *m2 )){
> > > }
> > > bool result = (m1 == first1 ) && !(* first1 < *m2 );
> > > if (! result ) {
> > > while ( first2 != m2 && !(* m1 < * first2 )) {
> > > ++ first2 ;
> > > }
> > > first1 = m1;
> > > std :: iter_swap (first1 , first2 );
> > > ++ first1 ;
> > > ++ first2 ;
> > > }
> > > if (( first1 != last1 ) && ( first2 != last2 )) {
> > > m1 = last1 ; m2 = first2 ;
> > > while (( m1 != first1 ) && (m2 != last2 )) {
> > > std :: iter_swap (--m1 , m2 );
> > > ++ m2;
> > > }
> > > std :: reverse (first1 , m1 );
> > > std :: reverse (first1 , last1 );
> > > std :: reverse (m2 , last2 );
> > > std :: reverse (first2 , last2 );
> > > }
> > > return ! result ;
> > > }
> > >
> > > template < class BidirectionalIterator >
> > > bool next_combination ( BidirectionalIterator first ,
> > > BidirectionalIterator middle ,
> > > BidirectionalIterator last )
> > > {
> > > return next_combination (first , middle , middle , last );
> > > }
> > > // Boost Code -->
> > >
> > > bool computeIncircle(
> > > const geos::geom::Polygon*
> poly,
> > > const geos::geom::Polygon*
> > negbuf,
> > > geos::geom::Coordinate& center,
> > > double& radius)
> > > {
> > > radius = 0.0;
> > >
> > > if (negbuf == NULL)
> > > {
> > > negbuf = poly;
> > > }
> > >
> > > try
> > > {
> > >
> > > if (!negbuf->isValid())
> > > {
> > > throw std::runtime_error("Invalid polygon.");
> > > }
> > >
> > > if (negbuf->getNumInteriorRing() > 0)
> > > {
> > > throw std::runtime_error("This sample only
> > supports
> > > non-holes polygon.");
> > > }
> > >
> > > const geos::geom::LineString* shell =
> > > negbuf->getExteriorRing();
> > > geos::geom::CoordinateArraySequence* coords =
> > >
> >
> static_cast<geos::geom::CoordinateArraySequence*>(shell->getCoordinates());
> > >
> > > // compute orientation of ring
> > > if (!geos::algorithm::CGAlgorithms::isCCW(coords))
> > > {
> > > coords->reverse(coords);
> > > }
> > >
> > > // compute angle bisector of every continuous three
> > points
> > > in polygon
> > > size_t cnt = coords->size() - 1;
> > > // <-- modified by sanak 2009.07.19
> > > //std::vector<geos::algorithm::HCoordinate*> bisects;
> > > std::vector<size_t> bisectids;
> > > std::map<size_t, geos::algorithm::HCoordinate*> bisects;
> > > // --> modified by sanak 2009.07.19
> > > for (size_t i = 0; i < cnt; i++)
> > > {
> > > size_t prev = (i == 0) ? cnt - 1 : i - 1;
> > > size_t curr = i;
> > > size_t next = (i == cnt - 1) ? 0 : i + 1;
> > >
> > > // <-- JTS 1.10 Code
> > > geos::geom::Coordinate a = coords->getAt(prev);
> > > geos::geom::Coordinate b = coords->getAt(curr);
> > > geos::geom::Coordinate c = coords->getAt(next);
> > >
> > > double len0 = b.distance(a);
> > > double len2 = b.distance(c);
> > > double frac = len0 / (len0 + len2);
> > > double dx = c.x - a.x;
> > > double dy = c.y - a.y;
> > >
> > > geos::geom::Coordinate splitPt =
> > > geos::geom::Coordinate(a.x + frac * dx, a.y + frac * dy);
> > > // --> JTS 1.10 Code
> > > // <-- modified by sanak 2009.07.19
> > > //bisects.push_back(new
> > > geos::algorithm::HCoordinate(b, splitPt));
> > > bisectids.push_back(i);
> > > bisects.insert(std::pair<size_t,
> > > geos::algorithm::HCoordinate*>(i, new geos::algorithm::HCoordinate(b,
> > > splitPt)));
> > > // --> modified by sanak 2009.07.19
> > > }
> > >
> > > // get intersection points of angle bisectors
> > > geos::geom::CoordinateArraySequence intsects;
> > > do
> > > {
> > > // <-- modified by sanak 2009.07.19
> > > //geos::algorithm::HCoordinate* hc0 =
> bisects[0];
> > > //geos::algorithm::HCoordinate* hc1 =
> bisects[1];
> > > geos::algorithm::HCoordinate* hc0 =
> > > bisects[bisectids[0]];
> > > geos::algorithm::HCoordinate* hc1 =
> > > bisects[bisectids[1]];
> > > // --> modified by sanak 2009.07.19
> > > geos::algorithm::HCoordinate* hcoord = new
> > > geos::algorithm::HCoordinate(*hc0, *hc1);
> > > try
> > > {
> > > // get intersection
> > > geos::geom::Coordinate coord =
> > > geos::geom::Coordinate(hcoord->getX(), hcoord->getY());
> > >
> > > // TODO:remove duplicate point
> > > intsects.add(coord);
> > > }
> > > catch (std::exception& ex)
> > > {
> > > // no intersection
> > > }
> > > if (hcoord)
> > > {
> > > delete hcoord;
> > > }
> > > }
> > > // <-- modified by sanak 2009.07.19
> > > //while (next_combination(bisects.begin(),
> bisects.begin()
> > +
> > > 2, bisects.end()));
> > > while (next_combination(bisectids.begin(),
> > bisectids.begin()
> > > + 2, bisectids.end()));
> > > // --> modified by sanak 2009.07.19
> > >
> > > // clean up
> > > if (bisects.size() > 0)
> > > {
> > > std::cout << "bissects ok." << std::endl;
> > > // <-- modified by sanak 2009.07.19
> > >
> > > //std::vector<geos::algorithm::HCoordinate*>::iterator it;
> > > //for (it = bisects.begin(); it !=
> bisects.end();
> > > it++)
> > > //{
> > > // delete *it;
> > > //}
> > > std::map<size_t,
> > > geos::algorithm::HCoordinate*>::iterator it;
> > > for (it = bisects.begin(); it != bisects.end();
> > > it++)
> > > {
> > > delete it->second;
> > > }
> > > // --> modified by sanak 2009.07.19
> > > }
> > > else
> > > std::cout << "No bissects." << std::endl;
> > >
> > > std::cout << "intsects.size()" << intsects.size() <<
> > > std::endl;
> > > if (intsects.size() > 0)
> > > {
> > > std::cout << "intsects ok." << std::endl;
> > >
> > > geos::geom::Point* incenter = NULL;
> > >
> > > geos::geom::GeometryFactory factory;
> > > const geos::geom::LineString* line =
> > > poly->getExteriorRing();
> > > //geos::geom::MultiPoint* rcpoints =
> > > factory.createMultiPoint(intsects); // for debug
> > >
> > > std::vector<geos::geom::Coordinate*>::iterator
> it;
> > > for (size_t i = 0; i < intsects.size(); i++)
> > > {
> > > geos::geom::Coordinate coord =
> > intsects[i];
> > > // exclude points that are outer of
> > polygon
> > > if
> > > (geos::algorithm::CGAlgorithms::isPointInRing(coord, coords))
> > > {
> > > geos::geom::Point* point =
> > > factory.createPoint(coord);
> > > // compute distance from every
> > > candidates of centres to polygon, and choose largest one
> > > double distance =
> > > line->distance(point);
> > > if (radius < distance)
> > > {
> > > radius = distance;
> > > incenter = point;
> > > }
> > > }
> > > //else std::cout << "Point not in Ring."
> > <<
> > > std::endl;
> > > }
> > >
> > > if (incenter != NULL && radius != 0.0)
> > > {
> > > geos::io::WKTWriter writer;
> > > std::cout << "wkt : " <<
> > > writer.write(incenter) << std::endl;
> > > std::cout << "radius : " <<
> radius
> > > << std::endl;
> > >
> > > center.x = incenter->getX();
> > > center.y = incenter->getY();
> > > //return true; // deleted by sanak
> > > 2009.07.19
> > > }
> > > // <-- added by sanak 2009.07.19
> > > else
> > > {
> > > geos::geom::Point* intPt =
> > > negbuf->getInteriorPoint();
> > > center.x = intPt->getX();
> > > center.y = intPt->getY();
> > > radius = line->distance(intPt);
> > > }
> > > return true;
> > > // --> added by sanak 2009.07.19
> > > }
> > > else
> > > {
> > > std::cout << "No intsects." << std::endl;
> > > throw std::runtime_error("Couldn't get
> > incenter.");
> > > }
> > > }
> > > catch (std::exception& ex)
> > > {
> > > std::cerr << "Exception : " << ex.what() << std::endl;
> > > throw ex;
> > > return false;
> > > }
> > > return false;
> > > }
> > >
> > > int main(int argc, char* argv[])
> > > {
> > > // check stdin
> > > #ifdef _MSC_VER
> > > if (::_isatty(fileno(stdin)))
> > > #else
> > > if (::isatty(fileno(stdin)))
> > > #endif // _MSC_VER
> > > {
> > > std::cerr << "Usage: %s < [wktfile(polygon)]" <<
> > std::endl;
> > > return -1;
> > > }
> > >
> > > std::string line;
> > > geos::io::WKTReader reader;
> > > std::vector<geos::geom::Polygon*> g;
> > > while (getline(std::cin, line, '\n'))
> > > {
> > > //std::cout << line << std::endl;
> > > geos::geom::Polygon* geom =
> > > static_cast<geos::geom::Polygon*>(reader.read(line));
> > >
> > > if (geom == NULL && geom->getGeometryType() !=
> "Polygon")
> > > {
> > > throw std::runtime_error("This sample only
> > supports
> > > polygon geometry");
> > > }
> > >
> > > geos::geom::Polygon* poly =
> > > static_cast<geos::geom::Polygon*>(geom);
> > > if (!poly->isValid())
> > > {
> > > throw std::runtime_error("Invalid polygon.");
> > > }
> > >
> > > if (poly->getNumInteriorRing() > 0)
> > > {
> > > throw std::runtime_error("This sample only
> > supports
> > > non-holes polygon.");
> > > }
> > >
> > > g.push_back(poly);
> > > }
> > >
> > >
> > > try
> > > {
> > >
> > > for (int i=0 ; i < g.size(); ++i)
> > > {
> > > geos::geom::Polygon* negbuf = NULL;
> > > double prevRadius = 0.0; // added by sanak
> > > 2009.07.19
> > > while (true)
> > > {
> > > geos::geom::Coordinate center;
> > > double radius = 0.0;
> > > if (computeIncircle(g[i], negbuf,
> center,
> > > radius))
> > > {
> > > // <-- modified by sanak
> > 2009.07.19
> > > //negbuf =
> > > static_cast<geos::geom::Polygon*>(g[i]->buffer(-radius));
> > > //if (negbuf == NULL &&
> > > negbuf->getGeometryType() != "Polygon")
> > > //{
> > > // // TODO:
> > > // throw
> > > std::runtime_error("Intercepted!!!!!!!");
> > > //}
> > > //double area =
> negbuf->getArea();
> > >
> > > //if (area < epsilonArea)
> > > geos::geom::Geometry* buf =
> > > g[i]->buffer(-(radius));
> > > std::string buftype =
> > > buf->getGeometryType();
> > > if (buftype == "Polygon")
> > > {
> > > negbuf =
> > > static_cast<geos::geom::Polygon*>(buf);
> > > }
> > > else if (buftype ==
> > "MultiPolygon")
> > > {
> > > size_t cnt =
> > > buf->getNumGeometries();
> > > double maxradius = 0.0;
> > > double tmpradius = 0.0;
> > > geos::geom::Coordinate
> > > tmpcenter;
> > > for (size_t j = 0; j <
> > cnt;
> > > j++)
> > > {
> > > const
> > > geos::geom::Polygon* item = static_cast<const
> > > geos::geom::Polygon*>(buf->getGeometryN(j));
> > > if
> > > (computeIncircle(g[i], item, tmpcenter, tmpradius))
> > > {
> > > if
> > > (tmpradius > maxradius)
> > > {
> > >
> > > negbuf = static_cast<geos::geom::Polygon*>(item->clone());
> > >
> > > center = tmpcenter;
> > >
> > > radius = tmpradius;
> > >
> > > maxradius = tmpradius;
> > > }
> > > }
> > > }
> > > }
> > >
> > > if ((radius - prevRadius) <
> > > epsilonRadius)
> > > // --> modified by sanak
> > 2009.07.19
> > > {
> > >
> > > std::cout << "center : "
> > <<
> > > center.toString() << std::endl;
> > > std::cout << "radius : "
> > <<
> > > radius << std::endl;
> > >
> > geos::geom::GeometryFactory
> > > factory;
> > > geos::geom::Point*
> > centerPt
> > > = factory.createPoint(center);
> > > geos::geom::Geometry*
> > > incircle = centerPt->buffer(radius);
> > > geos::io::WKTWriter
> > writer;
> > > std::cout << "wkt : " <<
> > > writer.write(incircle) << std::endl;
> > > // <-- modified by sanak
> > > 2009.07.19
> > > //return 0;
> > > break;
> > > // --> modified by sanak
> > > 2009.07.19
> > > }
> > > prevRadius = radius; // added by
> > > sanak 2009.07.19
> > > }
> > > else
> > > {
> > > std::cout << "Could not
> compute!"
> > <<
> > > std::endl;
> > > return -2;
> > > }
> > > }
> > > }
> > > }
> > > catch (std::exception& ex)
> > > {
> > > std::cerr << "Exception : " << ex.what() << std::endl;
> > > return -3;
> > > }
> > >
> > > return 0;
> > > }
> > >
> > > ------------------------------
> > >
> > > _______________________________________________
> > > geos-devel mailing list
> > > geos-devel at lists.osgeo.org
> > > http://lists.osgeo.org/mailman/listinfo/geos-devel
> > >
> > > End of geos-devel Digest, Vol 81, Issue 13
> > > ******************************************
> > >
> >
> >
> >
> > --
> > "#define QUESTION ((bb) || !(bb))" (Shakespeare)
> > -------------- next part --------------
> > An HTML attachment was scrubbed...
> > URL:
> >
> http://lists.osgeo.org/pipermail/geos-devel/attachments/20090719/f83b16e1/attachment.html
> >
> > ------------------------------
> >
> > _______________________________________________
> > geos-devel mailing list
> > geos-devel at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/geos-devel
> >
> > End of geos-devel Digest, Vol 81, Issue 14
> > ******************************************
> >
>
>
>
> --
> "#define QUESTION ((bb) || !(bb))" (Shakespeare)
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL:
> http://lists.osgeo.org/pipermail/geos-devel/attachments/20090719/d602a171/attachment.html
>
> ------------------------------
>
> _______________________________________________
> geos-devel mailing list
> geos-devel at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/geos-devel
>
> End of geos-devel Digest, Vol 81, Issue 15
> ******************************************
>
--
"#define QUESTION ((bb) || !(bb))" (Shakespeare)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/geos-devel/attachments/20090720/40ba7c0d/attachment-0001.html
More information about the geos-devel
mailing list