Hi Sanak,<br>Thanks a lot for the modified code: it worked on perfection!<br>Here is a screenshot of a run with my shapefile:<br><br><a href="http://ladybug.no-ip.org/files/inCircleFinal.png">http://ladybug.no-ip.org/files/inCircleFinal.png</a><br>
<br> Have a good rest of the weekend!<br> Jo<br><br><div class="gmail_quote">2009/7/19 <span dir="ltr"><<a href="mailto:geos-devel-request@lists.osgeo.org">geos-devel-request@lists.osgeo.org</a>></span><br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Send geos-devel mailing list submissions to<br>
<a href="mailto:geos-devel@lists.osgeo.org">geos-devel@lists.osgeo.org</a><br>
<br>
To subscribe or unsubscribe via the World Wide Web, visit<br>
<a href="http://lists.osgeo.org/mailman/listinfo/geos-devel" target="_blank">http://lists.osgeo.org/mailman/listinfo/geos-devel</a><br>
or, via email, send a message with subject or body 'help' to<br>
<a href="mailto:geos-devel-request@lists.osgeo.org">geos-devel-request@lists.osgeo.org</a><br>
<br>
You can reach the person managing the list at<br>
<a href="mailto:geos-devel-owner@lists.osgeo.org">geos-devel-owner@lists.osgeo.org</a><br>
<br>
When replying, please edit your Subject line so it is more specific<br>
than "Re: Contents of geos-devel digest..."<br>
<br>
<br>
Today's Topics:<br>
<br>
1. Re: Re: Computational Geometry Problem (Sanak Goe)<br>
<br>
<br>
----------------------------------------------------------------------<br>
<br>
Message: 1<br>
Date: Sun, 19 Jul 2009 21:24:34 +0900<br>
From: Sanak Goe <<a href="mailto:geosanak@gmail.com">geosanak@gmail.com</a>><br>
Subject: Re: [geos-devel] Re: Computational Geometry Problem<br>
To: GEOS Development List <<a href="mailto:geos-devel@lists.osgeo.org">geos-devel@lists.osgeo.org</a>><br>
Message-ID:<br>
<<a href="mailto:5f9be0a0907190524n5a8de434y38731cf04e7770d0@mail.gmail.com">5f9be0a0907190524n5a8de434y38731cf04e7770d0@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="iso-8859-1"<br>
<br>
Skipped content of type multipart/alternative-------------- next part --------------<br>
#include <locale.h><br>
#include <iostream><br>
#include <vector><br>
#include <algorithm><br>
#include <stdexcept><br>
#include <map> // added by sanak 2009.07.19<br>
#ifdef _MSC_VER<br>
#include <io.h><br>
#endif // _MSC_VER<br>
<br>
//Geos<br>
#include <geos_c.h><br>
<br>
#include <geos/algorithm/CGAlgorithms.h><br>
#include <geos/algorithm/HCoordinate.h><br>
#include <geos/geom/Coordinate.h><br>
#include <geos/geom/CoordinateArraySequence.h><br>
#include <geos/geom/LineString.h><br>
#include <geos/geom/MultiPoint.h><br>
#include <geos/geom/Point.h><br>
#include <geos/geom/Polygon.h><br>
#include <geos/io/WKTReader.h><br>
#include <geos/io/WKTWriter.h><br>
<br>
// <-- modified by sanak 2009.07.19<br>
//const double epsilonArea = 100.0;<br>
const double epsilonRadius = 0.1;<br>
// --> modified by sanak 2009.07.19<br>
<br>
// <-- Boost code<br>
template < class BidirectionalIterator ><br>
bool next_combination ( BidirectionalIterator first1 ,<br>
BidirectionalIterator last1 ,<br>
BidirectionalIterator first2 ,<br>
BidirectionalIterator last2 )<br>
{<br>
if (( first1 == last1 ) || ( first2 == last2 )) {<br>
return false ;<br>
}<br>
BidirectionalIterator m1 = last1 ;<br>
BidirectionalIterator m2 = last2 ; --m2;<br>
while (--m1 != first1 && !(* m1 < *m2 )){<br>
}<br>
bool result = (m1 == first1 ) && !(* first1 < *m2 );<br>
if (! result ) {<br>
while ( first2 != m2 && !(* m1 < * first2 )) {<br>
++ first2 ;<br>
}<br>
first1 = m1;<br>
std :: iter_swap (first1 , first2 );<br>
++ first1 ;<br>
++ first2 ;<br>
}<br>
if (( first1 != last1 ) && ( first2 != last2 )) {<br>
m1 = last1 ; m2 = first2 ;<br>
while (( m1 != first1 ) && (m2 != last2 )) {<br>
std :: iter_swap (--m1 , m2 );<br>
++ m2;<br>
}<br>
std :: reverse (first1 , m1 );<br>
std :: reverse (first1 , last1 );<br>
std :: reverse (m2 , last2 );<br>
std :: reverse (first2 , last2 );<br>
}<br>
return ! result ;<br>
}<br>
<br>
template < class BidirectionalIterator ><br>
bool next_combination ( BidirectionalIterator first ,<br>
BidirectionalIterator middle ,<br>
BidirectionalIterator last )<br>
{<br>
return next_combination (first , middle , middle , last );<br>
}<br>
// Boost Code --><br>
<br>
bool computeIncircle(<br>
const geos::geom::Polygon* poly,<br>
const geos::geom::Polygon* negbuf,<br>
geos::geom::Coordinate& center,<br>
double& radius)<br>
{<br>
radius = 0.0;<br>
<br>
if (negbuf == NULL)<br>
{<br>
negbuf = poly;<br>
}<br>
<br>
try<br>
{<br>
<br>
if (!negbuf->isValid())<br>
{<br>
throw std::runtime_error("Invalid polygon.");<br>
}<br>
<br>
if (negbuf->getNumInteriorRing() > 0)<br>
{<br>
throw std::runtime_error("This sample only supports non-holes polygon.");<br>
}<br>
<br>
const geos::geom::LineString* shell = negbuf->getExteriorRing();<br>
geos::geom::CoordinateArraySequence* coords = static_cast<geos::geom::CoordinateArraySequence*>(shell->getCoordinates());<br>
<br>
// compute orientation of ring<br>
if (!geos::algorithm::CGAlgorithms::isCCW(coords))<br>
{<br>
coords->reverse(coords);<br>
}<br>
<br>
// compute angle bisector of every continuous three points in polygon<br>
size_t cnt = coords->size() - 1;<br>
// <-- modified by sanak 2009.07.19<br>
//std::vector<geos::algorithm::HCoordinate*> bisects;<br>
std::vector<size_t> bisectids;<br>
std::map<size_t, geos::algorithm::HCoordinate*> bisects;<br>
// --> modified by sanak 2009.07.19<br>
for (size_t i = 0; i < cnt; i++)<br>
{<br>
size_t prev = (i == 0) ? cnt - 1 : i - 1;<br>
size_t curr = i;<br>
size_t next = (i == cnt - 1) ? 0 : i + 1;<br>
<br>
// <-- JTS 1.10 Code<br>
geos::geom::Coordinate a = coords->getAt(prev);<br>
geos::geom::Coordinate b = coords->getAt(curr);<br>
geos::geom::Coordinate c = coords->getAt(next);<br>
<br>
double len0 = b.distance(a);<br>
double len2 = b.distance(c);<br>
double frac = len0 / (len0 + len2);<br>
double dx = c.x - a.x;<br>
double dy = c.y - a.y;<br>
<br>
geos::geom::Coordinate splitPt = geos::geom::Coordinate(a.x + frac * dx, a.y + frac * dy);<br>
// --> JTS 1.10 Code<br>
// <-- modified by sanak 2009.07.19<br>
//bisects.push_back(new geos::algorithm::HCoordinate(b, splitPt));<br>
bisectids.push_back(i);<br>
bisects.insert(std::pair<size_t, geos::algorithm::HCoordinate*>(i, new geos::algorithm::HCoordinate(b, splitPt)));<br>
// --> modified by sanak 2009.07.19<br>
}<br>
<br>
// get intersection points of angle bisectors<br>
geos::geom::CoordinateArraySequence intsects;<br>
do<br>
{<br>
// <-- modified by sanak 2009.07.19<br>
//geos::algorithm::HCoordinate* hc0 = bisects[0];<br>
//geos::algorithm::HCoordinate* hc1 = bisects[1];<br>
geos::algorithm::HCoordinate* hc0 = bisects[bisectids[0]];<br>
geos::algorithm::HCoordinate* hc1 = bisects[bisectids[1]];<br>
// --> modified by sanak 2009.07.19<br>
geos::algorithm::HCoordinate* hcoord = new geos::algorithm::HCoordinate(*hc0, *hc1);<br>
try<br>
{<br>
// get intersection<br>
geos::geom::Coordinate coord = geos::geom::Coordinate(hcoord->getX(), hcoord->getY());<br>
<br>
// TODO:remove duplicate point<br>
intsects.add(coord);<br>
}<br>
catch (std::exception& ex)<br>
{<br>
// no intersection<br>
}<br>
if (hcoord)<br>
{<br>
delete hcoord;<br>
}<br>
}<br>
// <-- modified by sanak 2009.07.19<br>
//while (next_combination(bisects.begin(), bisects.begin() + 2, bisects.end()));<br>
while (next_combination(bisectids.begin(), bisectids.begin() + 2, bisectids.end()));<br>
// --> modified by sanak 2009.07.19<br>
<br>
// clean up<br>
if (bisects.size() > 0)<br>
{<br>
std::cout << "bissects ok." << std::endl;<br>
// <-- modified by sanak 2009.07.19<br>
//std::vector<geos::algorithm::HCoordinate*>::iterator it;<br>
//for (it = bisects.begin(); it != bisects.end(); it++)<br>
//{<br>
// delete *it;<br>
//}<br>
std::map<size_t, geos::algorithm::HCoordinate*>::iterator it;<br>
for (it = bisects.begin(); it != bisects.end(); it++)<br>
{<br>
delete it->second;<br>
}<br>
// --> modified by sanak 2009.07.19<br>
}<br>
else<br>
std::cout << "No bissects." << std::endl;<br>
<br>
std::cout << "intsects.size()" << intsects.size() << std::endl;<br>
if (intsects.size() > 0)<br>
{<br>
std::cout << "intsects ok." << std::endl;<br>
<br>
geos::geom::Point* incenter = NULL;<br>
<br>
geos::geom::GeometryFactory factory;<br>
const geos::geom::LineString* line = poly->getExteriorRing();<br>
//geos::geom::MultiPoint* rcpoints = factory.createMultiPoint(intsects); // for debug<br>
<br>
std::vector<geos::geom::Coordinate*>::iterator it;<br>
for (size_t i = 0; i < intsects.size(); i++)<br>
{<br>
geos::geom::Coordinate coord = intsects[i];<br>
// exclude points that are outer of polygon<br>
if (geos::algorithm::CGAlgorithms::isPointInRing(coord, coords))<br>
{<br>
geos::geom::Point* point = factory.createPoint(coord);<br>
// compute distance from every candidates of centres to polygon, and choose largest one<br>
double distance = line->distance(point);<br>
if (radius < distance)<br>
{<br>
radius = distance;<br>
incenter = point;<br>
}<br>
}<br>
//else std::cout << "Point not in Ring." << std::endl;<br>
}<br>
<br>
if (incenter != NULL && radius != 0.0)<br>
{<br>
geos::io::WKTWriter writer;<br>
std::cout << "wkt : " << writer.write(incenter) << std::endl;<br>
std::cout << "radius : " << radius << std::endl;<br>
<br>
center.x = incenter->getX();<br>
center.y = incenter->getY();<br>
//return true; // deleted by sanak 2009.07.19<br>
}<br>
// <-- added by sanak 2009.07.19<br>
else<br>
{<br>
geos::geom::Point* intPt = negbuf->getInteriorPoint();<br>
center.x = intPt->getX();<br>
center.y = intPt->getY();<br>
radius = line->distance(intPt);<br>
}<br>
return true;<br>
// --> added by sanak 2009.07.19<br>
}<br>
else<br>
{<br>
std::cout << "No intsects." << std::endl;<br>
throw std::runtime_error("Couldn't get incenter.");<br>
}<br>
}<br>
catch (std::exception& ex)<br>
{<br>
std::cerr << "Exception : " << ex.what() << std::endl;<br>
throw ex;<br>
return false;<br>
}<br>
return false;<br>
}<br>
<br>
int main(int argc, char* argv[])<br>
{<br>
// check stdin<br>
#ifdef _MSC_VER<br>
if (::_isatty(fileno(stdin)))<br>
#else<br>
if (::isatty(fileno(stdin)))<br>
#endif // _MSC_VER<br>
{<br>
std::cerr << "Usage: %s < [wktfile(polygon)]" << std::endl;<br>
return -1;<br>
}<br>
<br>
std::string line;<br>
geos::io::WKTReader reader;<br>
std::vector<geos::geom::Polygon*> g;<br>
while (getline(std::cin, line, '\n'))<br>
{<br>
//std::cout << line << std::endl;<br>
geos::geom::Polygon* geom =<br>
static_cast<geos::geom::Polygon*>(reader.read(line));<br>
<br>
if (geom == NULL && geom->getGeometryType() != "Polygon")<br>
{<br>
throw std::runtime_error("This sample only supports polygon geometry");<br>
}<br>
<br>
geos::geom::Polygon* poly = static_cast<geos::geom::Polygon*>(geom);<br>
if (!poly->isValid())<br>
{<br>
throw std::runtime_error("Invalid polygon.");<br>
}<br>
<br>
if (poly->getNumInteriorRing() > 0)<br>
{<br>
throw std::runtime_error("This sample only supports non-holes polygon.");<br>
}<br>
<br>
g.push_back(poly);<br>
}<br>
<br>
<br>
try<br>
{<br>
<br>
for (int i=0 ; i < g.size(); ++i)<br>
{<br>
geos::geom::Polygon* negbuf = NULL;<br>
double prevRadius = 0.0; // added by sanak 2009.07.19<br>
while (true)<br>
{<br>
geos::geom::Coordinate center;<br>
double radius = 0.0;<br>
if (computeIncircle(g[i], negbuf, center, radius))<br>
{<br>
// <-- modified by sanak 2009.07.19<br>
//negbuf = static_cast<geos::geom::Polygon*>(g[i]->buffer(-radius));<br>
//if (negbuf == NULL && negbuf->getGeometryType() != "Polygon")<br>
//{<br>
// // TODO:<br>
// throw std::runtime_error("Intercepted!!!!!!!");<br>
//}<br>
//double area = negbuf->getArea();<br>
<br>
//if (area < epsilonArea)<br>
geos::geom::Geometry* buf = g[i]->buffer(-(radius));<br>
std::string buftype = buf->getGeometryType();<br>
if (buftype == "Polygon")<br>
{<br>
negbuf = static_cast<geos::geom::Polygon*>(buf);<br>
}<br>
else if (buftype == "MultiPolygon")<br>
{<br>
size_t cnt = buf->getNumGeometries();<br>
double maxradius = 0.0;<br>
double tmpradius = 0.0;<br>
geos::geom::Coordinate tmpcenter;<br>
for (size_t j = 0; j < cnt; j++)<br>
{<br>
const geos::geom::Polygon* item = static_cast<const geos::geom::Polygon*>(buf->getGeometryN(j));<br>
if (computeIncircle(g[i], item, tmpcenter, tmpradius))<br>
{<br>
if (tmpradius > maxradius)<br>
{<br>
negbuf = static_cast<geos::geom::Polygon*>(item->clone());<br>
center = tmpcenter;<br>
radius = tmpradius;<br>
maxradius = tmpradius;<br>
}<br>
}<br>
}<br>
}<br>
<br>
if ((radius - prevRadius) < epsilonRadius)<br>
// --> modified by sanak 2009.07.19<br>
{<br>
<br>
std::cout << "center : " << center.toString() << std::endl;<br>
std::cout << "radius : " << radius << std::endl;<br>
geos::geom::GeometryFactory factory;<br>
geos::geom::Point* centerPt = factory.createPoint(center);<br>
geos::geom::Geometry* incircle = centerPt->buffer(radius);<br>
geos::io::WKTWriter writer;<br>
std::cout << "wkt : " << writer.write(incircle) << std::endl;<br>
// <-- modified by sanak 2009.07.19<br>
//return 0;<br>
break;<br>
// --> modified by sanak 2009.07.19<br>
}<br>
prevRadius = radius; // added by sanak 2009.07.19<br>
}<br>
else<br>
{<br>
std::cout << "Could not compute!" << std::endl;<br>
return -2;<br>
}<br>
}<br>
}<br>
}<br>
catch (std::exception& ex)<br>
{<br>
std::cerr << "Exception : " << ex.what() << std::endl;<br>
return -3;<br>
}<br>
<br>
return 0;<br>
}<br>
<br>
------------------------------<br>
<br>
_______________________________________________<br>
geos-devel mailing list<br>
<a href="mailto:geos-devel@lists.osgeo.org">geos-devel@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/geos-devel" target="_blank">http://lists.osgeo.org/mailman/listinfo/geos-devel</a><br>
<br>
End of geos-devel Digest, Vol 81, Issue 13<br>
******************************************<br>
</blockquote></div><br><br clear="all"><br>-- <br>"#define QUESTION ((bb) || !(bb))" (Shakespeare)<br><br>