/*
* This file is part of rasdaman community.
*
* Rasdaman community is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Rasdaman community is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with rasdaman community. If not, see .
*
* Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Peter Baumann /
rasdaman GmbH.
*
* For more information please see
* or contact Peter Baumann via .
/
/**
* SOURCE: polygon.cc
*
* MODULE: rasodmg
* CLASS: r_Polygon
*
* PURPOSE:
* This class maintains 2-D polygon sequences.
* It knows about them being closed or not.
* Input syntax (see constructor) is:
* polygon ::= point+
* point ::= '[' int ',' int ']'
* int ::= ASCII-integer
*
* COMMENTS:
* - see comment in shrinkPoly() about a potential error source
* - r_Error::r_Error_General should be replaced with more specific exception
*
*/
#include "rasodmg/polygon.hh"
#include
#include
#include
#if defined(SOLARIS)
#include
#endif
using std::sort;
//causes problems compiling on old red hat
//using std::iterator;
#include "raslib/miter.hh"
#include "rasodmg/marray.hh"
#include "debug/debug.hh"
static const char rcsid[] = "@(#)rasodmg, r_Polygon: $Header: /home/rasdev/CVS-repository/rasdaman/rasodmg/polygon.cc,v 1.28 2003/12/27 23:02:56 rasdev Exp $";
#ifndef LONG_MAX
const int LONG_MAX = (1<<31) - 1;
#endif
#ifndef LONG_MIN
const int LONG_MIN = (1<<31); // due to overflow
#endif
// ------------------------------------------------------------------
// r_Edge
// ------------------------------------------------------------------
r_Edge::r_Edge(const r_Point& newStart, const r_Point& newEnd) :
start(newStart), end(newEnd)
{
}
const r_Point&
r_Edge::getStart() const
{
return start;
}
const r_Point&
r_Edge::getEnd() const
{
return end;
}
double
r_Edge::getInvSlope() const
{
return (((double)end[0] - start[0]) / (end[1] - start[1]));
}
double
r_Edge::getSlope() const
{
return (((double)end[1] - start[1]) / (end[0] - start[0]));
}
double
r_Edge::getCurrX(r_Range y) const
{
double currX=0.0;
if(end[1]==start[1])
currX=end[0];
else
currX=getInvSlope()*(y - start[1]) + start[0];
return currX;
}
double
r_Edge::getCurrY(r_Range x) const
{
double currY=0.0;
if(end[0]==start[0])
currY=end[1];
else
currY=getSlope()*(x - start[0]) + start[1];
return currY;
}
void
r_Edge::print_status( std::ostream& s ) const
{
start.print_status(s);
s << "->";
end.print_status(s);
}
bool
r_Edge::isHorizontal() const
{
return start[1] == end[1] ? true:false;
}
// ------------------------------------------------------------------
// r_Polygon
// ------------------------------------------------------------------
const r_Dimension r_Polygon::polyPointDim=2;
r_Polygon::r_Polygon(const char* init) throw (r_Error)
: closed(false),
firstPointSet(false)
{
ENTER( "r_Polygon::r_Polygon, init=" << init );
if (init == NULL)
{
TALK( "r_Polygon::r_Polygon(" << (init?init: "NULL") << ")" );
RMInit::logOut << "r_Polygon::r_Polygon(" << (init?init: "NULL") << ")" << std::endl;
throw r_Error(POLYGONWRONGINITSTRING);
}
const int POINTBUFFERLEN=512;
const char* endPos = NULL;
size_t pointLen = 0;
char pointBuffer[POINTBUFFERLEN];
const char* startPos = index(init, '[');
if (startPos == NULL)
{
TALK( "r_Polygon::r_Polygon(" << init << ") the init string has to start with a '['" );
RMInit::logOut << "r_Polygon::r_Polygon(" << init << ") the init string has to start with a '['" << std::endl;
throw r_Error(POLYGONWRONGINITSTRING);
}
// while (true)
do
{
endPos = index(startPos, ']');
if (endPos == NULL)
{
TALK( "r_Polygon::r_Polygon(" << init << ") the init string has to contain valid r_Point definitions" );
RMInit::logOut << "r_Polygon::r_Polygon(" << init << ") the init string has to contain valid r_Point definitions" << std::endl;
throw r_Error(POLYGONWRONGINITSTRING);
}
pointLen = endPos - startPos;
if (pointLen >= POINTBUFFERLEN)
{
TALK( "r_Polygon::r_Polygon(" << init << ") the definition of one r_Point is too long, only 2 dimensions are allowed" );
RMInit::logOut << "r_Polygon::r_Polygon(" << init << ") the definition of one r_Point is too long, only 2 dimensions are allowed" << std::endl;
throw r_Error(POLYGONWRONGINITSTRING);
}
memset(pointBuffer, 0, POINTBUFFERLEN);
strncpy(pointBuffer, startPos, pointLen + 1);
addPoint(r_Point(pointBuffer));
startPos = index(endPos, '[');
// was an endless loop with break, changed it to a 'nice' loop -- PB 2003-sep-12
// if (startPos == NULL)
// break;
} while( startPos != NULL);
LEAVE( "r_Polygon::r_Polygon" );
}
r_Polygon::r_Polygon(r_Range x, r_Range y) : closed(false), firstPointSet(true)
{
firstPoint = r_Point(x, y);
currPoint = firstPoint;
}
r_Polygon::r_Polygon() : closed(false), firstPointSet(false)
{
}
r_Polygon::r_Polygon(const r_Polygon& old)
{
closed=old.closed;
firstPointSet=old.firstPointSet;
firstPoint=old.firstPoint;
currPoint=old.currPoint;
edges=old.edges;
}
r_Polygon&
r_Polygon::operator=(const r_Polygon& old)
{
if(this!=&old)
{
closed=old.closed;
firstPointSet=old.firstPointSet;
firstPoint=old.firstPoint;
currPoint=old.currPoint;
edges=old.edges;
}
return *this;
}
void
r_Polygon::addPoint(const r_Point& newPoint) throw(r_Error)
{
if (newPoint.dimension() != polyPointDim)
{
RMInit::logOut << "r_Polygon::addPoint(" << newPoint << ") only " << polyPointDim << " dimensional r_Points allowed" << std::endl;
throw r_Error(POLYGONWRONGPOINTDIMENSION);
}
if(closed)
{
RMInit::logOut << "r_Polygon::addPoint(" << newPoint << ") polygon closed" << std::endl;
throw r_Error(r_Error::r_Error_General);
}
if(firstPointSet)
{
// add an edge from currentPoint to newPoint
edges.push_back(r_Edge(currPoint, newPoint));
currPoint = newPoint;
//check if we have the first point
checkFistPoint();
}
else
{
firstPoint = newPoint;
currPoint = newPoint;
firstPointSet = true;
}
}
void
r_Polygon::addPointXY(r_Range x, r_Range y) throw(r_Error)
{
r_Point newPoint(x, y);
addPoint(newPoint);
}
void
r_Polygon::close()
{
if (closed)
return;
// add the final edge from currentPoint to firstPoint
edges.push_back(r_Edge(currPoint, firstPoint));
closed = true;
}
const std::vector&
r_Polygon::getEdges() const throw(r_Error)
{
// if the polygon is not closed we raise an exception.
if(!closed)
{
// TO DO: This should be an internal error sometimes.
RMInit::logOut << "r_Polygon::getEdges() polygon opened" << std::endl;
throw(r_Error(r_Error::r_Error_General));
}
return edges;
}
std::vector
r_Polygon::getPoints() const throw(r_Error)
{
if(!closed)
{
// TO DO: This should be an internal error sometimes.
RMInit::logOut << "r_Polygon::getPoints() polygon opened" << std::endl;
throw(r_Error(r_Error::r_Error_General));
}
std::vector retVal;
std::vector::const_iterator iter, iterEnd;
for(iter = edges.begin(), iterEnd=edges.end(); iter != iterEnd; ++iter)
{
retVal.push_back((*iter).getStart());
}
return retVal;
}
r_Polygon::r_Polygon_Type
r_Polygon::detectPolygonType() const throw(r_Error)
{
const unsigned int minNoEdges=3;
std::vector points=getPoints();
unsigned int i=0,j=0, k=0, n=0;
unsigned int flag = 0;
double z;
n=points.size();
//check if is at least a triangle
if (n < minNoEdges)
return r_Polygon::UNKNOWN;
//check if is a polygon in 2D
for (i=0; i 0)
flag |= 2;
if (flag == 3)
return r_Polygon::CONCAVE;
}
if (flag != 0)
return r_Polygon::CONVEX;
else
return r_Polygon::UNKNOWN;
}
void
r_Polygon::checkFistPoint()
{
if (firstPoint == currPoint)
closed = true;
else
closed = false;
}
void
r_Polygon::print_status( std::ostream& s ) const
{
std::vector::const_iterator iter = edges.begin();
std::vector::const_iterator iterEnd = edges.end();
s << "{";
while (iter!=iterEnd)
{
iter->print_status(s);
s << ", ";
++iter;
}
s << "} ";
if (closed)
s << "closed";
else //does not matter, because open will crash ...
s << "opened";
}
std::ostream& operator<<( std::ostream& s, const r_Polygon& d )
{
d.print_status( s );
return s;
}
void
r_Polygon::fillMArray( r_GMarray& myArray, bool fillInside, const std::string& bgr ) const throw(r_Error)
{
if(!closed)
{
// TO DO: This should be an internal error sometimes.
RMInit::logOut << "r_Polygon::fillMArray(...) polygon opened" << std::endl;
throw(r_Error(r_Error::r_Error_General));
}
// This code currently is not optimised. For every scanline all edges are checked.
// Normally you would keep a list of currently relevant edges and manage deletes
// and additions to this table by presorting. Anyway, focus was on correctnes
// for the time being and even that was not really easy.
// all edges of the poly except for horizontal ones.
std::set sortedEdges;
// the X values of the edges in the current scanline.
std::vector currXVals;
// just to save some typing.
r_Minterval myDom = myArray.spatial_domain();
unsigned long typeSize = myArray.get_type_length();
// build sortedEdges (sorted after startPoint y, startPoint x).
std::vector::const_iterator iter, iterEnd;
for(iter = edges.begin(), iterEnd=edges.end(); iter != iterEnd; ++iter)
sortedEdges.insert(*iter);
// now the actual filling is done. Remember that we only draw the outside!
// we iterate through the whole y range
for(r_Range y = myDom[1].low(); y <= myDom[1].high(); y++)
{
// update the currXVals by iterating through all edges.
for(std::multiset::const_iterator itera = sortedEdges.begin(); itera != sortedEdges.end(); itera++)
{
if( (*itera).getStart()[1] <= y && y <= (*itera).getEnd()[1] ||
(*itera).getEnd()[1] <= y && y <= (*itera).getStart()[1] )
{
currXVals.push_back((*itera).getCurrX(y));
}
}
// sort them.
sort(currXVals.begin(), currXVals.end());
// currently we can only draw concave polygons anyway (see below).
// So this is heavily simplified! Check version 1.3 for a
// blueprint of how it should look like.
if(fillInside)
{
if(currXVals.size() >= 1)
{
// currXVals is sorted, so just draw from the first to the last.
eraseLine(rint(currXVals.front()), rint(currXVals.back()), y, myArray, bgr);
}
}
else
{
// currXVals is sorted, so just draw from low to the first and
// from the last to high.
// This only works correctly if the polygon is clipped
// to the area of the image.
if(currXVals.size() >= 1)
{
eraseLine(myDom[0].low(), rint(currXVals.front())-1.0, y, myArray, bgr);
eraseLine(rint(currXVals.back())+1.0, myDom[0].high(), y, myArray, bgr);
}
else
{
eraseLine(myDom[0].low(), myDom[0].high(), y, myArray, bgr);
}
}
// Note:
// Couldn't get it working with an even/odd inside rule. Reason: If you
// want to do this correctly you have to classify two edges meeting into
// different categories to decide if they have to be put into currXVals
// once or twice. Otherwise you get problems. With this simplified version
// only convex polygons are filled correctly!
// delete the old currXVals
currXVals.clear();
}
}
r_Minterval
r_Polygon::getBoundingBox() const throw(r_Error)
{
if(!closed)
{
// TO DO: This should be an internal error sometimes.
TALK( "r_Polygon::getBoundingBox() polygon opened" );
RMInit::logOut << "r_Polygon::getBoundingBox() polygon opened" << std::endl;
throw(r_Error(r_Error::r_Error_General));
}
r_Minterval retVal(2);
r_Range minX = LONG_MAX;
r_Range maxX = LONG_MIN;
r_Range minY = LONG_MAX;
r_Range maxY = LONG_MIN;
r_Range currMinX=0, currMaxX=0, currMinY=0, currMaxY=0;
std::vector::const_iterator iter, iterEnd;
for(iter = edges.begin(), iterEnd=edges.end(); iter != iterEnd; ++iter)
{
currMinX = (*iter).getStart()[0] < (*iter).getEnd()[0] ? (*iter).getStart()[0] : (*iter).getEnd()[0];
currMaxX = (*iter).getStart()[0] > (*iter).getEnd()[0] ? (*iter).getStart()[0] : (*iter).getEnd()[0];
currMinY = (*iter).getStart()[1] < (*iter).getEnd()[1] ? (*iter).getStart()[1] : (*iter).getEnd()[1];
currMaxY = (*iter).getStart()[1] > (*iter).getEnd()[1] ? (*iter).getStart()[1] : (*iter).getEnd()[1];
minX = currMinX < minX ? currMinX : minX;
maxX = currMaxX > maxX ? currMaxX : maxX;
minY = currMinY < minY ? currMinY : minY;
maxY = currMaxY > maxY ? currMaxY : maxY;
}
retVal << r_Sinterval(minX, maxX) << r_Sinterval(minY, maxY);
return retVal;
}
void
r_Polygon::clip(const r_Minterval& clipDom) throw(r_Error)
{
if(!closed)
{
// TO DO: This should be an internal error sometimes.
TALK( "r_Polygon::getBoundingBox() polygon opened" );
RMInit::logOut << "r_Polygon::getBoundingBox() polygon opened" << std::endl;
throw(r_Error(r_Error::r_Error_General));
}
std::vector pointList;
// Disjunct polygons used to crash the program.
// check if the bounding box of the polygon overlaps the clipDom
if(clipDom.intersects_with(getBoundingBox()))
{
// We just clip all 4 edges
for(int s = r_Polygon::top; s <= r_Polygon::right; ++s)
{
pointList=clip1Side(clipDom, (r_Polygon::Side)s);
if(pointList.empty()) // do we have intersection points?
{
// return a polygon with one line only. This should delete everything.
pointList.push_back(r_Point(clipDom[0].low(), clipDom[1].low()));
pointList.push_back(r_Point(clipDom[0].low(), clipDom[1].high()));
}
fromPoints(pointList);
pointList.clear();
}
}
else
{
// return a polygon with one line only. This should delete everything.
pointList.push_back(r_Point(clipDom[0].low(), clipDom[1].low()));
pointList.push_back(r_Point(clipDom[0].low(), clipDom[1].high()));
fromPoints(pointList);
}
}
void
r_Polygon::scale(const r_Point& origin, const r_Minterval& mddDom,
const r_Minterval& clipDom, const double& scaleFactor) throw(r_Error)
{
r_Dimension dim = origin.dimension();
std::vector oldPoints = getPoints();
std::vector::const_iterator iter = oldPoints.begin();
std::vector::const_iterator iterEnd = oldPoints.end();
std::vector newPoints;
r_Point tmpPoint(dim);
r_Range coord=0;
// iterate through all points
while( iter != iterEnd )
{
// currently only 2-D, but who knows
for (int i=0; i::get_scaled_image().
coord = coord + mddDom[i].high() - clipDom[i].high();
tmpPoint[i] = coord;
}
newPoints.push_back(tmpPoint);
++iter;
}
fromPoints(newPoints);
}
void
r_Polygon::scale(const r_Point& origin, const double& scaleFactor) throw(r_Error)
{
r_Dimension dim = origin.dimension();
std::vector oldPoints = getPoints();
std::vector::const_iterator iter = oldPoints.begin();
std::vector::const_iterator iterEnd = oldPoints.end();
std::vector newPoints;
r_Point tmpPoint(dim);
r_Range coord=0;
//std::cout << "Polygon bounding box " << getBoundingBox() << std::endl;
// iterate through all points
while( iter != iterEnd )
{
// currently only 2-D, but who knows
for (int i=0; i oldPoints = getPoints();
std::vector::const_iterator iter = oldPoints.begin();
std::vector::const_iterator iterEnd = oldPoints.end();
std::vector newPoints;
r_Point tmpPoint(dim);
r_Range y=0;
// iterate through all points
while( iter != iterEnd )
{
y = mddDom[1].high() - (*iter)[1];
tmpPoint = (*iter);
tmpPoint[1] = y;
newPoints.push_back(tmpPoint);
++iter;
}
fromPoints(newPoints);
}
void
r_Polygon::fromPoints(const std::vector& newPoints) throw(r_Error)
{
std::vector::const_iterator iter, iterEnd;
if(newPoints.empty())
{
TALK( "r_Polygon::fromPoinst(....) newPoints is empty" );
RMInit::logOut << "r_Polygon::fromPoinst(....) newPoints is empty" << std::endl;
throw r_Error(r_Error::r_Error_General);
}
iter = newPoints.begin();
iterEnd = newPoints.end();
edges.clear();
firstPoint = *iter;
currPoint = *iter;
while( ++iter != iterEnd )
{
edges.push_back(r_Edge(currPoint, *iter));
currPoint = *iter;
}
edges.push_back(r_Edge(currPoint, firstPoint));
closed = true;
}
void
r_Polygon::eraseLine( r_Range x1, r_Range x2, r_Range y, r_GMarray& myArray, const std::string& bgr ) const throw(r_Error)
{
// Do nothing in that case (may happen due to rounding problems)
if(x2 < x1)
return;
r_Minterval eraseDom(2);
eraseDom << r_Sinterval(x1, x2) << r_Sinterval(y, y);
r_Minterval arrayDom = myArray.spatial_domain();
r_Bytes typeSize = myArray.get_type_length();
r_Bytes bgrSize = bgr.size();
const char *bgrContent=bgr.c_str();
char *currCell=NULL;
// Grrr. In RasDaMan the y are stored close together. So the whole fillPolygon
// should have been organised by x instead of y. Well, for now I just use an
// r_Miter here.
r_Miter eraseIter( &eraseDom, &arrayDom, typeSize, myArray.get_array() );
while( !eraseIter.isDone() )
{
currCell = eraseIter.nextCell();
// FIXME This potentially wont work for all types. I just set every byte to 0.
if(bgr.empty())
memset(currCell, 0, typeSize);
else
{
if( typeSize != bgrSize)
throw r_Error( r_Error::r_Error_TypeInvalid );
memmove(currCell, bgrContent, bgrSize);
}
}
}
std::vector
r_Polygon::clip1Side(const r_Minterval& b, r_Polygon::Side s)
{
// This routine clips a polygon against one (endless) edge of a bounding box.
// It is an implementation of the Sutherland-Hodgman algorithm geared more
// towards readability than efficiency. The algorithm classifies edges into
// four cases:
// Case 1: both ends are inside. Then add the end point to the list of points.
// Case 2: start inside, end outside. And only the intersection of the
// bounding box edge and the polygon edge.
// Case 3: completely outside. Add no points.
// Case 4: start outside, end inside. Add end and the intersection of the
// bounding box edge and the polygon edge.
std::vector out;
std::vector::const_iterator iter, iterEnd;
// iterate through all edges
for(iter = edges.begin(),iterEnd=edges.end(); iter != iterEnd; ++iter)
{
if(inside(b, (*iter).getEnd(), s))
{
if(inside(b, (*iter).getStart(), s))
{
// case 1
out.push_back((*iter).getEnd());
}
else
{
// case 4
out.push_back(intersect(b, *iter, s));
out.push_back((*iter).getEnd());
}
}
else
{
if(inside(b, (*iter).getStart(), s))
{
// case 2
out.push_back(intersect(b, *iter, s));
}
}
// do nothing for case 3
}
return out;
}
bool
r_Polygon::inside(const r_Minterval& b, const r_Point& p, r_Polygon::Side s)
{
switch(s)
{
case top:
return p[1] <= b[1].high();
case bottom:
return p[1] >= b[1].low();
case right:
return p[0] <= b[0].high();
case left:
return p[0] >= b[0].low();
default:
return false;
}
}
r_Point
r_Polygon::intersect(const r_Minterval& b, const r_Edge& e, r_Polygon::Side s)
{
switch(s)
{
case top:
return r_Point( e.getCurrX(b[1].high()), b[1].high() );
case bottom:
return r_Point( e.getCurrX(b[1].low()), b[1].low() );
case right:
return r_Point( b[0].high(), e.getCurrY(b[0].high()) );
default: //case left:
return r_Point( b[0].low(), e.getCurrY(b[0].low()) );
}
}
r_Point
r_Polygon::getMiddle() const throw(r_Error)
{
// Note that the summing up done here is a bit risky since overflows
// give incorrect results.
r_Point retVal(2);
double xSum = 0;
double ySum = 0;
int pointCount = 0;
std::vector myPoints = getPoints();
std::vector::const_iterator iter = myPoints.begin();
std::vector::const_iterator iterEnd = myPoints.end();
while( iter != iterEnd )
{
ySum += (*iter)[1];
xSum += (*iter)[0];
++pointCount;
++iter;
}
retVal[0] = rint((xSum / pointCount));
retVal[1] = rint((ySum / pointCount));
return retVal;
}
void
r_Polygon::shrinkPoly(int pixelCount) throw(r_Error)
{
// Ok now, we move all points towards the middle. Since we store edges we
// have to use this somewhat clumsy form with using points in between.
// Note that there is quite a bit of potential for error (e.g. points
// coinciding after moving them), Anyway, this was programmed as a quick
// hack for a problem at BLVA.
r_Point middle = getMiddle();
r_Dimension dim = middle.dimension();
std::vector oldPoints = getPoints();
std::vector::const_iterator iter = oldPoints.begin();
std::vector::const_iterator iterEnd = oldPoints.end();
std::vector newPoints;
r_Point tmpPoint(dim);
r_Range coord=0;
// iterate through all points
while( iter != iterEnd )
{
// currently only 2-D, but who knows
for (int i=0; i middle[i])
{
coord = coord - pixelCount;
}
else
{
if(coord < middle[i])
{
coord = coord + pixelCount;
}
}
tmpPoint[i] = coord;
}
newPoints.push_back(tmpPoint);
++iter;
}
fromPoints(newPoints);
}
int
r_Polygon::countEdges() const
{
return edges.size();
}
int
r_Polygon::countHorizontalEdges() const
{
int counter = 0;
std::vector::const_iterator iter = edges.begin();
for(int i=0;iisHorizontal() ? 1:0;
}
return counter;
}