/* ---- This file is part of SECONDO. Copyright (C) 2004, University in Hagen, Department of Computer Science, Database Systems for New Applications. SECONDO 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 2 of the License, or (at your option) any later version. SECONDO 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 SECONDO; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ---- //paragraph [1] Title: [{\Large \bf \begin {center}] [\end {center}}] //[TOC] [\tableofcontents] //[_] [\_] [1] Implementation of the Spatial Algebra February, 2003. Victor Teixeira de Almeida March-July, 2003. Zhiming Ding January, 2005 Leonardo Guerreiro Azevedo [TOC] 1 Overview This implementation file essentially contains the implementation of the classes ~Point~, ~Points~, ~Line~, and ~Region~ used in the Spatial Algebra. These classes respectively correspond to the memory representation for the type constructors ~point~, ~points~, ~line~, and ~region~. For more detailed information see SpatialAlgebra.h. 2 Defines and Includes */ #undef __TRACE__ //#define __TRACE__ cout << __FILE__ << "::" << __LINE__; #define __TRACE__ #include "Label.h" #include "../../Tools/Flob/Flob.h" #include "../../Tools/Flob/DbArray.h" #include "Algebra.h" #include "NestedList.h" #include "ListUtils.h" #include "GenericTC.h" #include "QueryProcessor.h" #include "StandardTypes.h" #include "SpatialAlgebra.h" #include "Algebras/FText/FTextAlgebra.h" #include "SecondoConfig.h" #include "AvlTree.h" #include "AVLSegment.h" #include "AlmostEqual.h" #include "Algebras/Relation-C++/RelationAlgebra.h" #include "RegionTools.h" #include "Symbols.h" #include "Stream.h" #include "RegionCreator.h" #include "StringUtils.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "../TopRel/Tree.h" #include "RobustSetOps.h" #include "DLine.h" #include "Segment.h" #include "DRM.h" #include "Disc.h" #include "Stack.h" #include "Algebras/Geoid/GeoDist.h" using namespace std; extern NestedList* nl; extern QueryProcessor* qp; typedef RegionCreator PRegionCreator; typedef robust::RealmChecker PRealmChecker; /* 3 Type investigation auxiliaries Within this algebra module, we have to handle with values of four different types: ~point~ and ~points~, ~line~ and ~region~. Later on we will examine nested list type descriptions. In particular, we are going to check whether they describe one of the four types just introduced. In order to simplify dealing with list expressions describing these types, we declare an enumeration, ~SpatialType~, containing the four types, and a function, ~SpatialTypeOfSymbol~, taking a nested list as argument and returning the corresponding ~SpatialType~ type name. */ enum SpatialType { stpoint, stpoints, stline, stregion, stbox, sterror, stsline, stcpoint }; SpatialType SpatialTypeOfSymbol( ListExpr symbol ) { if ( nl->AtomType( symbol ) == SymbolType ) { string s = nl->SymbolValue( symbol ); if ( s == Point::BasicType() ) return (stpoint); if ( s == Points::BasicType() ) return (stpoints); if ( s == Line::BasicType() ) return (stline); if ( s == Region::BasicType() ) return (stregion); if ( s == Rectangle<2>::BasicType() ) return (stbox); if ( s == SimpleLine::BasicType() ) return (stsline); if ( s == CPoint::BasicType() ) return (stcpoint); } return (sterror); } /* 9 Object Traversal functions These functions are utilities useful for traversing objects. They are basic functions to be called by the operations defined below. There are 6 combinations, pp, pl, pr, ll, lr, rr */ void SelectFirst_pl( const Points& P, const Line& L, object& obj, status& stat ) { P.SelectFirst(); L.SelectFirst(); Point p1(true); Point p2(true); HalfSegment hs; bool gotPt = P.GetPt( p1 ), gotHs = L.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); if( !gotPt && !gotHs ) { obj = none; stat = endboth; } else if( !gotPt ) { obj = second; stat = endfirst; } else if( !gotHs ) { obj = first; stat = endsecond; } else //both defined { stat = endnone; if( p1 < p2 ) obj = first; else if( p1 > p2 ) obj = second; else obj=both; } } void SelectNext_pl( const Points& P, const Line& L, object& obj, status& stat ) { // 1. get the current elements Point p1(true); Point p2(true); HalfSegment hs; bool gotPt = P.GetPt( p1 ), gotHs = L.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); //2. move the pointers if( !gotPt && !gotHs ) { //do nothing } else if( !gotPt ) { L.SelectNext(); gotHs = L.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); } else if( !gotHs ) { P.SelectNext(); gotPt = P.GetPt( p1 ); } else //both currently defined { if( p1 < p2 ) //then hs1 is the last output { P.SelectNext(); gotPt = P.GetPt( p1 ); } else if( p1 > p2 ) { L.SelectNext(); gotHs = L.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); } else { P.SelectNext(); gotPt = P.GetPt( p1 ); L.SelectNext(); gotHs = L.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); } } //3. generate the outputs if( !gotPt && !gotHs ) { obj = none; stat = endboth; } else if( !gotPt ) { obj = second; stat = endfirst; } else if( !gotHs ) { obj = first; stat = endsecond; } else //both defined { stat=endnone; if( p1 < p2 ) obj = first; else if( p1 > p2 ) obj = second; else obj = both; } } void SelectFirst_pr( const Points& P, const Region& R, object& obj, status& stat ) { P.SelectFirst(); R.SelectFirst(); Point p1(true); Point p2(true); HalfSegment hs; bool gotPt = P.GetPt( p1 ), gotHs = R.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); if( !gotPt && !gotHs ) { obj = none; stat = endboth; } else if( !gotPt ) { obj = second; stat = endfirst; } else if( !gotHs ) { obj = first; stat = endsecond; } else //both defined { stat = endnone; if( p1 < p2 ) obj = first; else if( p1 > p2 ) obj = second; else obj=both; } } void SelectNext_pr( const Points& P, const Region& R, object& obj, status& stat ) { // 1. get the current elements Point p1(true); Point p2(true); HalfSegment hs; bool gotPt = P.GetPt( p1 ), gotHs = R.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); //2. move the pointers if( !gotPt && !gotHs ) { //do nothing } else if( !gotPt ) { R.SelectNext(); gotHs = R.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); } else if( !gotHs ) { P.SelectNext(); gotPt = P.GetPt( p1 ); } else //both currently defined { if( p1 < p2 ) //then hs1 is the last output { P.SelectNext(); gotPt = P.GetPt( p1 ); } else if( p1 > p2 ) { R.SelectNext(); gotHs = R.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); } else { P.SelectNext(); gotPt = P.GetPt( p1 ); R.SelectNext(); gotHs = R.GetHs( hs ); if( gotHs ) p2 = hs.GetDomPoint(); } } //3. generate the outputs if( !gotPt && !gotHs ) { obj = none; stat = endboth; } else if( !gotPt ) { obj = second; stat = endfirst; } else if( !gotHs ) { obj = first; stat = endsecond; } else //both defined { stat = endnone; if( p1 < p2 ) obj = first; else if( p1 > p2 ) obj = second; else obj = both; } } void SelectFirst_lr( const Line& L, const Region& R, object& obj, status& stat ) { L.SelectFirst(); R.SelectFirst(); HalfSegment hs1, hs2; bool gotHs1 = L.GetHs( hs1 ), gotHs2 = R.GetHs( hs2 ); if( !gotHs1 && !gotHs2 ) { obj = none; stat = endboth; } else if( !gotHs1 ) { obj = second; stat = endfirst; } else if( !gotHs2 ) { obj = first; stat = endsecond; } else //both defined { stat = endnone; if( hs1 < hs2 ) obj = first; else if( hs1 > hs2 ) obj = second; else obj = both; } } void SelectNext_lr( const Line& L, const Region& R, object& obj, status& stat ) { // 1. get the current elements HalfSegment hs1, hs2; bool gotHs1 = L.GetHs( hs1 ), gotHs2 = R.GetHs( hs2 ); //2. move the pointers if( !gotHs1 && !gotHs2 ) { //do nothing } else if( !gotHs1 ) { R.SelectNext(); gotHs2 = R.GetHs( hs2 ); } else if( !gotHs2 ) { L.SelectNext(); gotHs1 = L.GetHs( hs1 ); } else //both currently defined { if( hs1 < hs2 ) //then hs1 is the last output { L.SelectNext(); gotHs1 = L.GetHs( hs1 ); } else if( hs1 > hs2 ) { R.SelectNext(); gotHs2 = R.GetHs( hs2 ); } else { L.SelectNext(); gotHs1 = L.GetHs( hs1 ); R.SelectNext(); gotHs2 = R.GetHs( hs2 ); } } //3. generate the outputs if( !gotHs1 && gotHs2 ) { obj = none; stat = endboth; } else if( !gotHs1 ) { obj = second; stat = endfirst; } else if( !gotHs2 ) { obj = first; stat = endsecond; } else //both defined { stat = endnone; if( hs1 < hs2 ) obj = first; else if( hs1 > hs2 ) obj = second; else obj=both; } } /* 3 Auxiliary classes */ // Generic print-function object for printing with STL::for_each template struct print : public unary_function { print(ostream& out) : os(out) {} void operator() (T x) { os << "\t" << x << "\n"; } ostream& os; }; /* 4 Type Constructor ~point~ A value of type ~point~ represents a point in the Euclidean plane or is undefined. 4.1 Implementation of the class ~Point~ */ ostream& operator<<( ostream& o, const Point& p ) { ios_base::fmtflags oldOptions = o.flags(); o.setf(ios_base::fixed,ios_base::floatfield); o.precision(8); if( p.IsDefined() ) o << "(" << p.GetX() << ", " << p.GetY() << ")"; else o << Symbol::UNDEFINED(); o.flags(oldOptions); return o; } const Rectangle<2> Point::BoundingBox(const Geoid* geoid /*=0*/) const { assert( IsDefined() ); if( IsDefined() ) { if( geoid && geoid->IsDefined() ){ // spherical case double minx = MAX(-180, x - ApplyFactor(x)); double maxx = MIN(180, x+ ApplyFactor(x)); double miny = MAX(-90, y-ApplyFactor(y)); double maxy = MIN(90,y+ApplyFactor(y)); double minMax[] = {minx,maxx,miny,maxy}; return Rectangle<2>( true, minMax); } else if(!geoid){ double minMax[] = { x - ApplyFactor(x), x + ApplyFactor(x), y - ApplyFactor(y), y + ApplyFactor(y) }; return Rectangle<2>( true, minMax); } } // else: return Rectangle<2>( false ); } ostream& Point::Print( ostream &os ) const { return os << *this; } string Point::toString(const Geoid* geoid /*=0*/) const { if(!IsDefined()){ return Symbol::UNDEFINED(); } stringstream s; s.setf(ios_base::fixed,ios_base::floatfield); s.precision(8); if(geoid){ // print as geographic coords double lat = GetY(); double lon = GetX(); string ns = (lat>=0)?"N":"S"; string ew = (lon>=0)?"E":"W"; lat = fabs(lat); int degLat = (int)lat; double resLat = (lat - degLat)*60.0; int minLat = ((int) resLat); double secLat = (resLat - minLat)*60.0; lon = fabs(lon); int degLon = (int)lon; double resLon = (lon - degLon)*60.0; int minLon = ((int) resLon); double secLon = (resLon - minLon)*60.0; s << "(" << degLat << "°" << minLat << "'" << secLat << "\"" << ns << ", " << degLon << "°" << minLon << "'" << secLon << "\"" << ew << ")"; } else { // print as euclidean coords s << *this; } return s.str(); } bool Point::Inside( const Rectangle<2>& r, const Geoid* geoid /*=0*/ ) const { assert( r.IsDefined() ); if( !IsDefined() || !r.IsDefined() || (geoid && !geoid->IsDefined()) ){ return false; } if( x < r.MinD(0) || x > r.MaxD(0) ) return false; else if( y < r.MinD(1) || y > r.MaxD(1) ) return false; return true; } void Point::ReadFromString(string value) { ListExpr list; if(!nl->ReadFromString(value,list)){ if(!nl->ReadFromString("("+value+")",list)){ SetDefined(false); return; } } if(listutils::isSymbolUndefined(list)){ SetDefined(false); return; } if( nl->HasLength(list,2) && listutils::isNumeric(nl->First(list)) && listutils::isNumeric(nl->Second(list))){ Set(listutils::getNumValue(nl->First(list)), listutils::getNumValue(nl->Second(list))); return; } if( nl->HasLength(list,3) && listutils::isNumeric(nl->First(list)) && listutils::isNumeric(nl->Third(list))){ Set(listutils::getNumValue(nl->First(list)), listutils::getNumValue(nl->Third(list))); return; } SetDefined(false); } templateclass Array> void Point::Intersection(const Point& p, PointsT& result, const Geoid* geoid /*=0*/ ) const{ result.Clear(); if(!IsDefined() || !p.IsDefined() || (geoid && !geoid->IsDefined()) ){ result.SetDefined(false); return; } result.SetDefined(true); if(AlmostEqual(*this, p)){ result += *this; } } templateclass Array> void Point::Intersection(const PointsT& ps, PointsT& result, const Geoid* geoid /*=0*/ ) const{ ps.Intersection(*this, result, geoid); } templateclass Array> void Point::Intersection(const LineT& l, PointsT& result, const Geoid* geoid /*=0*/) const{ l.Intersection(*this, result, geoid); } templateclass Array> void Point::Intersection(const RegionT& r, PointsT& result, const Geoid* geoid /*=0*/) const{ r.Intersection(*this, result, geoid); } templateclass Array> void Point::Intersection(const SimpleLineT& l, PointsT& result, const Geoid* geoid /*=0*/) const{ l.Intersection(*this, result, geoid); } templateclass Array> void Point::Minus(const Point& p, PointsT& result, const Geoid* geoid /*=0*/) const{ result.Clear(); if(!IsDefined() || !p.IsDefined() || (geoid && !geoid->IsDefined())){ result.SetDefined(false); return; } result.SetDefined(true); if(!AlmostEqual(*this, p)){ result += *this; } } templateclass Array> void Point::Minus(const PointsT& ps, PointsT& result, const Geoid* geoid /*=0*/) const{ result.Clear(); if(!IsDefined() || !ps.IsDefined() || (geoid && !geoid->IsDefined())){ result.SetDefined(false); return; } result.SetDefined(true); if(!ps.Contains(*this, geoid)){ result += *this; } } templateclass Array> void Point::Minus(const LineT& l, PointsT& result, const Geoid* geoid /*=0*/) const{ result.Clear(); if(!IsDefined() || !l.IsDefined() || (geoid && !geoid->IsDefined())){ result.SetDefined(false); return; } result.SetDefined(true); if(!l.Contains(*this, geoid)){ result += *this; } } templateclass Array> void Point::Minus(const RegionT& r, PointsT& result, const Geoid* geoid /*=0*/) const{ result.Clear(); if(!IsDefined() || !r.IsDefined() || (geoid && !geoid->IsDefined())){ result.SetDefined(false); return; } result.SetDefined(true); if(!r.Contains(*this, geoid)){ result += *this; } } templateclass Array> void Point::Minus(const SimpleLineT& l, PointsT& result, const Geoid* geoid /*=0*/) const{ result.Clear(); if(!IsDefined() || !l.IsDefined() || (geoid && !geoid->IsDefined())){ result.SetDefined(false); return; } result.SetDefined(true); if(!l.Contains(*this, geoid)){ result += *this; } } templateclass Array> void Point::Union(const Point& p, PointsT& result, const Geoid* geoid /*=0*/) const{ result.Clear(); if(!IsDefined() || !p.IsDefined() || (geoid && !geoid->IsDefined())){ result.SetDefined(false); return; } result.SetDefined(true); result.StartBulkLoad(); result += *this; result += p; result.EndBulkLoad(); } templateclass Array> void Point::Union(const PointsT& ps, PointsT& result, const Geoid* geoid /*=0*/) const{ ps.Union(*this, result, geoid); } templateclass Array> void Point::Union(const LineT& l, LineT& result, const Geoid* geoid /*=0*/) const{ l.Union(*this, result, geoid); } templateclass Array> void Point::Union(const RegionT& r, RegionT& result, const Geoid* geoid /*=0*/) const{ r.Union(*this, result, geoid); } templateclass Array> void Point::Union(const SimpleLineT& l, SimpleLineT& result, const Geoid* geoid /*=0*/) const{ l.Union(*this, result, geoid); } double Point::Distance( const Point& p, const Geoid* geoid /* = 0 */ ) const { assert( IsDefined() ); assert( p.IsDefined() ); assert( !geoid || geoid->IsDefined() ); if(geoid){ bool ok = false; double bearInitial = 0, bearFinal = 0; double d = DistanceOrthodromePrecise(p,*geoid,ok,bearInitial,bearFinal); if(ok){ return d; } else{ assert(false); return -1.0; } } else { double dx = p.x - x, dy = p.y - y; return sqrt( pow( dx, 2 ) + pow( dy, 2 ) ); } } double Point::Distance( const Rectangle<2>& r, const Geoid* geoid/*=0*/ ) const { assert( IsDefined() ); assert( r.IsDefined() ); assert( !geoid || geoid->IsDefined() ); if(geoid){ cout << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented!" << endl; // TODO: implement spherical gemetry case assert(false); } double rxmin = r.MinD(0), rxmax = r.MaxD(0), rymin = r.MinD(1), rymax = r.MaxD(1); double dx = ( (x > rxmin || AlmostEqual(x,rxmin)) && (x < rxmax || AlmostEqual(x,rxmax))) ? (0.0) : (min(abs(x-rxmin),abs(x-rxmax))); double dy = ( (y > rymin || AlmostEqual(y,rymin)) && (y < rymax || AlmostEqual(y,rymax))) ? (0.0) : (min(abs(y-rymin),abs(y-rymax))); return sqrt( pow( dx, 2 ) + pow( dy, 2 ) ); } bool Point::Intersects(const Rectangle<2>& r, const Geoid* geoid/*=0*/) const{ assert(IsDefined()); assert(r.IsDefined()); assert(!geoid); // not implemented yet return (x>=r.MinD(0) ) && (x<=r.MaxD(0)) && (y>=r.MinD(1) ) && (y<=r.MaxD(1)); } // calculate the enclosed angle between (a,b) and (b,c) in degrees double Point::calcEnclosedAngle( const Point &a, const Point &b, const Point &c, const Geoid* geoid /* = 0 */){ double beta = 0.0; errno = 0; if(geoid){ // use sperical trigonometry // very simplistic, works for short distances, but // should be improved! bool ok = true; double la = b.DistanceOrthodrome(c,*geoid, ok); // la = |(b,c)| assert(ok); double lb = a.DistanceOrthodrome(c,*geoid, ok); // lb = |(a,c)| assert(ok); double lc = a.DistanceOrthodrome(b,*geoid, ok); // lc = |(a,b)| assert(ok); assert(la != 0.0); assert(lc != 0.0); double cosb = (la*la + lc*lc - lb*lb) / (2*la*lc); cosb = max(cosb, -1.0); cosb = min(cosb, 1.0); beta = radToDeg( acos(cosb) ); assert(errno == 0); } else { // use euclidean geometry double la = b.Distance(c); double lb = a.Distance(c); double lc = a.Distance(b); double cosb = (la*la + lc*lc - lb*lb) / (2*la*lc); cosb = max(cosb, -1.0); cosb = min(cosb, 1.0); beta = radToDeg( acos(cosb) ); assert(errno == 0); } return beta; } // return the direction at starting point ~this~ double Point::Direction(const Point& p , const bool returnHeading, const Geoid* geoid, const bool atEndpoint /*=false*/ ) const{ if(!IsDefined() || !p.IsDefined() || (geoid && !geoid->IsDefined()) ){ throw( new SecondoException(string("invalid Arguments") + __FILE__ + " " + stringutils::int2str(__LINE__))); return -1.0; } if(geoid && (!checkGeographicCoord() || !p.checkGeographicCoord()) ){ cerr << __PRETTY_FUNCTION__ << ": Invalid geographic coordinate." << endl; throw( new SecondoException(string("invalid coordinate") + __FILE__ + " " + stringutils::int2str(__LINE__))); return -1.0; } if(AlmostEqual(*this, p)){ throw( new SecondoException(string("equal points") + __FILE__ + " " + stringutils::int2str(__LINE__))); return -1.0; } errno = 0; double direction = 0; if(!geoid){ //euclidean geometry Coord x1 = x, y1 = y, x2 = p.x, y2 = p.y; if( AlmostEqual( x1, x2 ) ){ if( y2 > y1 ){ direction = 90.0; } else { direction = 270.0; } } else if( AlmostEqual( y1, y2 ) ){ if( x2 > x1 ){ direction = 0.0; } else { direction = 180.0; } } else { direction = radToDeg( atan2( (y2 - y1) , (x2 - x1) ) ); direction = (direction <0.0)?direction+360.0:direction; } if(errno != 0) { cerr << __PRETTY_FUNCTION__ << ": Numerical error in euclidean." << endl; throw( new SecondoException(string("numerical error") + __FILE__ + " " + stringutils::int2str(__LINE__))); return -1.0; } // normalize return value: if(returnHeading){ // return result as standard heading direction = directionToHeading(direction); // already normalized } else { // return result as standard direction: Normalize to 0<=DIR<360 direction = AlmostEqual(direction,360.0)?0.0:direction; } return direction; } else {// spherical geometry double tc1=0.0, tc2=0.0; // true course at *this, p bool valid = true; double d = DistanceOrthodromePrecise(p, *geoid, valid, tc1, tc2); if(!valid || (d<0)){ throw( new SecondoException(string("error in distance computation") + __FILE__ + " " + stringutils::int2str(__LINE__))); return -1.0; } if(returnHeading){ return (atEndpoint)?tc2:tc1; } else { return headingToDirection((atEndpoint)?tc2:tc1); } } } // Returns the vertex, i.e. the southernmost/northernmost point // of the orthodrome *this -> p Point Point::orthodromeVertex(const double& rcourseDEG) const { if( !checkGeographicCoord() || (rcourseDEG<0) || (rcourseDEG>360) ){ return Point(false, 0.0, 0.0); } double Latv = 0; double Lov = 0; // Latitude Latv of Vertex double rcourse = directionToHeading(rcourseDEG); // orthodr. course in degree double cosLatv = cos(degToRad(this->GetY()))*sin(degToRad(rcourse)); Latv = radToDeg(acos(cosLatv)); if (AlmostEqual(rcourse,0) || AlmostEqual(rcourse,180)){ Latv = 90; } if (Latv > 90){ Latv = 180 - Latv; } if (this->GetY()<0){ Latv = - Latv; } // Longitude Lov of Vertex if (!AlmostEqual(Latv,0)){ double sinDlov = cos(degToRad(rcourse))/sin(degToRad(Latv)); double Dlov = radToDeg(asin(sinDlov)); if (rcourse > 180){ Dlov = -Dlov; } Lov = Dlov + this->GetX(); } else { Lov = 0; } return Point(true, Lov, Latv); } int Point::orthodromeAtLatitude( const Point &other, const double& latitudeDEG, double& lonMinDEG, double& lonMaxDEG) const{ if( !checkGeographicCoord() || !other.checkGeographicCoord() || AlmostEqual(*this,other) || (latitudeDEG<-90) || (latitudeDEG>90) ){ throw( new SecondoException(string("invalid parameter") + __FILE__ + " " + stringutils::int2str(__LINE__))); return 0; } double lon1 = degToRad(GetX()); double lat1 = degToRad(GetY()); double lon2 = degToRad(other.GetX()); double lat2 = degToRad(other.GetY()); double lat3 = degToRad(latitudeDEG); double l12 = lon1-lon2; double A = sin(lat1)*cos(lat2)*cos(lat3)*sin(l12); double B = sin(lat1)*cos(lat2)*cos(lat3)*cos(l12) - cos(lat1)*sin(lat2)*cos(lat3); double C = cos(lat1)*cos(lat2)*sin(lat3)*sin(l12); double lon = atan2(B,A); if(fabs(C) <= sqrt(A*A + B*B)){ double dlon = acos(C/sqrt(A*A+B*B)); double lon3_1=fmod2(lon1+dlon+lon+M_PI, 2*M_PI) - M_PI; double lon3_2=fmod2(lon1-dlon+lon+M_PI, 2*M_PI) - M_PI; lonMinDEG = radToDeg(MIN(lon3_1,lon3_2)); lonMaxDEG = radToDeg(MAX(lon3_1,lon3_2)); if(AlmostEqual(lon3_1,lon3_2)){ return 1; } return 2; } throw( new SecondoException(string("numerical error") + __FILE__ + " " + stringutils::int2str(__LINE__))); return 0; } bool isBetween(const double value, const double bound1, const double bound2){ return ( (MIN(bound1,bound2)<=value) && (value<=MAX(bound1,bound2)) ); } bool Point::orthodromeExtremeLatitudes(const Point &other, const Geoid &geoid, double &minLat, double &maxLat) const { if(!checkGeographicCoord() || !other.checkGeographicCoord() || !geoid.IsDefined() ){ minLat = 1000.0; maxLat = -1000.0; // cerr << __PRETTY_FUNCTION__ << ": Invalid parameter." << endl; throw( new SecondoException(string("invalid parameter") + __FILE__ + " " + stringutils::int2str(__LINE__))); return false; } if(AlmostEqual(*this,other)){ // cerr << __PRETTY_FUNCTION__ << ": *this == other." << endl; minLat=MIN(this->GetY(),other.GetY()); maxLat=MAX(this->GetY(),other.GetY()); // cerr << __PRETTY_FUNCTION__ << ": minLat=" << minLat << endl; // cerr << __PRETTY_FUNCTION__ << ": maxLat=" << maxLat << endl; return true; } double trueCourse = Direction(other,false,&geoid); if(AlmostEqual(trueCourse,0) || AlmostEqual(trueCourse,180)){ // cerr << __PRETTY_FUNCTION__ << ": Polar course." << endl; minLat=MIN(this->GetY(),other.GetY()); maxLat=MAX(this->GetY(),other.GetY()); // cerr << __PRETTY_FUNCTION__ << ": minLat=" << minLat << endl; // cerr << __PRETTY_FUNCTION__ << ": maxLat=" << maxLat << endl; return true; } Point vertex = orthodromeVertex(trueCourse); if(!vertex.IsDefined()){ // cerr << __PRETTY_FUNCTION__ << ": Cannot compute vertex." << endl; minLat=MIN(this->GetY(),other.GetY()); maxLat=MAX(this->GetY(),other.GetY()); // cerr << __PRETTY_FUNCTION__ << ": minLat=" << minLat << endl; // cerr << __PRETTY_FUNCTION__ << ": maxLat=" << maxLat << endl; throw( new SecondoException(string("invalid parameter") + __FILE__ + " " + stringutils::int2str(__LINE__))); return false; } // cerr << __PRETTY_FUNCTION__ << ": vertex=" << vertex << endl; double vertexReachedLonMin = 6666.6666; double vertexReachedLonMax = -6666.6666; int notransgressions = orthodromeAtLatitude(other,vertex.GetY(), vertexReachedLonMin, vertexReachedLonMax); if(notransgressions == 0){ // vertex not passed, should not happen cerr << __PRETTY_FUNCTION__ << ": Vertex not passed, should not happen." << endl; assert(false); return false; } else if(notransgressions == 1) { // vertex passed once. vertexReachedLonMax = vertexReachedLonMin; // } // else: vertex passed twice if( ( isBetween(directionToHeading(trueCourse),0.0,180.0) && ( (vertexReachedLonMinMAX(GetX(),other.GetX())) ) ) || ( isBetween(directionToHeading(trueCourse),180.0,360.0) && !( (vertexReachedLonMinMAX(GetX(),other.GetX())) ) ) ){ // vertex is not reached by this orthodrome // cerr << __PRETTY_FUNCTION__ << ": Vertex not passed." << endl; minLat = MIN(GetY(),other.GetY()); maxLat = MAX(GetY(),other.GetY()); } else if( (GetY()>=0) && (other.GetY()>=0) ){// both in northern hemissphere // vertex may be north of this and p // cerr << __PRETTY_FUNCTION__ << ": Both northern hemissphere." << endl; minLat = MIN(GetY(),other.GetY()); maxLat = MAX(vertex.GetY(),MAX(GetY(),other.GetY())); } else if( (GetY()<0) && (other.GetY()<0) ){// both in southern hemissphere // vertex may be south of this and p // cerr << __PRETTY_FUNCTION__ << ": Both southern hemissphere." << endl; minLat = MIN(vertex.GetY(),MIN(GetY(),other.GetY())); maxLat = MAX(GetY(),other.GetY()); } else { // both in different hemisspheres // vertex is either this or p // cerr << __PRETTY_FUNCTION__ << ": Both in different hemisspheres."<checkGeographicCoord() && p.checkGeographicCoord() && g.IsDefined(); if(!valid){ cout << __PRETTY_FUNCTION__ << ": Invalid coordinates." << endl; return-666.666; } initialBearingDEG = -666.666; finalBearingDEG = -666.666; if(AlmostEqual(*this,p)){ return 0; } errno = 0; double lat1 = GetY(); double lon1 = GetX(); double lat2 = p.GetY(); double lon2 = p.GetX(); double a = g.getR();// = 6378137 double f = g.getF();// = 1/298.257223563 double b = (1-f)*a; // = 6356752.314245 double L = degToRad(lon2-lon1); double U1 = atan( (1-f)*tan(degToRad(lat1)) ); double U2 = atan( (1-f)*tan(degToRad(lat2)) ); double sinU1 = sin(U1); double cosU1 = cos(U1); double sinU2 = sin(U2); double cosU2 = cos(U2); if(errno!=0){ cout << __PRETTY_FUNCTION__ << ": Error at (1)." << endl; valid = false; return -666.666; } double lambda = L; double lambdaP; int iterLimit = 100; int iter = 1; double sinLambda, cosLambda, sinSigma, cosSigma, sigma, sinAlpha, cosSqAlpha, cos2SigmaM, C; do{ sinLambda = sin(lambda); cosLambda = cos(lambda); sinSigma = sqrt( (cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda)); if (sinSigma==0){ // co-incident points valid = false; cout << __PRETTY_FUNCTION__ << ": Error: Co-incident points: *this=" << *this << ", p=" << p << endl; valid = false; return -666.666; } if(errno!=0){ cout << __PRETTY_FUNCTION__ << ": Error (1) in iteration " < epsilon) && (iter++iterLimit) { // formula failed to converge cout << __PRETTY_FUNCTION__ << ": Formula failed to converge after " << iter << " iterations. Delta=" << fabs(lambda-lambdaP) << " > epsilon=" << epsilon << "." << endl; valid = false; return -666.666; } double uSq = cosSqAlpha * (a*a - b*b) / (b*b); double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); double deltaSigma = B*sinSigma*(cos2SigmaM +(B/4)*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)- (B/6)*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))); double s = b*A*(sigma-deltaSigma); if(errno!=0){ cout << __PRETTY_FUNCTION__ << ": Error (3) " << iter << endl; valid = false; return -666.666; } // calculate initial/final bearings in addition to distance double fwdAz = atan2(cosU2*sinLambda, cosU1*sinU2-sinU1*cosU2*cosLambda); double revAz = atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda); initialBearingDEG = radToDeg(fwdAz); finalBearingDEG = radToDeg(revAz); // normalize heading results to 0 360.0){ initialBearingDEG = fmod(initialBearingDEG,360.0); } if( AlmostEqual(initialBearingDEG,0.0) ){ initialBearingDEG = 360.0; // Map NORTH to heading 360, not to 0 } while( finalBearingDEG<0.0 ){ finalBearingDEG += 360.0; } if(finalBearingDEG > 360.0){ finalBearingDEG = fmod(finalBearingDEG,360.0); } if( AlmostEqual(finalBearingDEG,0.0) ){ finalBearingDEG = 360.0; // Map NORTH to heading 360, not to 0 } if(errno!=0){ cout << __PRETTY_FUNCTION__ << ": Error (4) " << iter << endl; valid = false; return -666.666; } // cout << __PRETTY_FUNCTION__ << ": After " << iter // << " Iterations: Distance[m]=" << s // << ", initialBearingDEG=" << initialBearingDEG // << ", finalBearingDEG=" << finalBearingDEG << endl; valid = true; return s; } Point Point::MidpointTo(const Point& p, const Geoid* geoid /* = 0 */ ) const { if ( !this->IsDefined() || !p.IsDefined() ) { return Point( false ); } if( AlmostEqual(*this,p) ){ // no calculation required. return *this; } if(geoid){ // spherical geometry // TODO: Use ellipsoid ~geoid~ instead of ideal sphere! bool ok = this->checkGeographicCoord() && p.checkGeographicCoord(); double lon1 = degToRad( GetX() ); // X <-> geogr. Laenge (Longitude) double lat1 = degToRad( GetY() ); // Y <-> geogr. Breite (Latitude) double lat2 = degToRad( p.GetY() ); double dLon = degToRad( p.GetX() - GetX() ); errno = 0; double Bx = cos(lat2) * cos(dLon); double By = cos(lat2) * sin(dLon); double lat3 = atan2( sin(lat1)+sin(lat2), sqrt( (cos(lat1)+Bx)*(cos(lat1)+Bx) + By*By) ); double lon3 = lon1 + atan2(By, cos(lat1) + Bx); Point p( true, radToDeg(lon3), radToDeg(lat3) ); p.SetDefined( ok && (errno==0) && p.checkGeographicCoord() ); return p; } else { // euclidean geometry return Point( true, ( GetX() + p.GetX() )/2.0, ( GetY() + p.GetY() )/2.0 ); } } Point Point::MidpointTo(const Point& p, const Coord& f, const Geoid* geoid /*=0*/) const{ if( !IsDefined() || !p.IsDefined() || (geoid && !geoid->IsDefined()) || (f<0.0) || (f>1.0) ) { return Point(false); } if(AlmostEqual(f,0.0)) { return *this; } if(AlmostEqual(f,1.0)) { return p; } if(geoid){ // spherical TODO: use ellipsoid instead of ideal sphere double lon1 = degToRad(GetX()); // X <-> geogr. Laenge (Longitude) double lat1 = degToRad(GetY()); // Y <-> geogr. Breite (Latitude) double lon2 = degToRad(p.GetX()); double lat2 = degToRad(p.GetY()); if( AlmostEqual(lat1,lat2) && AlmostEqual(fabs(lon1-lon2),M_PI) ){ return Point(false); } errno = 0; double A = sin(1 - f) / sin(1); double B = sin(f) / sin(1); double C = A * cos(lat1) * cos(lon1) + B * cos(lat2)*cos(lon2); double D = A * cos(lat1) * sin(lon1) + B * cos(lat2) * sin(lon2); double lat = atan2(A * sin(lat1) + B * sin(lat2),sqrt(C*C + D*D)); double lon = atan2(D, C); Point p( true, radToDeg(lon), radToDeg(lat) ); p.SetDefined( (errno==0) && p.checkGeographicCoord() ); return p; } else { // euclidean double x1 = GetX(); double y1 = GetY(); double x2 = p.GetX(); double y2 = p.GetY(); double dx = x2-x1; double dy = y2-y1; double x = x1 + (f*dx); double y = y1 + (f*dy); return Point( true, x, y ); } } /* 4.2 List Representation The list representation of a point is ---- (x y) ---- 4.3 ~Out~-function */ ListExpr OutPoint( ListExpr typeInfo, Word value ) { Point* point = (Point*)(value.addr); if( point->IsDefined() ) return nl->TwoElemList( nl->RealAtom( point->GetX() ), nl->RealAtom( point->GetY() ) ); else return nl->SymbolAtom( Symbol::UNDEFINED() ); } /* 4.4 ~In~-function */ Word InPoint( const ListExpr typeInfo, const ListExpr instance, const int errorPos, ListExpr& errorInfo, bool& correct ) { correct = true; if( nl->ListLength( instance ) == 2 ) { ListExpr first = nl->First(instance); ListExpr second = nl->Second(instance); correct = listutils::isNumeric(first) && listutils::isNumeric(second); if(!correct){ return SetWord( Address(0) ); } else { return SetWord(new Point(true, listutils::getNumValue(first), listutils::getNumValue(second))); } } else if( listutils::isSymbolUndefined( instance ) ){ return SetWord(new Point(false)); } correct = false; return SetWord( Address(0) ); } /* 4.5 ~Create~-function */ Word CreatePoint( const ListExpr typeInfo ) { return SetWord( new Point( false ) ); } /* 4.6 ~Delete~-function */ void DeletePoint( const ListExpr typeInfo, Word& w ) { ((Point *)w.addr)->DeleteIfAllowed(); w.addr = 0; } /* 4.7 ~Close~-function */ void ClosePoint( const ListExpr typeInfo, Word& w ) { ((Point *)w.addr)->DeleteIfAllowed(); w.addr = 0; } /* 4.8 ~Clone~-function */ Word ClonePoint( const ListExpr typeInfo, const Word& w ) { return SetWord( new Point( *((Point *)w.addr) ) ); } /* 4.8 ~SizeOf~-function */ int SizeOfPoint() { return sizeof(Point); } /* 4.9 Function describing the signature of the type constructor */ ListExpr PointProperty() { return nl->TwoElemList( nl->FourElemList( nl->StringAtom("Signature"), nl->StringAtom("Example Type List"), nl->StringAtom("List Rep"), nl->StringAtom("Example List")), nl->FourElemList( nl->StringAtom("-> DATA"), nl->StringAtom(Point::BasicType()), nl->StringAtom("(x y)"), nl->StringAtom("(10 5)"))); } /* 4.10 Kind Checking Function This function checks whether the type constructor is applied correctly. Since type constructor ~point~ does not have arguments, this is trivial. */ bool CheckPoint( ListExpr type, ListExpr& errorInfo ) { return (listutils::isSymbol( type, Point::BasicType() )); } /* 4.12 Creation of the type constructor instance */ TypeConstructor point( Point::BasicType(), //name PointProperty, //property function describing signature OutPoint, InPoint, //Out and In functions 0, 0, //SaveToList and RestoreFromList functions CreatePoint, DeletePoint, //object creation and deletion OpenAttribute, SaveAttribute, // object open and save ClosePoint, ClonePoint, //object close, and clone Point::Cast, //cast function SizeOfPoint, //sizeof function CheckPoint); //kind checking function ostream& operator<<( ostream& o, const Points& ps ) { o << "<"; if( !ps.IsDefined() ) { o << " undef "; } else { for( int i = 0; i < ps.Size(); i++ ) { Point p(true); ps.Get( i, p ); o << " " << p; } } o << ">"; return o; } // use this when adding and sorting the DBArray int PointCompare( const void *a, const void *b ) { const Point *pa = (const Point*)a, *pb = (const Point*)b; assert(pa->IsDefined()); assert(pb->IsDefined()); if( *pa == *pb ) return 0; if( *pa < *pb ) return -1; return 1; } // use this when testing for containment or removing duplicates // in the DBArray int PointCompareAlmost( const void *a, const void *b ) { const Point *pa = (const Point*)a, *pb = (const Point*)b; assert(pa->IsDefined()); assert(pb->IsDefined()); if( AlmostEqual( *pa, *pb ) ) return 0; if( *pa < *pb ) return -1; return 1; } // Old implementation, should be replaced by the following one // to avoid problems when sorting and removing duplicates int PointHalfSegmentCompare( const void *a, const void *b ) { const Point *pa = (const Point *)a; const HalfSegment *hsb = (const HalfSegment *)b; assert(pa->IsDefined()); if( AlmostEqual( *pa, hsb->GetDomPoint() ) ) return 0; if( *pa < hsb->GetDomPoint() ) return -1; return 1; } // // use this when adding and sorting the DBArray // int PointHalfSegmentCompare( const void *a, const void *b ) // { // const Point *pa = (const Point *)a; // const HalfSegment *hsb = (const HalfSegment *)b; // // if( *pa == hsb->GetDomPoint() ) // return 0; // // if( *pa < hsb->GetDomPoint() ) // return -1; // // return 1; // } // use this when testing for containment or removing duplicates // in the DBArray int PointHalfSegmentCompareAlmost( const void *a, const void *b ) { const Point *pa = (const Point *)a; const HalfSegment *hsb = (const HalfSegment *)b; assert(pa->IsDefined()); if( AlmostEqual( *pa, hsb->GetDomPoint() ) ) return 0; if( *pa < hsb->GetDomPoint() ) return -1; return 1; } /* 5.2 List Representation The list representation of a point is ---- (x y) ---- 5.3 ~Out~-function */ ListExpr OutPoints( ListExpr typeInfo, Word value ) { Points* points = (Points*)(value.addr); if(!points->IsDefined()){ return nl->SymbolAtom(Symbol::UNDEFINED()); } if( points->IsEmpty() ) return nl->TheEmptyList(); Point p(true); assert( points->Get( 0, p ) ); ListExpr result = nl->OneElemList( OutPoint( nl->TheEmptyList(), SetWord( (void*)&p ) ) ); ListExpr last = result; for( int i = 1; i < points->Size(); i++ ) { assert( points->Get( i, p ) ); last = nl->Append( last, OutPoint( nl->TheEmptyList(), SetWord( (void*)&p ) ) ); } return result; } /* 5.4 ~In~-function */ Word InPoints( const ListExpr typeInfo, const ListExpr instance, const int errorPos, ListExpr& errorInfo, bool& correct ) { if(listutils::isSymbolUndefined(instance)) { Points* points = new Points(0); points->Clear(); points->SetDefined(false); correct=true; return SetWord( Address(points) ); } Points* points = new Points( max(0,nl->ListLength( instance) ) ); points->SetDefined( true ); if(nl->AtomType(instance)!=NoAtom) { points->DeleteIfAllowed(); correct = false; cout << __PRETTY_FUNCTION__ << ": Unexpected Atom!" << endl; return SetWord( Address(points) ); } ListExpr rest = instance; points->StartBulkLoad(); while( !nl->IsEmpty( rest ) ) { ListExpr first = nl->First( rest ); rest = nl->Rest( rest ); Point *p = (Point*)InPoint( nl->TheEmptyList(), first, 0, errorInfo, correct ).addr; if( correct && p->IsDefined() ) { (*points) += (*p); delete p; } else { if(p) { delete p; } cout << __PRETTY_FUNCTION__ << ": Incorrect or undefined point!" << endl; points->DeleteIfAllowed(); correct = false; return SetWord( Address(0) ); } } points->EndBulkLoad(); if( points->IsValid() ) { correct = true; return SetWord( points ); } points->DeleteIfAllowed(); correct = false; cout << __PRETTY_FUNCTION__ << ": Invalid points value!" << endl; return SetWord( Address(0) ); } /* 5.5 ~Create~-function */ Word CreatePoints( const ListExpr typeInfo ) { return SetWord( new Points( 0 ) ); } /* 5.6 ~Delete~-function */ void DeletePoints( const ListExpr typeInfo, Word& w ) { Points *ps = (Points *)w.addr; ps->Destroy(); ps->DeleteIfAllowed(false); w.addr = 0; } /* 5.7 ~Close~-function */ void ClosePoints( const ListExpr typeInfo, Word& w ) { ((Points *)w.addr)->DeleteIfAllowed(); w.addr = 0; } /* 5.8 ~Clone~-function */ Word ClonePoints( const ListExpr typeInfo, const Word& w ) { return SetWord( new Points( *((Points *)w.addr) ) ); } /* 7.8 ~Open~-function */ bool OpenPoints( SmiRecord& valueRecord, size_t& offset, const ListExpr typeInfo, Word& value ) { Points *ps = (Points*)Attribute::Open( valueRecord, offset, typeInfo ); value = SetWord( ps ); return true; } /* 7.8 ~Save~-function */ bool SavePoints( SmiRecord& valueRecord, size_t& offset, const ListExpr typeInfo, Word& value ) { Points *ps = (Points*)value.addr; Attribute::Save( valueRecord, offset, typeInfo, ps ); return true; } /* 5.8 ~SizeOf~-function */ int SizeOfPoints() { return sizeof(Points); } /* 5.11 Function describing the signature of the type constructor */ ListExpr PointsProperty() { return (nl->TwoElemList( nl->FourElemList(nl->StringAtom("Signature"), nl->StringAtom("Example Type List"), nl->StringAtom("List Rep"), nl->StringAtom("Example List")), nl->FourElemList(nl->StringAtom("-> DATA"), nl->StringAtom(Points::BasicType()), nl->StringAtom("(*) where point is ()"), nl->StringAtom("( (10 1)(4 5) )")))); } /* 5.12 Kind checking function This function checks whether the type constructor is applied correctly. Since type constructor ~point~ does not have arguments, this is trivial. */ bool CheckPoints( ListExpr type, ListExpr& errorInfo ) { return (nl->IsEqual( type, Points::BasicType() )); } /* 5.13 ~Cast~-function */ void* CastPoints(void* addr) { return (new (addr) Points()); } /* 5.14 Creation of the type constructor instance */ TypeConstructor points( Points::BasicType(), //name PointsProperty, //property function describing signature OutPoints, InPoints, //Out and In functions 0, 0, //SaveTo and RestoreFrom List functions CreatePoints, DeletePoints, //object creation and deletion OpenPoints, SavePoints, // object open and save ClosePoints, ClonePoints, //object close and clone CastPoints, //cast function SizeOfPoints, //sizeof function CheckPoints ); //kind checking function /* 6 Type Constructor ~halfsegment~ A ~halfsegment~ value is a pair of points, with a boolean flag indicating the dominating point . 6.1 Implementation of the class ~halfsegment~ */ void HalfSegment::Translate( const Coord& x, const Coord& y ) { lp.Translate( x, y ); rp.Translate( x, y ); } void HalfSegment::Set( bool ldp, const Point& lp, const Point& rp ) { assert( lp.IsDefined() ); assert( rp.IsDefined() ); assert( !AlmostEqual( lp, rp ) ); this->ldp = ldp; if( lp < rp ) { this->lp = lp; this->rp = rp; } else // rp > lp { this->lp = rp; this->rp = lp; } } int HalfSegment::Compare( const HalfSegment& hs ) const { const Point& dp = GetDomPoint(), sp = GetSecPoint(), DP = hs.GetDomPoint(), SP = hs.GetSecPoint(); if( dp < DP ) return -1; else if( dp > DP ) return 1; if( ldp != hs.ldp ) { if( ldp == false ) return -1; return 1; } else { bool v1 = IsVertical(); bool v2 = hs.IsVertical(); if( v1 && v2 ) // both are vertical { if( ( (CompareDouble(sp.GetY(),dp.GetY())>0) && ( CompareDouble(SP.GetY(),DP.GetY())>0) ) || ( (CompareDouble(sp.GetY(),dp.GetY())<0) && (CompareDouble(SP.GetY(),DP.GetY())<0) ) ) { if( sp < SP ) return -1; if( sp > SP ) return 1; return 0; } else if( CompareDouble(sp.GetY(),dp.GetY())>0) { if( ldp == true ) return 1; return -1; } else { if( ldp == true ) return -1; return 1; } } else if( AlmostEqual(dp.GetX(),sp.GetX()) ) { if( CompareDouble(sp.GetY(), dp.GetY())>0 ) { if( ldp == true ) return 1; return -1; } else if( CompareDouble(sp.GetY(),dp.GetY())<0 ) { if( ldp == true ) return -1; return 1; } } else if( AlmostEqual(DP.GetX(), SP.GetX()) ) { if( CompareDouble(SP.GetY() , DP.GetY())>0 ) { if( ldp == true ) return -1; return 1; } else if( CompareDouble(SP.GetY() , DP.GetY())<0 ) { if( ldp == true ) return 1; return -1; } } else { Coord xd = dp.GetX(), yd = dp.GetY(), xs = sp.GetX(), ys = sp.GetY(), Xd = DP.GetX(), Yd = DP.GetY(), Xs = SP.GetX(), Ys = SP.GetY(); double k = (yd - ys) / (xd - xs), K= (Yd -Ys) / (Xd - Xs); if( CompareDouble(k , K) <0 ) return -1; if( CompareDouble( k, K) > 0) return 1; if( sp < SP ) return -1; if( sp > SP ) return 1; return 0; } } assert( true ); // This code should never be reached return 0; } HalfSegment& HalfSegment::operator=( const HalfSegment& hs ) { ldp = hs.ldp; lp = hs.lp; rp = hs.rp; attr = hs.attr; return *this; } bool HalfSegment::operator==( const HalfSegment& hs ) const { return Compare(hs) == 0; } bool HalfSegment::operator!=( const HalfSegment& hs ) const { return !( *this == hs ); } bool HalfSegment::operator<( const HalfSegment& hs ) const { return Compare(hs) == -1; } bool HalfSegment::operator>( const HalfSegment& hs ) const { return Compare(hs) == 1; } int HalfSegment::LogicCompare( const HalfSegment& hs ) const { if( attr.faceno < hs.attr.faceno ) return -1; if( attr.faceno > hs.attr.faceno ) return 1; if( attr.cycleno < hs.attr.cycleno ) return -1; if( attr.cycleno > hs.attr.cycleno ) return 1; if( attr.edgeno < hs.attr.edgeno ) return -1; if( attr.edgeno > hs.attr.edgeno ) return 1; return 0; } int HalfSegmentCompare(const void *a, const void *b) { const HalfSegment *hsa = (const HalfSegment *)a, *hsb = (const HalfSegment *)b; return hsa->Compare( *hsb ); } int HalfSegmentLogicCompare(const void *a, const void *b) { const HalfSegment *hsa = (const HalfSegment *)a, *hsb = (const HalfSegment *)b; return hsa->LogicCompare( *hsb ); } int LRSCompare( const void *a, const void *b ) { const LRS *lrsa = (const LRS *)a, *lrsb = (const LRS *)b; if( lrsa->lrsPos < lrsb->lrsPos ) return -1; if( lrsa->lrsPos > lrsb->lrsPos ) return 1; return 0; } bool HalfSegment::insideLeft() const{ if(ldp){ return attr.insideAbove; } else { return !attr.insideAbove; } } ostream& operator<<(ostream &os, const HalfSegment& hs) { return os << "(" <<"F("<< hs.attr.faceno <<") C("<< hs.attr.cycleno <<") E(" << hs.attr.edgeno<<") DP(" << (hs.IsLeftDomPoint()? "L":"R") <<") IA("<< (hs.attr.insideAbove? "A":"U") <<") Co("< Yl && yl < Yr ) || ( yr > Yl && yr < Yr ) || ( Yl > yl && Yl < yr ) || ( Yr > yl && Yr < yr ) ) ) return true; return false; } if( !AlmostEqual( xl, xr ) ) // this segment is not vertical { k = (yr - yl) / (xr - xl); a = yl - k * xl; } if( !AlmostEqual( Xl, Xr ) ) // hs is not vertical { K = (Yr - Yl) / (Xr - Xl); A = Yl - K * Xl; } if( AlmostEqual( Xl, Xr ) ) //only hs is vertical { Coord y0 = k * Xl + a; if( ( Xl > xl || AlmostEqual( Xl, xl ) ) && ( Xl < xr || AlmostEqual( Xl, xr ) ) ) { if( ( ( y0 > Yl || AlmostEqual( y0, Yl ) ) && ( y0 < Yr || AlmostEqual( y0, Yr ) ) ) || ( ( y0 > Yr || AlmostEqual( y0, Yr ) ) && ( y0 < Yl || AlmostEqual( y0, Yl ) ) ) ) // (Xl, y0) is the intersection point return true; } return false; } if( AlmostEqual( xl, xr ) ) // only this segment is vertical { Coord Y0 = K * xl + A; if( ( xl > Xl || AlmostEqual( xl, Xl ) ) && ( xl < Xr || AlmostEqual( xl, Xr ) ) ) { if( ( ( Y0 > yl || AlmostEqual( Y0, yl ) ) && ( Y0 < yr || AlmostEqual( Y0, yr ) ) ) || ( ( Y0 > yr || AlmostEqual( Y0, yr ) ) && ( Y0 < yl || AlmostEqual( Y0, yl ) ) ) ) // (xl, Y0) is the intersection point return true; } return false; } // both segments are non-vertical if( AlmostEqual( k, K ) ) // both segments have the same inclination { /* if( ( AlmostEqual( A, a ) && (( ( xl > Xl || AlmostEqual( xl, Xl ) ) && ( xl < Xr || AlmostEqual( xl, Xr ) ) ))) || ( ( Xl > xl || AlmostEqual( xl, Xl ) ) && ( Xl < xr || AlmostEqual( xr, Xl ) ) ) ) // the segments are in the same straight line return true;*/ if( AlmostEqual( A, a ) && (( ( xl > Xl || AlmostEqual( xl, Xl ) ) && ( xl < Xr || AlmostEqual( xl, Xr ) ) ) || ( ( Xl > xl || AlmostEqual( xl, Xl ) ) && ( Xl < xr || AlmostEqual( xr, Xl ) ) ) )) // the segments are in the same straight line return true; } else { Coord x0 = (A - a) / (k - K); // y0 = x0 * k + a; if( ( x0 > xl || AlmostEqual( x0, xl ) ) && ( x0 < xr || AlmostEqual( x0, xr ) ) && ( x0 > Xl || AlmostEqual( x0, Xl ) ) && ( x0 < Xr || AlmostEqual( x0, Xr ) ) ) // the segments intersect at (x0, y0) return true; } return false; } bool HalfSegment::InnerIntersects( const HalfSegment& hs, const Geoid* geoid/*=0*/ ) const { if(geoid){ cerr << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." << endl; assert( false ); // TODO: implement spherical geometry case } double k = 0.0; double a = 0.0; double K = 0.0; double A = 0.0; Coord x0; //, y0; (x0, y0) is the intersection Coord xl = lp.GetX(), yl = lp.GetY(), xr = rp.GetX(), yr = rp.GetY(); if( xl != xr ) { k = (yr - yl) / (xr - xl); a = yl - k * xl; } Coord Xl = hs.GetLeftPoint().GetX(), Yl = hs.GetLeftPoint().GetY(), Xr = hs.GetRightPoint().GetX(), Yr = hs.GetRightPoint().GetY(); if( Xl != Xr ) { K = (Yr - Yl) / (Xr - Xl); A = Yl - K * Xl; } if( xl == xr && Xl == Xr ) //both l and L are vertical lines { if( xl != Xl ) return false; Coord ylow, yup, Ylow, Yup; if( yl < yr ) { ylow = yl; yup = yr; } else { ylow = yr; yup = yl; } if( Yl < Yr) { Ylow = Yl; Yup = Yr; } else { Ylow = Yr; Yup = Yl; } if( ylow >= Yup || yup <= Ylow ) return false; return true; } if( Xl == Xr ) //only L is vertical { double y0 = k * Xl + a; Coord yy = y0; //(Xl, y0) is the intersection of l and L if( Xl >= xl && Xl <= xr ) { if( ( yy > Yl && yy < Yr ) || ( yy > Yr && yy < Yl ) ) return true; return false; } else return false; } if( xl == xr ) //only l is vertical { double Y0 = K * xl + A; Coord YY = Y0; //(xl, Y0) is the intersection of l and L if( xl > Xl && xl < Xr ) { if( ( YY >= yl && YY <= yr ) || ( YY >= yr && YY <= yl ) ) return true; return false; } else return false; } //otherwise: both segments are non-vertical if( k == K ) { if( A != a ) //Parallel lines return false; //they are in the same straight line if( xr <= Xl || xl >= Xr ) return false; return true; } else { x0 = (A - a) / (k - K); // y0=x0*k+a; Coord xx = x0; if( xx >= xl && xx <= xr && xx > Xl && xx < Xr ) return true; return false; } } // corrected modulo function ( '%' has problems with negative values) double fmod2(const double &y, const double &x) { return y - x*floor(y/x); } bool HalfSegment::Intersection( const HalfSegment& hs, Point& resp, const Geoid* geoid /* = 0 */ ) const { if(geoid){ // spherical geometry not yet implemented! assert(false); if(!geoid->IsDefined()){ resp.SetDefined(false); } // see http://williams.best.vwh.net/avform.htm#Intersection Point res(false);//Destination point (UNDEFINED if no unique intersection) Point p1 = this->GetLeftPoint(); //1st segment's starting point Point p1end = this->GetRightPoint(); //1st segment's ending point Point p2 = hs.GetLeftPoint(); //2nd segment' s starting point Point p2end = hs.GetRightPoint(); //2nd segment's ending point if(!p1.checkGeographicCoord() || !p2.checkGeographicCoord() || !p1end.checkGeographicCoord() || !p2end.checkGeographicCoord() ) { return false; } double lon1 = degToRad(p1.GetX()); double lat1 = degToRad(p1.GetY()); double lon2 = degToRad(p2.GetX()); double lat2 = degToRad(p2.GetY()); double brng1 = p1.Direction(p1end,false,geoid); //Initial bearing from p1 double brng2 = p2.Direction(p2end,false,geoid); //Initial bearing from p2 if( (brng1<0) || (brng2<0) ){ return false; } double brng13 = degToRad(brng1); double brng23 = degToRad(brng2); double dLat = lat2-lat1; double dLon = lon2-lon1; errno = 0; double dist12 = 2*asin( sqrt( sin(dLat/2)*sin(dLat/2) + cos(lat1)*cos(lat2)*sin(dLon/2)*sin(dLon/2) ) ); if( (errno !=0) || (dist12 <= 0) ){ return false; } // initial/final bearings between points errno = 0; double brngA = acos( ( sin(lat2) - sin(lat1)*cos(dist12) ) / ( sin(dist12)*cos(lat1) ) ); if(errno != 0){ brngA = 0; // protect against rounding errno = 0; } double brngB = acos( ( sin(lat1) - sin(lat2)*cos(dist12) ) / ( sin(dist12)*cos(lat2) ) ); double brng12 = 0; double brng21 = 0; if (sin(lon2-lon1) > 0) { brng12 = brngA; brng21 = 2*M_PI - brngB; } else { brng12 = 2*M_PI - brngA; brng21 = brngB; } double alpha1 = fmod2((brng13 - brng12 +M_PI),(2*M_PI)) - M_PI;//angle 2-1-3 double alpha2 = fmod2((brng21 - brng23 +M_PI),(2*M_PI)) - M_PI;//angle 1-2-3 if( (sin(alpha1)==0) && (sin(alpha2)==0) ){// infinite intersections return false; } if( (errno != 0) || (sin(alpha1)*sin(alpha2) < 0) ){//ambiguous intersection return false; } //Ed Williams takes abs of alpha1/alpha2, but seems to break calculation? //alpha1 = fabs(alpha1); //alpha2 = fabs(alpha2); double alpha3 = acos( -cos(alpha1)*cos(alpha2) + sin(alpha1)*sin(alpha2)*cos(dist12) ); double dist13 = atan2( sin(dist12)*sin(alpha1)*sin(alpha2), cos(alpha2)+cos(alpha1)*cos(alpha3) ); double lat3 = asin( sin(lat1)*cos(dist13) + cos(lat1)*sin(dist13)*cos(brng13) ); double dLon13 = atan2( sin(brng13)*sin(dist13)*cos(lat1), cos(dist13)-sin(lat1)*sin(lat3) ); double lon3 = lon1+dLon13; lon3 = fmod2((lon3+M_PI),(2*M_PI)) - M_PI; // normalise to -180..180º res.SetDefined( (errno == 0) ); res.Set( radToDeg(lon3), radToDeg(lat3) ); return true; } else { // euclidean geometry double k = 0.0; double a = 0.0; double K = 0.0; double A = 0.0; if(!BoundingBox().Intersects(hs.BoundingBox())){ resp.SetDefined(false); return false; } resp.SetDefined( true ); Coord xl = lp.GetX(), yl = lp.GetY(), xr = rp.GetX(), yr = rp.GetY(), Xl = hs.GetLeftPoint().GetX(), Yl = hs.GetLeftPoint().GetY(), Xr = hs.GetRightPoint().GetX(), Yr = hs.GetRightPoint().GetY(); // Check for same endpoints if (AlmostEqual(lp,hs.GetLeftPoint())){ if (!hs.Contains(rp) && !this->Contains(hs.GetRightPoint())){ resp = lp; return true; } else { return false; //overlapping halfsegments } } else if (AlmostEqual(rp,hs.GetLeftPoint())){ if (!hs.Contains(lp) && !this->Contains(hs.GetRightPoint())){ resp = rp; return true; } else { return false; //overlapping halfsegments } } else if (AlmostEqual(lp,hs.GetRightPoint())){ if (!hs.Contains(rp) && !this->Contains(hs.GetLeftPoint())){ resp = lp; return true; } else { return false; //overlapping halfsegments } } else if (AlmostEqual(rp,hs.GetRightPoint())){ if (!hs.Contains(lp) && !this->Contains(hs.GetLeftPoint())){ resp = rp; return true; } else { return false; //overlapping halfsegments } } if( AlmostEqual( xl, xr ) && AlmostEqual( Xl, Xr ) ){ // both segments are vertical if( AlmostEqual( yr, Yl ) ){ resp.Set( xl, yr ); return true; } if( AlmostEqual( yl, Yr ) ) { resp.Set( xl, yl ); return true; } return false; } if( !AlmostEqual( xl, xr ) ) { // this segment is not vertical k = (yr - yl) / (xr - xl); a = yl - k * xl; } if( !AlmostEqual( Xl, Xr ) ){ // hs is not vertical K = (Yr - Yl) / (Xr - Xl); A = Yl - K * Xl; } if( AlmostEqual( Xl, Xr ) ) {//only hs is vertical Coord y0 = k * Xl + a; if( ( Xl > xl || AlmostEqual( Xl, xl ) ) && ( Xl < xr || AlmostEqual( Xl, xr ) ) ){ if( ( ( y0 > Yl || AlmostEqual( y0, Yl ) ) && ( y0 < Yr || AlmostEqual( y0, Yr ) ) ) || ( ( y0 > Yr || AlmostEqual( y0, Yr ) ) && ( y0 < Yl || AlmostEqual( y0, Yl ) ) ) ){ // (Xl, y0) is the intersection point resp.Set( Xl, y0 ); return true; } } return false; } if( AlmostEqual( xl, xr ) ){ // only this segment is vertical Coord Y0 = K * xl + A; if( ( xl > Xl || AlmostEqual( xl, Xl ) ) && ( xl < Xr || AlmostEqual( xl, Xr ) ) ){ if( ( ( Y0 > yl || AlmostEqual( Y0, yl ) ) && ( Y0 < yr || AlmostEqual( Y0, yr ) ) ) || ( ( Y0 > yr || AlmostEqual( Y0, yr ) ) && ( Y0 < yl || AlmostEqual( Y0, yl ) ) ) ){ // (xl, Y0) is the intersection point resp.Set( xl, Y0 ); return true; } } return false; } // both segments are non-vertical if( AlmostEqual( k, K ) ){ // both segments have the same inclination if( AlmostEqual( rp, hs.lp ) ){ resp = rp; return true; } if( AlmostEqual( lp, hs.rp ) ){ resp = lp; return true; } return false; } Coord x0 = (A - a) / (k - K); Coord y0 = x0 * k + a; if( ( x0 > xl || AlmostEqual( x0, xl ) ) && ( x0 < xr || AlmostEqual( x0, xr ) ) && ( x0 > Xl || AlmostEqual( x0, Xl ) ) && ( x0 < Xr || AlmostEqual( x0, Xr ) ) ){ // the segments intersect at (x0, y0) resp.Set( x0, y0 ); return true; } return false; } } bool HalfSegment::Intersection( const HalfSegment& hs, HalfSegment& reshs, const Geoid* geoid /* = 0 */) const { if(geoid) { Point p1s = GetLeftPoint(); Point p1e = GetRightPoint(); Point p2s = hs.GetLeftPoint(); Point p2e = hs.GetRightPoint(); bool ok = p1s.checkGeographicCoord() && p1e.checkGeographicCoord() && p2s.checkGeographicCoord() && p2e.checkGeographicCoord() && AlmostEqual(Distance( hs, geoid ),0.0); if(!ok){ // error OR no intersection (distance >=0) return false; } Point intersectpoint(false); if( Intersection( hs, intersectpoint, geoid) && intersectpoint.IsDefined()){ return false; // intersection is a point; }; if( (p1s.GetX() > 0) && (p1e.GetX() <0) ) { //crosses the +/-180° meridean! Point tmp(p1s); p1s=p1e; p1e=tmp; // swap start/endpoint for *this } if( (p2s.GetX() >0) && (p2e.GetX() <0) ) { // crosses the +/-180° meridean! Point tmp(p2s); p2s=p2e; p2e=tmp; // swap start/endpoint for hs } // no special case: crossing the equator (happens once at most!) Point leftborder(true); if(AlmostEqual(hs.Distance(p1s, geoid), 0.0)) {// LeftPoint on hs? leftborder = p1s; } else if(AlmostEqual(Distance(p2s, geoid), 0.0)) {// hs.LeftPoint on *this? leftborder = p2s; } else { return false; } Point rightborder(true); if(AlmostEqual(hs.Distance(p1e, geoid), 0.0)) {// RightPoint on hs? leftborder = p1e; } else if(AlmostEqual(Distance(p2e, geoid), 0.0)){// hs.RightPoint on *this? leftborder = p2e; } else { return false; } reshs.Set(true, leftborder, rightborder); return true; } // else: euclidean geometry double k, a, K, A; if( !BoundingBox().Intersects( hs.BoundingBox() ) ) return false; if( AlmostEqual( *this, hs ) ) { reshs = hs; return true; } Coord xl = lp.GetX(), yl = lp.GetY(), xr = rp.GetX(), yr = rp.GetY(), Xl = hs.GetLeftPoint().GetX(), Yl = hs.GetLeftPoint().GetY(), Xr = hs.GetRightPoint().GetX(), Yr = hs.GetRightPoint().GetY(); if( AlmostEqual( xl, xr ) && AlmostEqual( Xl, Xr ) ) //both l and L are vertical lines { Coord ylow, yup, Ylow, Yup; if( yl < yr ) { ylow = yl; yup = yr; } else { ylow = yr; yup = yl; } if( Yl < Yr ) { Ylow = Yl; Yup = Yr; } else { Ylow = Yr; Yup = Yl; } if( yup > Ylow && ylow < Yup && !AlmostEqual( yup, Ylow ) && !AlmostEqual( ylow, Yup ) ) { Point p1(true), p2(true); if( ylow > Ylow ) p1.Set( xl, ylow ); else p1.Set( xl, Ylow ); if( yup < Yup ) p2.Set( xl, yup ); else p2.Set( xl, Yup ); reshs.Set( true, p1, p2 ); return true; } else return false; } if( AlmostEqual( Xl, Xr ) || AlmostEqual( xl, xr ) ) //only L or l is vertical return false; //otherwise: both *this and *arg are non-vertical lines k = (yr - yl) / (xr - xl); a = yl - k * xl; K = (Yr - Yl) / (Xr - Xl); A = Yl - K * Xl; if( AlmostEqual( k, K ) && AlmostEqual( A, a ) ) { if( xr > Xl && xl < Xr && !AlmostEqual( xr, Xl ) && !AlmostEqual( xl, Xr ) ) { Point p1(true), p2(true); if( xl > Xl ) p1.Set( xl, yl ); else p1.Set( Xl, Yl ); if( xr < Xr ) p2.Set( xr, yr ); else p2.Set( Xr, Yr ); reshs.Set( true, p1, p2 ); return true; } } return false; } bool HalfSegment::Crosses( const HalfSegment& hs, const Geoid* geoid/*=0*/ ) const { if(geoid){ cerr << __PRETTY_FUNCTION__ << ": Spherial geometry not implemented." << endl; assert( false ); // TODO: implement spherical case. } double k = 0.0; double a = 0.0; double K = 0.0; double A = 0.0; if( !BoundingBox().Intersects( hs.BoundingBox(),geoid ) ){ return false; } Coord xl = lp.GetX(), yl = lp.GetY(), xr = rp.GetX(), yr = rp.GetY(), Xl = hs.GetLeftPoint().GetX(), Yl = hs.GetLeftPoint().GetY(), Xr = hs.GetRightPoint().GetX(), Yr = hs.GetRightPoint().GetY(); if( AlmostEqual( xl, xr ) && AlmostEqual( Xl, Xr ) ) // both segments are vertical return false; if( !AlmostEqual( xl, xr ) ) // this segment is not vertical { k = (yr - yl) / (xr - xl); a = yl - k * xl; } if( !AlmostEqual( Xl, Xr ) ) // hs is not vertical { K = (Yr - Yl) / (Xr - Xl); A = Yl - K * Xl; } if( AlmostEqual( Xl, Xr ) ) //only hs is vertical { double y0 = k * Xl + a; if( Xl > xl && !AlmostEqual( Xl, xl ) && Xl < xr && !AlmostEqual( Xl, xr ) ) { if( ( y0 > Yl && !AlmostEqual( y0, Yl ) && y0 < Yr && !AlmostEqual( y0, Yr ) ) || ( y0 > Yr && !AlmostEqual( y0, Yr ) && y0 < Yl && !AlmostEqual( y0, Yl ) ) ) // (Xl, y0) is the intersection point return true; } return false; } if( AlmostEqual( xl, xr ) ) // only this segment is vertical { double Y0 = K * xl + A; if( ( xl > Xl && !AlmostEqual( xl, Xl ) ) && ( xl < Xr && !AlmostEqual( xl, Xr ) ) ) { if( ( Y0 > yl && !AlmostEqual( Y0, yl ) && Y0 < yr && !AlmostEqual( Y0, yr ) ) || ( Y0 > yr && !AlmostEqual( Y0, yr ) && Y0 < yl && !AlmostEqual( Y0, yl ) ) ) // (xl, Y0) is the intersection point return true; } return false; } // both segments are non-vertical if( AlmostEqual( k, K ) ) // both segments have the same inclination return false; double x0 = (A - a) / (k - K); // y0 = x0 * k + a; if( x0 > xl && !AlmostEqual( x0, xl ) && x0 < xr && !AlmostEqual( x0, xr ) && x0 > Xl && !AlmostEqual( x0, Xl ) && x0 < Xr && !AlmostEqual( x0, Xr ) ) // the segments intersect at (x0, y0) return true; return false; } bool HalfSegment::Inside( const HalfSegment& hs, const Geoid* geoid/*=0*/ ) const { return hs.Contains( GetLeftPoint(),geoid ) && hs.Contains( GetRightPoint(),geoid ); } bool HalfSegment::Contains( const Point& p, const Geoid* geoid/*=0*/ ) const{ if( !p.IsDefined() || (geoid && !geoid->IsDefined()) ) { assert( p.IsDefined() ); return false; } if(geoid){ cerr << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." << endl; assert(!geoid); } if( AlmostEqual( p, lp ) || AlmostEqual( p, rp ) ){ return true; } Coord xl = lp.GetX(), yl = lp.GetY(), xr = rp.GetX(), yr = rp.GetY(), x = p.GetX(), y = p.GetY(); if( xl != xr && xl != x ) // the segment is not vertical { double k1 = (y - yl) / (x - xl), k2 = (yr - yl) / (xr - xl); if( AlmostEqual( k1, k2 ) ) { if( ( x > xl || AlmostEqual( x, xl ) ) && ( x < xr || AlmostEqual( x, xr ) ) ) // we check only this possibility because lp < rp and // therefore, in this case, xl < xr return true; } } else if( AlmostEqual( xl, xr ) && AlmostEqual( xl, x ) ) // the segment is vertical and the point is also in the // same x-position. In this case we just have to check // whether the point is inside the y-interval { if( (( y > yl || AlmostEqual( y, yl ) ) && ( y < yr || AlmostEqual( y, yr ) )) || (( y < yl || AlmostEqual( y, yl ) ) && ( y > yr || AlmostEqual( y, yr ) ) )) // Here we check both possibilities because we do not // know wheter yl < yr, given that we used the // AlmostEqual function in the previous condition return true; } return false; } double HalfSegment::Distance( const Point& p, const Geoid* geoid /* = 0 */) const { assert( p.IsDefined() ); assert( !geoid || geoid->IsDefined() ); if(!geoid){ // euclidean geometry Coord xl = GetLeftPoint().GetX(), yl = GetLeftPoint().GetY(), xr = GetRightPoint().GetX(), yr = GetRightPoint().GetY(), X = p.GetX(), Y = p.GetY(); double result, auxresult; if( xl == xr || yl == yr ) { if( xl == xr) //hs is vertical { if( (yl <= Y && Y <= yr) || (yr <= Y && Y <= yl) ) result = fabs( X - xl ); else { result = p.Distance( GetLeftPoint() ); auxresult = p.Distance( GetRightPoint() ); if( result > auxresult ) result = auxresult; } } else //hs is horizontal line: (yl==yr) { if( xl <= X && X <= xr ) result = fabs( Y - yl ); else { result = p.Distance( GetLeftPoint() ); auxresult = p.Distance( GetRightPoint() ); if( result > auxresult ) result = auxresult; } } } else { double k = (yr - yl) / (xr - xl), a = yl - k * xl, xx = (k * (Y - a) + X) / (k * k + 1), yy = k * xx + a; Coord XX = xx, YY = yy; Point PP( true, XX, YY ); if( xl <= XX && XX <= xr ) result = p.Distance( PP ); else { result = p.Distance( GetLeftPoint() ); auxresult = p.Distance( GetRightPoint() ); if( result > auxresult ) result = auxresult; } } return result; } Point p1 = GetLeftPoint(); Point p2 = GetRightPoint(); return geodist::getDist(p1.GetX(), p1.GetY(), p2.GetX(), p2.GetY(), p.GetX(), p.GetY(), geoid); /* // else: spherical geometry bool ok = true; double d13 = p.DistanceOrthodrome(GetLeftPoint(),*geoid,ok); // initial bearing from LeftPoint to p double theta13 = GetLeftPoint().Direction(p,false,geoid); // in deg // initial bearing from LeftPoint to RightPoint double theta12 = GetLeftPoint().Direction(GetRightPoint(),false,geoid);//[deg] double R = geoid->getR(); errno = 0; double res = asin(sin(d13/R)*sin((theta13-theta12)*M_PI/180.0))*R; ok = ok && (errno==0) && (theta13>=0.0) && (theta12>=0.0); if(ok){ return res; } // else: error return -666.666; */ } double HalfSegment::Distance( const HalfSegment& hs, const Geoid* geoid /* = 0 */) const { assert( !geoid || geoid->IsDefined() ); if(!geoid){ // euclidean geometry if( Intersects( hs ) ){ return 0.0; } double d1 = MIN( Distance( hs.GetLeftPoint(), geoid ), Distance( hs.GetRightPoint(), geoid ) ); double d2 = MIN( hs.Distance(this->GetLeftPoint(), geoid), hs.Distance(this->GetRightPoint(), geoid)); return MIN(d1,d2); } // else: spherical geometry cout << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." <& rect, const Geoid* geoid/*=0*/) const{ assert( rect.IsDefined() ); assert( !geoid || geoid->IsDefined() ); if(geoid){ cout << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." <Distance(p0,geoid); } else { hs.Set(true,p0,p1); dist = this->Distance(hs,geoid); } if(AlmostEqual(dist,0)){ return 0.0; } if(AlmostEqual(p1,p2)){ dist = MIN( dist, this->Distance(p1,geoid)); } else { hs.Set(true,p1,p2); dist = MIN( dist, this->Distance(hs,geoid)); } if(AlmostEqual(dist,0)){ return 0.0; } if(AlmostEqual(p2,p3)){ dist = MIN(dist, this->Distance(p2,geoid)); } else { hs.Set(true,p2,p3); dist = MIN( dist, this->Distance(hs,geoid)); } if(AlmostEqual(dist,0)){ return 0.0; } if(AlmostEqual(p3,p0)){ dist = MIN(dist, this->Distance(p3,geoid)); } else { hs.Set(true,p3,p0); dist = MIN( dist, this->Distance(hs,geoid)); } if(AlmostEqual(dist,0)){ return 0.0; } return dist; } bool HalfSegment::Intersects(const Rectangle<2>& rect, const Geoid* geoid/*=0*/) const{ if(!rect.IsDefined()){ return false; } if(rect.IsEmpty()){ return false; } if(lp.Intersects(rect,geoid)) return true; if(rp.Intersects(rect,geoid)) return true; // check for intersection of the 4 // segments of the rectangle Point p1(true, rect.MinD(0), rect.MinD(1)); Point p2(true, rect.MaxD(0), rect.MinD(1)); Point p3(true, rect.MaxD(0), rect.MaxD(1)); Point p4(true, rect.MinD(0), rect.MaxD(1)); HalfSegment hs1(true,p1,p2); if(Intersects(hs1)) return true; HalfSegment hs2(true,p2,p3); if(Intersects(hs2)) return true; HalfSegment hs3(true,p3,p4); if(Intersects(hs3)) return true; HalfSegment hs4(true,p4,p1); if(Intersects(hs4)) return true; return false; } double HalfSegment::MaxDistance(const Rectangle<2>& rect, const Geoid* geoid /*=0*/) const{ assert( rect.IsDefined() ); assert( !geoid || geoid->IsDefined() ); if(geoid){ cout << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." < y ) { abovey0 = yl; return true; } else if( xl < x && x <= xr ) { double k = (yr - yl) / (xr - xl), a = (yl - k * xl), y0 = k * x + a; Coord yy = y0; if( yy > y ) { abovey0 = y0; return true; } } } return false; } typedef unsigned int outcode; enum { TOP = 0x1, BOTTOM = 0x2, RIGHT = 0x4, LEFT = 0x8 }; outcode CompOutCode( double x, double y, double xmin, double xmax, double ymin, double ymax) { outcode code = 0; if (y > ymax) code |=TOP; else if (y < ymin) code |= BOTTOM; if ( x > xmax) code |= RIGHT; else if ( x < xmin) code |= LEFT; return code; } void HalfSegment::CohenSutherlandLineClipping( const Rectangle<2> &window, double &x0, double &y0, double &x1, double &y1, bool &accept, const Geoid* geoid/*=0*/ ) const { assert( window.IsDefined() ); assert( !geoid || geoid->IsDefined() ); if(geoid){ cerr << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." << endl; assert( false ); // TODO: implement spherical geometry case } // Outcodes for P0, P1, and whatever point lies outside the clip rectangle*/ outcode outcode0, outcode1, outcodeOut; double xmin = window.MinD(0) , xmax = window.MaxD(0), ymin = window.MinD(1), ymax = window.MaxD(1); bool done = false; accept = false; outcode0 = CompOutCode( x0, y0, xmin, xmax, ymin, ymax); outcode1 = CompOutCode( x1, y1, xmin, xmax, ymin, ymax); do{ if ( !(outcode0 | outcode1) ) { //"Trivial accept and exit"< &window, HalfSegment &hsInside, bool &inside, bool &isIntersectionPoint, Point &intersectionPoint, const Geoid* geoid/*=0*/) const { if( !window.IsDefined() ) { intersectionPoint.SetDefined( false ); assert( window.IsDefined() ); } assert( !geoid || geoid->IsDefined() ); if(geoid){ cerr << __PRETTY_FUNCTION__ << ": Spherical geometry not implemented." << endl; assert( false ); // TODO: implement spherical geometry case } double x0 = GetLeftPoint().GetX(), y0 = GetLeftPoint().GetY(), x1 = GetRightPoint().GetX(), y1 = GetRightPoint().GetY(); CohenSutherlandLineClipping(window, x0, y0, x1, y1, inside); isIntersectionPoint=false; intersectionPoint.SetDefined( false ); if (inside) { Point lp( true, x0, y0 ), rp(true, x1, y1 ); if (lp==rp) { intersectionPoint.SetDefined( true ); isIntersectionPoint = true; intersectionPoint=lp; } else { AttrType attr=this->GetAttr(); hsInside.Set(true, rp, lp); hsInside.SetAttr(attr); } } } /* 6.2 List Representation The list representation of a HalfSegment is ---- ( bool (lp rp)) ---- where the bool value indicate whether the dominating point is the left point. 6.3 ~In~ and ~Out~ Functions */ ListExpr OutHalfSegment( ListExpr typeInfo, Word value ) { HalfSegment* hs; hs = (HalfSegment*)(value.addr); Point lp = hs->GetLeftPoint(), rp = hs->GetRightPoint(); return nl->TwoElemList( nl-> BoolAtom(hs->IsLeftDomPoint() ), nl->TwoElemList( OutPoint( nl->TheEmptyList(), SetWord( &lp ) ), OutPoint( nl->TheEmptyList(), SetWord( &rp) ) ) ); } Word InHalfSegment( const ListExpr typeInfo, const ListExpr instance, const int errorPos, ListExpr& errorInfo, bool& correct ) { if ( nl->ListLength( instance ) == 2 ) { ListExpr first = nl->First(instance), second = nl->Second(instance), firstP = nl->TheEmptyList(), secondP = nl->TheEmptyList(); bool ldp; if( nl->IsAtom(first) && nl->AtomType(first) == BoolType ) { ldp = nl->BoolValue(first); if( nl->ListLength(second) == 2 ) { firstP = nl->First(second); secondP = nl->Second(second); } correct=true; Point *lp = (Point*)InPoint(nl->TheEmptyList(), firstP, 0, errorInfo, correct ).addr; if( correct && lp->IsDefined() ) { Point *rp = (Point*)InPoint(nl->TheEmptyList(), secondP, 0, errorInfo, correct ).addr; if( correct && rp->IsDefined() && *lp != *rp ) { HalfSegment *hs = new HalfSegment( ldp, *lp, *rp ); delete lp; delete rp; return SetWord( hs ); } delete rp; } delete lp; } } correct = false; return SetWord( Address(0) ); } /* ~getNext~ auxiliary function for Line::Transform */ templateclass Array> int getNext(const Array* hss, int pos, const char* usage){ HalfSegment hs; hss->Get(pos,hs); pos = hs.attr.partnerno; hss->Get(pos,hs); Point dp = hs.GetDomPoint(); int res = -1; // go backwards from pos, store a candidate as res int pos1 = pos-1; bool done = false; // search left of pos while( (pos1>=0) && ! done){ hss->Get(pos1,hs); Point dp1 = hs.GetDomPoint(); if(!AlmostEqual(dp,dp1)){ done = true; } else { if(usage[pos1] == 3){ return pos1; } else if(usage[pos1]==0){ res = pos1; } pos1--; } } // search right of pos pos1 = pos+1; done = false; while((pos1Size()) && ! done){ hss->Get(pos1,hs); Point dp1 = hs.GetDomPoint(); if(!AlmostEqual(dp,dp1)){ done = true; } else { if(usage[pos1] == 3){ return pos1; } else if(usage[pos1]==0){ res = pos1; } pos1++; } } return res; } /* ~markUsage~ This function marks all HalfSgements of a line with its usage. This means: 0 : Part of an unique cycle 1 : not Part of a cycle 2 : Part of an ambiguous cycle */ templateclass Array> void markUsage( LineT* line, char* usage, char* critical){ markUsage( (Array*) line->GetFLOB(0),usage, critical); } templateclass Array> void markUsage(const Array* hss, char* usage, char* critical ){ // step 1: mark halfsegments not belonging to a cycle memset(usage,0,hss->Size()); // meaning of elements of usage // 0 : not used => part of a cycle // 1 : not part of a cycle // 2 : part of a cycle // 3 : part of the current path (will be set to 1 or 2 later) // 4 : partner of a segment in the current path memset(critical,0,hss->Size()); // 1.1 mark critical points int count = 0; HalfSegment hs; Point lastPoint; for(int i=0;iSize();i++){ hss->Get(i,hs); Point currentPoint = hs.GetDomPoint(); if(i==0){ lastPoint = currentPoint; count = 1; } else { if(AlmostEqual(lastPoint,currentPoint)){ count++; } else { if(count != 2){ for(int r=1;r<=count-1;r++){ critical[i-r] = 1; } } lastPoint = currentPoint; count = 1; } } } if(count != 2){ for(int r=1;r<=count;r++){ critical[hss->Size()-r] = 1; } } // 1.2 mark segments for(int i=0;iSize();i++) { if(usage[i]==0 ){ // not used, but critical hss->Get(i,hs); vector path; usage[i] = 3; // part of current path usage[hs.attr.partnerno] = 4; // extend path int pos = i; path.push_back(pos); int next = getNext(hss,pos,usage); while(!path.empty()) { if(next < 0){ // no extension found => mark segments in path as non-cycle // go back until the next critical point is reachedA bool done = path.empty(); while(!done){ if(path.empty()){ done = true; } else { pos = path.back(); path.pop_back(); usage[pos] = 1; hss->Get(pos,hs); usage[hs.attr.partnerno] = 1; if(critical[pos]){ done = true; } else { } } } if(!path.empty()){ pos = path.back(); next = getNext(hss,pos,usage); } } else { if(usage[next] == 3){ // (sub) path found int p = path.back(); path.pop_back(); while(p!=next && !path.empty()){ usage[p] = 2; hss->Get(p,hs); usage[hs.attr.partnerno] = 2; p = path.back(); path.pop_back(); } usage[p] = 2; hss->Get(p,hs); usage[hs.attr.partnerno] = 2; if(!path.empty()){ pos = path.back(); // try to extend first part of path next = getNext(hss,pos,usage); } } else { // normal extension pos = next; path.push_back(pos); usage[pos] = 3; hss->Get(pos,hs); usage[hs.attr.partnerno] = 4; next = getNext(hss,pos,usage); } } } // while path not empty } } } /* Simple function marking some elements between min and max as used to avoid segments degenerated to single points. */ static double maxDist(const vector& orig, // points const int min, const int max, // range int& index, // resultindex const Geoid* geoid=0){ // search the point with the largest distance to the segment between // orig[min] and orig[max] double maxdist = 0.0; int maxindex = -1; try{ if(!AlmostEqual(orig.at(min),orig.at(max))){ HalfSegment hs(true,orig.at(min),orig.at(max)); for(int i=min+1; imaxdist){ maxdist = dist; maxindex = i; } } } else { // special case of a cycle Point p = orig[min]; for(int i=min+1; imaxdist){ maxdist = dist; maxindex = i; } } } index = maxindex; return maxdist; } catch (out_of_range&){ cerr << "min=" << min << " max=" << max << " size=" << orig.size() << endl; assert(false); return -1.0; } } /* Implementation of the Douglas Peucker algorithm. The return indoicates whether between min and max further points was used. */ bool douglas_peucker(const vector& orig, // original line const double epsilon, // maximum derivation bool* use, // result const int min, const int max, const bool force , const Geoid* geoid){ // current range // always use the endpoints use[min] = true; use[max] = true; if(min+1>=max){ // no inner points, nothing to do return false; } int index; double maxdist = maxDist(orig,min,max,index,geoid); bool cycle = AlmostEqual(orig[min], orig[max]); if( (maxdist<=epsilon) && // line closed enough !cycle && // no degenerated segment !force){ return false; // all ok, stop recursion } else { bool ins = douglas_peucker(orig,epsilon,use,min,index,cycle,geoid); if(index>=0){ douglas_peucker(orig,epsilon,use,index,max,cycle && !ins,geoid); } return true; } } void douglas_peucker(const vector& orig, // original chain of points const double epsilon, // maximum derivation bool* use, // result const bool forceThrow, const Geoid* geoid ){ for(unsigned int i=0;iIsDefined()){ return nl->SymbolAtom(Symbol::UNDEFINED()); } if( l->IsEmpty() ) return nl->TheEmptyList(); result = nl->TheEmptyList(); last = result; bool first = true; for( int i = 0; i < l->Size(); i++ ) { l->Get( i, hs ); if( hs.IsLeftDomPoint() == true ) { halfseg = OutHalfSegment( nl->TheEmptyList(), SetWord( (void*)&hs ) ); halfpoints = nl->Second( halfseg ); flatseg = nl->FourElemList( nl->First( nl->First( halfpoints ) ), nl->Second( nl->First( halfpoints ) ), nl->First( nl->Second( halfpoints ) ), nl->Second( nl->Second( halfpoints ) ) ); if( first == true ) { result = nl->OneElemList( flatseg ); last = result; first = false; } else last = nl->Append( last, flatseg ); } } return result; } /* 7.4 ~In~-function */ Word InLine( const ListExpr typeInfo, const ListExpr instance, const int errorPos, ListExpr& errorInfo, bool& correct ) { Line* l = new Line( 0 ); if(listutils::isSymbolUndefined(instance)){ l->SetDefined(false); correct=true; return SetWord(Address(l)); } HalfSegment * hs; l->StartBulkLoad(); ListExpr first, halfseg, halfpoint; ListExpr rest = instance; int edgeno = 0; if( !nl->IsAtom( instance ) ) { while( !nl->IsEmpty( rest ) ) { first = nl->First( rest ); rest = nl->Rest( rest ); if( nl->ListLength( first ) == 4 ) { halfpoint = nl->TwoElemList( nl->TwoElemList( nl->First(first), nl->Second(first)), nl->TwoElemList( nl->Third(first), nl->Fourth(first))); } else { // wrong list representation l->DeleteIfAllowed(); correct = false; return SetWord( Address(0) ); } halfseg = nl->TwoElemList(nl->BoolAtom(true), halfpoint); hs = (HalfSegment*)InHalfSegment( nl->TheEmptyList(), halfseg, 0, errorInfo, correct ).addr; if( correct ) { hs->attr.edgeno = edgeno++; *l += *hs; hs->SetLeftDomPoint( !hs->IsLeftDomPoint() ); *l += *hs; } delete hs; } l->EndBulkLoad(); correct = true; return SetWord( l ); } l->DeleteIfAllowed(); correct = false; return SetWord( Address(0) ); } /* 7.5 ~Create~-function */ Word CreateLine( const ListExpr typeInfo ) { return SetWord( new Line( 0 ) ); } /* 7.6 ~Delete~-function */ void DeleteLine( const ListExpr typeInfo, Word& w ) { Line *l = (Line *)w.addr; l->Destroy(); l->DeleteIfAllowed(false); w.addr = 0; } /* 7.7 ~Close~-function */ void CloseLine( const ListExpr typeInfo, Word& w ) { ((Line *)w.addr)->DeleteIfAllowed(); w.addr = 0; } /* 7.8 ~Clone~-function */ Word CloneLine( const ListExpr typeInfo, const Word& w ) { return SetWord( new Line( *((Line *)w.addr) ) ); } /* 7.8 ~Open~-function */ bool OpenLine( SmiRecord& valueRecord, size_t& offset, const ListExpr typeInfo, Word& value ) { Line *l = (Line*)Attribute::Open( valueRecord, offset, typeInfo ); value = SetWord( l ); return true; } /* 7.8 ~Save~-function */ bool SaveLine( SmiRecord& valueRecord, size_t& offset, const ListExpr typeInfo, Word& value ) { Line *l = (Line *)value.addr; Attribute::Save( valueRecord, offset, typeInfo, l ); return true; } /* 7.9 ~SizeOf~-function */ int SizeOfLine() { return sizeof(Line); } /* 7.11 Function describing the signature of the type constructor */ ListExpr LineProperty() { return nl->TwoElemList( nl->FourElemList( nl->StringAtom("Signature"), nl->StringAtom("Example Type List"), nl->StringAtom("List Rep"), nl->StringAtom("Example List")), nl->FourElemList( nl->StringAtom("-> DATA"), nl->StringAtom(Line::BasicType()), nl->StringAtom("(*) where segment is " "()"), nl->StringAtom("( (1 1 2 2)(3 3 4 4) )"))); } /* 7.12 Kind checking function This function checks whether the type constructor is applied correctly. Since type constructor ~line~ does not have arguments, this is trivial. */ bool CheckLine( ListExpr type, ListExpr& errorInfo ) { return nl->IsEqual( type, Line::BasicType() ); } /* 7.13 ~Cast~-function */ void* CastLine(void* addr) { return Line::Cast(addr); } /* 7.14 Creation of the type constructor instance */ TypeConstructor line( Line::BasicType(), //name LineProperty, //describing signature OutLine, InLine, //Out and In functions 0, 0, //SaveTo and RestoreFrom List functions CreateLine, DeleteLine, //object creation and deletion OpenLine, SaveLine, // object open and save CloseLine, CloneLine, //object close and clone CastLine, //cast function SizeOfLine, //sizeof function CheckLine ); //kind checking function Word InSimpleLine( const ListExpr typeInfo, const ListExpr instance1, const int errorPos, ListExpr& errorInfo, bool& correct ){ ListExpr instance = instance1; correct = true; if(listutils::isSymbolUndefined(instance)){ SimpleLine* line = new SimpleLine( 0 ); line->SetDefined(false); return SetWord(Address(line)); } if(nl->AtomType(instance)!=NoAtom){ correct=false; return SetWord(Address(0)); } bool startSmaller = true; if( nl->HasLength(instance,2) && (nl->AtomType(nl->Second(instance))==BoolType)){ startSmaller = nl->BoolValue(nl->Second(instance)); instance = nl->First(instance); } HalfSegment* hs; SimpleLine* line= new SimpleLine(10); int edgeno = 0; ListExpr rest = instance; line->StartBulkLoad(); while(!nl->IsEmpty(rest)){ ListExpr segment = nl->First(rest); if(!nl->HasLength(segment,4)){ correct=false; line->DeleteIfAllowed(); return SetWord(Address(0)); } ListExpr halfSegment = nl->TwoElemList( nl->BoolAtom(true), nl->TwoElemList( nl->TwoElemList(nl->First(segment), nl->Second(segment)), nl->TwoElemList(nl->Third(segment), nl->Fourth(segment)))); hs = static_cast( InHalfSegment( nl->TheEmptyList(), halfSegment, 0, errorInfo, correct ).addr); if(!correct){ delete hs; return SetWord(Address(0)); } hs->attr.edgeno = edgeno++; *line += *hs; hs->SetLeftDomPoint( !hs->IsLeftDomPoint() ); *line += *hs; delete hs; rest = nl->Rest(rest); } if(!line->EndBulkLoad()){ correct = false; line->DeleteIfAllowed(); return SetWord(Address(0)); }else{ correct = true; line->SetStartSmaller(startSmaller); return SetWord(line); } } ostream& operator<<(ostream& out, const LRS& lrs){ out << "lrsPos "<< lrs.lrsPos << ", hsPos " << lrs.hsPos; return out; } ListExpr OutSimpleLine( ListExpr typeInfo, Word value ) { ListExpr result, last; HalfSegment hs; ListExpr halfseg, halfpoints, flatseg; SimpleLine* l = static_cast(value.addr); if(!l->IsDefined()){ return nl->SymbolAtom(Symbol::UNDEFINED()); } if( l->IsEmpty() ){ return nl->TwoElemList( nl->TheEmptyList(), nl->BoolAtom( l->GetStartSmaller() ) ); } result = nl->TheEmptyList(); last = result; bool first = true; /* // version using halfsegment order for( int i = 0; i < l->Size(); i++ ) { l->Get( i, hs ); if( hs.IsLeftDomPoint() ){ halfseg = OutHalfSegment( nl->TheEmptyList(), SetWord( (void*)&hs ) ); halfpoints = nl->Second( halfseg ); flatseg = nl->FourElemList( nl->First( nl->First( halfpoints ) ), nl->Second( nl->First( halfpoints ) ), nl->First( nl->Second( halfpoints ) ), nl->Second( nl->Second( halfpoints ) ) ); if( first == true ) { result = nl->OneElemList( flatseg ); last = result; first = false; } else { last = nl->Append( last, flatseg ); } } } */ // version using lrs order LRS lrs; // for some curious reason, some halfsegments occur twice // in the lrs array set used; for(int i=0;ilrsSize(); i++){ l->Get(i,lrs); if(used.find(lrs.hsPos)==used.end()){ used.insert(lrs.hsPos); //cout << "LRS " << lrs << endl; l->Get(lrs.hsPos,hs); halfseg = OutHalfSegment( nl->TheEmptyList(), SetWord( (void*)&hs ) ); halfpoints = nl->Second( halfseg ); flatseg = nl->FourElemList( nl->First( nl->First( halfpoints ) ), nl->Second( nl->First( halfpoints ) ), nl->First( nl->Second( halfpoints ) ), nl->Second( nl->Second( halfpoints ) ) ); if( first == true ) { result = nl->OneElemList( flatseg ); last = result; first = false; } else { last = nl->Append( last, flatseg ); } } } return nl->TwoElemList( result, nl->BoolAtom( l->GetStartSmaller() ) ); } Word CreateSimpleLine( const ListExpr typeInfo ) { return SetWord( new SimpleLine( 0 ) ); } void DeleteSimpleLine( const ListExpr typeInfo, Word& w ) { SimpleLine *l = static_cast(w.addr); l->Destroy(); l->DeleteIfAllowed(false); w.addr = 0; } void CloseSimpleLine( const ListExpr typeInfo, Word& w ) { (static_cast(w.addr))->DeleteIfAllowed(); w.addr = 0; } Word CloneSimpleLine( const ListExpr typeInfo, const Word& w ) { return SetWord( new SimpleLine( *((SimpleLine *)w.addr) ) ); } int SizeOfSimpleLine() { return sizeof(SimpleLine); } ListExpr SimpleLineProperty() { return nl->TwoElemList( nl->FourElemList( nl->StringAtom("Signature"), nl->StringAtom("Example Type List"), nl->StringAtom("List Rep"), nl->StringAtom("Example List")), nl->FourElemList( nl->StringAtom("-> DATA"), nl->StringAtom(SimpleLine::BasicType()), nl->TextAtom("( (*) ) where segment is " "() and bool is TRUE or FALSE"), nl->StringAtom("( ( (1 1 2 2) (2 2 1 4) ) FALSE )"))); } bool CheckSimpleLine( ListExpr type, ListExpr& errorInfo ){ return nl->IsEqual( type, SimpleLine::BasicType() ); } void* CastSimpleLine(void* addr) { return SimpleLine::Cast(addr);; } TypeConstructor sline( SimpleLine::BasicType(), //name SimpleLineProperty, //describing signature OutSimpleLine, InSimpleLine, //Out and In functions 0, 0, //SaveTo and RestoreFrom List functions CreateSimpleLine, DeleteSimpleLine, //object creation and deletion OpenAttribute, SaveAttribute, // object open and save CloseSimpleLine, CloneSimpleLine, //object close and clone CastSimpleLine, //cast function SizeOfSimpleLine, //sizeof function CheckSimpleLine); bool IsInsideAbove(const HalfSegment& hs, const Point& thirdPoint){ if(AlmostEqual(hs.GetLeftPoint().GetX(),hs.GetRightPoint().GetX())){ // vertical segment return (MIN(hs.GetLeftPoint().GetX(),hs.GetRightPoint().GetX())IsDefined() ); if(geoid){ double size = p1.Distance(p1, geoid); assert( size > 0 ); return size; }else { double size = pow( (p1.GetX() - p2.GetX()),2) + pow( (p1.GetY() - p2.GetY()),2); size = sqrt(size); return size; } } //The angle function returns the angle of VP1P2 // P1 is the point on the window's edge double Angle(const Point &v, const Point &p1,const Point &p2) { assert( v.IsDefined() ); assert( p1.IsDefined() ); assert( p2.IsDefined() ); double coss; //If P1P2 is vertical and the window's edge // been tested is horizontal , then //the angle VP1P2 is equal to 90 degrees. On the //other hand, if P1P2 is vertical //and the window's edge been tested is vertical, then // the angle is 90 degrees. //Similar tests are applied when P1P2 is horizontal. if (p1.GetX() == p2.GetX()){ //the segment is vertical if (v.GetY()==p1.GetY()){ return M_PI/2; //horizontal edge } else { return 0; } } if (p1.GetY() == p2.GetY()){ //the segment is horizontal if (v.GetY()==p1.GetY()){ return 0; //horizontal edge } else { return M_PI/2; } } coss = double( ( (v.GetX() - p1.GetX()) * (p2.GetX() - p1.GetX()) ) + ( (v.GetY() - p1.GetY()) * (p2.GetY() - p1.GetY()) ) ) / (VectorSize(v,p1) * VectorSize(p2,p1)); //cout<0 ) //p2.GetX() is located to the left of p.GetX() direction = false; //RIGHT else direction = true; //LEFT } else { if ( (p.GetX()-p2.GetX())>0 ) //p2.GetX() is located to the right of p.GetX() direction = true; //LEFT else direction = false; //RIGHT } } return new EdgePoint(p,direction,reject); } void AddPointToEdgeArray( const Point &p, const HalfSegment &hs, const Rectangle<2> &window, vector pointsOnEdge[4] ) { EdgePoint *dp; Point v(true); AttrType attr; attr = hs.GetAttr(); Point p2(true); //If the left and right edges are been tested then //it is not need to check the angle //between the half segment and the edge. If the attribute //inside above is true, then //the direction is up (false), otherwise it is down (true). if (p.GetX() == window.MinD(0)) { dp = new EdgePoint(p,!attr.insideAbove,false); pointsOnEdge[WLEFT].push_back(*dp); } else if (p.GetX() == window.MaxD(0)) { dp = new EdgePoint(p,!attr.insideAbove,false); pointsOnEdge[WRIGHT].push_back(*dp); } if (p.GetY() == window.MinD(1)) { v.Set(window.MinD(0), window.MinD(1)); //In this case we don't know which point is outside the window, //so it is need to test both half segment's poinst. Moreover, //in order to use the same comparison that is used for //Top edge, it is need to choose the half segment point that //is over the bottom edge. if (hs.GetLeftPoint().GetY()>window.MinD(1)) dp = EdgePoint::GetEdgePoint(p,hs.GetLeftPoint(), attr.insideAbove,v,false); else dp = EdgePoint::GetEdgePoint(p,hs.GetRightPoint(), attr.insideAbove,v,false); pointsOnEdge[WBOTTOM].push_back(*dp); } else if (p.GetY() == window.MaxD(1)) { v.Set(window.MinD(0), window.MaxD(1)); //In this case we don't know which point is outside the window, //so it is need to test if (hs.GetLeftPoint().GetY()>window.MaxD(1)) dp = EdgePoint::GetEdgePoint(p,hs.GetLeftPoint(), attr.insideAbove,v,false); else dp = EdgePoint::GetEdgePoint(p,hs.GetRightPoint(), attr.insideAbove,v,false); pointsOnEdge[WTOP].push_back(*dp); } } bool GetAcceptedPoint( vector pointsOnEdge, int &i, const int end, EdgePoint &ep ) { //id is the indice of the current point in the scan //ep is the correct edge point that will be returned. ep = pointsOnEdge[i]; //discard all rejected points while (ep.rejected && i<=end) { i++; if (i>end) return false; EdgePoint epAux = pointsOnEdge[i]; //Discard all the points that was accepted but has a // corresponding rejection point. //In other words, point that has the same // coordinates and direction on the edge. if (!epAux.rejected && (epAux.direction==ep.direction) && (epAux.GetX() == ep.GetX()) && (epAux.GetY() == ep.GetY()) ) { while ( (i<=end) && (epAux.direction==ep.direction) && (epAux.GetX() == ep.GetX()) && (epAux.GetY() == ep.GetY()) ) { i++; if (i>end) return false; epAux = pointsOnEdge[i]; } } ep = epAux; } return true; } bool HalfSegment::RayDown( const Point& p, double &yIntersection ) const { if (this->IsVertical()) return false; const Coord& x = p.GetX(), y = p.GetY(), xl = GetLeftPoint().GetX(), yl = GetLeftPoint().GetY(), xr = GetRightPoint().GetX(), yr = GetRightPoint().GetY(); // between is true, iff xl <= x <= xr. const bool between = CompareDouble(x, xl) != -1 && CompareDouble(x, xr) != 1; if (!between) return false; const double k = (yr - yl) / (xr - xl); const double a = (yl - k * xl); const double y0 = k * x + a; if (CompareDouble(y0, y) == 1) // y0 > y: this is above p. return false; // y0 <= p: p is above or on this. yIntersection = y0; return true; } /* ~IsSpatialType~ This function checks whether the type given as a ListExpr is one of ~point~, ~points~, ~line~, or ~region~. */ bool IsSpatialType(ListExpr type){ if(!nl->IsAtom(type)){ return false; } if(nl->AtomType(type)!=SymbolType){ return false; } string t = nl->SymbolValue(type); if(t==Point::BasicType()) return true; if(t==Points::BasicType()) return true; if(t==Line::BasicType()) return true; if(t==Region::BasicType()) return true; return false; } /* This function check whether a region value is valid after the insertion of a new half segment. Whenever a half segment is about to be inserted, the state of the region is checked. A valid region must satisfy the following conditions: 1) any two cycles of the same region must be disconnect, which means that no edges of different cycles can intersect each other; 2) edges of the same cycle can only intersect with their endpoints, but no their middle points; 3) For a certain face, the holes must be inside the outer cycle; 4) For a certain face, any two holes can not contain each other; 5) Faces must have the outer cycle, but they can have no holes; 6) for a certain cycle, any two vertex can not be the same; 7) any cycle must be made up of at least 3 edges; 8) It is allowed that one face is inside another provided that their edges do not intersect. */ vector getCycle(const bool isHole, const vector& vhs){ // first extract the cycle bool used[vhs.size()]; for(unsigned int i=0;i vp; vp.push_back(sp); vp.push_back(cp); used[0] = true; for(unsigned int i=1; i< vhs.size();i++){ for(unsigned int j=1; j vp2; for(int i= vp.size()-1; i>=0; i--){ vp2.push_back(vp[i]); } return vp2; } else { return vp; } } /* Implementation of the Gauss Krueger Projection */ struct P3D{ double x; double y; double z; }; bool WGSGK::project(const Point& src, Point& result) const{ if(!src.IsDefined()){ cerr << __PRETTY_FUNCTION__ << ": Point Argument is undefined!" << endl; result.SetDefined(false); return false; } double x = src.GetX(); double y = src.GetY(); if(x<-180 || x>180 || y<-90 || y>90){ cerr << __PRETTY_FUNCTION__ << ": Point " << src << " is not a valid " << "geographic coordinate!" << endl; result.SetDefined(false); return false; } double a = x*Pi/180; double b = y*Pi/180; if(!useWGS){ BesselBLToGaussKrueger(b, a, result); return result.IsDefined(); } double l1 = a; double b1 = b; a=awgs; b=bwgs; double eq=eqwgs; double N=a/sqrt(1-eq*sin(b1)*sin(b1)); double Xq=(N+h1)*cos(b1)*cos(l1); double Yq=(N+h1)*cos(b1)*sin(l1); double Zq=((1-eq)*N+h1)*sin(b1); P3D p; HelmertTransformation(Xq, Yq, Zq, p); double X = p.x; double Y = p.y; double Z = p.z; a=abes; b=bbes; eq = eqbes; BLRauenberg(X, Y, Z, p); double b2 = p.x; double l2 = p.y; BesselBLToGaussKrueger(b2, l2, result); return result.IsDefined(); } bool WGSGK::project(const HalfSegment& src, HalfSegment& result) const{ result = src; Point p1(true),p2(true); if(!project(src.GetLeftPoint(),p1)) return false; if(!project(src.GetRightPoint(),p2)) return false; if(p2=0.000000000000001); result.SetDefined(true); result.Set((geoDezHeight/Pi)*180, (Latitude/Pi)*180); return result.IsDefined(); } /* 8.2 List Representation The list representation of a region is ---- (face1 face2 face3 ... ) where facei=(outercycle, holecycle1, holecycle2....) cyclei= (vertex1, vertex2, .....) where each vertex is a point. or undef ---- 8.3 ~Out~-function */ ListExpr OutRegion( ListExpr typeInfo, Word value ) { Region* cr = (Region*)(value.addr); if(!cr->IsDefined()){ return nl->SymbolAtom(Symbol::UNDEFINED()); } if( cr->IsEmpty() ) { return (nl->TheEmptyList()); } else { Region *RCopy=new Region(*cr, true); // in memory RCopy->LogicSort(); HalfSegment hs, hsnext; ListExpr regionNL = nl->TheEmptyList(); ListExpr regionNLLast = regionNL; ListExpr faceNL = nl->TheEmptyList(); ListExpr faceNLLast = faceNL; ListExpr cycleNL = nl->TheEmptyList(); ListExpr cycleNLLast = cycleNL; ListExpr pointNL; int currFace = -999999, currCycle= -999999; // avoid uninitialized use Point outputP(true), leftoverP(true); for( int i = 0; i < RCopy->Size(); i++ ) { RCopy->Get( i, hs ); if (i==0) { currFace = hs.attr.faceno; currCycle = hs.attr.cycleno; RCopy->Get( i+1, hsnext ); if ((hs.GetLeftPoint() == hsnext.GetLeftPoint()) || ((hs.GetLeftPoint() == hsnext.GetRightPoint()))) { outputP = hs.GetRightPoint(); leftoverP = hs.GetLeftPoint(); } else if ((hs.GetRightPoint() == hsnext.GetLeftPoint()) || ((hs.GetRightPoint() == hsnext.GetRightPoint()))) { outputP = hs.GetLeftPoint(); leftoverP = hs.GetRightPoint(); } else { cerr << "\n" << __PRETTY_FUNCTION__ << ": Wrong data format --- " << "discontiguous segments!" << endl << "\ths = " << hs << endl << "\thsnext = " << hsnext << endl; return nl->SymbolAtom(Symbol::UNDEFINED()); } pointNL = OutPoint( nl->TheEmptyList(), SetWord(&outputP) ); if (cycleNL == nl->TheEmptyList()) { cycleNL = nl->OneElemList(pointNL); cycleNLLast = cycleNL; } else { cycleNLLast = nl->Append( cycleNLLast, pointNL ); } } else { if (hs.attr.faceno == currFace) { if (hs.attr.cycleno == currCycle) { outputP=leftoverP; if (hs.GetLeftPoint() == leftoverP) leftoverP = hs.GetRightPoint(); else if (hs.GetRightPoint() == leftoverP) { leftoverP = hs.GetLeftPoint(); } else { cerr << "\n" << __PRETTY_FUNCTION__ << ": Wrong data format --- " << "discontiguous segment in cycle!" << endl << "\thh = " << hs << endl << "\tleftoverP = " << leftoverP << endl; return nl->SymbolAtom(Symbol::UNDEFINED()); } pointNL=OutPoint( nl->TheEmptyList(), SetWord( &outputP) ); if (cycleNL == nl->TheEmptyList()) { cycleNL=nl->OneElemList(pointNL); cycleNLLast = cycleNL; } else { cycleNLLast = nl->Append(cycleNLLast, pointNL); } } else { if (faceNL == nl->TheEmptyList()) { faceNL = nl->OneElemList(cycleNL); faceNLLast = faceNL; } else { faceNLLast = nl->Append(faceNLLast, cycleNL); } cycleNL = nl->TheEmptyList(); currCycle = hs.attr.cycleno; RCopy->Get( i+1, hsnext ); if ((hs.GetLeftPoint() == hsnext.GetLeftPoint()) || ((hs.GetLeftPoint() == hsnext.GetRightPoint()))) { outputP = hs.GetRightPoint(); leftoverP = hs.GetLeftPoint(); } else if ((hs.GetRightPoint() == hsnext.GetLeftPoint()) || ((hs.GetRightPoint() == hsnext.GetRightPoint()))) { outputP = hs.GetLeftPoint(); leftoverP = hs.GetRightPoint(); } else { cerr << "\n" << __PRETTY_FUNCTION__ << ": Wrong data format --- " << "discontiguous segments in cycle!" << endl << "\ths = " << hs << endl << "\thsnext = " << hsnext << endl; return nl->SymbolAtom(Symbol::UNDEFINED()); } pointNL = OutPoint( nl->TheEmptyList(), SetWord(&outputP) ); if (cycleNL == nl->TheEmptyList()) { cycleNL = nl->OneElemList(pointNL); cycleNLLast = cycleNL; } else { cycleNLLast = nl->Append(cycleNLLast, pointNL); } } } else { if (faceNL == nl->TheEmptyList()) { faceNL = nl->OneElemList(cycleNL); faceNLLast = faceNL; } else { faceNLLast = nl->Append(faceNLLast, cycleNL); } cycleNL = nl->TheEmptyList(); if (regionNL == nl->TheEmptyList()) { regionNL = nl->OneElemList(faceNL); regionNLLast = regionNL; } else { regionNLLast = nl->Append(regionNLLast, faceNL); } faceNL = nl->TheEmptyList(); currFace = hs.attr.faceno; currCycle = hs.attr.cycleno; RCopy->Get( i+1, hsnext ); if ((hs.GetLeftPoint() == hsnext.GetLeftPoint()) || ((hs.GetLeftPoint() == hsnext.GetRightPoint()))) { outputP = hs.GetRightPoint(); leftoverP = hs.GetLeftPoint(); } else if ((hs.GetRightPoint() == hsnext.GetLeftPoint()) || ((hs.GetRightPoint() == hsnext.GetRightPoint()))) { outputP = hs.GetLeftPoint(); leftoverP = hs.GetRightPoint(); } else { cerr << "\n" << __PRETTY_FUNCTION__ << ": Wrong data format --- " << "discontiguous segments in cycle!" << endl << "\ths = " << hs << endl << "\thsnext = " << hsnext << endl; return nl->SymbolAtom(Symbol::UNDEFINED()); } pointNL = OutPoint(nl->TheEmptyList(), SetWord(&outputP)); if (cycleNL == nl->TheEmptyList()) { cycleNL = nl->OneElemList(pointNL); cycleNLLast = cycleNL; } else { cycleNLLast = nl->Append(cycleNLLast, pointNL); } } } } if (faceNL == nl->TheEmptyList()) { faceNL = nl->OneElemList(cycleNL); faceNLLast = faceNL; } else { faceNLLast = nl->Append(faceNLLast, cycleNL); } cycleNL = nl->TheEmptyList(); if (regionNL == nl->TheEmptyList()) { regionNL = nl->OneElemList(faceNL); regionNLLast = regionNL; } else { regionNLLast = nl->Append(regionNLLast, faceNL); } faceNL = nl->TheEmptyList(); RCopy->DeleteIfAllowed(); return regionNL; } } /* 8.4 ~In~-function */ Word InRegion(const ListExpr typeInfo, const ListExpr instance, const int errorPos, ListExpr& errorInfo, bool& correct ){ if(listutils::isSymbol(instance,Symbol::UNDEFINED())){ Region* r = new Region(0); r->SetDefined(false); correct = true; return SetWord(Address(r)); } if(nl->AtomType(instance) != NoAtom){ correct = false; return SetWord(Address(0)); } ListExpr regNL = instance; vector > cycles; while(!nl->IsEmpty(regNL)){ ListExpr faceNL = nl->First(regNL); regNL = nl->Rest(regNL); if(nl->AtomType(faceNL)!=NoAtom){ correct = false; return SetWord(Address(0)); } bool firstCycle = true; Rectangle<2> faceRect; while(!nl->IsEmpty(faceNL)){ vector cycle; ListExpr cycleNL = nl->First(faceNL); faceNL = nl->Rest(faceNL); if(nl->AtomType(cycleNL)!=NoAtom){ correct=false; return SetWord(Address(0)); } Point firstPoint(false,0,0); bool fp = true; Rectangle<2> currentBB; while(!nl->IsEmpty(cycleNL)){ ListExpr pointNL = nl->First(cycleNL); cycleNL = nl->Rest(cycleNL); if(nl->ListLength(pointNL)!=2){ correct = false; return SetWord(Address(0)); } if(!listutils::isNumeric(nl->First(pointNL)) || !listutils::isNumeric(nl->Second(pointNL))){ correct = false; return SetWord(Address(0)); } Point p(true, listutils::getNumValue(nl->First(pointNL)), listutils::getNumValue(nl->Second(pointNL))); cycle.push_back(p); if(fp){ fp = false; firstPoint = p; currentBB = p.BoundingBox(); } else { currentBB = currentBB.Union(p.BoundingBox()); } } if(!AlmostEqual(firstPoint, cycle[cycle.size()-1])){ cycle.push_back(firstPoint); } if(firstCycle || !faceRect.Contains(currentBB)){ if(!getDir(cycle)){ reverseCycle(cycle); } firstCycle=false; faceRect = currentBB; } else { if(getDir(cycle) ){ reverseCycle(cycle); } } cycles.push_back(cycle); } } Region* res = buildRegion(cycles); correct = res!=0; return SetWord(res); } Word InRegion_old( const ListExpr typeInfo, const ListExpr instance, const int errorPos, ListExpr& errorInfo, bool& correct ) { Region* cr = new Region( 0 ); if(listutils::isSymbolUndefined(instance)){ cr->SetDefined(0); correct=true; return SetWord(Address(cr)); } cr->StartBulkLoad(); ListExpr RegionNL = instance; ListExpr FaceNL, CycleNL; int fcno=-1; int ccno=-1; int edno=-1; int partnerno = 0; if (!nl->IsAtom(instance)) { while( !nl->IsEmpty( RegionNL ) ) { FaceNL = nl->First( RegionNL ); RegionNL = nl->Rest( RegionNL); bool isCycle = true; //A face is composed by 1 cycle, and can have holes. //All the holes must be inside the face. (TO BE IMPLEMENTED0) //Region *faceCycle; fcno++; ccno=-1; edno=-1; if (nl->IsAtom( FaceNL )) { correct=false; return SetWord( Address(0) ); } while (!nl->IsEmpty( FaceNL) ) { CycleNL = nl->First( FaceNL ); FaceNL = nl->Rest( FaceNL ); ccno++; edno=-1; if (nl->IsAtom( CycleNL )) { correct=false; return SetWord( Address(0) ); } if (nl->ListLength( CycleNL) <3) { cerr << __PRETTY_FUNCTION__ << ": A cycle must have at least 3 edges!" << endl; correct=false; return SetWord( Address(0) ); } else { ListExpr firstPoint = nl->First( CycleNL ); ListExpr prevPoint = nl->First( CycleNL ); ListExpr flagedSeg = nl->TheEmptyList(); ListExpr currPoint = nl->TheEmptyList(); CycleNL = nl->Rest( CycleNL ); //Starting to compute a new cycle Points *cyclepoints= new Points( 8 ); // in memory Point *currvertex,p1(true),p2(true),firstP(true); //This function has the goal to store the half segments of //the cycle that is been treated. When the cycle's computation //is terminated the region rDir will be used to compute the //insideAbove //attribute of the half segments of this cycle. Region *rDir = new Region(32); rDir->StartBulkLoad(); currvertex = (Point*) InPoint ( nl->TheEmptyList(), firstPoint, 0, errorInfo, correct ).addr; if (!correct) { // todo: delete temp objects return SetWord( Address(0) ); } cyclepoints->StartBulkLoad(); (*cyclepoints) += (*currvertex); p1 = *currvertex; firstP = p1; cyclepoints->EndBulkLoad(); delete currvertex; while ( !nl->IsEmpty( CycleNL) ) { // cout<<"cycle "<First( CycleNL ); CycleNL = nl->Rest( CycleNL ); currvertex = (Point*) InPoint( nl->TheEmptyList(), currPoint, 0, errorInfo, correct ).addr; // cout<<"curvertex "<<*currvertex<Contains(*currvertex)) { cerr<< __PRETTY_FUNCTION__ << ": The same vertex: " <<(*currvertex) <<" appears repeatedly within the current cycle!"<StartBulkLoad(); (*cyclepoints) += (*currvertex); cyclepoints->EndBulkLoad(true,false,false); } delete currvertex; flagedSeg = nl->TwoElemList (nl-> BoolAtom(true), nl->TwoElemList(prevPoint, currPoint)); prevPoint=currPoint; edno++; //Create left dominating half segment HalfSegment * hs = (HalfSegment*)InHalfSegment ( nl->TheEmptyList(), flagedSeg, 0, errorInfo, correct ).addr; if(!correct){ if(hs){ cerr << __PRETTY_FUNCTION__ << ": Creation of left dominating " << "half segment (1) failed!" << endl; delete hs; } cr->DeleteIfAllowed(); return SetWord( Address(0) ); } hs->attr.faceno=fcno; hs->attr.cycleno=ccno; hs->attr.edgeno=edno; hs->attr.partnerno=partnerno; partnerno++; hs->attr.insideAbove = (hs->GetLeftPoint() == p1); //true (L-->R ),false (R--L) p1 = p2; if (( correct )&&( cr->InsertOk(*hs) )) { (*cr) += (*hs); // cout<<"cr+1 "<<*hs<IsLeftDomPoint() ) { (*rDir) += (*hs); // cout<<"rDr+1 "<<*hs<SetLeftDomPoint( false ); } else { hs->SetLeftDomPoint( true ); // cout<<"rDr+2 "<<*hs<TwoElemList (nl-> BoolAtom(true), nl->TwoElemList(firstPoint, currPoint)); HalfSegment * hs = (HalfSegment*)InHalfSegment ( nl->TheEmptyList(), flagedSeg, 0, errorInfo, correct ).addr; if(!correct){ if(hs){ cerr << __PRETTY_FUNCTION__ << ": Creation of " << "half segment (2) failed!" << endl; delete hs; } cr->DeleteIfAllowed(); return SetWord( Address(0) ); } hs->attr.faceno=fcno; hs->attr.cycleno=ccno; hs->attr.edgeno=edno; hs->attr.partnerno=partnerno; hs->attr.insideAbove = (hs->GetRightPoint() == firstP); //true (L-->R ),false (R--L), //the order of typing is last point than first point. partnerno++; //The last half segment of the region if (( correct )&&( cr->InsertOk(*hs) )) { (*cr) += (*hs); // cout<<"cr+3 "<<*hs<IsLeftDomPoint() ) { (*rDir) += (*hs); // cout<<"rDr+3 "<<*hs<SetLeftDomPoint( false ); } else { hs->SetLeftDomPoint( true ); // cout<<"rDr+4 "<<*hs<EndBulkLoad(true, false, false, false); //To calculate the inside above attribute bool direction = rDir->GetCycleDirection(); int h = cr->Size() - ( rDir->Size() * 2 ); while ( h < cr->Size()) { //after each left half segment of the region is its //correspondig right half segment HalfSegment hsIA; bool insideAbove; cr->Get(h,hsIA); /* The test for adjusting the inside above can be described as above, but was implemented in a different way that produces the same result. if ( (direction && hsIA->attr.insideAbove) || (!direction && !hsIA->attr.insideAbove) ) { //clockwise and l-->r or //counterclockwise and r-->l hsIA->attr.insideAbove=false; } else //clockwise and r-->r or //counterclockwise and l-->r true; */ if (direction == hsIA.attr.insideAbove) insideAbove = false; else insideAbove = true; if (!isCycle) insideAbove = !insideAbove; HalfSegment auxhsIA( hsIA ); auxhsIA.attr.insideAbove = insideAbove; cr->UpdateAttr(h,auxhsIA.attr); //Get right half segment cr->Get(h+1,hsIA); auxhsIA = hsIA; auxhsIA.attr.insideAbove = insideAbove; cr->UpdateAttr(h+1,auxhsIA.attr); h+=2; } //After the first face's cycle read the faceCycle variable is set. //Afterwards //it is tested if all the new cycles are inside the faceCycle. /* if (isCycle) faceCycle = new Region(rDir,false); else //To implement the test */ rDir->DeleteIfAllowed(); //After the end of the first cycle of the face, //all the following cycles are //holes, then isCycle is set to false. isCycle = false; } else { correct=false; return SetWord( Address(0) ); } } } } cr->SetNoComponents( fcno+1 ); cr->EndBulkLoad( true, true, true, false ); correct = true; return SetWord( cr ); } else { correct=false; return SetWord( Address(0) ); } } /* 8.5 ~Create~-function */ Word CreateRegion( const ListExpr typeInfo ) { //cout << "CreateRegion" << endl; return (SetWord( new Region( 0 ) )); } /* 8.6 ~Delete~-function */ void DeleteRegion( const ListExpr typeInfo, Word& w ) { //cout << "DeleteRegion" << endl; Region *cr = (Region *)w.addr; cr->Destroy(); cr->DeleteIfAllowed(false); w.addr = 0; } /* 8.7 ~Close~-function */ void CloseRegion( const ListExpr typeInfo, Word& w ) { //cout << "CloseRegion" << endl; ((Region *)w.addr)->DeleteIfAllowed(); w.addr = 0; } /* 8.8 ~Clone~-function */ Word CloneRegion( const ListExpr typeInfo, const Word& w ) { //cout << "CloneRegion" << endl; Region *cr = new Region( *((Region *)w.addr) ); return SetWord( cr ); } /* 7.8 ~Open~-function */ bool OpenRegion( SmiRecord& valueRecord, size_t& offset, const ListExpr typeInfo, Word& value ) { Region *r = (Region*)Attribute::Open( valueRecord, offset, typeInfo ); value = SetWord( r ); return true; } /* 7.8 ~Save~-function */ bool SaveRegion( SmiRecord& valueRecord, size_t& offset, const ListExpr typeInfo, Word& value ) { Region *r = (Region *)value.addr; Attribute::Save( valueRecord, offset, typeInfo, r ); return true; } /* 8.9 ~SizeOf~-function */ int SizeOfRegion() { return sizeof(Region); } /* 8.11 Function describing the signature of the type constructor */ ListExpr RegionProperty() { ListExpr listreplist = nl->TextAtom(); nl->AppendText(listreplist,"(*) where face is" " (*); " " and are *"); ListExpr examplelist = nl->TextAtom(); nl->AppendText(examplelist,"(((3 0)(10 1)(3 1))((3.1 0.1)" "(3.1 0.9)(6 0.8)))"); ListExpr remarkslist = nl->TextAtom(); nl->AppendText(remarkslist,"all must be completely within " "."); return (nl->TwoElemList( nl->FiveElemList(nl->StringAtom("Signature"), nl->StringAtom("Example Type List"), nl->StringAtom("List Rep"), nl->StringAtom("Example List"), nl->StringAtom("Remarks")), nl->FiveElemList(nl->StringAtom("-> DATA"), nl->StringAtom(Region::BasicType()), listreplist, examplelist, remarkslist))); } /* 8.12 Kind checking function This function checks whether the type constructor is applied correctly. Since type constructor ~point~ does not have arguments, this is trivial. */ bool CheckRegion( ListExpr type, ListExpr& errorInfo ) { return (nl->IsEqual( type, Region::BasicType() )); } /* 8.14 Creation of the type constructor instance */ TypeConstructor region( Region::BasicType(), //name RegionProperty, //describing signature OutRegion, InRegion, //Out and In functions 0, 0, //SaveTo and RestoreFrom List functions CreateRegion, DeleteRegion, //object creation and deletion OpenRegion, SaveRegion, // object open and save CloseRegion, CloneRegion, //object close and clone Region::Cast, //cast function SizeOfRegion, //sizeof function CheckRegion ); //kind checking function /* 8.15 type constructor for label data type */ GenTC