diff --git a/simgear/scene/tgdb/CMakeLists.txt b/simgear/scene/tgdb/CMakeLists.txt index 71edf2f6..42ca6548 100644 --- a/simgear/scene/tgdb/CMakeLists.txt +++ b/simgear/scene/tgdb/CMakeLists.txt @@ -20,6 +20,7 @@ set(HEADERS SGVertexArrayBin.hxx ShaderGeometry.hxx TreeBin.hxx + VPBElevationSlice.hxx VPBTechnique.hxx apt_signs.hxx obj.hxx @@ -38,6 +39,7 @@ set(SOURCES SGVasiDrawable.cxx ShaderGeometry.cxx TreeBin.cxx + VPBElevationSlice.cxx VPBTechnique.cxx apt_signs.cxx obj.cxx diff --git a/simgear/scene/tgdb/VPBElevationSlice.cxx b/simgear/scene/tgdb/VPBElevationSlice.cxx new file mode 100644 index 00000000..aed57bd9 --- /dev/null +++ b/simgear/scene/tgdb/VPBElevationSlice.cxx @@ -0,0 +1,976 @@ +// VPBVPBElevationSlice.cxx -- VirtualPlanetBuilder Elevation Slice +// +// Copyright (C) 2021 Stuart Buchanan +// +// This program 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. +// +// This program 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 this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +using namespace osgSim; + +namespace simgear { + +namespace VPBElevationSliceUtils +{ + +struct DistanceHeightXYZ +{ + + DistanceHeightXYZ(): + distance(0.0), + height(0.0) {} + + DistanceHeightXYZ(const DistanceHeightXYZ& dh): + distance(dh.distance), + height(dh.height), + position(dh.position) {} + + DistanceHeightXYZ(double d, double h, const osg::Vec3d& pos): + distance(d), + height(h), + position(pos) {} + + bool operator < (const DistanceHeightXYZ& rhs) const + { + // small distance values first + if (distance < rhs.distance) return true; + if (distance > rhs.distance) return false; + + // smallest heights first + return (height < rhs.height); + } + + bool operator == (const DistanceHeightXYZ& rhs) const + { + return distance==rhs.distance && height==rhs.height; + } + + bool operator != (const DistanceHeightXYZ& rhs) const + { + return distance!=rhs.distance || height!=rhs.height; + } + + bool equal_distance(const DistanceHeightXYZ& rhs, double epsilon=1e-6) const + { + return osg::absolute(rhs.distance - distance) <= epsilon; + } + + double distance; + double height; + osg::Vec3d position; +}; + +struct Point : public osg::Referenced, public DistanceHeightXYZ +{ + Point() {} + Point(double d, double h, const osg::Vec3d& pos): + DistanceHeightXYZ(d,h,pos) + { + } + + Point(const Point& point): + osg::Referenced(), + DistanceHeightXYZ(point) {} + + +}; + +struct Segment +{ + Segment(Point* p1, Point* p2) + { + if (*p1 < *p2) + { + _p1 = p1; + _p2 = p2; + } + else + { + _p1 = p2; + _p2 = p1; + } + } + + bool operator < ( const Segment& rhs) const + { + if (*_p1 < *rhs._p1) return true; + if (*rhs._p1 < *_p1) return false; + + return (*_p2 < *rhs._p2); + } + + enum Classification + { + UNCLASSIFIED, + IDENTICAL, + SEPARATE, + JOINED, + OVERLAPPING, + ENCLOSING, + ENCLOSED + }; + + Classification compare(const Segment& rhs) const + { + if (*_p1 == *rhs._p1 && *_p2==*rhs._p2) return IDENTICAL; + + const double epsilon = 1e-3; // 1mm + + double delta_distance = _p2->distance - rhs._p1->distance; + if (fabs(delta_distance) < epsilon) + { + if (fabs(_p2->height - rhs._p1->height) < epsilon) return JOINED; + } + + if (delta_distance==0.0) + { + return SEPARATE; + } + + if (rhs._p2->distance < _p1->distance || _p2->distance < rhs._p1->distance) return SEPARATE; + + bool rhs_p1_inside = (_p1->distance <= rhs._p1->distance) && (rhs._p1->distance <= _p2->distance); + bool rhs_p2_inside = (_p1->distance <= rhs._p2->distance) && (rhs._p2->distance <= _p2->distance); + + if (rhs_p1_inside && rhs_p2_inside) return ENCLOSING; + + bool p1_inside = (rhs._p1->distance <= _p1->distance) && (_p1->distance <= rhs._p2->distance); + bool p2_inside = (rhs._p1->distance <= _p2->distance) && (_p2->distance <= rhs._p2->distance); + + if (p1_inside && p2_inside) return ENCLOSED; + + if (rhs_p1_inside || rhs_p2_inside || p1_inside || p2_inside) return OVERLAPPING; + + return UNCLASSIFIED; + } + + double height(double d) const + { + double delta = (_p2->distance - _p1->distance); + return _p1->height + ((_p2->height - _p1->height) * (d - _p1->distance) / delta); + } + + double deltaHeight(Point& point) const + { + return point.height - height(point.distance); + } + + Point* createPoint(double d) const + { + if (d == _p1->distance) return _p1.get(); + if (d == _p2->distance) return _p2.get(); + + double delta = (_p2->distance - _p1->distance); + double r = (d - _p1->distance)/delta; + double one_minus_r = 1.0 - r; + return new Point(d, + _p1->height * one_minus_r + _p2->height * r, + _p1->position * one_minus_r + _p2->position * r); + } + + Point* createIntersectionPoint(const Segment& rhs) const + { + double A = _p1->distance; + double B = _p2->distance - _p1->distance; + double C = _p1->height; + double D = _p2->height - _p1->height; + + double E = rhs._p1->distance; + double F = rhs._p2->distance - rhs._p1->distance; + double G = rhs._p1->height; + double H = rhs._p2->height - rhs._p1->height; + + double div = D*F - B*H; + if (div==0.0) + { + return _p2.get(); + } + + double r = (G*F - E*H + A*H - C*F) / div; + + if (r<0.0) + { + return _p1.get(); + } + + if (r>1.0) + { + return _p2.get(); + } + + return new Point(A + B*r, C + D*r, _p1->position + (_p2->position - _p1->position)*r); + } + + + osg::ref_ptr _p1; + osg::ref_ptr _p2; + +}; + + +struct LineConstructor +{ + + typedef std::set SegmentSet; + + LineConstructor() {} + + + + void add(double d, double h, const osg::Vec3d& pos) + { + osg::ref_ptr newPoint = new Point(d,h,pos); + + + if (_previousPoint.valid() && newPoint->distance != _previousPoint->distance) + { + const double maxGradient = 100.0; + double gradient = fabs( (newPoint->height - _previousPoint->height) / (newPoint->distance - _previousPoint->distance) ); + + if (gradient < maxGradient) + { + _segments.insert( Segment(_previousPoint.get(), newPoint.get()) ); + } + } + + _previousPoint = newPoint; + } + + void endline() + { + _previousPoint = 0; + } + + void report() + { + } + + void pruneOverlappingSegments() + { + double epsilon = 0.001; + + for(SegmentSet::iterator itr = _segments.begin(); + itr != _segments.end(); + ) // itr increment is done manually at end of loop + { + SegmentSet::iterator nextItr = itr; + ++nextItr; + Segment::Classification classification = nextItr != _segments.end() ? itr->compare(*nextItr) : Segment::UNCLASSIFIED; + + while (classification>=Segment::OVERLAPPING) + { + + switch(classification) + { + case(Segment::OVERLAPPING): + { + // cases.... + // compute new end points for both segments + // need to work out which points are overlapping - lhs_p2 && rhs_p1 or lhs_p1 and rhs_p2 + // also need to check for cross cases. + + const Segment& lhs = *itr; + const Segment& rhs = *nextItr; + + bool rhs_p1_inside = (lhs._p1->distance <= rhs._p1->distance) && (rhs._p1->distance <= lhs._p2->distance); + bool lhs_p2_inside = (rhs._p1->distance <= lhs._p2->distance) && (lhs._p2->distance <= rhs._p2->distance); + + if (rhs_p1_inside && lhs_p2_inside) + { + double distance_between = osg::Vec2d(lhs._p2->distance - rhs._p1->distance, + lhs._p2->height - rhs._p1->height).length2(); + + if (distance_between < epsilon) + { + Segment newSeg(lhs._p2.get(), rhs._p2.get()); + _segments.insert(newSeg); + + _segments.erase(nextItr); + + nextItr = _segments.find(newSeg); + } + else + { + double dh1 = lhs.deltaHeight(*rhs._p1); + double dh2 = -rhs.deltaHeight(*lhs._p2); + + if (dh1 * dh2 < 0.0) + { + Point* cp = lhs.createIntersectionPoint(rhs); + + Segment seg1( lhs._p1.get(), lhs.createPoint(rhs._p2->distance) ); + Segment seg2( rhs._p1.get(), cp ); + Segment seg3( cp, lhs._p2.get() ); + Segment seg4( rhs.createPoint(lhs._p2->distance), lhs._p2.get() ); + + _segments.erase(nextItr); + _segments.erase(itr); + + _segments.insert(seg1); + _segments.insert(seg2); + _segments.insert(seg3); + _segments.insert(seg4); + + itr = _segments.find(seg1); + nextItr = itr; + ++nextItr; + + } + else if (dh1 <= 0.0 && dh2 <= 0.0) + { + Segment newSeg(rhs.createPoint(lhs._p2->distance), rhs._p2.get()); + + _segments.erase(nextItr); + _segments.insert(newSeg); + nextItr = itr; + ++nextItr; + + } + else if (dh1 >= 0.0 && dh2 >= 0.0) + { + Segment newSeg(lhs._p1.get(), lhs.createPoint(rhs._p1->distance)); + + _segments.erase(itr); + _segments.insert(newSeg); + itr = _segments.find(newSeg); + nextItr = itr; + ++nextItr; + } + else + { + ++nextItr; + } + } + } + else + { + ++nextItr; + } + + + break; + } + case(Segment::ENCLOSING): + { + // need to work out if rhs is below lhs or rhs is above lhs or crossing + + const Segment& enclosing = *itr; + const Segment& enclosed = *nextItr; + double dh1 = enclosing.deltaHeight(*enclosed._p1); + double dh2 = enclosing.deltaHeight(*enclosed._p2); + + if (dh1<=epsilon && dh2<=epsilon) + { + _segments.erase(nextItr); + + nextItr = itr; + ++nextItr; + + } + else if (dh1>=0.0 && dh2>=0.0) + { + + double d_left = enclosed._p1->distance - enclosing._p1->distance; + double d_right = enclosing._p2->distance - enclosed._p2->distance; + + if (d_left < epsilon && d_right < epsilon) + { + // treat ENCLOSED as ENCLOSING. + + nextItr = itr; + ++nextItr; + + _segments.erase(itr); + + itr = nextItr; + ++nextItr; + + } + else if (d_left < epsilon) + { + Segment newSeg(enclosing.createPoint(enclosed._p2->distance), enclosing._p2.get()); + nextItr = itr; + ++nextItr; + + _segments.erase(itr); + _segments.insert(newSeg); + + itr = nextItr; + ++nextItr; + } + else if (d_right < epsilon) + { + Segment newSeg(enclosing._p1.get(), enclosing.createPoint(enclosed._p1->distance) ); + + _segments.insert(newSeg); + _segments.erase(itr); + + itr = _segments.find(newSeg); + nextItr = itr; + ++nextItr; + } + else + { + Segment newSegLeft(enclosing._p1.get(), enclosing.createPoint(enclosed._p1->distance) ); + Segment newSegRight(enclosing.createPoint(enclosed._p2->distance), enclosing._p2.get()); + + _segments.erase(itr); + _segments.insert(newSegLeft); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + + } + else if (dh1 * dh2 < 0.0) + { + double d_left = enclosed._p1->distance - enclosing._p1->distance; + double d_right = enclosing._p2->distance - enclosed._p2->distance; + + if (d_left < epsilon && d_right < epsilon) + { + // treat ENCLOSED as ENCLOSING. + if (dh1 < 0.0) + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosing._p1.get(), cp); + Segment newSegRight(cp, enclosed._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + else + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosed._p1.get(), cp); + Segment newSegRight(cp, enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + + } + else if (d_left < epsilon) + { + if (dh1 < 0.0) + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosing._p1.get(), cp); + Segment newSegMid(cp, enclosed._p2.get()); + Segment newSegRight(enclosing.createPoint(enclosed._p2->distance), enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegMid); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + else + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosed._p1.get(), cp); + Segment newSegRight(cp, enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + } + else if (d_right < epsilon) + { + if (dh1 < 0.0) + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosing._p1.get(), cp); + Segment newSegRight(cp, enclosed._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + else + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosing._p1.get(), enclosing.createPoint(enclosed._p1->distance)); + Segment newSegMid(enclosed._p1.get(), cp); + Segment newSegRight(cp, enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegMid); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + } + else + { + if (dh1 < 0.0) + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosing._p1.get(), cp); + Segment newSegMid(cp, enclosed._p2.get()); + Segment newSegRight(enclosing.createPoint(enclosed._p2->distance), enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegMid); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + else + { + Point* cp = enclosing.createIntersectionPoint(enclosed); + Segment newSegLeft(enclosing._p1.get(), enclosing.createPoint(enclosed._p1->distance)); + Segment newSegMid(enclosed._p1.get(), cp); + Segment newSegRight(cp, enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(newSegLeft); + _segments.insert(newSegMid); + _segments.insert(newSegRight); + + itr = _segments.find(newSegLeft); + nextItr = itr; + ++nextItr; + } + + } + } + else + { + ++nextItr; + } + + break; + } + case(Segment::ENCLOSED): + { + // need to work out if lhs is below rhs or lhs is above rhs or crossing + const Segment& enclosing = *nextItr; + const Segment& enclosed = *itr; + double dh1 = enclosing.deltaHeight(*enclosed._p1); + double dh2 = enclosing.deltaHeight(*enclosed._p2); + + double d_left = enclosed._p1->distance - enclosing._p1->distance; + double d_right = enclosing._p2->distance - enclosed._p2->distance; + + if (d_left<=epsilon) + { + + if (dh1<=epsilon && dh2<=epsilon) + { + _segments.erase(itr); + + itr = nextItr; + ++nextItr; + } + else if (dh1>=0.0 && dh2>=0.0) + { + _segments.insert(Segment(enclosing.createPoint(enclosed._p2->distance), enclosed._p2.get())); + _segments.erase(nextItr); + + nextItr = itr; + ++nextItr; + } + else if (dh1 * dh2 < 0.0) + { + if (d_right<=epsilon) + { + // enclosed and enclosing effectively overlap + if (dh1 < 0.0) + { + Point* cp = enclosed.createIntersectionPoint(enclosing); + Segment segLeft(enclosed._p1.get(), cp); + Segment segRight(cp, enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(segLeft); + _segments.insert(segRight); + + itr = _segments.find(segLeft); + nextItr = itr; + ++nextItr; + } + else + { + Point* cp = enclosed.createIntersectionPoint(enclosing); + Segment segLeft(enclosing._p1.get(), cp); + Segment segRight(cp, enclosed._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(segLeft); + _segments.insert(segRight); + + itr = _segments.find(segLeft); + nextItr = itr; + ++nextItr; + } + } + else + { + // right hand side needs to be created + if (dh1 < 0.0) + { + Point* cp = enclosed.createIntersectionPoint(enclosing); + Segment segLeft(enclosed._p1.get(), cp); + Segment segRight(cp, enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(segLeft); + _segments.insert(segRight); + + itr = _segments.find(segLeft); + nextItr = itr; + ++nextItr; + } + else + { + Point* cp = enclosed.createIntersectionPoint(enclosing); + Segment segLeft(enclosing._p1.get(), cp); + Segment segMid(cp, enclosed._p2.get()); + Segment segRight(enclosing.createPoint(enclosed._p2->distance), enclosing._p2.get()); + + _segments.erase(itr); + _segments.erase(nextItr); + + _segments.insert(segLeft); + _segments.insert(segMid); + _segments.insert(segRight); + + itr = _segments.find(segLeft); + nextItr = itr; + ++nextItr; + } + } + } + else + { + ++nextItr; + } + } + else + { + SG_LOG(SG_TERRAIN, SG_ALERT, "*** ENCLOSED: is not coincendet with left handside of ENCLOSING, case not handled, advancing."); + } + + break; + } + default: + SG_LOG(SG_TERRAIN, SG_ALERT, "** Not handled, advancing"); + ++nextItr; + break; + } + + classification = ((itr != _segments.end()) && (nextItr != _segments.end())) ? itr->compare(*nextItr) : Segment::UNCLASSIFIED; + } + + // increment iterator it it's not already at end of container. + if (itr!=_segments.end()) ++itr; + } + } + + unsigned int numOverlapping(SegmentSet::const_iterator startItr) const + { + if (startItr==_segments.end()) return 0; + + SegmentSet::const_iterator nextItr = startItr; + ++nextItr; + + unsigned int num = 0; + while (nextItr!=_segments.end() && startItr->compare(*nextItr)>=Segment::OVERLAPPING) + { + ++num; + ++nextItr; + } + return num; + } + + unsigned int totalNumOverlapping() const + { + unsigned int total = 0; + for(SegmentSet::const_iterator itr = _segments.begin(); + itr != _segments.end(); + ++itr) + { + total += numOverlapping(itr); + } + return total; + } + + void copyPoints(VPBElevationSlice::Vec3dList& intersections, VPBElevationSlice::DistanceHeightList& distanceHeightIntersections) + { + if (_segments.empty()) return; + + SegmentSet::iterator prevItr = _segments.begin(); + SegmentSet::iterator nextItr = prevItr; + ++nextItr; + + intersections.push_back( prevItr->_p1->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(prevItr->_p1->distance, prevItr->_p1->height) ); + + intersections.push_back( prevItr->_p2->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(prevItr->_p2->distance, prevItr->_p2->height) ); + + for(; + nextItr != _segments.end(); + ++nextItr,++prevItr) + { + Segment::Classification classification = prevItr->compare(*nextItr); + switch(classification) + { + case(Segment::SEPARATE): + { + intersections.push_back( nextItr->_p1->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p1->distance, nextItr->_p1->height) ); + + intersections.push_back( nextItr->_p2->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p2->distance, nextItr->_p2->height) ); + break; + } + case(Segment::JOINED): + { +#if 1 + intersections.push_back( nextItr->_p2->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p2->distance, nextItr->_p2->height) ); +#else + intersections.push_back( nextItr->_p1->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p1->distance, nextItr->_p1->height) ); + + intersections.push_back( nextItr->_p2->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p2->distance, nextItr->_p2->height) ); +#endif + break; + } + default: + { + intersections.push_back( nextItr->_p1->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p1->distance, nextItr->_p1->height) ); + + intersections.push_back( nextItr->_p2->position ); + distanceHeightIntersections.push_back( VPBElevationSlice::DistanceHeight(nextItr->_p2->distance, nextItr->_p2->height) ); + break; + } + } + + } + + } + + SegmentSet _segments; + osg::ref_ptr _previousPoint; + osg::Plane _plane; + osg::ref_ptr _em; + +}; + +} + +VPBElevationSlice::VPBElevationSlice() +{ + /** Remove the osgSim::DatabaseCacheReadCallback that enables automatic loading + * of external PagedLOD tiles to ensure that the highest level of detail is used in intersections. + * This automatic loading of tiles is done by the intersection traversal that is done within + * the computeIntersections(..) method, so can result in long intersection times when external + * tiles have to be loaded. + * The external loading of tiles can be disabled by removing the read callback, this is done by + * calling the setDatabaseCacheReadCallback(DatabaseCacheReadCallback*) method with a value of 0.*/ + _intersectionVisitor.setReadCallback(0); +} + +void VPBElevationSlice::computeIntersections(osg::Node* scene, osg::Node::NodeMask traversalMask) +{ + osg::Plane plane; + osg::Polytope boundingPolytope; + + + // set up the main intersection plane + osg::Vec3d planeNormal = (_endPoint - _startPoint) ^ _upVector; + planeNormal.normalize(); + plane.set( planeNormal, _startPoint ); + + // set up the start cut off plane + osg::Vec3d startPlaneNormal = _upVector ^ planeNormal; + startPlaneNormal.normalize(); + boundingPolytope.add( osg::Plane(startPlaneNormal, _startPoint) ); + + // set up the end cut off plane + osg::Vec3d endPlaneNormal = planeNormal ^ _upVector; + endPlaneNormal.normalize(); + boundingPolytope.add( osg::Plane(endPlaneNormal, _endPoint) ); + + osg::ref_ptr intersector = new osgUtil::PlaneIntersector(plane, boundingPolytope); + + + intersector->setRecordHeightsAsAttributes(true); + + _intersectionVisitor.reset(); + _intersectionVisitor.setTraversalMask(traversalMask); + _intersectionVisitor.setIntersector( intersector.get() ); + + scene->accept(_intersectionVisitor); + + osgUtil::PlaneIntersector::Intersections& intersections = intersector->getIntersections(); + + typedef osgUtil::PlaneIntersector::Intersection::Polyline Polyline; + typedef osgUtil::PlaneIntersector::Intersection::Attributes Attributes; + + if (!intersections.empty()) + { + + osgUtil::PlaneIntersector::Intersections::iterator itr; + for(itr = intersections.begin(); + itr != intersections.end(); + ++itr) + { + osgUtil::PlaneIntersector::Intersection& intersection = *itr; + + if (intersection.matrix.valid()) + { + // transform points on polyline + for(Polyline::iterator pitr = intersection.polyline.begin(); + pitr != intersection.polyline.end(); + ++pitr) + { + *pitr = (*pitr) * (*intersection.matrix); + } + + // matrix no longer needed. + intersection.matrix = 0; + } + } + + VPBElevationSliceUtils::LineConstructor constructor; + constructor._plane = plane; + + // convert into distance/height + for(itr = intersections.begin(); + itr != intersections.end(); + ++itr) + { + osgUtil::PlaneIntersector::Intersection& intersection = *itr; + for(Polyline::iterator pitr = intersection.polyline.begin(); + pitr != intersection.polyline.end(); + ++pitr) + { + const osg::Vec3d& v = *pitr; + osg::Vec2d delta_xy( v.x() - _startPoint.x(), v.y() - _startPoint.y()); + double distance = delta_xy.length(); + + constructor.add( distance, v.z(), v); + } + constructor.endline(); + } + + // copy final results + _intersections.clear(); + _distanceHeightIntersections.clear(); + + // constructor.report(); + + unsigned int numOverlapping = constructor.totalNumOverlapping(); + + while(numOverlapping>0) + { + unsigned int previousNumOverlapping = numOverlapping; + + constructor.pruneOverlappingSegments(); + // constructor.report(); + + numOverlapping = constructor.totalNumOverlapping(); + if (previousNumOverlapping == numOverlapping) break; + } + + constructor.copyPoints(_intersections, _distanceHeightIntersections); + } + else + { + SG_LOG(SG_TERRAIN, SG_DEBUG, "No intersections found."); + } + +} + +VPBElevationSlice::Vec3dList VPBElevationSlice::computeVPBElevationSlice(osg::Node* scene, const osg::Vec3d& startPoint, const osg::Vec3d& endPoint, const osg::Vec3d& upVector, osg::Node::NodeMask traversalMask) +{ + VPBElevationSlice es; + es.setStartPoint(startPoint); + es.setEndPoint(endPoint); + es.setUpVector(upVector); + es.computeIntersections(scene, traversalMask); + return es.getIntersections(); +} + +} \ No newline at end of file diff --git a/simgear/scene/tgdb/VPBElevationSlice.hxx b/simgear/scene/tgdb/VPBElevationSlice.hxx new file mode 100644 index 00000000..0f2bf976 --- /dev/null +++ b/simgear/scene/tgdb/VPBElevationSlice.hxx @@ -0,0 +1,88 @@ +// VPBVPBElevationSlice.hxx -- VirtualPlanetBuilder Elevation Slice +// +// Copyright (C) 2021 Stuart Buchanan +// +// This program 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. +// +// This program 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 this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + +#ifndef VPBELEVATIONSLICE +#define VPBELEVATIONSLICE 1 + +#include + +// include so we can get access to the DatabaseCacheReadCallback +#include + +namespace simgear { + +class VPBElevationSlice +{ + public : + + VPBElevationSlice(); + + /** Set the start point of the slice.*/ + void setStartPoint(const osg::Vec3d& startPoint) { _startPoint = startPoint; } + + /** Get the start point of the slice.*/ + const osg::Vec3d& getStartPoint() const { return _startPoint; } + + /** Set the end point of the slice.*/ + void setEndPoint(const osg::Vec3d& endPoint) { _endPoint = endPoint; } + + /** Get the end point of the slice.*/ + const osg::Vec3d& getEndPoint() const { return _endPoint; } + + /** Set the up Vector of the slice.*/ + void setUpVector(const osg::Vec3d& upVector) { _upVector = upVector; } + + /** Get the start point of the slice.*/ + const osg::Vec3d& getUpVector() const { return _upVector; } + + + typedef std::vector Vec3dList; + + /** Get the intersections in the form of a vector of Vec3d. */ + const Vec3dList& getIntersections() const { return _intersections; } + + typedef std::pair DistanceHeight; + typedef std::vector DistanceHeightList; + + /** Get the intersections in the form a vector of pair representing distance along the slice and height. */ + const DistanceHeightList& getDistanceHeightIntersections() const { return _distanceHeightIntersections; } + + + /** Compute the intersections with the specified scene graph, the results are stored in vectors of Vec3d. + * Note, if the topmost node is a CoordinateSystemNode then the input points are assumed to be geocentric, + * with the up vector defined by the EllipsoidModel attached to the CoordinateSystemNode. + * If the topmost node is not a CoordinateSystemNode then a local coordinates frame is assumed, with a local up vector. */ + void computeIntersections(osg::Node* scene, osg::Node::NodeMask traversalMask=0xffffffff); + + /** Compute the vertical distance between the specified scene graph and a single HAT point.*/ + static Vec3dList computeVPBElevationSlice(osg::Node* scene, const osg::Vec3d& startPoint, const osg::Vec3d& endPoint, const osg::Vec3d& upVector, osg::Node::NodeMask traversalMask=0xffffffff); + + protected : + osg::Vec3d _startPoint; + osg::Vec3d _endPoint; + osg::Vec3d _upVector; + Vec3dList _intersections; + DistanceHeightList _distanceHeightIntersections; + + osgUtil::IntersectionVisitor _intersectionVisitor; +}; + +} + +#endif diff --git a/simgear/scene/tgdb/VPBTechnique.cxx b/simgear/scene/tgdb/VPBTechnique.cxx index d8078261..e82d7789 100644 --- a/simgear/scene/tgdb/VPBTechnique.cxx +++ b/simgear/scene/tgdb/VPBTechnique.cxx @@ -18,7 +18,6 @@ #include -#include #include #include @@ -44,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -134,8 +134,6 @@ void VPBTechnique::setFilterMatrixAs(FilterType filterType) void VPBTechnique::init(int dirtyMask, bool assumeMultiThreaded) { - OSG_INFO<<"Doing VPBTechnique::init()"< lock(_writeBufferMutex); @@ -144,26 +142,15 @@ void VPBTechnique::init(int dirtyMask, bool assumeMultiThreaded) if (dirtyMask==0) return; + osgTerrain::TileID tileID = tile->getTileID(); + SG_LOG(SG_TERRAIN, SG_DEBUG, "Init of tile " << tileID.x << "," << tileID.y << " level " << tileID.level << " " << dirtyMask); + osg::ref_ptr buffer = new BufferData; Locator* masterLocator = computeMasterLocator(); osg::Vec3d centerModel = computeCenterModel(*buffer, masterLocator); - // We use TerrainTile::IMAGERY_DIRTY to (re-)generate features from a .STG file, such as line features - // We use TerrainTile::ELEVATION_DIRTY to regenerate the geometry - if (dirtyMask & TerrainTile::ELEVATION_DIRTY) { - generateGeometry(*buffer, masterLocator, centerModel); - applyColorLayers(*buffer, masterLocator); - applyTrees(*buffer, masterLocator); - } - - if (dirtyMask & TerrainTile::IMAGERY_DIRTY) { - applyLineFeatures(*buffer, masterLocator); - } - - -/* if ((dirtyMask & TerrainTile::IMAGERY_DIRTY)==0) { generateGeometry(*buffer, masterLocator, centerModel); @@ -189,7 +176,7 @@ void VPBTechnique::init(int dirtyMask, bool assumeMultiThreaded) applyTrees(*buffer, masterLocator); applyLineFeatures(*buffer, masterLocator); } -*/ + if (buffer->_transform.valid()) buffer->_transform->setThreadSafeRefUnref(true); if (!_currentBufferData || !assumeMultiThreaded) @@ -837,8 +824,6 @@ void VPBTechnique::generateGeometry(BufferData& buffer, Locator* masterLocator, if (mat) { effectProp = new SGPropertyNode(); makeChild(effectProp, "inherits-from")->setStringValue(mat->get_effect_name()); - //effectPropeffect->getChild("texture", 0, true)->getChild("image", 0, true)->setStringValue("Terrain/forest1a.png"); - //SG_LOG(SG_TERRAIN, SG_ALERT, "Created VPB effect property " << mat->get_effect_name()); } else { SG_LOG( SG_TERRAIN, SG_ALERT, "Unable to get effect for VPB - no matching material in library"); effectProp = new SGPropertyNode; @@ -1763,39 +1748,39 @@ void VPBTechnique::applyTrees(BufferData& buffer, Locator* masterLocator) void VPBTechnique::applyLineFeatures(BufferData& buffer, Locator* masterLocator) { - int line_features_lod_range = 6; + unsigned int line_features_lod_range = 6; + float minWidth = 9999.9; - SGPropertyNode* propertyNode = _options->getPropertyNode().get(); + const unsigned int tileLevel = _terrainTile->getTileID().level; + const SGPropertyNode* propertyNode = _options->getPropertyNode().get(); if (propertyNode) { - line_features_lod_range = propertyNode->getIntValue("/sim/rendering/static-lod/line-features-lod-range", line_features_lod_range); + const SGPropertyNode* static_lod = propertyNode->getNode("/sim/rendering/static-lod"); + line_features_lod_range = static_lod->getIntValue("line-features-lod-level", line_features_lod_range); + if (static_lod->getChildren("lod-level").size() > tileLevel) { + minWidth = static_lod->getChildren("lod-level")[tileLevel]->getFloatValue("line-features-min-width", minWidth); + } } - // Do not generate vegetation for tiles too far away - if (_terrainTile->getTileID().level < line_features_lod_range) { + if (tileLevel < line_features_lod_range) { + // Do not generate vegetation for tiles too far away return; } - SGMaterialLibPtr matlib = _options->getMaterialLib(); + SG_LOG(SG_TERRAIN, SG_DEBUG, "Generating roads of width > " << minWidth << " for tile LoD level " << tileLevel); + + const SGMaterialLibPtr matlib = _options->getMaterialLib(); SGMaterial* mat = 0; - osg::Vec3d world; - SGGeod loc; + const osg::Vec3d world = buffer._transform->getMatrix().getTrans(); if (! matlib) { SG_LOG(SG_TERRAIN, SG_ALERT, "Unable to get materials library to generate roads"); return; } - // Determine the center of the tile, sadly non-trivial. - osg::Vec3d tileloc = computeCenter(buffer, masterLocator); - masterLocator->convertLocalToModel(tileloc, world); - const SGVec3d world2 = SGVec3d(world.x(), world.y(), world.z()); - loc = SGGeod::fromCart(world2); - - const osg::Matrixd R_vert = makeZUpFrameRelative(loc); - // Get all appropriate roads. We assume that the VPB terrain tile is smaller than a Bucket size. - SGBucket bucket = SGBucket(loc); + const SGGeod loc = SGGeod::fromCart(toSG(world)); + const SGBucket bucket = SGBucket(loc); auto roads = std::find_if(_lineFeatureLists.begin(), _lineFeatureLists.end(), [bucket](BucketLineFeatureBinList b){return (b.first == bucket);}); if (roads == _lineFeatureLists.end()) return; @@ -1803,7 +1788,7 @@ void VPBTechnique::applyLineFeatures(BufferData& buffer, Locator* masterLocator) SGMaterialCache* matcache = _options->getMaterialLib()->generateMatCache(loc, _options); for (; roads != _lineFeatureLists.end(); ++roads) { - LineFeatureBinList roadBins = roads->second; + const LineFeatureBinList roadBins = roads->second; for (auto rb = roadBins.begin(); rb != roadBins.end(); ++rb) { @@ -1826,7 +1811,7 @@ void VPBTechnique::applyLineFeatures(BufferData& buffer, Locator* masterLocator) auto lineFeatures = rb->getLineFeatures(); for (auto r = lineFeatures.begin(); r != lineFeatures.end(); ++r) { - generateLineFeature(buffer, masterLocator, *r, world, R_vert, v, t, n, xsize, ysize); + if (r->_width > minWidth) generateLineFeature(buffer, masterLocator, *r, world, v, t, n, xsize, ysize); } if (v->size() == 0) continue; @@ -1854,7 +1839,7 @@ void VPBTechnique::applyLineFeatures(BufferData& buffer, Locator* masterLocator) } } -void VPBTechnique::generateLineFeature(BufferData& buffer, Locator* masterLocator, LineFeatureBin::LineFeature road, osg::Vec3d modelCenter, const osg::Matrixd R_vert, osg::Vec3Array* v, osg::Vec2Array* t, osg::Vec3Array* n, unsigned int xsize, unsigned int ysize) +void VPBTechnique::generateLineFeature(BufferData& buffer, Locator* masterLocator, LineFeatureBin::LineFeature road, osg::Vec3d modelCenter, osg::Vec3Array* v, osg::Vec2Array* t, osg::Vec3Array* n, unsigned int xsize, unsigned int ysize) { if (road._nodes.size() < 2) { SG_LOG(SG_TERRAIN, SG_ALERT, "Coding error - LineFeatureBin::LineFeature with fewer than two nodes"); @@ -1867,7 +1852,11 @@ void VPBTechnique::generateLineFeature(BufferData& buffer, Locator* masterLocato auto road_iter = road._nodes.begin(); ma = *road_iter - modelCenter; - osg::Vec3d intersect = getMeshIntersection(buffer, masterLocator, ma); + // We're in Earth-centered coordinates, so "up" is simply directly away from (0,0,0) + osg::Vec3d up = modelCenter; + up.normalize(); + + osg::Vec3d intersect = getMeshIntersection(buffer, masterLocator, ma, up); // If we've got an intersection, add it to the list if (intersect != ma) roadPoints.push_back(intersect); @@ -1877,28 +1866,22 @@ void VPBTechnique::generateLineFeature(BufferData& buffer, Locator* masterLocato for (; road_iter != road._nodes.end(); road_iter++) { mb = *road_iter - modelCenter; - intersect = getMeshIntersection(buffer, masterLocator, mb); + intersect = getMeshIntersection(buffer, masterLocator, mb, up); if (intersect != mb) roadPoints.push_back(intersect); - osgSim::ElevationSlice es; - es.setStartPoint(ma); - es.setEndPoint(mb); - es.computeIntersections(buffer._geometry); + auto esl = VPBElevationSlice::computeVPBElevationSlice(buffer._geometry, ma, mb, up); - auto esl = es.getIntersections(); - for(auto eslitr = esl.begin(); eslitr != esl.end(); ++eslitr) - { - roadPoints.push_back(*eslitr); + for(auto eslitr = esl.begin(); eslitr != esl.end(); ++eslitr) { + roadPoints.push_back(*eslitr); } // Now traverse the next segment - ma = mb; + ma = intersect; } // We now have a series of points following the topography of the elevation mesh. float xtex, ytex; osg::Vec3d local, origin, top, start, end; - osg::Vec3d up = osg::Matrixd::inverse(R_vert) * osg::Vec3d(0,0,1); auto iter = roadPoints.begin(); start = *iter + up; @@ -1951,17 +1934,10 @@ void VPBTechnique::generateLineFeature(BufferData& buffer, Locator* masterLocato } // Find the intersection of a given SGGeod with the terrain mesh -osg::Vec3d VPBTechnique::getMeshIntersection(BufferData& buffer, Locator* masterLocator, osg::Vec3d pt) +osg::Vec3d VPBTechnique::getMeshIntersection(BufferData& buffer, Locator* masterLocator, osg::Vec3d pt, osg::Vec3d up) { - osg::Vec3d local, origin, top; - - // Get a vertical line segment. - masterLocator->convertModelToLocal(pt, local); - masterLocator->convertLocalToModel(osg::Vec3d(local.x(), local.y(), -1000.0), origin); - masterLocator->convertLocalToModel(osg::Vec3d(local.x(), local.y(), 100000.0), top); - osg::ref_ptr intersector; - intersector = new osgUtil::LineSegmentIntersector(origin, top); + intersector = new osgUtil::LineSegmentIntersector(pt - up*100.0, pt + up*8000.0); osgUtil::IntersectionVisitor visitor(intersector.get()); buffer._geometry->accept(visitor); @@ -1969,7 +1945,7 @@ osg::Vec3d VPBTechnique::getMeshIntersection(BufferData& buffer, Locator* master // We have an intersection with the terrain model, so return it return intersector->getFirstIntersection().getWorldIntersectPoint(); } else { - // No intersection. Likely this point is outside our geometry. So just return the cartesian element. + // No intersection. Likely this point is outside our geometry. So just return the original element. return pt; } } @@ -2039,15 +2015,17 @@ void VPBTechnique::releaseGLObjects(osg::State* state) const } // Simple vistor to check for any underlying terrain meshes that intersect with a given constraint therefore may need to be modified -// (e.g elevation lowered to ensure the terrain doesn't poke through an airport mesh, orn roads generated) +// (e.g elevation lowered to ensure the terrain doesn't poke through an airport mesh, or line features generated) class TerrainVisitor : public osg::NodeVisitor { public: - osg::ref_ptr _constraint; // Object to flatten the terrain under. + osg::ref_ptr _constraint; // Object describing the volume to be modified. int _dirtyMask; // Dirty mask to apply. - TerrainVisitor( osg::ref_ptr node, int mask) : + int _minLevel; // Minimum LoD level to modify. + TerrainVisitor( osg::ref_ptr node, int mask, int minLevel) : osg::NodeVisitor(osg::NodeVisitor::TRAVERSE_ALL_CHILDREN), _constraint(node), - _dirtyMask(mask) + _dirtyMask(mask), + _minLevel(minLevel) { } virtual ~TerrainVisitor() { } @@ -2055,13 +2033,17 @@ class TerrainVisitor : public osg::NodeVisitor { void apply(osg::Node& node) { osgTerrain::TerrainTile* tile = dynamic_cast(&node); - if (tile) - { + if (tile) { // Determine if the constraint should affect this tile. + const int level = tile->getTileID().level; const osg::BoundingSphere tileBB = tile->getBound(); - if (tileBB.intersects(_constraint->getBound())) { + if ((level >= _minLevel) && tileBB.intersects(_constraint->getBound())) { // Dirty any existing terrain tiles containing this constraint, which will force regeneration + osgTerrain::TileID tileID = tile->getTileID(); + SG_LOG(SG_TERRAIN, SG_DEBUG, "Setting dirty mask for tile " << tileID.x << "," << tileID.y << " level " << tileID.level << " " << _dirtyMask); tile->setDirtyMask(_dirtyMask); + } else if (tileBB.intersects(_constraint->getBound())) { + traverse(node); } } else { traverse(node); @@ -2077,7 +2059,7 @@ void VPBTechnique::addElevationConstraint(osg::ref_ptr constraint, os const std::lock_guard lock(VPBTechnique::_constraint_mutex); // Lock the _constraintGroup for this scope _constraintGroup->addChild(constraint.get()); - TerrainVisitor ftv(constraint, TerrainTile::ALL_DIRTY); + TerrainVisitor ftv(constraint, TerrainTile::ALL_DIRTY, 0); terrain->accept(ftv); } @@ -2120,24 +2102,25 @@ void VPBTechnique::addLineFeatureList(SGBucket bucket, LineFeatureBinList roadLi // We need to trigger a re-build of the appropriate Terrain tile, so create a pretend node and run the TerrainVisitor to // "dirty" the TerrainTile that it intersects with. osg::ref_ptr n = new osg::Node(); - SGVec3d coord1, coord2; + + //SGVec3d coord1, coord2; + //SGGeodesy::SGGeodToCart(SGGeod::fromDegM(bucket.get_center_lon() -0.5*bucket.get_width(), bucket.get_center_lat() -0.5*bucket.get_height(), 0.0), coord1); + //SGGeodesy::SGGeodToCart(SGGeod::fromDegM(bucket.get_center_lon() +0.5*bucket.get_width(), bucket.get_center_lat() +0.5*bucket.get_height(), 0.0), coord2); + //osg::BoundingBox bbox = osg::BoundingBox(toOsg(coord1), toOsg(coord2)); + //n->setInitialBound(bbox); - SGGeodesy::SGGeodToCart(SGGeod::fromDegM(bucket.get_center_lon() -0.5*bucket.get_width(), bucket.get_center_lat() -0.5*bucket.get_height(), 0.0), coord1); - SGGeodesy::SGGeodToCart(SGGeod::fromDegM(bucket.get_center_lon() +0.5*bucket.get_width(), bucket.get_center_lat() +0.5*bucket.get_height(), 0.0), coord2); - //SGGeodesy::SGGeodToCart(SGGeod::fromDegM(bucket.get_center_lon(), bucket.get_center_lat(), 0.0), coord); - //n->setInitialBound(osg::BoundingSphere(toOsg(coord), max(bucket.get_width_m(), bucket.get_height_m()))); - - osg::BoundingBox bbox = osg::BoundingBox(toOsg(coord1), toOsg(coord2)); - n->setInitialBound(bbox); + SGVec3d coord; + SGGeodesy::SGGeodToCart(SGGeod::fromDegM(bucket.get_center_lon(), bucket.get_center_lat(), 0.0), coord); + n->setInitialBound(osg::BoundingSphere(toOsg(coord), max(bucket.get_width_m(), bucket.get_height_m()))); SG_LOG(SG_TERRAIN, SG_DEBUG, "Adding line features to " << bucket.gen_index_str()); - TerrainVisitor ftv(n, TerrainTile::IMAGERY_DIRTY); + TerrainVisitor ftv(n, TerrainTile::ALL_DIRTY, 0); terrainNode->accept(ftv); } void VPBTechnique::unloadLineFeatures(SGBucket bucket) { - SG_LOG(SG_TERRAIN, SG_ALERT, "Erasing all roads with entry " << bucket); + SG_LOG(SG_TERRAIN, SG_DEBUG, "Erasing all roads with entry " << bucket); const std::lock_guard lock(VPBTechnique::_lineFeatureLists_mutex); // Lock the _lineFeatureLists for this scope // C++ 20... //std::erase_if(_lineFeatureLists, [bucket](BucketLineFeatureBinList p) { return p.first == bucket; } ); diff --git a/simgear/scene/tgdb/VPBTechnique.hxx b/simgear/scene/tgdb/VPBTechnique.hxx index fd27aa70..6dd8aa81 100644 --- a/simgear/scene/tgdb/VPBTechnique.hxx +++ b/simgear/scene/tgdb/VPBTechnique.hxx @@ -132,13 +132,12 @@ class VPBTechnique : public TerrainTechnique virtual void generateLineFeature(BufferData& buffer, Locator* masterLocator, LineFeatureBin::LineFeature road, osg::Vec3d modelCenter, - const osg::Matrixd R_vert, osg::Vec3Array* v, osg::Vec2Array* t, osg::Vec3Array* n, unsigned int xsize, unsigned int ysize); - virtual osg::Vec3d getMeshIntersection(BufferData& buffer, Locator* masterLocator, osg::Vec3d pt); + virtual osg::Vec3d getMeshIntersection(BufferData& buffer, Locator* masterLocator, osg::Vec3d pt, osg::Vec3d up); OpenThreads::Mutex _writeBufferMutex; osg::ref_ptr _currentBufferData;