FGFS-4.1/OpenSceneGraph
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
132b355d4f
OpenSceneGraph/src/osgSim/SphereSegment.cpp
Robert Osfield 0e16b64665 Coverted tabs to space in core libraries.
2005-11-17 13:35:53 +00:00

3042 lines
103 KiB
C++
Raw Blame History

/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2005 Robert Osfield
*
* This library is open source and may be redistributed and/or modified under
* the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or
* (at your option) any later version. The full license is in LICENSE file
* included with this distribution, and on the openscenegraph.org website.
*
* This library 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
* OpenSceneGraph Public License for more details.
*/
#include <osgSim/SphereSegment>
#include <osg/Notify>
#include <osg/CullFace>
#include <osg/LineWidth>
#include <osg/Transform>
#include <osg/Geometry>
#include <osg/TriangleIndexFunctor>
#include <osg/ShapeDrawable>
#include <osg/io_utils>
#include <algorithm>
#include <list>
using namespace osgSim;
// Define the collection of nested classes, all Drawables, which make
// up the parts of the sphere segment.
/**
SphereSegment::Surface is the Drawable which represents the specified area of the
sphere's surface.
*/
class SphereSegment::Surface: public osg::Drawable
{
public:
Surface(SphereSegment* ss): osg::Drawable(), _ss(ss) {}
Surface():_ss(0)
{
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::Surface() default constructor"<<std::endl;
}
Surface(const Surface& rhs, const osg::CopyOp& co=osg::CopyOp::SHALLOW_COPY):osg::Drawable(rhs,co), _ss(0)
{
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::Surface() copy constructor"<<std::endl;
}
META_Object(osgSim,Surface)
void drawImplementation(osg::State& state) const;
virtual osg::BoundingBox computeBound() const;
protected:
private:
SphereSegment* _ss;
};
void SphereSegment::Surface::drawImplementation(osg::State& state) const
{
_ss->Surface_drawImplementation(state);
}
osg:: BoundingBox SphereSegment::Surface::computeBound() const
{
osg:: BoundingBox bbox;
_ss->Surface_computeBound(bbox);
return bbox;
}
/**
SphereSegment::EdgeLine is the Drawable which represents the line around the edge
of the specified area of the sphere's EdgeLine.
*/
class SphereSegment::EdgeLine: public osg::Drawable
{
public:
EdgeLine(SphereSegment* ss): osg::Drawable(), _ss(ss) { init(); }
EdgeLine():_ss(0)
{
init();
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::EdgeLine() default constructor"<<std::endl;
}
EdgeLine(const EdgeLine& rhs, const osg::CopyOp& co=osg::CopyOp::SHALLOW_COPY):osg::Drawable(rhs,co), _ss(0)
{
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::EdgeLine() copy constructor"<<std::endl;
}
META_Object(osgSim,EdgeLine)
void drawImplementation(osg::State& state) const;
protected:
void init()
{
// switch off lighting.
getOrCreateStateSet()->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
//getOrCreateStateSet()->setAttributeAndModes(new osg::LineWidth(2.0),osg::StateAttribute::OFF);
}
virtual osg::BoundingBox computeBound() const;
private:
SphereSegment* _ss;
};
void SphereSegment::EdgeLine::drawImplementation(osg::State& state) const
{
_ss->EdgeLine_drawImplementation(state);
}
osg::BoundingBox SphereSegment::EdgeLine::computeBound() const
{
osg::BoundingBox bbox;
_ss->EdgeLine_computeBound(bbox);
return bbox;
}
/**
SphereSegment::Side is a Drawable which represents one of the
planar areas, at either the minimum or maxium azimuth.
*/
class SphereSegment::Side: public osg::Drawable
{
public:
Side(SphereSegment* ss, SphereSegment::SideOrientation po, SphereSegment::BoundaryAngle pa):
osg::Drawable(), _ss(ss), _planeOrientation(po), _BoundaryAngle(pa) {}
Side():_ss(0)
{
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::Side() default constructor"<<std::endl;
}
Side(const Side& rhs, const osg::CopyOp& co=osg:: CopyOp::SHALLOW_COPY): osg::Drawable(rhs,co), _ss(0)
{
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::Side() copy constructor"<<std::endl;
}
META_Object(osgSim,Side)
void drawImplementation(osg::State& state) const;
protected:
virtual osg::BoundingBox computeBound() const;
private:
SphereSegment* _ss;
SphereSegment::SideOrientation _planeOrientation;
SphereSegment::BoundaryAngle _BoundaryAngle;
};
void SphereSegment::Side::drawImplementation(osg::State& state) const
{
_ss->Side_drawImplementation(state, _planeOrientation, _BoundaryAngle);
}
osg::BoundingBox SphereSegment::Side::computeBound() const
{
osg::BoundingBox bbox;
_ss->Side_computeBound(bbox, _planeOrientation, _BoundaryAngle);
return bbox;
}
/**
SphereSegment::Spoke is a Drawable which represents a spoke.
*/
class SphereSegment::Spoke: public osg::Drawable
{
public:
Spoke(SphereSegment* ss, SphereSegment::BoundaryAngle azAngle, SphereSegment::BoundaryAngle elevAngle):
osg::Drawable(), _ss(ss), _azAngle(azAngle), _elevAngle(elevAngle) { init(); }
Spoke():_ss(0)
{
init();
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::Spoke() default constructor"<<std::endl;
}
Spoke(const Spoke& rhs, const osg::CopyOp& co=osg:: CopyOp::SHALLOW_COPY): osg::Drawable(rhs,co), _ss(0)
{
osg::notify(osg::WARN)<<
"Warning: unexpected call to osgSim::SphereSegment::Spoke() copy constructor"<<std::endl;
}
META_Object(osgSim,Spoke)
void drawImplementation(osg::State& state) const;
protected:
void init()
{
// switch off lighting.
getOrCreateStateSet()->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
//getOrCreateStateSet()->setAttributeAndModes(new osg::LineWidth(2.0),osg::StateAttribute::OFF);
}
virtual osg::BoundingBox computeBound() const;
private:
SphereSegment* _ss;
SphereSegment::BoundaryAngle _azAngle, _elevAngle;
};
void SphereSegment::Spoke::drawImplementation(osg::State& state) const
{
_ss->Spoke_drawImplementation(state, _azAngle, _elevAngle);
}
osg::BoundingBox SphereSegment::Spoke::computeBound() const
{
osg::BoundingBox bbox;
_ss->Spoke_computeBound(bbox, _azAngle, _elevAngle);
return bbox;
}
SphereSegment::SphereSegment(const osg::Vec3& centre, float radius, const osg::Vec3& vec, float azRange,
float elevRange, int density):
osg::Geode(),
_centre(centre), _radius(radius),
_density(density),
_drawMask(DrawMask(ALL))
{
// Rather than store the vector, we'll work out the azimuth boundaries and elev
// boundaries now, rather than at draw time.
setArea(vec, azRange, elevRange);
init();
}
void SphereSegment::setCentre(const osg::Vec3& c)
{
_centre = c;
dirtyAllDrawableDisplayLists();
dirtyAllDrawableBounds();
dirtyBound();
}
const osg::Vec3& SphereSegment::getCentre() const
{
return _centre;
}
void SphereSegment::setRadius(float r)
{
_radius = r;
dirtyAllDrawableDisplayLists();
dirtyAllDrawableBounds();
dirtyBound();
}
float SphereSegment::getRadius() const
{
return _radius;
}
void SphereSegment::setArea(const osg::Vec3& v, float azRange, float elevRange)
{
osg::Vec3 vec(v);
vec.normalize(); // Make sure we're unit length
// Calculate the elevation range
float elev = asin(vec.z()); // Elevation angle
elevRange /= 2.0f;
_elevMin = elev - elevRange;
_elevMax = elev + elevRange;
// Calculate the azimuth range, cater for trig ambiguities
float xyLen = cos(elev);
float az;
if(vec.x() != 0.0f) az = asin(vec.x()/xyLen);
else az = acos(vec.y()/xyLen);
azRange /= 2.0f;
_azMin = az - azRange;
_azMax = az + azRange;
dirtyAllDrawableDisplayLists();
dirtyAllDrawableBounds();
dirtyBound();
}
void SphereSegment::getArea(osg::Vec3& vec, float& azRange, float& elevRange) const
{
azRange = _azMax - _azMin;
elevRange = _elevMax - _elevMin;
float az = azRange/2.0f;
float elev = elevRange/2.0f;
vec.set(cos(elev)*sin(az), cos(elev)*cos(az), sin(elev));
}
void SphereSegment::setArea(float azMin, float azMax,
float elevMin, float elevMax)
{
_azMin=azMin;
_azMax=azMax;
_elevMin=elevMin;
_elevMax=elevMax;
dirtyAllDrawableDisplayLists();
dirtyAllDrawableBounds();
dirtyBound();
}
void SphereSegment::getArea(float &azMin, float &azMax,
float &elevMin, float &elevMax) const
{
azMin=_azMin;
azMax=_azMax;
elevMin=_elevMin;
elevMax=_elevMax;
}
void SphereSegment::setDensity(int density)
{
_density = density;
dirtyAllDrawableDisplayLists();
}
int SphereSegment::getDensity() const
{
return _density;
}
void SphereSegment::init()
{
addDrawable(new Surface(this));
addDrawable(new EdgeLine(this));
addDrawable(new Side(this,AZIM,MIN));
addDrawable(new Side(this,AZIM,MAX));
addDrawable(new Side(this,ELEV,MIN));
addDrawable(new Side(this,ELEV,MAX));
addDrawable(new Spoke(this,MIN,MIN));
addDrawable(new Spoke(this,MIN,MAX));
addDrawable(new Spoke(this,MAX,MIN));
addDrawable(new Spoke(this,MAX,MAX));
}
void SphereSegment::dirtyAllDrawableDisplayLists()
{
for(DrawableList::iterator itr = _drawables.begin();
itr != _drawables.end();
++itr)
{
(*itr)->dirtyDisplayList();
}
}
void SphereSegment::dirtyAllDrawableBounds()
{
for(DrawableList::iterator itr = _drawables.begin();
itr != _drawables.end();
++itr)
{
(*itr)->dirtyBound();
}
}
void SphereSegment::Surface_drawImplementation(osg::State& /* state */) const
{
const float azIncr = (_azMax - _azMin)/_density;
const float elevIncr = (_elevMax - _elevMin)/_density;
// Draw the area on the sphere surface if needed
// ---------------------------------------------
if(_drawMask & SURFACE)
{
glColor4fv(_surfaceColor.ptr());
bool drawBackSide = true;
bool drawFrontSide = true;
// draw back side.
if (drawBackSide)
{
for(int i=0; i+1<=_density; i++)
{
// Because we're drawing quad strips, we need to work out
// two azimuth values, to form each edge of the (z-vertical)
// strips
float az1 = _azMin + (i*azIncr);
float az2 = _azMin + ((i+1)*azIncr);
glBegin(GL_QUAD_STRIP);
for (int j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
// QuadStrip Edge formed at az1
// ----------------------------
// Work out the sphere normal
float x = cos(elev)*sin(az1);
float y = cos(elev)*cos(az1);
float z = sin(elev);
glNormal3f(-x, -y, -z);
glVertex3f(_centre.x() + _radius*x,
_centre.y() + _radius*y,
_centre.z() + _radius*z);
// QuadStrip Edge formed at az2
// ----------------------------
// Work out the sphere normal
x = cos(elev)*sin(az2);
y = cos(elev)*cos(az2);
// z = sin(elev); z doesn't change
glNormal3f(-x, -y, -z);
glVertex3f(_centre.x() + _radius*x,
_centre.y() + _radius*y,
_centre.z() + _radius*z);
}
glEnd();
}
}
// draw front side
if (drawFrontSide)
{
for(int i=0; i+1<=_density; i++)
{
// Because we're drawing quad strips, we need to work out
// two azimuth values, to form each edge of the (z-vertical)
// strips
float az1 = _azMin + (i*azIncr);
float az2 = _azMin + ((i+1)*azIncr);
glBegin(GL_QUAD_STRIP);
for (int j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
// QuadStrip Edge formed at az1
// ----------------------------
// Work out the sphere normal
float x = cos(elev)*sin(az2);
float y = cos(elev)*cos(az2);
float z = sin(elev);
glNormal3f(x, y, z);
glVertex3f(_centre.x() + _radius*x,
_centre.y() + _radius*y,
_centre.z() + _radius*z);
// QuadStrip Edge formed at az2
// ----------------------------
// Work out the sphere normal
// z = sin(elev); z doesn't change
x = cos(elev)*sin(az1);
y = cos(elev)*cos(az1);
glNormal3f(x, y, z);
glVertex3f(_centre.x() + _radius*x,
_centre.y() + _radius*y,
_centre.z() + _radius*z);
}
glEnd();
}
}
}
}
bool SphereSegment::Surface_computeBound(osg::BoundingBox& bbox) const
{
bbox.init();
float azIncr = (_azMax - _azMin)/_density;
float elevIncr = (_elevMax - _elevMin)/_density;
for(int i=0; i<=_density; i++){
float az = _azMin + (i*azIncr);
for(int j=0; j<=_density; j++){
float elev = _elevMin + (j*elevIncr);
bbox.expandBy(
osg::Vec3(_centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev))
);
}
}
return true;
}
void SphereSegment::EdgeLine_drawImplementation(osg::State& /* state */) const
{
const float azIncr = (_azMax - _azMin)/_density;
const float elevIncr = (_elevMax - _elevMin)/_density;
// Draw the edgeline if necessary
// ------------------------------
if(_drawMask & EDGELINE)
{
glColor4fv(_edgeLineColor.ptr());
// Top edge
glBegin(GL_LINE_STRIP);
int i;
for(i=0; i<=_density; i++)
{
float az = _azMin + (i*azIncr);
glVertex3f(
_centre.x() + _radius*cos(_elevMax)*sin(az),
_centre.y() + _radius*cos(_elevMax)*cos(az),
_centre.z() + _radius*sin(_elevMax));
}
glEnd();
// Bottom edge
glBegin(GL_LINE_STRIP);
for(i=0; i<=_density; i++)
{
float az = _azMin + (i*azIncr);
glVertex3f(
_centre.x() + _radius*cos(_elevMin)*sin(az),
_centre.y() + _radius*cos(_elevMin)*cos(az),
_centre.z() + _radius*sin(_elevMin));
}
glEnd();
// Left edge
glBegin(GL_LINE_STRIP);
int j;
for(j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
glVertex3f(
_centre.x() + _radius*cos(elev)*sin(_azMin),
_centre.y() + _radius*cos(elev)*cos(_azMin),
_centre.z() + _radius*sin(elev));
}
glEnd();
// Right edge
glBegin(GL_LINE_STRIP);
for(j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
glVertex3f(
_centre.x() + _radius*cos(elev)*sin(_azMax),
_centre.y() + _radius*cos(elev)*cos(_azMax),
_centre.z() + _radius*sin(elev));
}
glEnd();
#if 0
// Split right
glBegin(GL_LINE_STRIP);
glVertex3f(
_centre.x() + _radius*cos(_elevMin)*sin(_azMax),
_centre.y() + _radius*cos(_elevMin)*cos(_azMax),
_centre.z() + _radius*sin(_elevMin));
glVertex3f(_centre.x(), _centre.y(), _centre.z());
glVertex3f(
_centre.x() + _radius*cos(_elevMax)*sin(_azMax),
_centre.y() + _radius*cos(_elevMax)*cos(_azMax),
_centre.z() + _radius*sin(_elevMax));
glEnd();
// Split left
glBegin(GL_LINE_STRIP);
glVertex3f(
_centre.x() + _radius*cos(_elevMin)*sin(_azMin),
_centre.y() + _radius*cos(_elevMin)*cos(_azMin),
_centre.z() + _radius*sin(_elevMin));
glVertex3f(_centre.x(), _centre.y(), _centre.z());
glVertex3f(
_centre.x() + _radius*cos(_elevMax)*sin(_azMin),
_centre.y() + _radius*cos(_elevMax)*cos(_azMin),
_centre.z() + _radius*sin(_elevMax));
glEnd();
#endif
}
}
bool SphereSegment::EdgeLine_computeBound(osg::BoundingBox& bbox) const
{
bbox.init();
float azIncr = (_azMax - _azMin)/_density;
float elevIncr = (_elevMax - _elevMin)/_density;
// Top edge
int i;
for(i=0; i<=_density; i++)
{
float az = _azMin + (i*azIncr);
bbox.expandBy(
_centre.x() + _radius*cos(_elevMax)*sin(az),
_centre.y() + _radius*cos(_elevMax)*cos(az),
_centre.z() + _radius*sin(_elevMax));
}
// Bottom edge
for(i=0; i<=_density; i++)
{
float az = _azMin + (i*azIncr);
bbox.expandBy(
_centre.x() + _radius*cos(_elevMin)*sin(az),
_centre.y() + _radius*cos(_elevMin)*cos(az),
_centre.z() + _radius*sin(_elevMin));
}
// Left edge
int j;
for(j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
bbox.expandBy(
_centre.x() + _radius*cos(elev)*sin(_azMin),
_centre.y() + _radius*cos(elev)*cos(_azMin),
_centre.z() + _radius*sin(elev));
}
// Right edge
for(j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
bbox.expandBy(
_centre.x() + _radius*cos(elev)*sin(_azMax),
_centre.y() + _radius*cos(elev)*cos(_azMax),
_centre.z() + _radius*sin(elev));
}
return true;
}
void SphereSegment::Side_drawImplementation(osg::State& /* state */,
SphereSegment::SideOrientation orientation,
SphereSegment::BoundaryAngle boundaryAngle) const
{
// Draw the planes if necessary
// ----------------------------
if(_drawMask & SIDES)
{
bool drawBackSide = true;
bool drawFrontSide = true;
int start, end, delta;
glColor4fv(_planeColor.ptr());
// draw back side.
if (drawBackSide)
{
if(orientation == AZIM) // This is a plane at a given azimuth
{
const float az = (boundaryAngle==MIN?_azMin:_azMax);
const float elevIncr = (_elevMax - _elevMin)/_density;
// Normal
osg::Vec3 normal = osg::Vec3(cos(_elevMin)*sin(az), cos(_elevMin)*cos(az), sin(_elevMin))
^ osg::Vec3(cos(_elevMax)*sin(az), cos(_elevMax)*cos(az), sin(_elevMax));
if (boundaryAngle==MIN)
{
start = _density;
end = 0;
}
else
{
start = 0;
end = _density;
normal = -normal; // Make sure normals orientationint 'outwards'
}
delta = end>start?1:-1;
if (drawBackSide)
{
// Tri fan
glNormal3f(-normal.x(),-normal.y(),-normal.z());
glBegin(GL_TRIANGLE_FAN);
glVertex3fv(_centre.ptr());
for (int j=start; j!=end+delta; j+=delta)
{
float elev = _elevMin + (j*elevIncr);
glVertex3f( _centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
}
glEnd();
}
if (boundaryAngle==MIN)
{
start = 0;
end = _density;
}
else
{
start = _density;
end = 0;
}
delta = end>start?1:-1;
if (drawFrontSide)
{
glNormal3fv(normal.ptr());
glBegin(GL_TRIANGLE_FAN);
glVertex3fv(_centre.ptr());
for (int j=start; j!=end+delta; j+=delta)
{
float elev = _elevMin + (j*elevIncr);
glVertex3f( _centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
}
glEnd();
}
}
else if(orientation == ELEV) // This is a plane at a given elevation
{
const float elev = (boundaryAngle==MIN?_elevMin:_elevMax);
const float azIncr = (_azMax - _azMin)/_density;
// Normal
osg::Vec3 normal = osg::Vec3(cos(elev)*sin(_azMax), cos(elev)*cos(_azMax), sin(elev))
^ osg::Vec3(cos(elev)*sin(_azMin), cos(elev)*cos(_azMin), sin(elev));
if (boundaryAngle==MIN)
{
start = _density;
end = 0;
normal = -normal; // Make sure normals orientationint 'outwards'
}
else
{
start = 0;
end = _density;
}
delta = end>start?1:-1;
if (drawBackSide)
{
glNormal3f(-normal.x(),-normal.y(),-normal.z());
// Tri fan
glBegin(GL_TRIANGLE_FAN);
glVertex3fv(_centre.ptr());
for (int j=start; j!=end+delta; j+=delta)
{
float az = _azMin + (j*azIncr);
glVertex3f( _centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
}
glEnd();
}
if (boundaryAngle==MIN)
{
start = 0;
end = _density;
}
else
{
start = _density;
end = 0;
}
delta = end>start?1:-1;
if (drawFrontSide)
{
glNormal3fv(normal.ptr());
// Tri fan
glBegin(GL_TRIANGLE_FAN);
glVertex3fv(_centre.ptr());
for (int j=start; j!=end+delta; j+=delta)
{
float az = _azMin + (j*azIncr);
glVertex3f( _centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
}
glEnd();
}
}
}
}
}
bool SphereSegment::Side_computeBound(osg::BoundingBox& bbox,
SphereSegment::SideOrientation orientation,
SphereSegment::BoundaryAngle boundaryAngle) const
{
bbox.init();
bbox.expandBy(_centre);
if(orientation == AZIM) // This is a plane at a given azimuth
{
const float az = (boundaryAngle==MIN?_azMin:_azMax);
const float elevIncr = (_elevMax - _elevMin)/_density;
for (int j=0; j<=_density; j++)
{
float elev = _elevMin + (j*elevIncr);
bbox.expandBy(
_centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
}
}
else if(orientation == ELEV) // This is a plane at a given elevation
{
const float elev = (boundaryAngle==MIN?_elevMin:_elevMax);
const float azIncr = (_azMax - _azMin)/_density;
for(int i=0; i<=_density; i++)
{
float az = _azMin + (i*azIncr);
bbox.expandBy(
_centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
}
}
return true;
}
void SphereSegment::Spoke_drawImplementation(osg::State&, BoundaryAngle azAngle, BoundaryAngle elevAngle) const
{
if(_drawMask & SPOKES){
glColor4fv(_spokeColor.ptr());
const float az = (azAngle==MIN?_azMin:_azMax);
const float elev = (elevAngle==MIN?_elevMin:_elevMax);
glBegin(GL_LINES);
glVertex3fv(_centre.ptr());
glVertex3f( _centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
glEnd();
}
}
bool SphereSegment::Spoke_computeBound(osg::BoundingBox& bbox, BoundaryAngle azAngle, BoundaryAngle elevAngle) const
{
const float az = (azAngle==MIN?_azMin:_azMax);
const float elev = (elevAngle==MIN?_elevMin:_elevMax);
bbox.expandBy(_centre);
bbox.expandBy( _centre.x() + _radius*cos(elev)*sin(az),
_centre.y() + _radius*cos(elev)*cos(az),
_centre.z() + _radius*sin(elev));
return true;
}
void SphereSegment::setDrawMask(DrawMask dm)
{
_drawMask=dm;
dirtyAllDrawableDisplayLists();
dirtyAllDrawableBounds();
dirtyBound();
}
struct ActivateTransparencyOnType
{
ActivateTransparencyOnType(const std::type_info& t): _t(t) {}
void operator()(osg::ref_ptr<osg::Drawable>& dptr) const
{
if(typeid(*dptr)==_t)
{
osg::StateSet* ss = dptr->getOrCreateStateSet();
ss->setRenderingHint(osg::StateSet::TRANSPARENT_BIN);
ss->setAttributeAndModes(new osg::CullFace(osg::CullFace::BACK),osg::StateAttribute::ON);
ss->setMode(GL_BLEND,osg::StateAttribute::ON);
dptr->dirtyDisplayList();
}
}
const std::type_info& _t;
};
struct DeactivateTransparencyOnType
{
DeactivateTransparencyOnType(const std::type_info& t): _t(t) {}
void operator()(osg::ref_ptr<osg::Drawable>& dptr) const
{
if(typeid(*dptr)==_t)
{
osg::StateSet* ss = dptr->getStateSet();
if(ss) ss->setRenderingHint(osg::StateSet::OPAQUE_BIN);
dptr->dirtyDisplayList();
}
}
const std::type_info& _t;
};
void SphereSegment::setSurfaceColor(const osg::Vec4& c)
{
_surfaceColor=c;
if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(Surface)));
else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(Surface)));
}
void SphereSegment::setSpokeColor(const osg::Vec4& c)
{
_spokeColor=c;
if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(Spoke)));
else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(Spoke)));
}
void SphereSegment::setEdgeLineColor(const osg::Vec4& c)
{
_edgeLineColor=c;
if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(EdgeLine)));
else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(EdgeLine)));
}
void SphereSegment::setSideColor(const osg::Vec4& c)
{
_planeColor=c;
if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(Side)));
else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(Side)));
}
void SphereSegment::setAllColors(const osg::Vec4& c)
{
setSurfaceColor(c);
setSpokeColor(c);
setEdgeLineColor(c);
setSideColor(c);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// SphereSegment interesection code.
class PolytopeVisitor : public osg::NodeVisitor
{
public:
typedef std::pair<osg::Matrix, osg::Polytope> MatrixPolytopePair;
typedef std::vector<MatrixPolytopePair> PolytopeStack;
struct Hit
{
Hit(const osg::Matrix& matrix, osg::NodePath& nodePath, osg::Drawable* drawable):
_matrix(matrix),
_nodePath(nodePath),
_drawable(drawable) {}
osg::Matrix _matrix;
osg::NodePath _nodePath;
osg::ref_ptr<osg::Drawable> _drawable;
};
typedef std::vector<Hit> HitList;
PolytopeVisitor(const osg::Matrix& matrix, const osg::Polytope& polytope):
osg::NodeVisitor(osg::NodeVisitor::TRAVERSE_ACTIVE_CHILDREN)
{
_polytopeStack.push_back(MatrixPolytopePair());
_polytopeStack.back().first = matrix;
_polytopeStack.back().second.setAndTransformProvidingInverse(polytope, _polytopeStack.back().first);
}
void reset()
{
_polytopeStack.clear();
_hits.clear();
}
void apply(osg::Node& node)
{
if (_polytopeStack.back().second.contains(node.getBound()))
{
traverse(node);
}
}
void apply(osg::Transform& transform)
{
if (_polytopeStack.back().second.contains(transform.getBound()))
{
const osg::Polytope& polytope = _polytopeStack.front().second;
osg::Matrix matrix = _polytopeStack.back().first;
transform.computeLocalToWorldMatrix(matrix, this);
_polytopeStack.push_back(MatrixPolytopePair());
_polytopeStack.back().first = matrix;
_polytopeStack.back().second.setAndTransformProvidingInverse(polytope, matrix);
traverse(transform);
_polytopeStack.back();
}
}
void apply(osg::Geode& node)
{
if (_polytopeStack.back().second.contains(node.getBound()))
{
for(unsigned int i=0; i<node.getNumDrawables(); ++i)
{
if (_polytopeStack.back().second.contains(node.getDrawable(i)->getBound()))
{
_hits.push_back(Hit(_polytopeStack.back().first,getNodePath(),node.getDrawable(i)));
}
}
traverse(node);
}
}
HitList& getHits() { return _hits; }
protected:
PolytopeStack _polytopeStack;
HitList _hits;
};
SphereSegment::LineList SphereSegment::computeIntersection(const osg::Matrixd& transform, osg::Node* subgraph)
{
osg::notify(osg::INFO)<<"Creating line intersection between sphere segment and subgraph."<<std::endl;
osg::BoundingBox bb = getBoundingBox();
osg::Polytope polytope;
polytope.add(osg::Plane(1.0,0.0,0.0,-bb.xMin()));
polytope.add(osg::Plane(-1.0,0.0,0.0,bb.xMax()));
polytope.add(osg::Plane(0.0,1.0,0.0,-bb.yMin()));
polytope.add(osg::Plane(0.0,-1.0,0.0,bb.yMax()));
polytope.add(osg::Plane(0.0,0.0,1.0,-bb.zMin()));
polytope.add(osg::Plane(0.0,0.0,-1.0,bb.zMax()));
osg::Plane pl;
pl.set(osg::Vec3(-1.0,0.0,0.0), bb.corner(0));
PolytopeVisitor polytopeVisitor(transform, polytope);
subgraph->accept(polytopeVisitor);
if (polytopeVisitor.getHits().empty())
{
osg::notify(osg::INFO)<<"No hits found."<<std::endl;
return LineList();
}
// create a LineList to store all the compute line segments
LineList all_lines;
// compute the line intersections with each of the hit drawables
osg::notify(osg::INFO)<<"Hits found. "<<polytopeVisitor.getHits().size()<<std::endl;
PolytopeVisitor::HitList& hits = polytopeVisitor.getHits();
for(PolytopeVisitor::HitList::iterator itr = hits.begin();
itr != hits.end();
++itr)
{
SphereSegment::LineList lines = computeIntersection(itr->_matrix, itr->_drawable.get());
all_lines.insert(all_lines.end(), lines.begin(), lines.end());
}
// join all the lines that have ends that are close together..
return all_lines;
}
osg::Node* SphereSegment::computeIntersectionSubgraph(const osg::Matrixd& transform, osg::Node* subgraph)
{
osg::notify(osg::INFO)<<"Creating line intersection between sphere segment and subgraph."<<std::endl;
osg::BoundingBox bb = getBoundingBox();
osg::Polytope polytope;
polytope.add(osg::Plane(1.0,0.0,0.0,-bb.xMin()));
polytope.add(osg::Plane(-1.0,0.0,0.0,bb.xMax()));
polytope.add(osg::Plane(0.0,1.0,0.0,-bb.yMin()));
polytope.add(osg::Plane(0.0,-1.0,0.0,bb.yMax()));
polytope.add(osg::Plane(0.0,0.0,1.0,-bb.zMin()));
polytope.add(osg::Plane(0.0,0.0,-1.0,bb.zMax()));
osg::Plane pl;
pl.set(osg::Vec3(-1.0,0.0,0.0), bb.corner(0));
PolytopeVisitor polytopeVisitor(transform, polytope);
subgraph->accept(polytopeVisitor);
if (polytopeVisitor.getHits().empty())
{
osg::notify(osg::INFO)<<"No hits found."<<std::endl;
return 0;
}
// create a LineList to store all the compute line segments
osg::Group* group = new osg::Group;
// compute the line intersections with each of the hit drawables
osg::notify(osg::INFO)<<"Hits found. "<<polytopeVisitor.getHits().size()<<std::endl;
PolytopeVisitor::HitList& hits = polytopeVisitor.getHits();
for(PolytopeVisitor::HitList::iterator itr = hits.begin();
itr != hits.end();
++itr)
{
group->addChild(computeIntersectionSubgraph(itr->_matrix, itr->_drawable.get()));
}
// join all the lines that have ends that are close together..
return group;
}
namespace SphereSegmentIntersector
{
struct dereference_less
{
template<class T, class U>
inline bool operator() (const T& lhs,const U& rhs) const
{
return *lhs < *rhs;
}
};
struct SortFunctor
{
typedef std::vector< osg::Vec3 > VertexArray;
SortFunctor(VertexArray& vertices):
_vertices(vertices) {}
bool operator() (unsigned int p1, unsigned int p2) const
{
return _vertices[p1]<_vertices[p2];
}
VertexArray& _vertices;
};
struct TriangleIntersectOperator
{
TriangleIntersectOperator():
_radius(-1.0f),
_azMin(0.0f),
_azMax(0.0f),
_elevMin(0.0f),
_elevMax(0.0f),
_numOutside(0),
_numInside(0),
_numIntersecting(0) {}
enum SurfaceType
{
NO_SURFACE = 0,
RADIUS_SURFACE,
LEFT_SURFACE,
RIGHT_SURFACE,
BOTTOM_SURFACE,
TOP_SURFACE
};
struct Triangle;
struct Edge : public osg::Referenced
{
typedef std::vector<Triangle*> TriangleList;
enum IntersectionType
{
NO_INTERSECTION,
POINT_1,
POINT_2,
MID_POINT,
BOTH_ENDS
};
Edge(unsigned int p1, unsigned int p2, SurfaceType intersectionEdge=NO_SURFACE):
_intersectionType(NO_INTERSECTION),
_intersectionVertexIndex(0),
_p1Outside(false),
_p2Outside(false),
_intersectionEdge(intersectionEdge)
{
if (p1>p2)
{
_p1 = p2;
_p2 = p1;
}
else
{
_p1 = p1;
_p2 = p2;
}
}
bool operator < (const Edge& edge) const
{
if (_p1<edge._p1) return true;
else if (_p1>edge._p1) return false;
else return _p2<edge._p2;
}
inline void addTriangle(Triangle* tri)
{
TriangleList::iterator itr = std::find(_triangles.begin(), _triangles.end(), tri);
if (itr==_triangles.end()) _triangles.push_back(tri);
}
void removeFromToTraverseList(Triangle* tri)
{
TriangleList::iterator itr = std::find(_toTraverse.begin(), _toTraverse.end(), tri);
if (itr!=_toTraverse.end()) _toTraverse.erase(itr);
}
bool completlyOutside() const { return _p1Outside && _p2Outside; }
unsigned int _p1;
unsigned int _p2;
TriangleList _triangles;
// intersection information
IntersectionType _intersectionType;
osg::Vec3 _intersectionVertex;
unsigned int _intersectionVertexIndex;
bool _p1Outside;
bool _p2Outside;
TriangleList _toTraverse;
SurfaceType _intersectionEdge;
};
typedef std::list< osg::ref_ptr<Edge> > EdgeList;
struct Triangle : public osg::Referenced
{
Triangle(unsigned int p1, unsigned int p2, unsigned int p3):
_p1(p1), _p2(p2), _p3(p3),
_e1(0), _e2(0), _e3(0)
{
sort();
}
bool operator < (const Triangle& rhs) const
{
if (_p1 < rhs._p1) return true;
else if (_p1 > rhs._p1) return false;
else if (_p2 < rhs._p2) return true;
else if (_p2 > rhs._p2) return false;
else return (_p3 < rhs._p3);
}
bool operator == (const Triangle& rhs) const
{
return (_p1 == rhs._p1) && (_p2 != rhs._p2) && (_p3 != rhs._p3);
}
bool operator != (const Triangle& rhs) const
{
return (_p1 != rhs._p1) || (_p2 != rhs._p2) || (_p3 != rhs._p3);
}
void sort()
{
if (_p1>_p2) std::swap(_p1,_p2);
if (_p1>_p3) std::swap(_p1,_p3);
if (_p2>_p3) std::swap(_p2,_p3);
}
Edge* oppositeActiveEdge(Edge* edge)
{
if (edge!=_e1 && edge!=_e2 && edge!=_e3)
{
osg::notify(osg::INFO)<<"Edge problem"<<std::endl;
return 0;
}
if (edge!=_e1 && _e1 && _e1->_intersectionType!=Edge::NO_INTERSECTION) return _e1;
if (edge!=_e2 && _e2 && _e2->_intersectionType!=Edge::NO_INTERSECTION) return _e2;
if (edge!=_e3 && _e3 && _e3->_intersectionType!=Edge::NO_INTERSECTION) return _e3;
return 0;
}
unsigned int _p1;
unsigned int _p2;
unsigned int _p3;
Edge* _e1;
Edge* _e2;
Edge* _e3;
};
struct Region
{
enum Classification
{
INSIDE = -1,
INTERSECTS = 0,
OUTSIDE = 1
};
Region():
_radiusSurface(OUTSIDE),
_leftSurface(OUTSIDE),
_rightSurface(OUTSIDE),
_bottomSurface(OUTSIDE),
_topSurface(OUTSIDE) {}
void classify(const osg::Vec3& vertex, double radius2, double azimMin, double azimMax, double elevMin, double elevMax)
{
double azimCenter = (azimMax+azimMin)*0.5;
double azimRange = (azimMax-azimMin)*0.5;
double rad2 = vertex.length2();
double length_xy = sqrtf(vertex.x()*vertex.x() + vertex.y()*vertex.y());
double elevation = atan2((double)vertex.z(),length_xy);
// radius surface
if (rad2 > radius2) _radiusSurface = OUTSIDE;
else if (rad2 < radius2) _radiusSurface = INSIDE;
else _radiusSurface = INTERSECTS;
// bottom surface
if (elevation<elevMin) _bottomSurface = OUTSIDE;
else if (elevation>elevMin) _bottomSurface = INSIDE;
else _bottomSurface = INTERSECTS;
// top surface
if (elevation>elevMax) _topSurface = OUTSIDE;
else if (elevation<elevMax) _topSurface = INSIDE;
else _topSurface = INTERSECTS;
double dot_alphaMin = cos(azimMin)*vertex.x() - sin(azimMin)*vertex.y();
if (dot_alphaMin<0.0) _leftSurface = OUTSIDE;
else if (dot_alphaMin>0.0) _leftSurface = INSIDE;
else _leftSurface = INTERSECTS;
double dot_alphaMax = cos(azimMax)*vertex.x() - sin(azimMax)*vertex.y();
if (dot_alphaMax>0.0) _rightSurface = OUTSIDE;
else if (dot_alphaMax<0.0) _rightSurface = INSIDE;
else _rightSurface = INTERSECTS;
double azim = atan2(vertex.x(),vertex.y());
double azimDelta1 = azim-azimCenter;
double azimDelta2 = 2.0*osg::PI + azim-azimCenter;
double azimDelta = std::min(fabs(azimDelta1), fabs(azimDelta2));
if (azimDelta>azimRange)
{
_leftRightSurfaces = OUTSIDE;
}
else if (azimDelta<azimRange)
{
_leftRightSurfaces = INSIDE;
}
else if (azimDelta==azimRange)
{
_leftRightSurfaces = INTERSECTS;
}
}
Classification _radiusSurface;
Classification _leftRightSurfaces;
Classification _leftSurface;
Classification _rightSurface;
Classification _bottomSurface;
Classification _topSurface;
};
struct RegionCounter
{
RegionCounter():
_numVertices(0),
_outside_radiusSurface(0),
_inside_radiusSurface(0),
_intersects_radiusSurface(0),
_outside_leftRightSurfaces(0),
_inside_leftRightSurfaces(0),
_intersects_leftRightSurfaces(0),
_outside_leftSurface(0),
_inside_leftSurface(0),
_intersects_leftSurface(0),
_outside_rightSurface(0),
_inside_rightSurface(0),
_intersects_rightSurface(0),
_outside_bottomSurface(0),
_inside_bottomSurface(0),
_intersects_bottomSurface(0),
_outside_topSurface(0),
_inside_topSurface(0),
_intersects_topSurface(0) {}
void add(const Region& region)
{
++_numVertices;
if (region._radiusSurface == Region::OUTSIDE) ++_outside_radiusSurface;
if (region._radiusSurface == Region::INSIDE) ++_inside_radiusSurface;
if (region._radiusSurface == Region::INTERSECTS) ++_intersects_radiusSurface;
if (region._leftRightSurfaces == Region::OUTSIDE) ++_outside_leftRightSurfaces;
if (region._leftRightSurfaces == Region::INSIDE) ++_inside_leftRightSurfaces;
if (region._leftRightSurfaces == Region::INTERSECTS) ++_intersects_leftRightSurfaces;
if (region._leftSurface == Region::OUTSIDE) ++_outside_leftSurface;
if (region._leftSurface == Region::INSIDE) ++_inside_leftSurface;
if (region._leftSurface == Region::INTERSECTS) ++_intersects_leftSurface;
if (region._rightSurface == Region::OUTSIDE) ++_outside_rightSurface;
if (region._rightSurface == Region::INSIDE) ++_inside_rightSurface;
if (region._rightSurface == Region::INTERSECTS) ++_intersects_rightSurface;
if (region._bottomSurface == Region::OUTSIDE) ++_outside_bottomSurface;
if (region._bottomSurface == Region::INSIDE) ++_inside_bottomSurface;
if (region._bottomSurface == Region::INTERSECTS) ++_intersects_bottomSurface;
if (region._topSurface == Region::OUTSIDE) ++_outside_topSurface;
if (region._topSurface == Region::INSIDE) ++_inside_topSurface;
if (region._topSurface == Region::INTERSECTS) ++_intersects_topSurface;
}
Region::Classification overallClassification() const
{
// if all vertices are outside any of the surfaces then we are completely outside
if (_outside_radiusSurface==_numVertices ||
_outside_leftRightSurfaces==_numVertices ||
_outside_topSurface==_numVertices ||
_outside_bottomSurface==_numVertices) return Region::OUTSIDE;
// if all the vertices on all the sides and inside then we are completely inside
if (_inside_radiusSurface==_numVertices &&
_inside_leftRightSurfaces==_numVertices &&
_inside_topSurface==_numVertices &&
_inside_bottomSurface==_numVertices) return Region::INSIDE;
return Region::INTERSECTS;
}
bool intersecting(SurfaceType surfaceType) const
{
// if all vertices are outside any of the surfaces then we are completely outside
if ((surfaceType!=RADIUS_SURFACE && _outside_radiusSurface!=0) ||
(surfaceType!=LEFT_SURFACE && _outside_leftSurface!=0) ||
(surfaceType!=RIGHT_SURFACE && _outside_rightSurface!=0) ||
(surfaceType!=TOP_SURFACE && _outside_topSurface!=0) ||
(surfaceType!=BOTTOM_SURFACE && _outside_bottomSurface!=0)) return false;
// if all the vertices on all the sides and inside then we are completely inside
if ((surfaceType!=RADIUS_SURFACE && _inside_radiusSurface!=0) &&
(surfaceType!=LEFT_SURFACE && _inside_leftSurface!=0) &&
(surfaceType!=RIGHT_SURFACE && _inside_rightSurface!=0) &&
(surfaceType!=TOP_SURFACE && _inside_topSurface!=0) &&
(surfaceType!=BOTTOM_SURFACE && _inside_bottomSurface!=0)) return false;
return true;
}
int numberOfIntersectingSurfaces() const
{
int sidesThatIntersect = 0;
if (_outside_radiusSurface!=_numVertices && _inside_radiusSurface!=_numVertices) ++sidesThatIntersect;
if (_outside_leftSurface!=_numVertices && _inside_leftSurface!=_numVertices) ++sidesThatIntersect;
if (_outside_rightSurface!=_numVertices && _inside_rightSurface!=_numVertices) ++sidesThatIntersect;
if (_outside_topSurface!=_numVertices && _inside_topSurface!=_numVertices) ++sidesThatIntersect;
if (_outside_bottomSurface!=_numVertices && _inside_bottomSurface!=_numVertices) ++sidesThatIntersect;
return sidesThatIntersect;
}
unsigned int _numVertices;
unsigned int _outside_radiusSurface;
unsigned int _inside_radiusSurface;
unsigned int _intersects_radiusSurface;
unsigned int _outside_leftRightSurfaces;
unsigned int _inside_leftRightSurfaces;
unsigned int _intersects_leftRightSurfaces;
unsigned int _outside_leftSurface;
unsigned int _inside_leftSurface;
unsigned int _intersects_leftSurface;
unsigned int _outside_rightSurface;
unsigned int _inside_rightSurface;
unsigned int _intersects_rightSurface;
unsigned int _outside_bottomSurface;
unsigned int _inside_bottomSurface;
unsigned int _intersects_bottomSurface;
unsigned int _outside_topSurface;
unsigned int _inside_topSurface;
unsigned int _intersects_topSurface;
};
typedef std::vector< osg::Vec3 > VertexArray;
typedef std::vector< Region > RegionArray;
typedef std::vector< bool > BoolArray;
typedef std::vector< unsigned int > IndexArray;
typedef std::vector< osg::ref_ptr<Triangle> > TriangleArray;
typedef std::set< osg::ref_ptr<Edge>, dereference_less > EdgeSet;
VertexArray _originalVertices;
RegionArray _regions;
BoolArray _vertexInIntersectionSet;
IndexArray _candidateVertexIndices;
IndexArray _remapIndices;
TriangleArray _triangles;
EdgeSet _edges;
osg::Vec3 _centre;
double _radius;
double _azMin, _azMax, _elevMin, _elevMax;
unsigned int _numOutside;
unsigned int _numInside;
unsigned int _numIntersecting;
SphereSegment::LineList _generatedLines;
void computePositionAndRegions(const osg::Matrixd& matrix, osg::Vec3Array& array)
{
_originalVertices.resize(array.size());
_regions.resize(array.size());
_vertexInIntersectionSet.resize(array.size(), false);
_candidateVertexIndices.clear();
double radius2 = _radius*_radius;
for(unsigned int i=0; i<array.size(); ++i)
{
osg::Vec3 vertex = array[i]*matrix - _centre;
_originalVertices[i] = vertex;
_regions[i].classify(vertex, radius2, _azMin, _azMax, _elevMin, _elevMax);
}
}
inline void operator()(unsigned int p1, unsigned int p2, unsigned int p3)
{
RegionCounter rc;
rc.add(_regions[p1]);
rc.add(_regions[p2]);
rc.add(_regions[p3]);
Region::Classification classification = rc.overallClassification();
// reject if outside.
if (classification==Region::OUTSIDE)
{
++_numOutside;
return;
}
if (rc.numberOfIntersectingSurfaces()==0)
{
++_numInside;
return;
}
++_numIntersecting;
_triangles.push_back(new Triangle(p1,p2,p3));
if (!_vertexInIntersectionSet[p1])
{
_vertexInIntersectionSet[p1] = true;
_candidateVertexIndices.push_back(p1);
}
if (!_vertexInIntersectionSet[p2])
{
_vertexInIntersectionSet[p2] = true;
_candidateVertexIndices.push_back(p2);
}
if (!_vertexInIntersectionSet[p3])
{
_vertexInIntersectionSet[p3] = true;
_candidateVertexIndices.push_back(p3);
}
}
void removeDuplicateVertices()
{
osg::notify(osg::INFO)<<"Removing duplicates : num vertices in "<<_candidateVertexIndices.size()<<std::endl;
if (_candidateVertexIndices.size()<2) return;
std::sort(_candidateVertexIndices.begin(), _candidateVertexIndices.end(), SortFunctor(_originalVertices));
_remapIndices.resize(_originalVertices.size());
for(unsigned int i=0; i< _originalVertices.size(); ++i)
{
_remapIndices[i] = i;
}
bool verticesRemapped = false;
IndexArray::iterator itr = _candidateVertexIndices.begin();
unsigned int lastUniqueIndex = *(itr++);
for(; itr != _candidateVertexIndices.end(); ++itr)
{
//unsigned int i = *itr;
// osg::notify(osg::INFO)<<" i="<<i<<" lastUniqueIndex="<<lastUniqueIndex<<std::endl;
if (_originalVertices[*itr]==_originalVertices[lastUniqueIndex])
{
osg::notify(osg::INFO)<<"Combining vertex "<<*itr<<" with "<<lastUniqueIndex<<std::endl;
_remapIndices[*itr] = lastUniqueIndex;
verticesRemapped = true;
}
else
{
lastUniqueIndex = *itr;
}
}
if (verticesRemapped)
{
osg::notify(osg::INFO)<<"Remapping triangle vertices "<<std::endl;
for(TriangleArray::iterator titr = _triangles.begin();
titr != _triangles.end();
++titr)
{
(*titr)->_p1 = _remapIndices[(*titr)->_p1];
(*titr)->_p2 = _remapIndices[(*titr)->_p2];
(*titr)->_p3 = _remapIndices[(*titr)->_p3];
(*titr)->sort();
}
}
}
void removeDuplicateTriangles()
{
osg::notify(osg::INFO)<<"Removing duplicate triangles : num triangles in "<<_triangles.size()<<std::endl;
if (_triangles.size()<2) return;
std::sort(_triangles.begin(), _triangles.end(), dereference_less());
unsigned int lastUniqueTriangle = 0;
unsigned int numDuplicates = 0;
for(unsigned int i=1; i<_triangles.size(); ++i)
{
if ( *(_triangles[lastUniqueTriangle]) != *(_triangles[i]) )
{
++lastUniqueTriangle;
if (lastUniqueTriangle!=i)
{
_triangles[lastUniqueTriangle] = _triangles[i];
}
}
else
{
++numDuplicates;
}
}
if (lastUniqueTriangle<_triangles.size()-1)
{
_triangles.erase(_triangles.begin()+lastUniqueTriangle+1, _triangles.end());
}
osg::notify(osg::INFO)<<"Removed duplicate triangles : num duplicates found "<<numDuplicates<<std::endl;
osg::notify(osg::INFO)<<"Removed duplicate triangles : num triangles out "<<_triangles.size()<<std::endl;
}
void buildEdges(Triangle* tri)
{
tri->_e1 = addEdge(tri->_p1, tri->_p2, tri);
tri->_e2 = addEdge(tri->_p2, tri->_p3, tri);
tri->_e3 = addEdge(tri->_p1, tri->_p3, tri);
}
void buildEdges()
{
_edges.clear();
for(TriangleArray::iterator itr = _triangles.begin();
itr != _triangles.end();
++itr)
{
Triangle* tri = itr->get();
RegionCounter rc;
rc.add(_regions[tri->_p1]);
rc.add(_regions[tri->_p2]);
rc.add(_regions[tri->_p3]);
int numIntersections = rc.numberOfIntersectingSurfaces();
if (numIntersections>=1)
{
buildEdges(tri);
}
}
osg::notify(osg::INFO)<<"Number of edges "<<_edges.size()<<std::endl;
unsigned int numZeroConnections = 0;
unsigned int numSingleConnections = 0;
unsigned int numDoubleConnections = 0;
unsigned int numMultiConnections = 0;
osg::notify(osg::INFO)<<"Number of edges "<<_edges.size()<<std::endl;
for(EdgeSet::iterator eitr = _edges.begin();
eitr != _edges.end();
++eitr)
{
const Edge* edge = eitr->get();
unsigned int numConnections = edge->_triangles.size();
if (numConnections==0) ++numZeroConnections;
else if (numConnections==1) ++numSingleConnections;
else if (numConnections==2) ++numDoubleConnections;
else ++numMultiConnections;
}
osg::notify(osg::INFO)<<"Number of numZeroConnections "<<numZeroConnections<<std::endl;
osg::notify(osg::INFO)<<"Number of numSingleConnections "<<numSingleConnections<<std::endl;
osg::notify(osg::INFO)<<"Number of numDoubleConnections "<<numDoubleConnections<<std::endl;
osg::notify(osg::INFO)<<"Number of numMultiConnections "<<numMultiConnections<<std::endl;
}
Edge* addEdge(unsigned int p1, unsigned int p2, Triangle* tri)
{
osg::ref_ptr<Edge> edge = new Edge(p1, p2);
EdgeSet::iterator itr = _edges.find(edge);
if (itr==_edges.end())
{
edge->addTriangle(tri);
_edges.insert(edge);
return edge.get();
}
else
{
Edge* edge = const_cast<Edge*>(itr->get());
edge->addTriangle(tri);
return edge;
}
}
SphereSegment::LineList connectIntersections(EdgeList& hitEdges)
{
SphereSegment::LineList lineList;
osg::notify(osg::INFO)<<"Number of edge intersections "<<hitEdges.size()<<std::endl;
if (hitEdges.empty()) return lineList;
// now need to build the toTraverse list for each hit edge,
// but should only contain traingles that actually hit the intersection surface
EdgeList::iterator hitr;
for(hitr = hitEdges.begin();
hitr != hitEdges.end();
++hitr)
{
Edge* edge = hitr->get();
edge->_toTraverse.clear();
//osg::notify(osg::INFO)<<"edge= "<<edge<<std::endl;
for(Edge::TriangleList::iterator titr = edge->_triangles.begin();
titr != edge->_triangles.end();
++titr)
{
Triangle* tri = *titr;
// count how many active edges there are on this triangle
unsigned int numActiveEdges = 0;
unsigned int numEdges = 0;
if (tri->_e1 && tri->_e1->_intersectionType!=Edge::NO_INTERSECTION) ++numActiveEdges;
if (tri->_e2 && tri->_e2->_intersectionType!=Edge::NO_INTERSECTION) ++numActiveEdges;
if (tri->_e3 && tri->_e3->_intersectionType!=Edge::NO_INTERSECTION) ++numActiveEdges;
if (tri->_e1) ++numEdges;
if (tri->_e2) ++numEdges;
if (tri->_e3) ++numEdges;
// if we have one or more then add it into the edges to traverse list
if (numActiveEdges>1)
{
//osg::notify(osg::INFO)<<" adding tri="<<tri<<std::endl;
edge->_toTraverse.push_back(tri);
}
// osg::notify(osg::INFO)<<"Number active edges "<<numActiveEdges<<" num orignal edges "<<numEdges<<std::endl;
}
}
while(!hitEdges.empty())
{
// find the an open edge
for(hitr = hitEdges.begin();
hitr != hitEdges.end();
++hitr)
{
Edge* edge = hitr->get();
if (edge->_toTraverse.size()==1) break;
}
if (hitr == hitEdges.end())
{
hitr = hitEdges.begin();
}
// osg::notify(osg::INFO)<<"New line "<<std::endl;
osg::Vec3Array* newLine = new osg::Vec3Array;
lineList.push_back(newLine);
Edge* edge = hitr->get();
while (edge)
{
// osg::notify(osg::INFO)<<" vertex "<<edge->_intersectionVertex<<std::endl;
newLine->push_back(edge->_intersectionVertex+_centre/*+osg::Vec3(0.0f,0.0f,200.0f)*/);
Edge* newEdge = 0;
Triangle* tri = !(edge->_toTraverse.empty()) ? edge->_toTraverse.back() : 0;
if (tri)
{
newEdge = tri->oppositeActiveEdge(edge);
edge->removeFromToTraverseList(tri);
newEdge->removeFromToTraverseList(tri);
// osg::notify(osg::INFO)<<" tri="<<tri<<" edge="<<edge<<" newEdge="<<newEdge<<std::endl;
if (edge==newEdge)
{
osg::notify(osg::INFO)<<" edge returned to itself problem "<<std::endl;
}
}
else
{
newEdge = 0;
}
if (edge->_toTraverse.empty())
{
edge->_intersectionType = Edge::NO_INTERSECTION;
// remove edge for the hitEdges.
hitr = std::find(hitEdges.begin(), hitEdges.end(), edge);
if (hitr!=hitEdges.end()) hitEdges.erase(hitr);
}
// move on to next edge in line.
edge = newEdge;
}
}
return lineList;
}
template<class I>
SphereSegment::LineList computeIntersections(I intersector)
{
// collect all the intersecting edges
EdgeList hitEdges;
for(EdgeSet::iterator itr = _edges.begin();
itr != _edges.end();
++itr)
{
Edge* edge = const_cast<Edge*>(itr->get());
if (intersector(edge))
{
hitEdges.push_back(edge);
}
}
return connectIntersections(hitEdges);
}
template<class I>
void trim(SphereSegment::LineList& lineList, osg::Vec3Array* sourceLine, I intersector)
{
if (sourceLine->empty()) return;
// osg::notify(osg::INFO)<<"Testing line of "<<sourceLine->size()<<std::endl;
unsigned int first=0;
while (first<sourceLine->size())
{
// find first valid vertex.
for(; first<sourceLine->size(); ++first)
{
if (intersector.distance((*sourceLine)[first]-_centre)>=0.0) break;
}
if (first==sourceLine->size())
{
// osg::notify(osg::INFO)<<"No valid points found"<<std::endl;
return;
}
// find last valid vertex.
unsigned int last = first+1;
for(; last<sourceLine->size(); ++last)
{
if (intersector.distance((*sourceLine)[last]-_centre)<0.0) break;
}
if (first==0 && last==sourceLine->size())
{
// osg::notify(osg::INFO)<<"Copying complete line"<<std::endl;
lineList.push_back(sourceLine);
}
else
{
// osg::notify(osg::INFO)<<"Copying partial line line"<<first<<" "<<last<<std::endl;
osg::Vec3Array* newLine = new osg::Vec3Array;
if (first>0)
{
newLine->push_back(intersector.intersectionPoint((*sourceLine)[first-1]-_centre, (*sourceLine)[first]-_centre)+_centre);
}
for(unsigned int i=first; i<last; ++i)
{
newLine->push_back((*sourceLine)[i]);
}
if (last<sourceLine->size())
{
newLine->push_back(intersector.intersectionPoint((*sourceLine)[last-1]-_centre, (*sourceLine)[last]-_centre)+_centre);
}
lineList.push_back(newLine);
}
first = last;
}
}
// handle a paird of surfaces that work to enclose a convex region, which means that
// points can be inside either surface to be valid, and be outside both surfaces to be invalid.
template<class I>
void trim(SphereSegment::LineList& lineList, osg::Vec3Array* sourceLine, I intersector1, I intersector2)
{
if (sourceLine->empty()) return;
// osg::notify(osg::INFO)<<"Testing line of "<<sourceLine->size()<<std::endl;
unsigned int first=0;
while (first<sourceLine->size())
{
// find first valid vertex.
for(; first<sourceLine->size(); ++first)
{
osg::Vec3 v = (*sourceLine)[first]-_centre;
if (intersector1.distance(v)>=0.0 || intersector2.distance(v)>=0.0 ) break;
}
if (first==sourceLine->size())
{
// osg::notify(osg::INFO)<<"No valid points found"<<std::endl;
return;
}
// find last valid vertex.
unsigned int last = first+1;
for(; last<sourceLine->size(); ++last)
{
osg::Vec3 v = (*sourceLine)[last]-_centre;
if (intersector1.distance(v)<0.0 && intersector2.distance(v)<0.0 ) break;
}
if (first==0 && last==sourceLine->size())
{
// osg::notify(osg::INFO)<<"Copying complete line"<<std::endl;
lineList.push_back(sourceLine);
}
else
{
osg::notify(osg::INFO)<<"Copying partial line line"<<first<<" "<<last<<std::endl;
osg::Vec3Array* newLine = new osg::Vec3Array;
if (first>0)
{
osg::Vec3 start = (*sourceLine)[first-1]-_centre;
osg::Vec3 end = (*sourceLine)[first]-_centre;
float end1 = intersector1.distance(end);
float end2 = intersector2.distance(end);
// work out which intersector to use by discounting the one that
// isn't a plausible candiate.
bool possible1 = end1>=0.0;
bool possible2 = end2>=0.0;
if (possible1 && possible2)
{
double start1 = intersector1.distance(start);
double start2 = intersector2.distance(start);
// need to check which intersection is latest.
double end1 = intersector1.distance(end);
double delta1 = (end1-start1);
double end2 = intersector2.distance(end);
double delta2 = (end2-start2);
double r1 = fabs(delta1)>0.0 ? start1/delta1 : 0.0;
double r2 = fabs(delta2)>0.0 ? start2/delta2 : 0.0;
// choose intersection which is nearest the end point.
if (r1<r2)
{
osg::notify(osg::INFO)<<"start point, 1 near to end than 2"<<r1<<" "<<r2<<std::endl;
possible1 = true;
possible2 = false;
}
else
{
osg::notify(osg::INFO)<<"start point, 2 near to end than 1"<<std::endl;
possible1 = false;
possible2 = true;
}
}
if (possible1)
{
newLine->push_back(intersector1.intersectionPoint(start, end)+_centre);
}
else
{
newLine->push_back(intersector2.intersectionPoint(start, end)+_centre);
}
}
for(unsigned int i=first; i<last; ++i)
{
newLine->push_back((*sourceLine)[i]);
}
if (last<sourceLine->size())
{
osg::Vec3 start = (*sourceLine)[last-1]-_centre;
osg::Vec3 end = (*sourceLine)[last]-_centre;
double start1 = intersector1.distance(start);
double start2 = intersector2.distance(start);
// work out which intersector to use by discounting the one that
// isn't a plausible candiate.
bool possible1 = start1>=0.0;
bool possible2 = start2>=0.0;
if (possible1 && possible2)
{
double end1 = intersector1.distance(end);
double end2 = intersector2.distance(end);
possible1 = end1<0.0;
possible2 = end2<0.0;
if (possible1 && possible2)
{
// need to check which intersection is latest.
double end1 = intersector1.distance(end);
double delta1 = (end1-start1);
double end2 = intersector2.distance(end);
double delta2 = (end2-start2);
double r1 = fabs(delta1)>0.0 ? start1/delta1 : 0.0;
double r2 = fabs(delta2)>0.0 ? start2/delta2 : 0.0;
// choose intersection which is nearest the end point.
if (r1>r2)
{
osg::notify(osg::INFO)<<"end point, 1 near to end than 2"<<r1<<" "<<r2<<std::endl;
possible1 = true;
possible2 = false;
}
else
{
osg::notify(osg::INFO)<<"end point, 2 near to end than 1"<<std::endl;
possible1 = false;
possible2 = true;
}
}
}
if (possible1)
{
newLine->push_back(intersector1.intersectionPoint(start, end)+_centre);
}
else
{
newLine->push_back(intersector2.intersectionPoint(start, end)+_centre);
}
}
lineList.push_back(newLine);
}
first = last;
}
}
template<class I>
void trim(SphereSegment::LineList& lineList, I intersector)
{
SphereSegment::LineList newLines;
// collect all the intersecting edges
for(SphereSegment::LineList::iterator itr = lineList.begin();
itr != lineList.end();
++itr)
{
osg::Vec3Array* line = itr->get();
trim(newLines, line, intersector);
}
lineList.swap(newLines);
}
template<class I>
void trim(SphereSegment::LineList& lineList, I intersector1, I intersector2)
{
SphereSegment::LineList newLines;
// collect all the intersecting edges
for(SphereSegment::LineList::iterator itr = lineList.begin();
itr != lineList.end();
++itr)
{
osg::Vec3Array* line = itr->get();
trim(newLines, line, intersector1, intersector2);
}
lineList.swap(newLines);
}
struct LinePair
{
LinePair(osg::Vec3Array* line):
_line(line),
_lineEnd(0),
_neighbourLine(0),
_neighbourLineEnd(0) {}
bool operator < (const LinePair& linePair) const
{
return _distance < linePair._distance;
}
void consider(osg::Vec3Array* testline)
{
if (_neighbourLine.valid())
{
float distance = ((*_line)[0]-(*testline)[0]).length();
if (distance<_distance)
{
_lineEnd = 0;
_neighbourLine = testline;
_neighbourLineEnd = 0;
_distance = distance;
}
distance = ((*_line)[0]-(*testline)[testline->size()-1]).length();
if (distance<_distance)
{
_lineEnd = 0;
_neighbourLine = testline;
_neighbourLineEnd = testline->size()-1;
_distance = distance;
}
distance = ((*_line)[_line->size()-1]-(*testline)[0]).length();
if (distance<_distance)
{
_lineEnd = _line->size()-1;
_neighbourLine = testline;
_neighbourLineEnd = 0;
_distance = distance;
}
distance = ((*_line)[_line->size()-1]-(*testline)[testline->size()-1]).length();
if (distance<_distance)
{
_lineEnd = _line->size()-1;
_neighbourLine = testline;
_neighbourLineEnd = testline->size()-1;
_distance = distance;
}
}
else
{
_neighbourLine = testline;
if (_neighbourLine==_line)
{
_lineEnd = 0;
_neighbourLineEnd = _neighbourLine->size()-1;
_distance = ((*_line)[_lineEnd]-(*_neighbourLine)[_neighbourLineEnd]).length();
}
else
{
_distance = ((*_line)[0]-(*_neighbourLine)[0]).length();
_lineEnd = 0;
_neighbourLineEnd = 0;
float distance = ((*_line)[0]-(*_neighbourLine)[_neighbourLine->size()-1]).length();
if (distance<_distance)
{
_lineEnd = 0;
_neighbourLineEnd = _neighbourLine->size()-1;
_distance = distance;
}
distance = ((*_line)[_line->size()-1]-(*_neighbourLine)[0]).length();
if (distance<_distance)
{
_lineEnd = _line->size()-1;
_neighbourLineEnd = 0;
_distance = distance;
}
distance = ((*_line)[_line->size()-1]-(*_neighbourLine)[_neighbourLine->size()-1]).length();
if (distance<_distance)
{
_lineEnd = _line->size()-1;
_neighbourLineEnd = _neighbourLine->size()-1;
_distance = distance;
}
}
}
};
bool contains(osg::Vec3Array* line) const
{
return _line==line || _neighbourLine==line;
}
osg::ref_ptr<osg::Vec3Array> _line;
unsigned int _lineEnd;
osg::ref_ptr<osg::Vec3Array> _neighbourLine;
unsigned int _neighbourLineEnd;
float _distance;
};
void joinEnds(float fuseDistance, bool doFuse, bool allowJoinToSelf)
{
SphereSegment::LineList fusedLines;
SphereSegment::LineList unfusedLines;
// first seperat the already fused lines from the unfused ones.
for(SphereSegment::LineList::iterator itr = _generatedLines.begin();
itr != _generatedLines.end();
++itr)
{
osg::Vec3Array* line = itr->get();
if (line->size()>=2)
{
if ((*line)[0]==(*line)[line->size()-1])
{
fusedLines.push_back(line);
}
else
{
unfusedLines.push_back(line);
}
}
}
while (unfusedLines.size()>=1)
{
// generate a set of line pairs to establish which
// line pair has the minimum distance.
typedef std::multiset<LinePair> LinePairSet;
LinePairSet linePairs;
for(unsigned int i=0; i<unfusedLines.size(); ++i)
{
unsigned int j = allowJoinToSelf ? i : i+1;
if (j<unfusedLines.size())
{
LinePair linePair(unfusedLines[i].get());
for(; j<unfusedLines.size(); ++j)
{
linePair.consider(unfusedLines[j].get());
}
linePairs.insert(linePair);
}
}
if (linePairs.empty())
{
osg::notify(osg::INFO)<<"Line Pairs empty"<<std::endl;
break;
}
for(LinePairSet::iterator itr = linePairs.begin();
itr != linePairs.end();
++itr)
{
osg::notify(osg::INFO)<<"Line "<<itr->_line.get()<<" "<<itr->_lineEnd<<" neighbour "<<itr->_neighbourLine.get()<<" "<<itr->_neighbourLineEnd<<" distance="<<itr->_distance<<std::endl;
}
LinePair linePair = *linePairs.begin();
if (linePair._distance > fuseDistance)
{
osg::notify(osg::INFO)<<"Completed work, shortest distance left is "<<linePair._distance<<std::endl;
break;
}
if (linePair._line == linePair._neighbourLine)
{
osg::notify(osg::INFO)<<"Fusing line to itself"<<std::endl;
osg::Vec3Array* line = linePair._line.get();
osg::Vec3 average = ((*line)[0]+(*line)[line->size()-1])*0.5f;
if (doFuse)
{
(*line)[0] = average;
(*line)[line->size()-1] = average;
}
else
{
// add start of line to end.
line->push_back((*line)[0]);
}
fusedLines.push_back(line);
SphereSegment::LineList::iterator fitr = std::find(unfusedLines.begin(), unfusedLines.end(), line);
if (fitr != unfusedLines.end())
{
unfusedLines.erase(fitr);
}
else
{
osg::notify(osg::INFO)<<"Error couldn't find line in unfused list, exiting fusing loop."<<std::endl;
break;
}
}
else
{
osg::Vec3Array* line1 = linePair._line.get();
int fuseEnd1 = linePair._lineEnd;
int openEnd1 = fuseEnd1==0 ? line1->size()-1 : 0;
int direction1 = openEnd1<fuseEnd1 ? 1 : -1;
osg::Vec3Array* line2 = linePair._neighbourLine.get();
int fuseEnd2 = linePair._neighbourLineEnd;
int openEnd2 = fuseEnd2==0 ? line2->size()-1 : 0;
int direction2 = fuseEnd2<openEnd2 ? 1 : -1;
osg::Vec3Array* newline = new osg::Vec3Array;
// copy across all but fuse end of line1
for(int i=openEnd1;
i != fuseEnd1;
i += direction1)
{
newline->push_back((*line1)[i]);
}
// add the average of the two fused ends
if (doFuse)
{
osg::Vec3 average = ((*line1)[fuseEnd1] + (*line2)[fuseEnd2])*0.5f;
newline->push_back(average);
}
else
{
newline->push_back((*line1)[fuseEnd1]);
newline->push_back((*line2)[fuseEnd2]);
}
// copy across from the next point in from fuseEnd2 to the openEnd2.
for(int j=fuseEnd2 + direction2;
j != openEnd2 + direction2;
j += direction2)
{
newline->push_back((*line2)[j]);
}
// remove line1 from unfused list.
SphereSegment::LineList::iterator fitr = std::find(unfusedLines.begin(), unfusedLines.end(), line1);
if (fitr != unfusedLines.end())
{
unfusedLines.erase(fitr);
}
// remove line2 from unfused list.
fitr = std::find(unfusedLines.begin(), unfusedLines.end(), line2);
if (fitr != unfusedLines.end())
{
unfusedLines.erase(fitr);
}
// add the newline into the unfused for further processing.
unfusedLines.push_back(newline);
osg::notify(osg::INFO)<<"Fusing two seperate lines "<<newline<<std::endl;
}
_generatedLines = fusedLines;
_generatedLines.insert(_generatedLines.end(), unfusedLines.begin(), unfusedLines.end());
}
}
};
bool computeQuadraticSolution(double a, double b, double c, double& s1, double& s2)
{
// avoid division by zero.
if (a==0.0)
{
s1 = 0.0;
s2 = 0.0;
return false;
}
double inside_sqrt = b*b - 4.0*a*c;
// avoid sqrt of negative number
if (inside_sqrt<0.0)
{
s1 = 0.0;
s2 = 0.0;
return false;
}
double rhs = sqrt(inside_sqrt);
s1 = (-b + rhs)/(2.0*a);
s2 = (-b - rhs)/(2.0*a);
return true;
}
struct AzimPlaneIntersector
{
AzimPlaneIntersector(TriangleIntersectOperator& tif, double azim, bool lowerOutside):
_tif(tif),
_lowerOutside(lowerOutside)
{
_plane.set(cos(azim),-sin(azim),0.0,0.0);
_endPlane.set(sin(azim),cos(azim),0.0,0.0);
}
TriangleIntersectOperator& _tif;
osg::Plane _plane;
osg::Plane _endPlane;
bool _lowerOutside;
inline bool operator() (TriangleIntersectOperator::Edge* edge)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
osg::Vec3& v1 = _tif._originalVertices[edge->_p1];
osg::Vec3& v2 = _tif._originalVertices[edge->_p2];
double d1 = _plane.distance(v1);
double d2 = _plane.distance(v2);
edge->_p1Outside = _lowerOutside ? (d1<0.0) : (d1>0.0);
edge->_p2Outside = _lowerOutside ? (d2<0.0) : (d2>0.0);
// if both points inside then disgard
if (d1<0.0 && d2<0.0) return false;
// if both points outside then disgard
if (d1>0.0 && d2>0.0) return false;
if (d1==0.0)
{
if (d2==0.0)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::BOTH_ENDS;
}
else
{
edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_1;
}
}
else if (d2==0.0)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_2;
}
else
{
double div = d2-d1;
if (div==0.0)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
return false;
}
double r = -d1 / div;
if (r<0.0 || r>1.0)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
return false;
}
// osg::notify(osg::INFO)<<"r = "<<r<<std::endl;
double one_minus_r = 1.0-r;
edge->_intersectionType = TriangleIntersectOperator::Edge::MID_POINT;
edge->_intersectionVertex = v1*one_minus_r + v2*r;
}
return true;
}
// compute the intersection between line segment and surface
osg::Vec3 intersectionPoint(const osg::Vec3& v1, const osg::Vec3& v2)
{
double d1 = _plane.distance(v1);
double d2 = _plane.distance(v2);
double div = d2-d1;
if (div==0.0)
{
return v1;
}
double r = -d1 / div;
double one_minus_r = 1.0-r;
return v1*one_minus_r + v2*r;
}
// positive distance to the inside.
double distance(const osg::Vec3& v)
{
return _lowerOutside ? _plane.distance(v) : -_plane.distance(v) ;
}
};
struct ElevationIntersector
{
ElevationIntersector(TriangleIntersectOperator& tif, double elev, bool lowerOutside):
_tif(tif),
_elev(elev),
_lowerOutside(lowerOutside) {}
TriangleIntersectOperator& _tif;
double _elev;
bool _lowerOutside;
inline bool operator() (TriangleIntersectOperator::Edge* edge)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
osg::Vec3& v1 = _tif._originalVertices[edge->_p1];
osg::Vec3& v2 = _tif._originalVertices[edge->_p2];
double length_xy1 = sqrt(v1.x()*v1.x() + v1.y()*v1.y());
double elev1 = atan2((double)v1.z(),length_xy1);
double length_xy2 = sqrt(v2.x()*v2.x() + v2.y()*v2.y());
double elev2 = atan2((double)v2.z(),length_xy2);
edge->_p1Outside = _lowerOutside ? (elev1<_elev) : (elev1>_elev);
edge->_p2Outside = _lowerOutside ? (elev2<_elev) : (elev2>_elev);
// if both points inside then disgard
if (elev1<_elev && elev2<_elev) return false;
// if both points outside then disgard
if (elev1>_elev && elev2>_elev) return false;
if (elev1==_elev)
{
if (elev2==_elev)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::BOTH_ENDS;
}
else
{
edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_1;
}
}
else if (elev2==_elev)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_2;
}
else
{
double dx = v2.x()-v1.x();
double dy = v2.y()-v1.y();
double dz = v2.z()-v1.z();
double t = tan(_elev);
double tt = t*t;
double a = dz*dz-tt*(dx*dx + dy*dy);
double b = 2.0*(v1.z()*dz - tt*(v1.x()*dx + v1.y()*dy));
double c = v1.z()*v1.z() - tt*(v1.x()*v1.x() + v1.y()*v1.y());
double s1, s2;
if (!computeQuadraticSolution(a,b,c,s1,s2))
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
return false;
}
double r = 0.0;
if (s1>=0.0 && s1<=1.0)
{
r = s1;
}
else if (s2>=0.0 && s2<=1.0)
{
r = s2;
}
else
{
osg::notify(osg::INFO)<<"neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
return false;
}
double one_minus_r = 1.0-r;
edge->_intersectionType = TriangleIntersectOperator::Edge::MID_POINT;
edge->_intersectionVertex = v1*one_minus_r + v2*r;
}
return true;
}
// compute the intersection between line segment and surface
osg::Vec3 intersectionPoint(const osg::Vec3& v1, const osg::Vec3& v2)
{
double dx = v2.x()-v1.x();
double dy = v2.y()-v1.y();
double dz = v2.z()-v1.z();
double t = tan(_elev);
double tt = t*t;
double a = dz*dz-tt*(dx*dx + dy*dy);
double b = 2.0*(v1.z()*dz - tt*(v1.x()*dx + v1.y()*dy));
double c = v1.z()*v1.z() - tt*(v1.x()*v1.x() + v1.y()*v1.y());
double s1, s2;
if (!computeQuadraticSolution(a,b,c,s1,s2))
{
osg::notify(osg::INFO)<<"Warning::neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
return v1;
}
double r = 0.0;
if (s1>=0.0 && s1<=1.0)
{
r = s1;
}
else if (s2>=0.0 && s2<=1.0)
{
r = s2;
}
else
{
osg::notify(osg::INFO)<<"Warning::neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
return v1;
}
double one_minus_r = 1.0-r;
return v1*one_minus_r + v2*r;
}
// positive distance to the inside.
double distance(const osg::Vec3& v)
{
double length_xy = sqrt(v.x()*v.x() + v.y()*v.y());
double computedElev = atan2((double)v.z(),length_xy);
return _lowerOutside ? computedElev-_elev : _elev-computedElev ;
}
};
struct RadiusIntersector
{
RadiusIntersector(TriangleIntersectOperator& tif):
_tif(tif) {}
TriangleIntersectOperator& _tif;
inline bool operator() (TriangleIntersectOperator::Edge* edge)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
osg::Vec3& v1 = _tif._originalVertices[edge->_p1];
osg::Vec3& v2 = _tif._originalVertices[edge->_p2];
double radius1 = v1.length();
double radius2 = v2.length();
edge->_p1Outside = radius1>_tif._radius;
edge->_p2Outside = radius2>_tif._radius;
// if both points inside then disgard
if (radius1<_tif._radius && radius2<_tif._radius) return false;
// if both points outside then disgard
if (radius1>_tif._radius && radius2>_tif._radius) return false;
if (radius1==_tif._radius)
{
if (radius2==_tif._radius)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::BOTH_ENDS;
}
else
{
edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_1;
}
}
else if (radius2==_tif._radius)
{
edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_2;
}
else
{
double dx = v2.x()-v1.x();
double dy = v2.y()-v1.y();
double dz = v2.z()-v1.z();
double a = dx*dx + dy*dy + dz*dz;
double b = 2.0*(v1.x()*dx + v1.y()*dy + v1.z()*dz);
double c = v1.x()*v1.x() + v1.y()*v1.y() + v1.z()*v1.z() - _tif._radius*_tif._radius;
double s1, s2;
if (!computeQuadraticSolution(a,b,c,s1,s2))
{
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
return false;
}
double r = 0.0;
if (s1>=0.0 && s1<=1.0)
{
r = s1;
}
else if (s2>=0.0 && s2<=1.0)
{
r = s2;
}
else
{
osg::notify(osg::INFO)<<"neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
return false;
}
double one_minus_r = 1.0-r;
edge->_intersectionType = TriangleIntersectOperator::Edge::MID_POINT;
edge->_intersectionVertex = v1*one_minus_r + v2*r;
}
return true;
}
// compute the intersection between line segment and surface
osg::Vec3 intersectionPoint(const osg::Vec3& v1, const osg::Vec3& v2)
{
double dx = v2.x()-v1.x();
double dy = v2.y()-v1.y();
double dz = v2.z()-v1.z();
double a = dx*dx + dy*dy + dz*dz;
double b = 2.0*(v1.x()*dx + v1.y()*dy + v1.z()*dz);
double c = v1.x()*v1.x() + v1.y()*v1.y() + v1.z()*v1.z() - _tif._radius*_tif._radius;
double s1, s2;
if (!computeQuadraticSolution(a,b,c,s1,s2))
{
osg::notify(osg::INFO)<<"Warning: neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
return v1;
}
double r = 0.0;
if (s1>=0.0 && s1<=1.0)
{
r = s1;
}
else if (s2>=0.0 && s2<=1.0)
{
r = s2;
}
else
{
osg::notify(osg::INFO)<<"Warning: neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
return v1;
}
double one_minus_r = 1.0-r;
return v1*one_minus_r + v2*r;
}
// positive distance to the inside.
double distance(const osg::Vec3& v)
{
return _tif._radius-v.length();
}
};
}
using namespace SphereSegmentIntersector;
SphereSegment::LineList SphereSegment::computeIntersection(const osg::Matrixd& matrix, osg::Drawable* drawable)
{
// cast to Geometry, return empty handed if Drawable not a Geometry.
osg::Geometry* geometry = dynamic_cast<osg::Geometry*>(drawable);
if (!geometry) return SphereSegment::LineList();
// get vertices from geometry, return empty handed if a Vec3Array not present.
osg::Vec3Array* vertices = dynamic_cast<osg::Vec3Array*>(geometry->getVertexArray());
if (!vertices) return SphereSegment::LineList();
typedef osg::TriangleIndexFunctor<TriangleIntersectOperator> TriangleIntersectFunctor;
TriangleIntersectFunctor tif;
tif._centre = _centre;
tif._radius = _radius;
tif._azMin = _azMin;
tif._azMax = _azMax;
tif._elevMin = _elevMin;
tif._elevMax = _elevMax;
tif.computePositionAndRegions(matrix, *vertices);
// traverse the triangles in the Geometry dedicating intersections
geometry->accept(tif);
osg::notify(osg::INFO)<<"_numOutside = "<<tif._numOutside<<std::endl;
osg::notify(osg::INFO)<<"_numInside = "<<tif._numInside<<std::endl;
osg::notify(osg::INFO)<<"_numIntersecting = "<<tif._numIntersecting<<std::endl;
tif.removeDuplicateVertices();
tif.removeDuplicateTriangles();
tif.buildEdges();
RadiusIntersector radiusIntersector(tif);
AzimPlaneIntersector azMinIntersector(tif,_azMin, true);
AzimPlaneIntersector azMinEndIntersector(tif,_azMin-osg::PI*0.5, true);
AzimPlaneIntersector azMaxIntersector(tif,_azMax, false);
AzimPlaneIntersector azMaxEndIntersector(tif,_azMax-osg::PI*0.5, true);
ElevationIntersector elevMinIntersector(tif,_elevMin, true);
ElevationIntersector elevMaxIntersector(tif,_elevMax, false);
// create the line intersections with the terrain
SphereSegment::LineList radiusLines = tif.computeIntersections(radiusIntersector);
SphereSegment::LineList elevMinLines = tif.computeIntersections(elevMinIntersector);
SphereSegment::LineList elevMaxLines = tif.computeIntersections(elevMaxIntersector);
SphereSegment::LineList azMinLines;
SphereSegment::LineList azMaxLines;
double azimRange = _azMax-_azMin;
if (azimRange<2.0*osg::PI)
{
azMinLines = tif.computeIntersections(azMinIntersector);
azMaxLines = tif.computeIntersections(azMaxIntersector);
// trim the azimuth intersection lines by the radius
tif.trim(azMinLines,radiusIntersector);
tif.trim(azMaxLines,radiusIntersector);
// trim the azim intersection lines by the elevation
tif.trim(azMinLines, elevMinIntersector);
tif.trim(azMaxLines, elevMinIntersector);
// trim the centeral ends of the azim lines
tif.trim(azMinLines,azMinEndIntersector);
tif.trim(azMaxLines,azMaxEndIntersector);
if (azimRange<=osg::PI)
{
// trim the radius and elevation intersection lines by the azimMin
tif.trim(radiusLines, azMinIntersector);
tif.trim(elevMinLines, azMinIntersector);
tif.trim(elevMaxLines, azMinIntersector);
// trim the radius and elevation intersection lines by the azimMax
tif.trim(radiusLines, azMaxIntersector);
tif.trim(elevMinLines, azMaxIntersector);
tif.trim(elevMaxLines, azMaxIntersector);
}
else
{
// need to have new intersector which handles convex azim planes
tif.trim(radiusLines, azMinIntersector, azMaxIntersector);
tif.trim(elevMinLines, azMinIntersector, azMaxIntersector);
tif.trim(elevMaxLines, azMinIntersector, azMaxIntersector);
}
}
// trim elevation intersection lines by radius
tif.trim(elevMinLines,radiusIntersector);
tif.trim(elevMaxLines,radiusIntersector);
// trim the radius and elevation intersection lines by the elevMin
tif.trim(radiusLines, elevMinIntersector);
// trim the radius and elevation intersection lines by the elevMax
tif.trim(radiusLines, elevMaxIntersector);
// collect all lines together.
tif._generatedLines.insert(tif._generatedLines.end(), radiusLines.begin(), radiusLines.end());
tif._generatedLines.insert(tif._generatedLines.end(), azMinLines.begin(), azMinLines.end());
tif._generatedLines.insert(tif._generatedLines.end(), azMaxLines.begin(), azMaxLines.end());
tif._generatedLines.insert(tif._generatedLines.end(), elevMinLines.begin(), elevMinLines.end());
tif._generatedLines.insert(tif._generatedLines.end(), elevMaxLines.begin(), elevMaxLines.end());
osg::notify(osg::INFO)<<"number of seperate lines = "<<tif._generatedLines.size()<<std::endl;
float fuseDistance = 1.0;
tif.joinEnds(fuseDistance, true, true);
osg::notify(osg::INFO)<<"number of seperate lines after fuse = "<<tif._generatedLines.size()<<std::endl;
float joinDistance = 1e8;
tif.joinEnds(joinDistance, false, false);
osg::notify(osg::INFO)<<"number of seperate lines after join = "<<tif._generatedLines.size()<<std::endl;
tif.joinEnds(joinDistance, false, true);
osg::notify(osg::INFO)<<"number of seperate lines after second join = "<<tif._generatedLines.size()<<std::endl;
return tif._generatedLines;
}
osg::Node* SphereSegment::computeIntersectionSubgraph(const osg::Matrixd& matrix, osg::Drawable* drawable)
{
SphereSegment::LineList generatedLines = computeIntersection(matrix, drawable);
osg::Geode* geode = new osg::Geode;
geode->getOrCreateStateSet()->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
for(osgSim::SphereSegment::LineList::iterator itr = generatedLines.begin();
itr != generatedLines.end();
++itr)
{
osg::Geometry* geom = new osg::Geometry;
geode->addDrawable(geom);
osg::Vec3Array* vertices = itr->get();
geom->setVertexArray(vertices);
geom->addPrimitiveSet(new osg::DrawArrays(GL_LINE_STRIP, 0, vertices->getNumElements()));
}
#if 0
float radius = 0.1f;
for(unsigned int i=0; i<tif._originalVertices.size(); ++i)
{
osg::ShapeDrawable* sd = new osg::ShapeDrawable(new osg::Sphere(tif._originalVertices[i]+tif._centre,radius));
TriangleIntersectFunctor::RegionCounter rc;
rc.add(tif._regions[i]);
TriangleIntersectFunctor::Region::Classification region = rc.overallClassification();
if (region==TriangleIntersectFunctor::Region::OUTSIDE)
{
sd->setColor(osg::Vec4(1.0,0.0,0.0,1.0));
}
else if (region==TriangleIntersectFunctor::Region::INSIDE)
{
sd->setColor(osg::Vec4(1.0,1.0,0.0,1.0));
}
else if (region==TriangleIntersectFunctor::Region::INTERSECTS)
{
sd->setColor(osg::Vec4(1.0,1.0,1.0,1.0));
}
geode->addDrawable(sd);
}
#endif
return geode;
}
Reference in New Issue View Git Blame Copy Permalink