[geos-devel] Computational Geometry Problem
Jo
doublebyte at gmail.com
Sun Jul 19 17:56:53 EDT 2009
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-0001.html
More information about the geos-devel
mailing list