diff options
author | Constantin Jucovschi <cj@ubuntu.localdomain> | 2009-04-24 07:20:22 -0400 |
---|---|---|
committer | Constantin Jucovschi <cj@ubuntu.localdomain> | 2009-04-24 07:20:22 -0400 |
commit | 8f27e65bddd7d4b8515ce620fb485fdd78fcdf89 (patch) | |
tree | bd328a4dd4f92d32202241b5e3a7f36177792c5f /rasodmg/polycutout.cc | |
download | rasdaman-upstream-8f27e65bddd7d4b8515ce620fb485fdd78fcdf89.tar.gz rasdaman-upstream-8f27e65bddd7d4b8515ce620fb485fdd78fcdf89.tar.xz rasdaman-upstream-8f27e65bddd7d4b8515ce620fb485fdd78fcdf89.zip |
Initial commitv8.0
Diffstat (limited to 'rasodmg/polycutout.cc')
-rw-r--r-- | rasodmg/polycutout.cc | 709 |
1 files changed, 709 insertions, 0 deletions
diff --git a/rasodmg/polycutout.cc b/rasodmg/polycutout.cc new file mode 100644 index 0000000..78f3d75 --- /dev/null +++ b/rasodmg/polycutout.cc @@ -0,0 +1,709 @@ +/* +* 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 <http://www.gnu.org/licenses/>. +* +* Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Peter Baumann / +rasdaman GmbH. +* +* For more information please see <http://www.rasdaman.org> +* or contact Peter Baumann via <baumann@rasdaman.com>. +/ +/** + * SOURCE: polycutout.cc + * + * MODULE: rasodmg + * CLASS: r_PolygonCutOut + * + * COMMENTS: + * + * This class is intended to perform the polygon cut out operation + * Except for the r_PolygonCutOut class, everything else is just + * for internal use + * Attention: this version works good (I hope :-) but has just a tiny geometrical problem: + * as usual in discrete space, lines are build up by little horiontal segments. + * Edges are considered to be inside, except the first little segment of a line. + * This is only for avoiding a bigger problem, I will develop a better algorythm +*/ + +/* +This module is intended to perform polygon cut out operation. It supports the exact ESRI specifications +regarding "clean" polygons, this means: not self intersecting, closed, if you walk from an point to the +next of edge the inside is on the right. It also supports multiple rings, but the outer one has to have +the "inside" inside. +Don't forget that the rasdaman coordinates are with x left->right and y top->down, or at least +I have considered them so because of the usual coordinates of tiffs +*/ + +static const char rcsid[] = "@(#)rasodmg, r_PolygonCutOut: $Id: polycutout.cc,v 1.11 2002/10/09 09:58:05 hoefner Exp $"; + +#include "rasodmg/polycutout.hh" +#include "rasodmg/marray.hh" +#include "raslib/miter.hh" + +#include <math.h> + +r_SegmentIterator::r_SegmentIterator(r_Point &s,r_Point &e) + { + start = s; + end = e; + reset(); + } + +void r_SegmentIterator::reset() + { + dx = end[0] - start[0]; + dy = end[1] - start[1]; + + if( dx >= 0) + { if( dy >= 0 ) + { if( dx >= dy ) { cadran = 0;} + else { cadran = 1; swap(dx,dy);} + } + else + { dy = -dy; + if( dx >= dy ) { cadran = 7;} + else { cadran = 6; swap(dx,dy);} + } + } + else + { dx = -dx; + if( dy >= 0 ) + { if( dx >= dy ) { cadran = 3;} + else { cadran = 2; swap(dx,dy);} + } + else + { dy = -dy; + if( dx >= dy ) { cadran = 4;} + else { cadran = 5; swap(dx,dy);} + } + } + cx = 0; + cy = 0; + beta = 0; +// std::cout<<"reset: ["<<start<<','<<end<<"] cadran="<<cadran<<" dx="<<dx<<" dy="<<dy<<std::endl; + } + +r_Point r_SegmentIterator::next() + { + r_Point resultat = createCurrentPoint(); + //prepare next + cx++; + beta += dy; + if( (beta<<1) > dx ) { beta -= dx; cy++;} + return resultat; + } + +bool r_SegmentIterator::hasMore() + { +// std::cout<<"hasMore: cx="<<cx<<" lim="<<dx+1<<std::endl; + return (cx <= dx ) ? true : false; + } + +int r_SegmentIterator::cosFunc() // limited use,1000 * cos(alfa) + { // use local variables + double ldx = end[0] - start[0]; + double ldy = end[1] - start[1]; + + double num = ldx*ldx + ldy*ldy; + if(num == 0) return 0; + + return (1000 * ldx)/sqrt(num); + } + +void r_SegmentIterator::swap(r_Range &x, r_Range &y) + { + r_Range temp = x; + x=y; + y=temp; + } + +r_Point r_SegmentIterator::createCurrentPoint() + { +// std::cout<<" cCP: "<<cx<<':'<<cy<<std::endl; + r_Range rcx = cx; + r_Range rcy = cy; + switch(cadran) + { case 0: break; + case 1: swap(rcx,rcy); break; + case 2: swap(rcx,rcy); rcx=-rcx; break; + case 3: rcx=-rcx; break; + case 4: rcx=-rcx; rcy=-rcy; break; + case 5: swap(rcx,rcy); rcx=-rcx; rcy=-rcy; break; + case 6: swap(rcx,rcy); rcy=-rcy; break; + case 7: rcy=-rcy; break; + } + //std::cout<<"cCP: cx="<<cx<<" cy="<<cy<<" rx="<<rcx<<" ry="<<rcy<<" cadran="<<cadran<<std::endl; + return r_Point(start[0] + rcx, start[1] + rcy); + } + + + +r_Line::r_Line() + { + a = b = c = 0; + } + +r_Line::r_Line(double nA,double nB,double nC) + { + a = nA; b = nB; c = nC; + } + +r_Line::r_Line(r_Point &s,r_Point &e) + { + a = e[1] - s[1]; + b = s[0] - e[0]; + c = (double)e[0]*(double)s[1] - (double)s[0]*(double)e[1]; + } + +double r_Line::getA() { return a; } +double r_Line::getB() { return b; } +double r_Line::getC() { return c; } + +float r_Line::ecuatia(r_Point &p) + { + return a * p[0] + b * p[1] + c; + } + +ostream& operator<<(ostream& os,r_Line &l) + { + return os<<"{L"<<l.a<<','<<l.b<<','<<l.c<<'}'; + } + + + + +r_PolygonCutOut::r_PolygonCutOut() + { + table = 0; + usedCount = 0; + mArray = 0; + + imgWidth = 0; + imgHeight = 0; + imgX = 0; + imgY = 0; + tableWidth = 0; + tableHeight = 0; + } +r_PolygonCutOut::~r_PolygonCutOut() + { + clearTables(); + } + +void r_PolygonCutOut::setMArray(r_GMarray& myArray) + { + mArray = &myArray; + + r_Minterval currDomain; + currDomain = mArray->spatial_domain(); + + imgX = currDomain[0].low(); + imgY = currDomain[1].low(); + + imgWidth = currDomain[0].get_extent(); + imgHeight = currDomain[1].get_extent(); + } + +void r_PolygonCutOut::addPolygon(const r_Polygon &p) + { + polygons.push_back(p); + } + +bool r_PolygonCutOut::compute() + { +//std::cout<<"PCO: compute in"<<std::endl; + if(!imgWidth || !imgHeight || polygons.size()==0) return false; + +//std::cout<<"PCO: initTable - in"<<std::endl; + initTable(); +//std::cout<<"PCO: initTable - out"<<std::endl; + + list<r_Polygon>::iterator polyIter = polygons.begin(); + + for(int i = 0; i < polygons.size(); i++, polyIter++) + { + r_Polygon &currPolygon = *polyIter; + + const std::vector<r_Edge>& edges = currPolygon.getEdges(); + std::vector<r_Edge>::const_iterator edgeIterator = edges.begin(); + + for(int j = 0; j < edges.size(); j++,edgeIterator++) + { + r_Point start = edgeIterator->getStart(); + r_Point end = edgeIterator->getEnd(); + + if(start == end) continue; //avoid such edges, causes problems + + int inside = -computeInside(start,end); // minus because rasdaman coord are flipped + // we have to see if this is correct like this or not + + //std::cout<<"compute edge "<<j<<" ["<<start<<"-"<<end<<"] inside="<<inside<<std::endl; + if(inside == 0 ) computeOneHorSegment(start,end); + else computeOneSegment(start,end,inside); + } + } + for(int line=0;line<tableHeight;line++) + { ordonateLine(line); + minimizeLine(line); + } +//print(); +// std::cout<<"PCO: compute out"<<std::endl; + return true; + } + +bool r_PolygonCutOut::initTable() + { + clearTables(); + tableWidth = computeTableWidth(); + tableHeight= imgHeight; + +//std::cout<<"PCO: initTable: "<<tableWidth<<','<<tableHeight<<std::endl; + table = new TablePoint[tableWidth*tableHeight]; + usedCount = new int[tableHeight]; + + // if(table==NULL || usedCount==NULL) std::cout<<"sunt nule"<<std::endl; + + for(int i=0;i<tableWidth*tableHeight;i++) + { + table[i].x =-1; + table[i].inside= 0; + table[i].cosFunc= 0; + } + + for(int j=0;j<tableHeight;j++) + usedCount[j]=0; + + return true; + } + +void r_PolygonCutOut::clearTables() + { + if(table) delete[] table; + table = NULL; + if(usedCount) delete[] usedCount; + usedCount =NULL; + } +int r_PolygonCutOut::computeTableWidth() + { + // once it was: tableWidth = poly->getCountNodes() + poly->countHorizontalEdges() + 2; + + list<r_Polygon>::iterator iter=polygons.begin(); + + int nodes = 0; + int horizEdges = 0; + + for(int i=0;i<polygons.size();i++, iter++) + { + r_Polygon &polygon = *iter; + + nodes += polygon.countEdges(); + horizEdges += polygon.countHorizontalEdges(); + } + return nodes + horizEdges; // niet + 2; + } + +void r_PolygonCutOut::ordonateLine(int line) + { + int count=usedCount[line]; + TablePoint *tLine = &getTP(line,0); + bool change=false; + do{ + change=false; + for(int i=0;i<count-1;i++) + { bool swap=false; + + if(tLine[i].x > tLine[i+1].x) swap=true; + if(tLine[i].x == tLine[i+1].x && + tLine[i].inside != tLine[i+1].inside && + tLine[i].cosFunc > tLine[i+1].cosFunc) + swap=true; + if(swap) + { TablePoint temp = tLine[i]; + tLine[i] = tLine[i+1]; + tLine[i+1] = temp; + change = true; + } + } + }while(change); + + for(int j=0;j<count-1;j++) + { if(tLine[j] == tLine[j+1]) + tLine[j].inside=0; + } + + } + +void r_PolygonCutOut::minimizeLine(int line) + { + int count=usedCount[line]; + if(count <= 2) return; //nothing to do + + int inside = 0; + TablePoint *tLine = &getTP(line,0); + + int level=0; + for(int i=0;i<count-1;i++) + { + if(level==0) { inside=tLine[i].inside; + level++; + continue; + } + + if(tLine[i].inside==inside) level++; + else level--; + + if(level!=0) tLine[i].inside=0; // means disabled + } + + int s,d; + for(s=0,d=0;s<count;s++) + { + if(tLine[s].inside==0) continue; + if(s!=d) tLine[d]=tLine[s]; + d++; + } + usedCount[line]=d; + } + +void r_PolygonCutOut::computeOneSegment(r_Point start, r_Point end, int inside) + { + int lastLine = -1; + int lastX = -1; + r_SegmentIterator segiter(start,end); + + int cosFunc=segiter.cosFunc(); + int cosNow = cosFunc; // first point: cosFunc, lastPoint -cosFunc, rest 0 + + bool firstLine = true; + bool lastPoint = false; + // this firstLine-lastPoint are here only top be sure that the first and the last point are in the table + // and not replaced by neighbours + while(segiter.hasMore()) + { r_Point cp = segiter.next(); + if(segiter.hasMore()==false) // last point + { + cosNow=-cosFunc; + lastPoint = true; + } + + r_Range coordX = cp[0]; + r_Range coordY = cp[1]; + + int tableLine = coordY - imgY; + + if(tableLine>=0 && tableLine < tableHeight) + { + if(tableLine != lastLine) + { addPoint(tableLine,coordX,inside,cosNow); + + if(lastLine != -1) firstLine = false; + } + else if(firstLine == false) + { + if(inside > 0 && coordX < lastX) replacePoint(tableLine,coordX,inside,cosNow); + else if(inside < 0 && coordX > lastX) replacePoint(tableLine,coordX,inside,cosNow); + else if(lastPoint) replacePoint(tableLine,coordX,inside,cosNow); + } + } + + cosNow = 0; + lastLine = tableLine; + lastX = coordX; + } + } + +void r_PolygonCutOut::computeOneHorSegment(r_Point start, r_Point end) + { + r_SegmentIterator segiter(start,end); + int cosFunc=segiter.cosFunc(); + + r_Range startX = start[0]; + r_Range endX = end[0]; + r_Range commonY = start[1]; + + int tableLine = commonY - imgY; + + if(tableLine>=0 && tableLine < tableHeight) + { + if(startX <= endX) + { addPoint(tableLine,startX, 1, cosFunc); + addPoint(tableLine,endX, -1,-cosFunc); + } + else { // in this case cosFunc < 0 + addPoint(tableLine,endX, 1,-cosFunc); + addPoint(tableLine,startX,-1, cosFunc); + } + } + } + +int r_PolygonCutOut::computeInside(r_Point start, r_Point end) + { + r_Line line(start,end); + r_Point control(start[0]+1,start[1]); + float f = line.ecuatia(control); + if(f==0) return 0; + if(f<0 ) return -1; + return 1; + } + + +r_PolygonCutOut::TablePoint& r_PolygonCutOut::getTP(int line, int column) + { + TablePoint *tp = table + (line * tableWidth + column); + return *tp; + } + +void r_PolygonCutOut::print(int onlyLine) + { + std::cout<<"r_PolygonCutOut::print: "<<tableWidth<<':'<<tableHeight<<std::endl; + + if(onlyLine==-1) // print all + { + for(int i=0;i<tableHeight;i++) + printLine(i); + } + else + { + if(0<=onlyLine && onlyLine <tableHeight) printLine(onlyLine); + else std::cout<<"Line "<<onlyLine<<" doesn't exist"<<std::endl; + } + } + +void r_PolygonCutOut::printLine(r_Range line) + { + int count=usedCount[line]; + TablePoint *tp = &getTP(line,0); + + std::cout<<"line "<<line<<'('<<count<<") "; + for(int i=0;i<count;i++) + { char inside='0'; + if(tp[i].inside==-1) inside='-'; + if(tp[i].inside==1 ) inside='+'; + + std::cout<<'['<<tp[i].x<<','<<inside<<','<<tp[i].cosFunc<<']'; + } + std::cout<<std::endl; + } + +void r_PolygonCutOut::addPoint(int tableLine,r_Range coordX, int inside,int cosFunc) + { + int which = usedCount[tableLine]; + TablePoint *tp = &getTP(tableLine,which); + + tp->x = coordX; + tp->inside = inside; + tp->cosFunc = cosFunc; + usedCount[tableLine]++; + } + +void r_PolygonCutOut::replacePoint(int tableLine,r_Range coordX,int inside, int cosFunc) + { + int which = usedCount[tableLine]-1; + TablePoint *tp = &getTP(tableLine,which); + tp->x = coordX; + tp->inside = inside; + tp->cosFunc = cosFunc; + } + +bool r_PolygonCutOut::TablePoint::operator==(r_PolygonCutOut::TablePoint &tp) + { return (x==tp.x && inside==tp.inside) ? true: false; + } + + +// debug only int lastline = -1; +void +r_PolygonCutOut::eraseLine( r_Range x1, r_Range x2, r_Range y, const string& bgr ) throw(r_Error) +{ + // can heapen if we are outside the domain of the array + if(x2 < x1) return; + + //if(lastline != y && x1 != 0) + // std::cout<<"Erasing line: y="<<y<<" ["<<x1<<","<<x2<<"]"<<std::endl; + //lastline = y; + + r_Minterval eraseDom(2); + eraseDom << r_Sinterval(x1, x2) << r_Sinterval(y, y); + r_Minterval arrayDom = mArray->spatial_domain(); + r_Bytes typeSize = mArray->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, mArray->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); + } + } +} + + +bool +r_PolygonCutOut::fillMArrayOutside(const string& bgr ) throw(r_Error) +{ + if( compute() == false ) return false; + + r_Minterval currDomain = mArray->spatial_domain(); + + r_Range imgMinX = imgX; + r_Range imgMaxX = currDomain[0].high(); + + for(int lineY = 0; lineY < tableHeight ; lineY++) + { + r_Range coordY = imgY + lineY; + + int countUsed = usedCount[lineY]; + + if(countUsed==0) + { // polygon didn't touch this line, so all line is outside + eraseLine(imgMinX, imgMaxX, coordY, bgr); + continue; + } + + TablePoint *tp = &getTP(lineY, 0); + + r_Range xLeft = imgMinX-1; + r_Range xRight = imgMaxX+1; + bool erase = false; + + if( tp[countUsed-1].x < imgMinX) + { // the whole polygon is on the left of this line so we have to erase the whole line + eraseLine( imgMinX, imgMaxX, coordY, bgr); + continue; + } + + if( tp[0].x > imgMaxX) + { // the whole polygon is on the right of this line so we have to erase the whole line + eraseLine( imgMinX, imgMaxX, coordY, bgr); + continue; + } + + for(int i=0;i<countUsed; i++, tp++) + { + if(tp->x < imgMinX) continue; + + if(tp->inside == -1) + { // inside is <-, so outside is -> + xLeft = tp->x; + if(xLeft > imgMaxX) break; + + if(i == countUsed-1) + { xRight = imgMaxX+1; erase = true; + } + } + + if(tp->inside == 1) + { // inside is ->, so outside is <- + xRight = tp->x > imgMaxX ? imgMaxX+1 : tp->x; + erase = true; + } + + if(erase) + { eraseLine( xLeft+1, xRight-1, coordY, bgr); + erase = false; + } + } + + } + clearTables(); + return true; + } + +/* + Attention, this two functions look very similar, but with have significant differences + One of them is the absence of some -1 or +1 in the next function. This is important + because we consider the edges to be "inside", and the edges are in the table, so + erasing the inside is erase [start,end] and erasing the outside is erase [start+1,end-1] +*/ + +bool +r_PolygonCutOut::fillMArrayInside(const string& bgr ) throw(r_Error) +{ + if( compute() == false ) return false; + + r_Minterval currDomain = mArray->spatial_domain(); + + r_Range imgMinX = imgX; + r_Range imgMaxX = currDomain[0].high(); + + for(int lineY = 0; lineY < tableHeight ; lineY++) + { + r_Range coordY = imgY + lineY; + + int countUsed = usedCount[lineY]; + + if(countUsed==0) + { // polygon didn't touch this line, so all line is outside + continue; + } + + TablePoint *tp = &getTP(lineY, 0); + + r_Range xLeft = imgMinX; + r_Range xRight = imgMaxX; + bool erase = false; + + if( tp[countUsed-1].x < imgMinX) + { // the whole polygon is on the left of this line so we have to erase the whole line + continue; + } + + if( tp[0].x > imgMaxX) + { // the whole polygon is on the right of this line so we have to erase the whole line + continue; + } + + for(int i=0;i<countUsed; i++, tp++) + { + if(tp->x < imgMinX) continue; + + if(tp->inside == 1) + { // inside is ->, so outside is <- + xLeft = tp->x; + if(xLeft > imgMaxX) break; + + if(i == countUsed-1) + { xRight = imgMaxX; erase = true; + } + } + + if(tp->inside == -1) + { // inside is <-, so outside is -> + xRight = tp->x > imgMaxX ? imgMaxX : tp->x; + erase = true; + } + + if(erase) + { eraseLine( xLeft, xRight, coordY, bgr); + erase = false; + } + } + + } + clearTables(); + return true; + } |