/* * 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: 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 r_SegmentIterator::r_SegmentIterator(r_Point &s,r_Point &e) { start = s; end = e; reset(); } void r_SegmentIterator::reset() { dx = end - start; dy = end - start; 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: ["<spatial_domain(); imgX = currDomain.low(); imgY = currDomain.low(); imgWidth = currDomain.get_extent(); imgHeight = currDomain.get_extent(); } void r_PolygonCutOut::addPolygon(const r_Polygon &p) { polygons.push_back(p); } bool r_PolygonCutOut::compute() { //std::cout<<"PCO: compute in"<::iterator polyIter = polygons.begin(); for(int i = 0; i < polygons.size(); i++, polyIter++) { r_Polygon &currPolygon = *polyIter; const std::vector& edges = currPolygon.getEdges(); std::vector::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 "<getCountNodes() + poly->countHorizontalEdges() + 2; list::iterator iter=polygons.begin(); int nodes = 0; int horizEdges = 0; for(int i=0;i 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=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; r_Range endX = end; r_Range commonY = start; 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+1,start); 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: "<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="<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.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.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;ix < 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.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.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;ix < 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; }