[geos-devel] Re: Computational Geometry Problem
Sanak Goe
geosanak at gmail.com
Sun Jul 19 08:24:34 EDT 2009
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;
}
More information about the geos-devel
mailing list