From 598dda59f5bcfea2e6c58a60b758a94fc1938438 Mon Sep 17 00:00:00 2001 From: Robert Osfield Date: Tue, 25 Oct 2005 13:28:48 +0000 Subject: [PATCH] From Geoff Michel, added support for constrain delaunay triangultion, and osgdelaunay example. --- Make/makedirdefs | 1 + VisualStudio/VisualStudio.dsw | 30 + .../examples/osgdelaunay/osgdelaunay.dsp | 101 ++ examples/osgdelaunay/GNUmakefile | 19 + examples/osgdelaunay/GNUmakefile.inst | 13 + examples/osgdelaunay/osgdelaunay.cpp | 1378 +++++++++++++++++ include/osgUtil/DelaunayTriangulator | 110 +- src/osgUtil/DelaunayTriangulator.cpp | 906 ++++++++++- 8 files changed, 2531 insertions(+), 27 deletions(-) create mode 100644 VisualStudio/examples/osgdelaunay/osgdelaunay.dsp create mode 100644 examples/osgdelaunay/GNUmakefile create mode 100644 examples/osgdelaunay/GNUmakefile.inst create mode 100644 examples/osgdelaunay/osgdelaunay.cpp diff --git a/Make/makedirdefs b/Make/makedirdefs index c5e25b046..2aa13b0fd 100644 --- a/Make/makedirdefs +++ b/Make/makedirdefs @@ -199,6 +199,7 @@ ifeq ($(PRODUCER_INSTALLED),yes) osgcluster\ osgcopy\ osgcubemap\ + osgdelaunay\ osgdepthshadow\ osgdepthpartition\ osgdistortion\ diff --git a/VisualStudio/VisualStudio.dsw b/VisualStudio/VisualStudio.dsw index 5a3cea76d..c469d0f6d 100644 --- a/VisualStudio/VisualStudio.dsw +++ b/VisualStudio/VisualStudio.dsw @@ -2175,6 +2175,36 @@ Package=<4> ############################################################################### +Project: "Example osgdelaunay"=.\examples\osgdelaunay\osgdelaunay.dsp - Package Owner=<4> + +Package=<5> +{{{ +}}} + +Package=<4> +{{{ + Begin Project Dependency + Project_Dep_Name Core osg + End Project Dependency + Begin Project Dependency + Project_Dep_Name Core osgDB + End Project Dependency + Begin Project Dependency + Project_Dep_Name Core osgGA + End Project Dependency + Begin Project Dependency + Project_Dep_Name Core osgProducer + End Project Dependency + Begin Project Dependency + Project_Dep_Name Core osgUtil + End Project Dependency + Begin Project Dependency + Project_Dep_Name Core osgText + End Project Dependency +}}} + +############################################################################### + Project: "Example osgtesselate"=.\examples\osgtesselate\osgtesselate.dsp - Package Owner=<4> Package=<5> diff --git a/VisualStudio/examples/osgdelaunay/osgdelaunay.dsp b/VisualStudio/examples/osgdelaunay/osgdelaunay.dsp new file mode 100644 index 000000000..3c3047571 --- /dev/null +++ b/VisualStudio/examples/osgdelaunay/osgdelaunay.dsp @@ -0,0 +1,101 @@ +# Microsoft Developer Studio Project File - Name="Example osgdelaunay" - Package Owner=<4> +# Microsoft Developer Studio Generated Build File, Format Version 6.00 +# ** DO NOT EDIT ** + +# TARGTYPE "Win32 (x86) Console Application" 0x0103 + +CFG=Example osgdelaunay - Win32 Release +!MESSAGE This is not a valid makefile. To build this project using NMAKE, +!MESSAGE use the Export Makefile command and run +!MESSAGE +!MESSAGE NMAKE /f "osgdelaunay.mak". +!MESSAGE +!MESSAGE You can specify a configuration when running NMAKE +!MESSAGE by defining the macro CFG on the command line. For example: +!MESSAGE +!MESSAGE NMAKE /f "osgdelaunay.mak" CFG="Example osgdelaunay - Win32 Release" +!MESSAGE +!MESSAGE Possible choices for configuration are: +!MESSAGE +!MESSAGE "Example osgdelaunay - Win32 Release" (based on "Win32 (x86) Console Application") +!MESSAGE "Example osgdelaunay - Win32 Debug" (based on "Win32 (x86) Console Application") +!MESSAGE + +# Begin Project +# PROP AllowPerConfigDependencies 0 +# PROP Scc_ProjName "" +# PROP Scc_LocalPath "" +CPP=cl.exe +RSC=rc.exe + +!IF "$(CFG)" == "Example osgdelaunay - Win32 Release" + +# PROP BASE Use_MFC 0 +# PROP BASE Use_Debug_Libraries 0 +# PROP BASE Output_Dir "Release" +# PROP BASE Intermediate_Dir "Release" +# PROP BASE Target_Dir "" +# PROP Use_MFC 0 +# PROP Use_Debug_Libraries 0 +# PROP Output_Dir "Release" +# PROP Intermediate_Dir "Release" +# PROP Ignore_Export_Lib 0 +# PROP Target_Dir "" +MTL=midl.exe +# ADD BASE CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c +# ADD CPP /nologo /MD /W3 /GR /GX /O2 /I "../../../include" /I "../../../../OpenThreads/include" /I "../../../../Producer/include" /I "../../../../3rdParty/include" /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /Zm200 /c +# ADD BASE RSC /l 0x809 /d "NDEBUG" +# ADD RSC /l 0x809 /d "NDEBUG" +BSC32=bscmake.exe +# ADD BASE BSC32 /nologo +# ADD BSC32 /nologo +LINK32=link.exe +# ADD BASE LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /machine:I386 +# ADD LINK32 OpenThreadsWin32.lib opengl32.lib /nologo /subsystem:console /pdb:none /machine:I386 /out:"../../../bin/osgdelaunay.exe" /libpath:"../../../lib" /libpath:"../../../../OpenThreads/lib/win32" /libpath:"../../../../Producer/lib" /libpath:"../../../../3rdParty/lib" + +!ELSEIF "$(CFG)" == "Example osgdelaunay - Win32 Debug" + +# PROP BASE Use_MFC 0 +# PROP BASE Use_Debug_Libraries 1 +# PROP BASE Output_Dir "Debug" +# PROP BASE Intermediate_Dir "Debug" +# PROP BASE Target_Dir "" +# PROP Use_MFC 0 +# PROP Use_Debug_Libraries 1 +# PROP Output_Dir "Debug" +# PROP Intermediate_Dir "Debug" +# PROP Ignore_Export_Lib 0 +# PROP Target_Dir "" +MTL=midl.exe +# ADD BASE CPP /nologo /W3 /Gm /GX /Zi /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c +# ADD CPP /nologo /MDd /W3 /Gm /GR /GX /Zi /Od /I "../../../include" /I "../../../../OpenThreads/include" /I "../../../../Producer/include" /I "../../../../3rdParty/include" /D "_CONSOLE" /D "_MBCS" /D "FL_DLL" /D "WIN32" /D "_DEBUG" /FR /YX /FD /Zm200 /c +# ADD BASE RSC /l 0x809 /d "_DEBUG" +# ADD RSC /l 0x809 /d "_DEBUG" +BSC32=bscmake.exe +# ADD BASE BSC32 /nologo +# ADD BSC32 /nologo +LINK32=link.exe +# ADD BASE LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept +# ADD LINK32 OpenThreadsWin32d.lib opengl32.lib /nologo /subsystem:console /debug /machine:I386 /nodefaultlib:"libcmt" /out:"../../../bin/osgdelaunayd.exe" /pdbtype:sept /libpath:"../../../lib" /libpath:"../../../../OpenThreads/lib/win32" /libpath:"../../../../Producer/lib" /libpath:"../../../../3rdParty/lib" +# SUBTRACT LINK32 /incremental:no + +!ENDIF + +# Begin Target + +# Name "Example osgdelaunay - Win32 Release" +# Name "Example osgdelaunay - Win32 Debug" +# Begin Source File + +SOURCE=..\..\..\examples\osgdelaunay\osgdelaunay.cpp +# End Source File +# Begin Source File + +SOURCE=..\..\icons\osg_icon.rc +# End Source File +# End Target +# Begin Group "Resource Files" + +# PROP Default_Filter "ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe" +# End Group +# End Project diff --git a/examples/osgdelaunay/GNUmakefile b/examples/osgdelaunay/GNUmakefile new file mode 100644 index 000000000..8ec42afe2 --- /dev/null +++ b/examples/osgdelaunay/GNUmakefile @@ -0,0 +1,19 @@ +TOPDIR = ../.. +include $(TOPDIR)/Make/makedefs + +CXXFILES =\ + osgdelaunay.cpp\ + +LIBS += -losgProducer -lProducer -losgText -losgGA -losgDB -losgUtil -losg $(GL_LIBS) $(X_LIBS) $(OTHER_LIBS) + +INSTFILES = \ + $(CXXFILES)\ + GNUmakefile.inst=GNUmakefile + +EXEC = osgdelaunay + +INC += $(X_INC) + +include $(TOPDIR)/Make/makerules + + diff --git a/examples/osgdelaunay/GNUmakefile.inst b/examples/osgdelaunay/GNUmakefile.inst new file mode 100644 index 000000000..1ae517488 --- /dev/null +++ b/examples/osgdelaunay/GNUmakefile.inst @@ -0,0 +1,13 @@ +TOPDIR = ../.. +include $(TOPDIR)/Make/makedefs + +CXXFILES =\ + osgdelaunay.cpp\ + +LIBS += -losgProducer -lProducer -losgDB -losgText -losgUtil -losg $(GL_LIBS) $(X_LIBS) $(OTHER_LIBS) + +EXEC = osgdelaunay + +INC += $(X_INC) + +include $(TOPDIR)/Make/makerules diff --git a/examples/osgdelaunay/osgdelaunay.cpp b/examples/osgdelaunay/osgdelaunay.cpp new file mode 100644 index 000000000..83b53fca7 --- /dev/null +++ b/examples/osgdelaunay/osgdelaunay.cpp @@ -0,0 +1,1378 @@ +/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2003 Robert Osfield +* +* This application is open source and may be redistributed and/or modified +* freely and without restriction, both in commericial and non commericial applications, +* as long as this copyright notice is maintained. +* +* This application 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. +*/ + +/** Example of use of delaunay triangulator with constraints. +* this could be a method of generating terrains, a constraint forces certain edges to +* exist in the triangulation. +*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // tesselator triangulates the constrained triangles + +#include +#include + + +/** here are 2 common types of constraint +* Area - forces an area to be filled; replacement geometry is a canopy and optional wall +* Linear - constructs a closed loop of constant width around a line. +*/ +class WallConstraint: public osgUtil::DelaunayConstraint { // forces lines to eb edge + // wall constraint - can generate a wall at the coordinates of the constraint +public: +/** if you derive a class from DelaunayConstraint then you can create +* a specific geometry creation routine. + */ + WallConstraint() : height(0), txxrepWall(10), txyrepWall(10) { } + + /** or create a wall around the constraint area: */ + virtual osg::Geometry * makeWallGeometry(const osg::Vec3Array *points) const; + + /** for basic purposes, you can call these routines to make simple fill in geometries */ + virtual osg::DrawArrays* makeWall(void ) const { // build a wall height high around the constraint + const osg::Vec3Array *_line= dynamic_cast(getVertexArray()); + return (new osg::DrawArrays(osg::PrimitiveSet::QUAD_STRIP,0,2*_line->size())); + } + + + virtual osg::Vec3Array *getWall(const osg::Vec3Array *points,const float height) const; + virtual osg::Vec2Array *getWallTexcoords(const osg::Vec3Array *points,const float height) const; + virtual osg::Vec3Array *getWallNormals(const osg::Vec3Array *points) const { + osg::ref_ptr nrms=new osg::Vec3Array; + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprgetMode()==osg::PrimitiveSet::LINE_LOOP || + prset->getMode()==osg::PrimitiveSet::LINE_STRIP) { // loops and walls + // start with the last point on the loop + osg::Vec3 prevp=(*vertices)[prset->index (prset->getNumIndices()-1)]; + for (unsigned int i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + osg::Vec3 nrm=(curp-prevp)^osg::Vec3(0,0,1); + nrm.normalize(); + nrms->push_back(nrm); + nrms->push_back(nrm); + prevp=curp; + } + const osg::Vec3 curp=(*vertices)[prset->index (0)]; + osg::Vec3 nrm=(curp-prevp)^osg::Vec3(0,0,1); + nrm.normalize(); + nrms->push_back(nrm); + nrms->push_back(nrm); + } + } + return nrms.release(); + } + + + + // geometry creation parameters + void setWallTexrep(const float w,const float h) { txxrepWall=w;txyrepWall=h;} + + /** Wall Geometry will return with this texture applied: */ + void setTexture(const char *tx) { texture=tx;} + /** fence/wall height */ + void setHeight(const float h) { height=h;} +protected: + float height; + std::string texture; + float txxrepWall, txyrepWall; +}; +class ArealConstraint: public osgUtil::DelaunayConstraint { // forces edges of an area to fit triangles + // areal constraint - general nonuniform field, forest, lake etc. +public: +/** if you derive a class from DelaunayConstraint then you can create +* a specific geometry creation routine. + */ + ArealConstraint() : txxrepArea(10), txyrepArea(10),txxrepWall(10), txyrepWall(10) { } + + /** return a geometry that fills the constraint. + */ + virtual osg::Geometry * makeAreal( osg::Vec3Array *points); + + /** or create a wall around the constraint area: */ + virtual osg::Geometry * makeWallGeometry( osg::Vec3Array *points) ; + + /** for basic purposes, you can call these routines to make simple fill in geometries */ + virtual osg::DrawArrays* makeWall(void ) const; + virtual osg::Vec3Array *getWall(const osg::Vec3Array *points,const float height) const; + virtual osg::Vec2Array *getWallTexcoords(const osg::Vec3Array *points,const float height) const; + virtual osg::Vec3Array *getWallNormals(const osg::Vec3Array *points) const; + /** Canopies are the same triangles as the terrain but offset by height above + * (height might be 0). */ + virtual osg::DrawArrays* makeCanopy(void ) const; + virtual osg::Vec3Array *getCanopy(const osg::Vec3Array *points,const float height) const; + virtual osg::Vec2Array *getCanopyTexcoords(const osg::Vec3Array *points) const; + virtual osg::Vec3Array *getCanopyNormals(const osg::Vec3Array *points) const; + + // geometry creation parameters + void setTexrep(const float w,const float h) { txxrepArea=w;txyrepArea=h;} + void setWallTexrep(const float w,const float h) { txxrepWall=w;txyrepWall=h;} + /** Geometry will return with this texture applied: */ + void setWallTexture(const char *tx) { walltexture=tx;} + /** Geometry will return with this texture applied: */ + void setTexture(const char *tx) { texture=tx;} + /** fence/wall height */ + void setHeight(const float h) { height=h;} + std::string walltexture; +protected: + float height; + std::string texture; + float txxrepArea, txyrepArea; + float txxrepWall, txyrepWall; +}; + +class LinearConstraint: public osgUtil::DelaunayConstraint { +/** forces edges of a "road" to fit triangles +* if 2 roads cross, then the overlap will be replaced by a 'cross road' + * and the roads built up to the cross roads with a texture along its length. */ +public: + LinearConstraint() : DelaunayConstraint(), txxrepAlong(10), txyrepAcross(10), width(2) { } + + /** geometry creation parameters */ + /* Width of linear feature (eg road, railway) */ + void setWidth(const float w) { width=w;} + + /** Texture repeat distance across linear (often equal to width) and along its length */ + virtual void setTexrep(const float w,const float h) { txyrepAcross=h;txxrepAlong=w; } + + /** generate constant width around line - creates the area to be cut into the terrain. */ + virtual void setVertices( osg::Vec3Array *lp, const float width); + + /** return a geometry that fills the constraint. + */ + virtual osg::Geometry *makeGeometry(const osg::Vec3Array *points) ; + + /** return normals array - flat shaded */ + osg::Vec3Array* getNormals(const osg::Vec3Array *points); + + /** Roads apply a texture proportional to length along the road line. */ + virtual osg::DrawArrays* makeRoad( ) const; + virtual osg::Vec3Array *getRoadVertices(const osg::Vec3Array *points) const; + virtual osg::Vec2Array *getRoadTexcoords(const osg::Vec3Array *points) ; + + virtual osg::Vec3Array *getRoadNormals(const osg::Vec3Array *points) const; + /** Geometry will return with this texture applied: */ + void setTexture(const char *tx) { texture=tx;} + +protected: + osg::ref_ptr _tcoords; + osg::ref_ptr _edgecoords; + float txxrepAlong, txyrepAcross; + std::string texture; + float width; // width of a linear feature + osg::ref_ptr _midline; // defines the midline of a road, rail etc. +}; + +/** a specific type of constaint - that replaces an area with a pyramid */ + +class pyramid : public osgUtil::DelaunayConstraint { +/** sample user constriant - creates hole in terrain to fit base of pyramid, and + * geometry of an Egyptian pyramid to fit the hole. */ +public: + pyramid() : _side(100.) {} + + void setpos(const osg::Vec3 p, const float size) { _pos=p;_side=size;} + + virtual osg::Geometry * makeGeometry(const osg::Vec3Array *points) const + { + // create pyramid geometry. Centre plus points around base + const osg::Vec3Array *_line= dynamic_cast(getVertexArray()); + osg::Geometry *gm=new osg::Geometry; + osg::Vec3Array *pts=new osg::Vec3Array; + osg::Vec3Array *norms=new osg::Vec3Array; + osg::Vec2Array *tcoords=new osg::Vec2Array; + int ip; + + pts->push_back(_pos+osg::Vec3(0,0,_side)*0.5); + for (ip=0; ip<4; ip++) { + pts->push_back((*_line)[ip]); + } + for (ip=1; ip<5; ip++) { + osg::Vec3 nrm=((*pts)[ip]-(*pts)[0])^((*pts)[ip==4?0:ip+1]-(*pts)[ip]); + nrm.normalize( ); + norms->push_back(nrm); + } + + gm->setNormalBinding(osg::Geometry::BIND_PER_PRIMITIVE); + gm->setVertexArray(pts); + osg::StateSet *dstate= gm->getOrCreateStateSet( ); + dstate->setMode( GL_LIGHTING, osg::StateAttribute::ON ); + + osg::Image* image = osgDB::readImageFile("stoneWall.jpg"); + if (image) + { + osg::Texture2D* txt = new osg::Texture2D; + txt->setImage(image); + txt->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT ); + txt->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::REPEAT ); + dstate->setTextureAttributeAndModes(0,txt,osg::StateAttribute::ON); + } + gm->setNormalArray(norms); + //// gm->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::TRIANGLE_FAN,0,6)); + osg::DrawElementsUInt *dui=new osg::DrawElementsUInt(GL_TRIANGLES); + for (ip=0; ip<4; ip++) { + dui->push_back(0); + dui->push_back(ip+1); + dui->push_back(ip==3?1:ip+2); + } + tcoords->push_back(osg::Vec2(2,4)); + tcoords->push_back(osg::Vec2(0,0)); + tcoords->push_back(osg::Vec2(4,0)); + tcoords->push_back(osg::Vec2(0,0)); + tcoords->push_back(osg::Vec2(4,0)); + gm->setTexCoordArray(0,tcoords); + gm->addPrimitiveSet(dui); + return gm; + } + virtual void calcVertices( void) { // must have a position first + osg::Vec3Array *edges=new osg::Vec3Array; + osg::Vec3 valong; + edges->push_back(_pos+osg::Vec3(0.5,0.5,0)*_side); + edges->push_back(_pos+osg::Vec3(-0.5,0.5,0)*_side); + edges->push_back(_pos+osg::Vec3(-0.5,-0.5,0)*_side); + edges->push_back(_pos+osg::Vec3(0.5,-0.5,0)*_side); + setVertexArray(edges); + } +private: + osg::Vec3 _pos; // where the pyramid is + float _side ; // length of side +}; + +float getheight(const float x, const float y) +{ // returns the x,y,height of terrain + return 150*sin(x*.0020)*cos(y*.0020); +} +osg::Vec3d getpt(const int np) +{ // returns the x,y,height of terrain up to maxp^2 points + static int maxp =40; + int i=np/maxp; + int j=np%maxp; + // make the random scale 0.00 if you want an equispaced XY grid. + float x=3000.0/(maxp-1)*i+0.00052*rand(); + float y=3000.0/(maxp-1)*j+0.0005*rand(); + float z=getheight(x,y); + if (np>=maxp*maxp) z=-1.e32; + return osg::Vec3d(x,y,z); +} +osg::Node* createHUD(const int ndcs,std::string what) +{ // add a string reporting the type of winding rule tesselation applied + osg::Geode* geode = new osg::Geode(); + + std::string timesFont("fonts/arial.ttf"); + + // turn lighting off for the text and disable depth test to ensure its always ontop. + osg::StateSet* stateset = geode->getOrCreateStateSet(); + stateset->setMode(GL_LIGHTING,osg::StateAttribute::OFF); + + // Disable depth test, and make sure that the hud is drawn after everything + // else so that it always appears ontop. + stateset->setMode(GL_DEPTH_TEST,osg::StateAttribute::OFF); + stateset->setRenderBinDetails(11,"RenderBin"); + + osg::Vec3 position(50.0f,900.0f,0.0f); + osg::Vec3 delta(0.0f,-35.0f,0.0f); + + { + osgText::Text* text = new osgText::Text; + geode->addDrawable( text ); + std::ostringstream cue; + cue<<"Delaunay triangulation with constraints level "<setFont(timesFont); + text->setPosition(position); + text->setText(cue.str()); + text->setColor(osg::Vec4(1.0,1.0,0.8,1.0)); + position += delta*(ndcs+2); + + text = new osgText::Text; + geode->addDrawable( text ); + + text->setFont(timesFont); + text->setPosition(position); + text->setText("(use 'W' wireframe & 'T' texture to visualise mesh)"); + text->setColor(osg::Vec4(1.0,1.0,0.8,1.0)); + position += delta; + + } + { + osgText::Text* text = new osgText::Text; + geode->addDrawable( text ); + + text->setFont(timesFont); + text->setPosition(position); + text->setText("Press 'n' to add another constraint."); + + } + + // create the hud. + osg::MatrixTransform* modelview_abs = new osg::MatrixTransform; + modelview_abs->setReferenceFrame(osg::Transform::ABSOLUTE_RF); + modelview_abs->setMatrix(osg::Matrix::identity()); + modelview_abs->addChild(geode); + + osg::Projection* projection = new osg::Projection; + projection->setMatrix(osg::Matrix::ortho2D(0,1280,0,1024)); + projection->addChild(modelview_abs); + + return projection; + +} +osg::Group *makedelaunay(const int ndcs) +{ // create a terrain tile. This is just an example! + // ndcs is the number of delaunay constraints to be applied + osg::ref_ptr grp=new osg::Group; + osg::ref_ptr geode=new osg::Geode; + osg::ref_ptr trig=new osgUtil::DelaunayTriangulator(); + osg::StateSet *stateset=geode->getOrCreateStateSet(); + + osg::Vec3Array *points=new osg::Vec3Array; + + osg::Image* image = osgDB::readImageFile("Lands-Flowers-Small_Purple.png"); + if (image) + { + osg::Texture2D* texture = new osg::Texture2D; + texture->setImage(image); + texture->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT ); + texture->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::REPEAT ); + stateset->setTextureAttributeAndModes(0,texture,osg::StateAttribute::ON); + } + + geode->setStateSet( stateset ); + unsigned int i; + + int eod=0; + while (eod>=0) { + osg::Vec3d pos=getpt(eod); + if (pos.z()>-10000) { + points->push_back(pos); + eod++; + } else { + eod=-9999; + } + } + std::vector < pyramid* > pyrlist; + osg::ref_ptr wc; // This example does not remove the interior + osg::ref_ptr dc2; + osg::ref_ptr forest; + osg::ref_ptr dc3; + osg::ref_ptr dc6; + osg::ref_ptr dc6a; + osg::ref_ptr dc8; + osg::ref_ptr forestroad; + osg::ref_ptr forestroad2; + osg::ref_ptr forestroad3; + osg::ref_ptr dc; + std::ostringstream what; + if (1==1) { // add a simple constraint of few points + osg::ref_ptr dc=new osgUtil::DelaunayConstraint; + osg::Vec3Array *bounds=new osg::Vec3Array; + unsigned int nmax=4; + for (i=0 ; ipush_back(osg::Vec3(x,y,getheight(x,y))); + } + dc->setVertexArray(bounds); + dc->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP,0,nmax) ); + + trig->addInputConstraint(dc.get()); + what << nmax << " point simple constraint\n"; + } + if (ndcs>0) { // add 5 pyramids + for (unsigned int ipy=0; ipy<5/*5*/; ipy++) { + osg::ref_ptr pyr=new pyramid; + float x=2210+ipy*120, y=1120+ipy*220; + pyr->setpos(osg::Vec3(x,y,getheight(x,y)),125.0+10*ipy); + pyr->calcVertices(); // make vertices + pyr->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,0,4) ); + trig->addInputConstraint(pyr.get()); + pyrlist.push_back(pyr.get()); + } + what << 5 << " pyramids\n"; + if (ndcs>1) { + // add a simple constraint feature - this can cut holes in the terrain or just leave the triangles + // with edges forced to the constraint. + dc=new osgUtil::DelaunayConstraint; + osg::Vec3Array *bounds=new osg::Vec3Array; + for (i=0 ; i<12; i++) { + float x=610.0+420*sin(i/3.0),y=610.0+420*cos(i/3.0); + bounds->push_back(osg::Vec3(x,y,getheight(x,y))); + } + dc->setVertexArray(bounds); + dc->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,0,12) ); + + trig->addInputConstraint(dc.get()); + what << 12 << " point closed loop"; + + if (ndcs>2) { + wc=new WallConstraint; // This example does not remove the interior + // eg to force terrain edges that are on ridges in the terrain etc. + // use wireframe to see the constrained edges. + // NB this is not necessarily a closed loop of edges. + // we do however build a wall at the coordinates. + bounds=new osg::Vec3Array; + for (i=0 ; i<5; i++) { + float x=1610.0+420*sin(i/1.0),y=1610.0+420*cos(i/1.0); + bounds->push_back(osg::Vec3(x,y,getheight(x,y))); + } + wc->setVertexArray(bounds); + wc->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP,0,5) ); + wc->setHeight(12.0); + trig->addInputConstraint(wc.get()); + what << " with interior removed\n"; + what << 5 << " point wall derived constraint\n"; + + if (ndcs>3) { + // add a removed area and replace it with a different texture + dc2=new ArealConstraint; + bounds=new osg::Vec3Array; + for (i=0 ; i<18; i++) { + float x=1610.0+420*sin(i/3.0),y=610.0+220*cos(i/3.0); + bounds->push_back(osg::Vec3(x,y,getheight(x,y))); + } + dc2->setVertexArray(bounds); + dc2->setTexrep(100,100); // texture is repeated at this frequency + dc2->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,0,18) ); + trig->addInputConstraint(dc2.get()); + what << 18 << " point area replaced\n"; + + if (ndcs>4) { + dc3=new LinearConstraint; + // a linear feature or 'road' + osg::Vec3Array *verts=new osg::Vec3Array; + for (i=0 ; i<32; i++) { + float x=610.0+50*i+90*sin(i/5.0),y=1110.0+90*cos(i/5.0); + verts->push_back(osg::Vec3(x,y,getheight(x,y))); + } + dc3->setVertices(verts,9.5); // width of road + for (osg::Vec3Array::iterator vit=points->begin(); vit!=points->end(); ) { + if (dc3->contains(*vit)) { + vit=points->erase(vit); + } else { + vit++; + } + } + trig->addInputConstraint(dc3.get()); + what << 32 << " point road constraint\n"; + if (ndcs>5) { + // add a removed area and replace it with a 'forest' with textured roof and walls + forest=new ArealConstraint; + bounds=new osg::Vec3Array; + for (i=0 ; i<12; i++) { + float x=610.0+420*sin(i/2.0),y=1810.0+420*cos(i/2.0); + bounds->push_back(osg::Vec3(x,y,getheight(x,y))); + } + forest->setVertexArray(bounds); + forest->setHeight(50); + forest->setWallTexrep(100,50); + forest->setTexrep(100,100); // texture is repeated at this frequency + forest->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,0,12) ); + if (ndcs==6) trig->addInputConstraint(forest.get()); + what << 12 << " point forest constraint\n"; + + if (ndcs>6) { // add roads that intersect forest + osg::ref_ptr forestplus=new osgUtil::DelaunayConstraint; + forestroad=new LinearConstraint; + verts=new osg::Vec3Array; + for (i=0 ; i<12; i++) { + int ip=(i-6)*(i-6); + float xp=410.0+20.0*ip; + float y=1210.0+150*i; + verts->push_back(osg::Vec3(xp,y,getheight(xp,y))); + } + forestroad->setVertices(verts,22); // add road + forestplus->merge(forestroad.get()); + forestroad2=new LinearConstraint; + verts=new osg::Vec3Array; + for (i=0 ; i<12; i++) { + int ip=(i-6)*(i-6); + float xp=810.0-10.0*ip; + float y=1010.0+150*i; + verts->push_back(osg::Vec3(xp,y,getheight(xp,y))); + } + forestroad2->setVertices(verts,22); // add road + forestplus->merge(forestroad2.get()); + forestroad3=new LinearConstraint; + verts=new osg::Vec3Array; + for (i=0 ; i<6; i++) { + int ip=(i-6)*(i-6); + float xp=210.0+140.0*i+ip*10.0; + float y=1510.0+150*i; + verts->push_back(osg::Vec3(xp,y,getheight(xp,y))); + } + forestroad3->setVertices(verts,22); // add road + forestplus->merge(forestroad3.get()); + forestplus->merge(forest.get()); + forestplus->handleOverlaps(); + for (osg::Vec3Array::iterator vit=points->begin(); vit!=points->end(); ) { + if (forestroad->contains(*vit)) { + vit=points->erase(vit); + } else if (forestroad2->contains(*vit)) { + vit=points->erase(vit); + } else if (forestroad3->contains(*vit)) { + vit=points->erase(vit); + } else { + vit++; + } + } + trig->addInputConstraint(forestplus.get()); + what << " roads intersect forest constraint\n"; + if (ndcs>7) { + // this option adds a more complex DC + // made of several (ok 2 - extend your own way) overlapping DC's + osg::ref_ptr dcoverlap=new osgUtil::DelaunayConstraint; + float x=1200; float y=1900; + { + verts=new osg::Vec3Array; + dc6=new LinearConstraint; + verts->push_back(osg::Vec3(x-180,y,getheight(x-180,y))); + verts->push_back(osg::Vec3(x+180,y,getheight(x+180,y))); + dc6->setVertices(verts,22); // width of road + dcoverlap->merge(dc6.get()); + } + { + dc6a= new LinearConstraint; + verts=new osg::Vec3Array; + verts->push_back(osg::Vec3(x,y-180,getheight(x,y-180))); + verts->push_back(osg::Vec3(x-20,y,getheight(x,y))); + verts->push_back(osg::Vec3(x,y+180,getheight(x,y+180))); + dc6a->setVertices(verts,22); // width of road + dcoverlap->merge(dc6a.get()); + } + what << "2 intersecting roads, with added points\n"; + if (ndcs>9) { + // add yet more roads + dc8= new LinearConstraint; + verts=new osg::Vec3Array; + float rad=60.0; + for (float theta=0; theta<4*osg::PI; theta+=0.1*osg::PI) { + float xp=x+rad*cos(theta), yp=y+rad*sin(theta); + verts->push_back(osg::Vec3(xp,yp,getheight(xp,yp))); + rad+=2.5; + } + dc8->setVertices(verts,16); // width of road + dcoverlap->merge(dc8.get()); + what << "Spiral road crosses several other constraints."; + } + dcoverlap->handleOverlaps(); + if (ndcs>8) { + // remove vertices cleans up the texturing at the intersection. + dcoverlap->removeVerticesInside(dc6.get()); + dcoverlap->removeVerticesInside(dc6a.get()); + if (dc8.valid()) dcoverlap->removeVerticesInside(dc8.get()); + what << " remove internal vertices to improve texturing."; + } + for (osg::Vec3Array::iterator vit=points->begin(); vit!=points->end(); ) { + if (dcoverlap->contains(*vit)) { + vit=points->erase(vit); + } else { + vit++; + } + } + trig->addInputConstraint(dcoverlap.get()); + } + } + } + } + } + } + } + } // ndcs>0 + trig->setInputPointArray(points); + + /** NB you need to supply a vec3 array for the triangulator to calculate normals into */ + osg::Vec3Array *norms=new osg::Vec3Array; + trig->setOutputNormalArray(norms); + + trig->triangulate(); + osg::notify(osg::WARN) << " End of trig\n " < gm=new osg::Geometry; + gm->setVertexArray(points); // points may have been modified in order by triangulation. + /** calculate texture coords for terrain points */ + if (image) { + float repeat=150.0, ry=150.0; // how often to repeat texture + osg::Vec2Array *tcoords=new osg::Vec2Array; + for (osg::Vec3Array::iterator itr=points->begin(); itr!=points->end(); itr++) { + osg::Vec2 tcatxy((*itr).x()/repeat,(*itr).y()/ry); + tcoords->push_back(tcatxy); + } + gm->setTexCoordArray(0,tcoords); + } + gm->addPrimitiveSet(trig->getTriangles()); + gm->setNormalArray(trig->getOutputNormalArray()); + gm->setNormalBinding(osg::Geometry::BIND_PER_PRIMITIVE); + geode->addDrawable(gm.get()); + if (ndcs>0) { + for ( std::vector < pyramid* >::iterator itr=pyrlist.begin(); itr!=pyrlist.end(); itr++) { + trig->removeInternalTriangles(*itr); + geode->addDrawable((*itr)->makeGeometry(points)); // this fills the holes of each pyramid with geometry + } + + if (ndcs>2) { + trig->removeInternalTriangles(dc.get()); + + wc->setTexture("Brick-Norman-Brown.TGA"); // wall looks like brick + geode->addDrawable(wc->makeWallGeometry(points)); // this creates wall at wc drawarrays + if (ndcs>3) { + trig->removeInternalTriangles(dc2.get()); + osg::ref_ptr arpts=dc2->getPoints(points); + dc2->setTexture("Lands-Needles_2.png"); + geode->addDrawable(dc2->makeAreal(arpts.get())); // this creates fill in geometry + + if (ndcs>4) { // a simple "road" + trig->removeInternalTriangles(dc3.get()); + dc3->setTexture ("road.png"); + dc3->setTexrep(40,9.5); // texture is repeated at this frequency + geode->addDrawable(dc3->makeGeometry(points)); // this creates road geometry + + if (ndcs>5) { + if (ndcs>6) { // road & forest overlap - order of removal is important + trig->removeInternalTriangles(forestroad.get()); + trig->removeInternalTriangles(forestroad2.get()); + trig->removeInternalTriangles(forestroad3.get()); + } + trig->removeInternalTriangles(forest.get()); + forest->setTexture("forestRoof.png"); + osg::ref_ptr locpts=forest->getPoints(points); + geode->addDrawable(forest->makeAreal(locpts.get())); + + forest->setWallTexture("forestEdge.png"); + geode->addDrawable(forest->makeWallGeometry(locpts.get()) ); + for (osg::Vec3Array::iterator vit=(*locpts).begin(); vit!=(*locpts).end(); vit++) { + (*vit)+=osg::Vec3(0,0,30); + } + + if (ndcs>6) {// road & forest overlap + forestroad->setTexture ("road.png"); + forestroad->setTexrep(40,22); // texture is repeated at this frequency + geode->addDrawable(forestroad->makeGeometry(points)); // this creates road geometry + forestroad2->setTexture ("road.png"); + forestroad2->setTexrep(40,22); // texture is repeated at this frequency + geode->addDrawable(forestroad2->makeGeometry(points)); // this creates road geometry + forestroad3->setTexture ("road.png"); + forestroad3->setTexrep(40,22); // texture is repeated at this frequency + geode->addDrawable(forestroad3->makeGeometry(points)); // this creates road geometry + if (ndcs>7) {// several overlapping DC's - add geom + trig->removeInternalTriangles(dc6.get()); + // dc6->makeDrawable(); + // dc6a->makeDrawable(); + dc6->setTexture ("road.png"); + dc6->setTexrep(40,22); // texture is repeated at this frequency + geode->addDrawable(dc6->makeGeometry(points)); // this creates road geometry + trig->removeInternalTriangles(dc6a.get()); + dc6a->setTexture ("road.png"); + dc6a->setTexrep(40,22); // texture is repeated at this frequency + geode->addDrawable(dc6a->makeGeometry(points)); // this creates road geometry + if (dc8.valid()) { + trig->removeInternalTriangles(dc8.get()); + dc8->setTexture ("road.png"); + dc8->setTexrep(40,16); // texture is repeated at this frequency + geode->addDrawable(dc8->makeGeometry(points)); // this creates road geometry + } + } + } + } + } + } + } + } + grp->addChild(geode.get()); + grp->addChild(createHUD(ndcs,what.str())); + return grp.release(); +} + +class KeyboardEventHandler : public osgGA::GUIEventHandler +{ // extra event handler traps 'n' key to re-triangulate the basic terrain. +public: + + KeyboardEventHandler(osg::Node *nd,osgProducer::Viewer &vr): + _scene(nd), viewer(vr), iview(0) {} + + virtual bool handle(const osgGA::GUIEventAdapter& ea,osgGA::GUIActionAdapter&) + { + switch(ea.getEventType()) + { + case(osgGA::GUIEventAdapter::KEYDOWN): + { + if (_scene && ea.getKey()=='n') + { + // re-tesselate the scene graph. + // the same contours are re-tesselated using a new method. Old contours + // & tesselation type are held internally in the derived Geode class tesselateDemoGeometry. + // cxTesselateVisitor tsv; + // _scene->accept(tsv); + iview++; + if (iview>10) iview=0; + osg::ref_ptr loadedModel = makedelaunay(iview); + viewer.setSceneData(loadedModel.get()); + return true; + } + break; + } + default: + break; + } + return false; + } + + virtual void accept(osgGA::GUIEventHandlerVisitor& v) + { + v.visit(*this); + } + + osg::Node *_scene; + osgProducer::Viewer &viewer; + int iview; +}; +int main( int argc, char **argv ) +{ + + // use an ArgumentParser object to manage the program arguments. + osg::ArgumentParser arguments(&argc,argv); + + // set up the usage document, in case we need to print out how to use this program. + arguments.getApplicationUsage()->setApplicationName(arguments.getApplicationName()); + arguments.getApplicationUsage()->setDescription(arguments.getApplicationName()+" is the standard OpenSceneGraph example which loads and visualises 3d models."); + arguments.getApplicationUsage()->setCommandLineUsage(arguments.getApplicationName()+" [options] filename ..."); + arguments.getApplicationUsage()->addCommandLineOption("--image ","Load an image and render it on a quad"); + arguments.getApplicationUsage()->addCommandLineOption("--dem ","Load an image/DEM and render it on a HeightField"); + arguments.getApplicationUsage()->addCommandLineOption("-h or --help","Display command line paramters"); + arguments.getApplicationUsage()->addCommandLineOption("--help-env","Display environmental variables available"); + arguments.getApplicationUsage()->addCommandLineOption("--help-keys","Display keyboard & mouse bindings available"); + arguments.getApplicationUsage()->addCommandLineOption("--help-all","Display all command line, env vars and keyboard & mouse bindigs."); + + + // construct the viewer. + osgProducer::Viewer viewer(arguments); + + // set up the value with sensible default event handlers. + viewer.setUpViewer(osgProducer::Viewer::STANDARD_SETTINGS); + + // get details on keyboard and mouse bindings used by the viewer. + viewer.getUsage(*arguments.getApplicationUsage()); + + // if user request help write it out to cout. + bool helpAll = arguments.read("--help-all"); + unsigned int helpType = ((helpAll || arguments.read("-h") || arguments.read("--help"))? osg::ApplicationUsage::COMMAND_LINE_OPTION : 0 ) | + ((helpAll || arguments.read("--help-env"))? osg::ApplicationUsage::ENVIRONMENTAL_VARIABLE : 0 ) | + ((helpAll || arguments.read("--help-keys"))? osg::ApplicationUsage::KEYBOARD_MOUSE_BINDING : 0 ); + if (helpType) + { + arguments.getApplicationUsage()->write(std::cout, helpType); + return 1; + } + + // report any errors if they have occured when parsing the program aguments. + if (arguments.errors()) + { + arguments.writeErrorMessages(std::cout); + return 1; + } + + if (arguments.argc()<1) + { + arguments.getApplicationUsage()->write(std::cout,osg::ApplicationUsage::COMMAND_LINE_OPTION); + return 1; + } + + osg::Timer_t start_tick = osg::Timer::instance()->tick(); + + // create the scene from internal specified terrain/constraints. + osg::ref_ptr loadedModel = makedelaunay(0); + + // if no model has been successfully loaded report failure. + if (!loadedModel) + { + std::cout << arguments.getApplicationName() <<": No data loaded" << std::endl; + return 1; + } + + + // any option left unread are converted into errors to write out later. + arguments.reportRemainingOptionsAsUnrecognized(); + + // report any errors if they have occured when parsing the program aguments. + if (arguments.errors()) + { + arguments.writeErrorMessages(std::cout); + } + + osg::Timer_t end_tick = osg::Timer::instance()->tick(); + + std::cout << "Time to load = "<delta_s(start_tick,end_tick)<0.0) { + osg::Vec3 off(0,0,height); + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprgetMode()==osg::PrimitiveSet::LINE_LOOP || + prset->getMode()==osg::PrimitiveSet::LINE_STRIP) { // nothing else loops + // start with the last point on the loop + for (unsigned int i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + wall->push_back(curp); + wall->push_back(curp+off); + } + const osg::Vec3 curp=(*vertices)[prset->index (0)]; + wall->push_back(curp); + wall->push_back(curp+off); + } + } + } + return wall; +} +osg::Vec2Array * WallConstraint::getWallTexcoords(const osg::Vec3Array *points,const float height) const +{ // return array of points for a wall height high around the constraint + osg::Vec2Array *tcoords= NULL; + if (height>0.0) { + float texrepRound=txxrepWall; + tcoords= new osg::Vec2Array; + float circumference=0; // distance around wall to get exact number of repeats of texture + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprindex (prset->getNumIndices()-1)]; + unsigned int i; + for (i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + circumference+=(curp-prevp).length(); + prevp=curp; + } + const osg::Vec3 curp=(*vertices)[prset->index (0)]; + circumference+=(curp-prevp).length(); + + int nround=(int)(circumference/txxrepWall); + if (nround<1) nround=1; // at least one repeat. + texrepRound=circumference/nround; + + float ds=0; + prevp=(*vertices)[prset->index (prset->getNumIndices()-1)]; + if (tcoords) { + for (i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + osg::Vec2 tci=osg::Vec2f(ds/texrepRound,0/txyrepWall); + tcoords->push_back(tci); + tci=osg::Vec2f(ds/texrepRound,height/txyrepWall); + tcoords->push_back(tci); + ds+=(curp-prevp).length(); + prevp=curp; + } + osg::Vec2 tci=osg::Vec2f(ds/texrepRound,0/txyrepWall); + tcoords->push_back(tci); + tci=osg::Vec2f(ds/texrepRound,height/txyrepWall); + tcoords->push_back(tci); + } + } // per primitiveset + + } + return tcoords; +} +osg::Geometry *WallConstraint::makeWallGeometry(const osg::Vec3Array *points) const +{ + osg::ref_ptr gm=new osg::Geometry; // the wall + if (texture!="") { + osg::Image* image = osgDB::readImageFile(texture.c_str()); + if (image) + { + osg::Texture2D* txt = new osg::Texture2D; + osg::StateSet* stateset = gm->getOrCreateStateSet(); + txt->setImage(image); + txt->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT ); + txt->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP ); + stateset->setTextureAttributeAndModes(0,txt,osg::StateAttribute::ON); + osg::Material* material = new osg::Material; + material->setAmbient(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,0.0f,1.0f)); + material->setDiffuse(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,1.0f,1.0f)); + stateset->setAttribute(material,osg::StateAttribute::ON); + stateset->setMode( GL_LIGHTING, osg::StateAttribute::ON ); + } + } + gm->setVertexArray(getWall(points,height)); + gm->addPrimitiveSet(makeWall()); + gm->setTexCoordArray(0,getWallTexcoords(points,height)); + gm->setNormalBinding(osg::Geometry::BIND_PER_VERTEX); + gm->setNormalArray(getWallNormals(points)); // this creates normals to walls + + return gm.release(); +} + +osg::Vec3Array *ArealConstraint::getWallNormals(const osg::Vec3Array *points) const +{ + osg::Vec3Array *nrms=new osg::Vec3Array; + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprgetMode()==osg::PrimitiveSet::LINE_LOOP) { // nothing else loops + // start with the last point on the loop + osg::Vec3 prevp=(*vertices)[prset->index (prset->getNumIndices()-1)]; + for (unsigned int i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + osg::Vec3 nrm=(curp-prevp)^osg::Vec3(0,0,1); + nrm.normalize(); + nrms->push_back(nrm); + nrms->push_back(nrm); + prevp=curp; + } + const osg::Vec3 curp=(*vertices)[prset->index (0)]; + osg::Vec3 nrm=(curp-prevp)^osg::Vec3(0,0,1); + nrm.normalize(); + nrms->push_back(nrm); + nrms->push_back(nrm); + } + } + return nrms; +} + + +osg::Vec3Array * ArealConstraint::getWall(const osg::Vec3Array* points,const float height) const +{ // return array of points for a wall height high around the constraint + osg::Vec3Array *wall=new osg::Vec3Array; + if (height>0.0) { + osg::Vec3 off(0,0,height); + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprgetMode()==osg::PrimitiveSet::LINE_LOOP) { // nothing else loops + // start with the last point on the loop + for (unsigned int i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + wall->push_back(curp); + wall->push_back(curp+off); + } + const osg::Vec3 curp=(*vertices)[prset->index (0)]; + wall->push_back(curp); + wall->push_back(curp+off); + } + } + } + return wall; +} + +osg::Vec2Array * ArealConstraint::getWallTexcoords(const osg::Vec3Array *points,const float height) const +{ // return array of points for a wall height high around the constraint + osg::Vec2Array *tcoords= NULL; + if (height>0.0) { + float texrepRound=txxrepWall; + tcoords= new osg::Vec2Array; + float circumference=0; // distance around wall to get exact number of repeats of texture + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprindex (prset->getNumIndices()-1)]; + unsigned int i; + for (i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + circumference+=(curp-prevp).length(); + prevp=curp; + } + const osg::Vec3 curp=(*vertices)[prset->index (0)]; + circumference+=(curp-prevp).length(); + + int nround=(int)(circumference/txxrepWall); + if (nround<1) nround=1; // at least one repeat. + texrepRound=circumference/nround; + + float ds=0; + prevp=(*vertices)[prset->index (prset->getNumIndices()-1)]; + if (tcoords) { + for (i=0; igetNumIndices(); i++) { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + osg::Vec2 tci=osg::Vec2f(ds/texrepRound,0/txyrepWall); + tcoords->push_back(tci); + tci=osg::Vec2f(ds/texrepRound,height/txyrepWall); + tcoords->push_back(tci); + ds+=(curp-prevp).length(); + prevp=curp; + } + osg::Vec2 tci=osg::Vec2f(ds/texrepRound,0/txyrepWall); + tcoords->push_back(tci); + tci=osg::Vec2f(ds/texrepRound,height/txyrepWall); + tcoords->push_back(tci); + } + } // per primitiveset + } + return tcoords; +} +osg::DrawArrays* ArealConstraint::makeCanopy( void ) const +{ + return (new osg::DrawArrays(osg::PrimitiveSet::TRIANGLES,0,3*_interiorTris.size())); +} +osg::Vec3Array *ArealConstraint::getCanopy(const osg::Vec3Array *points,const float height) const +{ // returns the array of vertices in the canopy + osg::Vec3 off(0,0,height); + osg::Vec3Array *internals=new osg::Vec3Array; + trilist::const_iterator tritr; + for (tritr=_interiorTris.begin(); tritr!=_interiorTris.end();tritr++) { + for (int i=0; i<3; i++) { + int index=(*tritr)[i]; + internals->push_back((*points)[index]+off); + } + } + return internals; +} +osg::Vec3Array *ArealConstraint::getCanopyNormals(const osg::Vec3Array *points) const +{ + osg::Vec3Array *nrms=new osg::Vec3Array; + trilist::const_iterator tritr; + for (tritr=_interiorTris.begin(); tritr!=_interiorTris.end();tritr++) { + osg::Vec3 e1=(*points)[(*tritr)[1]]-(*points)[(*tritr)[0]]; + osg::Vec3 e2=(*points)[(*tritr)[2]]-(*points)[(*tritr)[0]]; + osg::Vec3 nrm=e1^e2; + nrm.normalize(); + nrms->push_back(nrm); + } + return nrms; +} + +osg::Vec2Array *ArealConstraint::getCanopyTexcoords(const osg::Vec3Array *points) const +{ + osg::Vec3Array::const_iterator tritr; + osg::ref_ptr tcoords= new osg::Vec2Array ; + for (tritr=points->begin(); tritr!=points->end();tritr++) { + // calculate tcoords for terrain from xy drape. + osg::Vec2 tci=osg::Vec2f(tritr->x()/txxrepArea, tritr->y()/txyrepArea); + tcoords->push_back(tci); + } + return tcoords.release(); +} + +osg::DrawArrays * ArealConstraint::makeWall(void) const +{ // build a wall height high around the constraint + const osg::Vec3Array *_line= dynamic_cast(getVertexArray()); + return (new osg::DrawArrays(osg::PrimitiveSet::QUAD_STRIP,0,2+2*_line->size())); +} + +osg::Geometry *ArealConstraint::makeWallGeometry( osg::Vec3Array *pt) +{ + osg::ref_ptr gm=new osg::Geometry; // the wall + osg::ref_ptr edges=new osg::Geometry; // edges of bounds + edges->setVertexArray(pt); + osg::DrawElementsUInt *trgeom=getTriangles(); + edges->addPrimitiveSet(trgeom); + + osg::ref_ptr tscx=new osgUtil::Tesselator; // this assembles all the constraints + tscx->setTesselationType(osgUtil::Tesselator::TESS_TYPE_GEOMETRY); + tscx->setBoundaryOnly(true); + tscx->setWindingType( osgUtil::Tesselator::TESS_WINDING_NONZERO); + // find all edges. + const osg::Vec3Array *points=dynamic_cast(getVertexArray()); + + tscx->retesselatePolygons(*(edges)); // find all edges + + if (walltexture!="") { + osg::Image* image = osgDB::readImageFile(walltexture.c_str()); + if (image) + { + osg::Texture2D* txt = new osg::Texture2D; + osg::StateSet* stateset = gm->getOrCreateStateSet(); + txt->setImage(image); + txt->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT ); + txt->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::CLAMP ); + stateset->setTextureAttributeAndModes(0,txt,osg::StateAttribute::ON); + } + } + points=dynamic_cast(edges->getVertexArray()); + int nstart=0; + osg::ref_ptr coords=new osg::Vec3Array; + osg::ref_ptr tcoords=new osg::Vec2Array; + for (unsigned int i=0; igetNumPrimitiveSets(); i++) { + osg::PrimitiveSet *pr=edges->getPrimitiveSet(i); + if (pr->getMode() == osg::PrimitiveSet::LINE_LOOP) { + float ds=0; + for (unsigned int icon=0; icongetNumIndices(); icon++) { + unsigned int ithis=pr->index(icon); + osg::Vec3 pt= (*points)[ithis]; + coords->push_back(pt); + coords->push_back(pt+osg::Vec3(0,0,height)); + tcoords->push_back(osg::Vec2(ds/txxrepWall,0)); + tcoords->push_back(osg::Vec2(ds/txxrepWall,1.0)); + if (icongetNumIndices()-1) ds+=((*points)[pr->index(icon+1)]-(*points)[ithis]).length(); + else ds+=((*points)[pr->index(0)]-(*points)[ithis]).length(); + } + // repeat first point + unsigned int ithis=pr->index(0); + coords->push_back((*points)[ithis]); + coords->push_back((*points)[ithis]+osg::Vec3(0,0,height)); + tcoords->push_back(osg::Vec2(ds/txxrepWall,0)); + tcoords->push_back(osg::Vec2(ds/txxrepWall,1.0)); + gm->setVertexArray(coords.get()); + gm->setTexCoordArray(0,tcoords.get()); + gm->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::QUAD_STRIP,nstart,2+2*pr->getNumIndices())); + nstart+=2+2*pr->getNumIndices(); + } + } + + return gm.release(); +} + + +osg::Geometry * ArealConstraint::makeAreal( osg::Vec3Array *points) +{ + osg::ref_ptr gm; // the fill in area + if (_interiorTris.size()>0) { + gm =new osg::Geometry; // the forest roof + gm->setVertexArray(points); + osg::DrawElementsUInt *trgeom=getTriangles(); + gm->addPrimitiveSet(trgeom); + gm->setTexCoordArray(0,getCanopyTexcoords(points)); + osg::Image* image = osgDB::readImageFile(texture); + if (image) + { + osg::Texture2D* txt = new osg::Texture2D; + osg::StateSet* stateset = gm->getOrCreateStateSet(); + txt->setImage(image); + txt->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT ); + txt->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::REPEAT ); + stateset->setTextureAttributeAndModes(0,txt,osg::StateAttribute::ON); + osg::Material* material = new osg::Material; + material->setAmbient(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,1.0f,1.0f)); + material->setDiffuse(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,1.0f,1.0f)); + stateset->setAttribute(material,osg::StateAttribute::ON); + stateset->setMode( GL_LIGHTING, osg::StateAttribute::ON ); + } + } + return gm.release(); +} + + +void LinearConstraint::setVertices( osg::Vec3Array *lp, const float w) +{ // generate constant width around line (calls setvertices(edges)) + osg::ref_ptr edges=new osg::Vec3Array; + _tcoords=new osg::Vec2Array; // texture coordinates for replacement geometry + _edgecoords=new osg::Vec3Array; // posiiton coordinates for replacement geometry + width=w; + _midline=lp; + float ds=0; + for(unsigned int i=0;isize();i++) { + osg::Vec3 valong; + osg::Vec3 pos[2]; + + if (i==0) { + valong=(*lp)[i+1]-(*lp)[i]; + } else if (i==lp->size()-1) { + valong=(*lp)[i]-(*lp)[i-1]; + } else { + valong=(*lp)[i+1]-(*lp)[i-1]; + } + valong.normalize(); + osg::Vec3 vperp=valong^osg::Vec3(0,0,1); + pos[0]=(*lp)[i]-vperp*.5*width; + pos[1]=(*lp)[i]+vperp*.5*width; + edges->push_back(pos[0]); + _edgecoords->push_back(pos[0]); + _tcoords->push_back(osg::Vec2(0/txyrepAcross,ds/txxrepAlong)); + edges->insert(edges->begin() ,pos[1]); + _edgecoords->insert(_edgecoords->begin() ,pos[1]); + _tcoords->insert(_tcoords->begin() ,osg::Vec2(width/txyrepAcross,ds/txxrepAlong)); + if (isize()-1) ds+=((*lp)[i+1]-(*lp)[i]).length(); + } + setVertexArray(edges.get()); + addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,0,edges->size()) ); +} + +osg::DrawArrays* LinearConstraint::makeRoad(void ) const +{ + return new osg::DrawArrays(osg::PrimitiveSet::QUAD_STRIP,0,2*_midline->size()); + +} + +osg::Vec3Array *LinearConstraint::getRoadNormals(const osg::Vec3Array *points) const +{ + osg::Vec3Array *nrms=new osg::Vec3Array; + for(unsigned int i=0;i<_midline->size();i++) { + osg::Vec3 valong; // vector along midline of road + if (i==0) { + valong=(*_midline)[i+1]-(*_midline)[i]; + } else if (i==_midline->size()-1) { + valong=(*_midline)[i]-(*_midline)[i-1]; + } else { + valong=(*_midline)[i+1]-(*_midline)[i-1]; + } + osg::Vec3 vperp=valong^osg::Vec3(0,0,1); + osg::Vec3 nrm=vperp^valong; // normal to linear + nrm.normalize(); + nrms->push_back(nrm); // repeated for each vertex of linear. + nrms->push_back(nrm); + } + return nrms; +} +osg::Vec3Array *LinearConstraint::getRoadVertices(const osg::Vec3Array *points) const +{ + osg::Vec3Array *linearEdges=new osg::Vec3Array; + for(unsigned int i=0;i<_midline->size();i++) { + osg::Vec3 valong; // vector along midline of road + if (i==0) { + valong=(*_midline)[i+1]-(*_midline)[i]; + } else if (i==_midline->size()-1) { + valong=(*_midline)[i]-(*_midline)[i-1]; + } else { + valong=(*_midline)[i+1]-(*_midline)[i-1]; + } + valong.normalize(); + osg::Vec3 vperp=valong^osg::Vec3(0,0,1); // vector across road + // sides of linear + linearEdges->push_back((*_midline)[i]-vperp*.5*width); + linearEdges->push_back((*_midline)[i]+vperp*.5*width); + } + return linearEdges; +} + +osg::Vec2Array *LinearConstraint::getRoadTexcoords(const osg::Vec3Array *points) { + // need to create a vec2 array from the coordinates that fits the road + osg::Vec3Array::const_iterator tritr; + osg::ref_ptr tcoords= new osg::Vec2Array ; + for (tritr=points->begin(); tritr!=points->end();tritr++) { + osg::Vec2 tci(-1.,-1.); + int ib=0; + // osg::Vec3Array *varr=dynamic_cast(getVertexArray()); + bool ptfound=false; + for (osg::Vec3Array::iterator vit=_edgecoords->begin(); vit!= _edgecoords->end() && !ptfound; vit++) { + if ((*vit)==(*tritr)) { + tci=_tcoords->at(ib); + ptfound=true; + } + ib++; + } + if (!ptfound) { // search for surrounding points and interpolate + ib=0; + osg::Vec3 pminus=(_edgecoords->back()); // need pminus for interpolation + int ibm1=_edgecoords->size()-1; + for (osg::Vec3Array::iterator vit=_edgecoords->begin(); vit!= _edgecoords->end() /*&& !ptfound*/; vit++) { + osg::Vec3 pplus=(*vit)-(*tritr); + osg::Vec3 dpm=pminus-(*tritr); + pplus.set (pplus.x(),pplus.y(),0); + dpm.set (dpm.x(),dpm.y(),0); + float dprod=pplus*dpm/(pplus.length() * dpm.length()); + if (dprod<-0.9999) { // *tritr lies between.... + osg::Vec2 tminus=_tcoords->at(ibm1); + osg::Vec2 tplus=_tcoords->at(ib); + float frac=(dpm.length()/(dpm.length()+pplus.length())); + tci=tminus+((tplus-tminus)*frac); + ptfound=true; + } + ibm1=ib; + ib++; + pminus=(*vit); + } + } + tcoords->push_back(tci); + } + // some extra points are not interpolated as they lie between 2 interpolated vertices + for (tritr=points->begin(); tritr!=points->end();tritr++) { + int ib=tritr-points->begin(); + osg::Vec2 tci=tcoords->at(ib); + if (tci.x()<-.99 && tci.y()<-.99) { + // search through each of the primitivesets + osg::Vec3Array::const_iterator ptitr; + // osg::notify(osg::WARN) << "Not calculated " << (*tritr).x() <<"," << (*tritr).y() << std::endl; + for (ptitr=points->begin(); ptitr!=points->end();ptitr++) { + } + } + } + return tcoords.release(); +} +osg::Vec3Array * LinearConstraint::getNormals(const osg::Vec3Array *points) +{ + osg::ref_ptr norms=new osg::Vec3Array; + for (osg::DrawElementsUInt::iterator uiitr=prim_tris_->begin(); uiitr!=prim_tris_->end();uiitr+=3) { + osg::Vec3 e1=(*points)[*(uiitr+1)]-(*points)[(*uiitr)]; + osg::Vec3 e2=(*points)[*(uiitr+2)]-(*points)[*(uiitr+1)]; + osg::Vec3 n=e1^e2; + n.normalize(); + // if (n.z()<0) n=-n; + norms->push_back(n); + } + return norms.release(); +} + +osg::Geometry * LinearConstraint::makeGeometry(const osg::Vec3Array *points) +{ + osg::ref_ptr gm=new osg::Geometry; // the fill in road/railway + if (_midline->size()>0) { + osg::ref_ptr locpts=getPoints(points); + if (texture!="") { + osg::Image* image = osgDB::readImageFile(texture.c_str()); + if (image) + { + osg::Texture2D* txt = new osg::Texture2D; + osg::StateSet* stateset = gm->getOrCreateStateSet(); + txt->setImage(image); + txt->setWrap( osg::Texture2D::WRAP_S, osg::Texture2D::REPEAT ); + txt->setWrap( osg::Texture2D::WRAP_T, osg::Texture2D::REPEAT ); + stateset->setTextureAttributeAndModes(0,txt,osg::StateAttribute::ON); + osg::Material* material = new osg::Material; + material->setAmbient(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,1.0f,1.0f)); + material->setDiffuse(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,1.0f,1.0f)); + stateset->setAttribute(material,osg::StateAttribute::ON); + stateset->setMode( GL_LIGHTING, osg::StateAttribute::ON ); + } + gm->setTexCoordArray(0,getRoadTexcoords(locpts.get())); + } + gm->setVertexArray(locpts.get()); + gm->setNormalArray(getNormals(locpts.get())); + gm->setNormalBinding(osg::Geometry::BIND_PER_PRIMITIVE); + gm->addPrimitiveSet(getTriangles()); + } + + return gm.release(); + +} + + diff --git a/include/osgUtil/DelaunayTriangulator b/include/osgUtil/DelaunayTriangulator index da2a7858e..809a5fef4 100644 --- a/include/osgUtil/DelaunayTriangulator +++ b/include/osgUtil/DelaunayTriangulator @@ -14,22 +14,96 @@ #ifndef OSGUTIL_DELAUNAYTRIANGULATOR_ #define OSGUTIL_DELAUNAYTRIANGULATOR_ +#include + #include #include #include #include #include +#include #include namespace osgUtil { -/** Utility class that triangulates an irregular network of sample points. +/** DelaunayTriangulator: Utility class that triangulates an irregular network of sample points. Just create a DelaunayTriangulator, assign it the sample point array and call its triangulate() method to start the triangulation. Then you can obtain the generated primitive by calling the getTriangles() method. + + Add DelaunayConstraints (or derived class) to control the triangulation edges. */ +class OSGUTIL_EXPORT DelaunayConstraint: public osg::Geometry { + // controls the edges in a Delaunay triangulation. + // constraints can be linear (with width), areal (contains an area) + // uses: to replace part of a terrain with an alternative textured model (roads, lakes). + // the primitive sets in this are either LINE_LOOP or LINE_STRIP +public: + DelaunayConstraint() { } + + /** Each primitiveset is a list of vertices which may be closed by joining up to its start + * to make a loop. Constraints should be simple lines, not crossing themselves. + * Constraints which cross other constraints can cause difficulties - see the example + * for methods of dealing with them. */ + + /** collect up indices of triangle from delaunay triangles. + * The delaunay triangles inside the DelaunayConstraint area can be used to fill + * the area or generate geometry that terrain follows the area in some way. + * These triangles can form a canopy or a field. */ + void addtriangle(const int i1,const int i2, const int i3); + + /** Get the filling primitive. One: + * triangulate must have bneen called and + * two: triangle list is filled when + * DelaunayTriangulator::removeInternalTriangles is called. + * These return the triangles removed from the delaunay triangulation by + * DelaunayTriangulator::removeInternalTriangles. */ + inline const osg::DrawElementsUInt *getTriangles() const; + inline osg::DrawElementsUInt *getTriangles(); + + /** Call BEFORE makeDrawable to reorder points to make optimised set + */ + osg::Vec3Array *getPoints(const osg::Vec3Array *points); + + /** converts simple list of triangles into a drawarray. + */ + osg::DrawElementsUInt *makeDrawable(); + + /** Add vertices and constraint loops from dco + * Can be used to generate extra vertices where dco crosses 'this' using + * osgUtil::tesselator to insert overlap vertices. + */ + void merge(DelaunayConstraint *dco); + + /** remove from line the vertices that are inside dco + */ + void removeVerticesInside(const DelaunayConstraint *dco); + + /** return winding number as a float of loop around testpoint; may use multiple loops + * does not reject points on the edge or very very close to the edge */ + float windingNumber(const osg::Vec3 testpoint) const ; + + /** true if testpoint is internal (or external) to constraint. */ + virtual bool contains(const osg::Vec3 testpoint) const; + virtual bool outside(const osg::Vec3 testpoint) const; + + /** Tesselate the constraint loops so that the crossing points are interpolated + * and added to the contraints for the triangulation. */ + void DelaunayConstraint::handleOverlaps(void); + +protected: + virtual ~DelaunayConstraint() {} + + typedef std::vector< int* > trilist; // array of indices in points array defining triangles + + trilist _interiorTris; // list of triangles that fits the area. + + osg::ref_ptr prim_tris_; // returns a PrimitiveSet to draw the interior of this DC +}; + + class OSGUTIL_EXPORT DelaunayTriangulator: public osg::Referenced { public: @@ -37,6 +111,8 @@ public: explicit DelaunayTriangulator(osg::Vec3Array *points, osg::Vec3Array *normals = 0); DelaunayTriangulator(const DelaunayTriangulator ©, const osg::CopyOp ©op = osg::CopyOp::SHALLOW_COPY); + typedef std::vector< osg::ref_ptr > linelist; + /** Get the const input point array. */ inline const osg::Vec3Array *getInputPointArray() const; @@ -46,6 +122,17 @@ public: /** Set the input point array. */ inline void setInputPointArray(osg::Vec3Array *points); + /** Add an input constraint loop. + ** the edges of the loop will constrain the triangulation. + ** if remove!=0, the internal triangles of the constraint will be removed; + ** the user may the replace the constraint line with an equivalent geometry. + ** GWM July 2005 */ + void addInputConstraint(DelaunayConstraint *dc) { + constraint_lines.push_back(dc); + return; + } + + /** Get the const output normal array (optional). */ inline const osg::Vec3Array *getOutputNormalArray() const; @@ -63,15 +150,24 @@ public: /** Get the generated primitive (call triangulate() first). */ inline osg::DrawElementsUInt *getTriangles(); + + /** remove the triangles internal to the constraint loops. + * (Line strips cannot remove any internal triangles). */ + void removeInternalTriangles(DelaunayConstraint *constraint); protected: virtual ~DelaunayTriangulator(); DelaunayTriangulator &operator=(const DelaunayTriangulator &) { return *this; } + int getindex(const osg::Vec3 pt,const osg::Vec3Array *points); private: osg::ref_ptr points_; osg::ref_ptr normals_; osg::ref_ptr prim_tris_; + + // GWM these lines provide required edges in the triangulated shape. + linelist constraint_lines; + }; // INLINE METHODS @@ -116,6 +212,18 @@ inline osg::DrawElementsUInt *DelaunayTriangulator::getTriangles() return prim_tris_.get(); } +inline const osg::DrawElementsUInt *DelaunayConstraint::getTriangles() const +{ + return prim_tris_.get(); +} + +inline osg::DrawElementsUInt *DelaunayConstraint::getTriangles() +{ + return prim_tris_.get(); +} + + } #endif + diff --git a/src/osgUtil/DelaunayTriangulator.cpp b/src/osgUtil/DelaunayTriangulator.cpp index 9064251cf..178cb4acd 100644 --- a/src/osgUtil/DelaunayTriangulator.cpp +++ b/src/osgUtil/DelaunayTriangulator.cpp @@ -1,4 +1,19 @@ +/* -*-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 +// NB this algorithm makes heavy use of the osgUtil::Tesselator for constrained triangulation. +// truly it is built on the shoulders of giants. #include #include @@ -6,8 +21,9 @@ #include #include -#include #include +#include //GWM July 2005 map is used in constraints. +#include // tesselator triangulates the constrained triangles namespace osgUtil { @@ -53,20 +69,24 @@ inline osg::Vec3 compute_circumcircle( ((a.x() - c.x()) * (a.x() + c.x()) + (a.y() - c.y()) * (a.y() + c.y())) / 2 * (b.x() - c.x())) / D; - r2 = (c.x() - cx) * (c.x() - cx) + (c.y() - cy) * (c.y() - cy); + // r2 = (c.x() - cx) * (c.x() - cx) + (c.y() - cy) * (c.y() - cy); + // the return r square is compared with r*r many times in an inner loop + // so for efficiency use the inefficient sqrt once rather than 30* multiplies later. + r2 = sqrt((c.x() - cx) * (c.x() - cx) + (c.y() - cy) * (c.y() - cy)); } return osg::Vec3(cx, cy, r2); } // Test whether a point (only the x and y coordinates are used) lies inside -// a circle; the circle is passed as a vector: (Cx, Cy, r^2). +// a circle; the circle is passed as a vector: (Cx, Cy, r). inline bool point_in_circle(const osg::Vec3 &point, const osg::Vec3 &circle) { float r2 = (point.x() - circle.x()) * (point.x() - circle.x()) + (point.y() - circle.y()) * (point.y() - circle.y()); - return r2 <= circle.z(); + return r2 <= circle.z()*circle.z(); +// return r2 <= circle.z(); } @@ -86,7 +106,8 @@ class Edge { public: // Comparison object (for sorting) - struct Less { + struct Less + { inline bool operator()(const Edge &e1, const Edge &e2) const { if (e1.ibs() < e2.ibs()) return true; @@ -126,7 +147,8 @@ private: // CLASS: Triangle -class Triangle { +class Triangle +{ public: Triangle(Vertex_index a, Vertex_index b, Vertex_index c, osg::Vec3Array *points) : a_(a), @@ -146,14 +168,175 @@ public: inline const Edge &get_edge(int idx) const { return edge_[idx]; } inline const osg::Vec3 &get_circumcircle() const { return cc_; } + inline osg::Vec3 compute_centroid(const osg::Vec3Array *points) const + { + return ((*points)[a_] +(*points)[b_] + (*points)[c_])/3; + } + inline osg::Vec3 compute_normal(osg::Vec3Array *points) const { osg::Vec3 N = ((*points)[b_] - (*points)[a_]) ^ ((*points)[c_] - (*points)[a_]); return N / N.length(); } + bool isedge(const unsigned int ip1,const unsigned int ip2) const + { // is one of the edges of this triangle from ip1-ip2 + bool isedge=ip1==a() && ip2==b(); + if (!isedge) + { + isedge=ip1==b() && ip2==c(); + if (!isedge) + { + isedge=ip1==c() && ip2==a(); + } + } + return isedge; + } + // GWM July 2005 add test for triangle intersected by p1-p2. + // return true for unused edge + const bool intersected(const unsigned int ip1,const unsigned int ip2,const osg::Vec2 p1 ,const osg::Vec2 p2,const int iedge, osg::Vec3Array *points) const + { + // return true if edge iedge of triangle is intersected by ip1,ip2 + Vertex_index ie1,ie2; + if (iedge==0) + { + ie1=a(); + ie2=b(); + } + else if (iedge==1) + { + ie1=b(); + ie2=c(); + } + else if (iedge==2) + { + ie1=c(); + ie2=a(); + } + if (ip1==ie1 || ip2==ie1) return false; + if (ip1==ie2 || ip2==ie2) return false; + + osg::Vec2 tp1((*points)[ie1].x(),(*points)[ie1].y()); + osg::Vec2 tp2((*points)[ie2].x(),(*points)[ie2].y()); + return intersect(tp1,tp2,p1,p2); + } + + bool intersectedby(const osg::Vec2 p1,const osg::Vec2 p2,osg::Vec3Array *points) const { + // true if line [p1,p2] cuts at least one edge of this triangle + osg::Vec2 tp1((*points)[a()].x(),(*points)[a()].y()); + osg::Vec2 tp2((*points)[b()].x(),(*points)[b()].y()); + osg::Vec2 tp3((*points)[c()].x(),(*points)[c()].y()); + bool ip=intersect(tp1,tp2,p1,p2); + if (!ip) + { + ip=intersect(tp2,tp3,p1,p2); + if (!ip) + { + ip=intersect(tp3,tp1,p1,p2); + } + } + return ip; + } + int whichEdge(osg::Vec3Array *points,const osg::Vec2 p1, const osg::Vec2 p2, + const unsigned int e1,const unsigned int e2) const + { + int icut=0; + // find which edge of triangle is cut by line (p1-p2) and is NOT e1-e2 indices. + // return 1 - cut is on edge b-c; 2==c-a + osg::Vec2 tp1((*points)[a()].x(),(*points)[a()].y()); // triangle vertices + osg::Vec2 tp2((*points)[b()].x(),(*points)[b()].y()); + osg::Vec2 tp3((*points)[c()].x(),(*points)[c()].y()); + bool ip=intersect(tp2,tp3,p1,p2); + if (ip && (a()==e1 || a()==e2)) { return 1;} + ip=intersect(tp3,tp1,p1,p2); + if (ip && (b()==e1 || b()==e2)) { return 2;} + ip=intersect(tp1,tp2,p1,p2); + if (ip && (c()==e1 || c()==e2)) { return 3;} + return icut; + } + + bool usesVertex(const unsigned int ip) const + { + return ip==a_ || ip==b_ || ip==c_; + } + + int lineBisectTest(const osg::Vec3 apt,const osg::Vec3 bpt,const osg::Vec3 cpt, const unsigned int ip1, const osg::Vec2 p2) const + { + osg::Vec2 vp2tp=p2-osg::Vec2(apt.x(), apt.y()); // vector from p1 to a. + // test is: cross product (z component) with ab,ac is opposite signs + // AND dot product with ab,ac has at least one positive value. + osg::Vec2 vecba=osg::Vec2(bpt.x(), bpt.y())-osg::Vec2(apt.x(), apt.y()); + osg::Vec2 vecca=osg::Vec2(cpt.x(), cpt.y())-osg::Vec2(apt.x(), apt.y()); + float cprodzba=vp2tp.x()*vecba.y() - vp2tp.y()*vecba.x(); + float cprodzca=vp2tp.x()*vecca.y() - vp2tp.y()*vecca.x(); + // osg::notify(osg::WARN) << "linebisect test " << ip1 << " tri " << a_<<","<< b_<<","<< c_< and lies between the 2 edges which meet at vertex + // is that which uses index ip1. + // line is to p2 + // return value is 0 - no crossing; 1,2,3 - which edge of the triangle is cut. + if (a_==ip1) + { + // first vertex is the vertex - test that a_ to p2 lies beteen edges a,b and a,c + osg::Vec3 apt=(*points)[a_]; + osg::Vec3 bpt=(*points)[b_]; + osg::Vec3 cpt=(*points)[c_]; + return lineBisectTest(apt,bpt,cpt,ip1,p2)?1:0; + } + else if (b_==ip1) + { + // second vertex is the vertex - test that b_ to p2 lies beteen edges a,b and a,c + osg::Vec3 apt=(*points)[b_]; + osg::Vec3 bpt=(*points)[c_]; + osg::Vec3 cpt=(*points)[a_]; + return lineBisectTest(apt,bpt,cpt,ip1,p2)?2:0; + } + else if (c_==ip1) + { + // 3rd vertex is the vertex - test that c_ to p2 lies beteen edges a,b and a,c + osg::Vec3 apt=(*points)[c_]; + osg::Vec3 bpt=(*points)[a_]; + osg::Vec3 cpt=(*points)[b_]; + return lineBisectTest(apt,bpt,cpt,ip1,p2)?3:0; + } + return 0; + } + private: + + bool intersect(const osg::Vec2 p1,const osg::Vec2 p2,const osg::Vec2 p3,const osg::Vec2 p4) const + { + // intersection point of p1,p2 and p3,p4 + // test from http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/ + // the intersection must be internal to the lines, not an end point. + float det=((p4.y()-p3.y())*(p2.x()-p1.x())-(p4.x()-p3.x())*(p2.y()-p1.y())); + if (det!=0) + { + // point on line is P=p1+ua.(p2-p1) and Q=p3+ub.(p4-p3) + float ua=((p4.x()-p3.x())*(p1.y()-p3.y())-(p4.y()-p3.y())*(p1.x()-p3.x()))/det; + float ub=((p2.x()-p1.x())*(p1.y()-p3.y())-(p2.y()-p1.y())*(p1.x()-p3.x()))/det; + if (ua> 0.00 && ua< 1 && ub> 0.0000 && ub< 1) + { + return true; + } + } + return false; + } + Vertex_index a_; Vertex_index b_; Vertex_index c_; @@ -161,18 +344,22 @@ private: Edge edge_[3]; }; +typedef std::list Triangle_list; // comparison function for sorting sample points by the X coordinate bool Sample_point_compare(const osg::Vec3 &p1, const osg::Vec3 &p2) { - return p1.x() < p2.x(); + // replace pure sort by X coordinate with X then Y. + // errors can occur if the delaunay triangulation specifies 2 points at same XY and different Z + if (p1.x() != p2.x()) return p1.x() < p2.x(); + if (p1.y() != p2.y()) return p1.y() < p2.y(); // GWM 30.06.05 - further rule if x coords are same. + osg::notify(osg::INFO) << "Two points are coincident at "< Edge_set; -typedef std::list Triangle_list; - DelaunayTriangulator::DelaunayTriangulator(): @@ -199,17 +386,289 @@ DelaunayTriangulator::~DelaunayTriangulator() { } +const Triangle * getTriangleWithEdge(const unsigned int ip1,const unsigned int ip2, const Triangle_list *triangles) +{ // find triangle in list with edge from ip1 to ip2... + Triangle_list::const_iterator trconnitr; // connecting triangle + int idx=0; + for (trconnitr=triangles->begin(); trconnitr!=triangles->end(); ) + { + if (trconnitr->isedge (ip1,ip2)) + { + // this is the triangle. + return &(*trconnitr); + } + ++trconnitr; + idx++; + } + return NULL; //-1; +} + +int DelaunayTriangulator::getindex(const osg::Vec3 pt,const osg::Vec3Array *points) +{ + // return index of pt in points (or -1) + for (unsigned int i=0; isize(); i++) + { + if (pt.x()==(*points)[i].x() &&pt.y()==(*points)[i].y() ) + { + return i; + } + } + return -1; +} + +Triangle_list fillHole(osg::Vec3Array *points, std::vector vindexlist) +{ + // eg clockwise vertex neighbours around the hole made by the constraint + Triangle_list triangles; // returned list + osg::ref_ptr gtess=new osg::Geometry; // add all the contours to this for analysis + osg::ref_ptr constraintverts=new osg::Vec3Array; + osg::ref_ptr tscx=new osgUtil::Tesselator; // this assembles all the constraints + + for (std::vector::iterator itint=vindexlist.begin(); itint!=vindexlist.end(); itint++) + { + // osg::notify(osg::WARN)<< "add point " << (*itint) << " at " << (*points)[*itint].x() << ","<< (*points)[*itint].y() <push_back((*points)[*itint]); + } + + unsigned int npts=vindexlist.size(); + + gtess->setVertexArray(constraintverts.get()); + gtess->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::POLYGON,0,npts)); + tscx->setTesselationNormal(osg::Vec3(0.0,0.0,1.0)); + tscx->setTesselationType(osgUtil::Tesselator::TESS_TYPE_GEOMETRY); + tscx->setBoundaryOnly(false); + tscx->setWindingType( osgUtil::Tesselator::TESS_WINDING_ODD); // the commonest tesselation is default, ODD. GE2 allows intersections of constraints to be found. + tscx->retesselatePolygons(*(gtess.get())); // this should insert extra vertices where constraints overlap + + // extract triangles from gtess + + unsigned int ipr; + for (ipr=0; iprgetNumPrimitiveSets(); ipr++) + { + unsigned int ic; + osg::PrimitiveSet* prset=gtess->getPrimitiveSet(ipr); + // osg::notify(osg::WARN)<< "gtess set " << ipr << " nprims " << prset->getNumPrimitives() << + // " type " << prset->getMode() << std::endl; + unsigned int pidx,pidx1,pidx2; + switch (prset->getMode()) { + case osg::PrimitiveSet::TRIANGLES: + for (ic=0; icgetNumIndices()-2; ic+=3) + { + if (prset->index(ic)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic)]); + pidx=points->size(); + } + else + { + pidx=vindexlist[prset->index(ic)]; + } + + if (prset->index(ic+1)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic+1)]); + pidx1=points->size(); + } + else + { + pidx1=vindexlist[prset->index(ic+1)]; + } + + if (prset->index(ic+2)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic+2)]); + pidx2=points->size(); + } + else + { + pidx2=vindexlist[prset->index(ic+2)]; + } + triangles.push_back(Triangle(pidx, pidx1, pidx2, points)); + // osg::notify(osg::WARN)<< "vert " << prset->index(ic) << " in array"<getNumIndices()-2; ic++) + { + if (prset->index(ic)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic)]); + pidx=points->size(); + } else { + pidx=vindexlist[prset->index(ic)]; + } + if (prset->index(ic+1)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic+1)]); + pidx1=points->size(); + } + else + { + pidx1=vindexlist[prset->index(ic+1)]; + } + + if (prset->index(ic+2)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic+2)]); + pidx2=points->size(); + } + else + { + pidx2=vindexlist[prset->index(ic+2)]; + } + + if (ic%2==0) + { + triangles.push_back(Triangle(pidx, pidx1, pidx2, points)); + } + else + { + triangles.push_back(Triangle(pidx1, pidx, pidx2, points)); + } + // osg::notify(osg::WARN)<< "vert " << prset->index(ic) << " in array"<index(0)]; + if (prset->index(0)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(0)]); + pidx=points->size(); + } + else + { + pidx=vindexlist[prset->index(0)]; + } + // osg::notify(osg::WARN)<< "tfan has " << prset->getNumIndices() << " indices"<getNumIndices()-1; ic++) + { + if (prset->index(ic)>=npts) + { + // this is an added point. + points->push_back((*constraintverts)[prset->index(ic)]); + pidx1=points->size(); + } + else + { + pidx1=vindexlist[prset->index(ic)]; + } + + if (prset->index(ic+1)>=npts) + { // this is an added point. + points->push_back((*constraintverts)[prset->index(ic+1)]); + pidx2=points->size(); + } + else + { + pidx2=vindexlist[prset->index(ic+1)]; + } + triangles.push_back(Triangle(pidx, pidx1, pidx2, points)); + } + } + break; + default: + osg::notify(osg::WARN)<< "WARNING set " << ipr << " nprims " << prset->getNumPrimitives() << + " type " << prset->getMode() << " Type not triangle, tfan or strip"<< std::endl; + break; + } + } + return triangles; +} + +void DelaunayConstraint::removeVerticesInside(const DelaunayConstraint *dco) +{ /** remove vertices from this which are internal to dco. + * retains potins that are extremely close to edge of dco + * defined as edge of dco subtends>acs(0.999999) + */ + int nrem=0; + osg::Vec3Array *vertices= dynamic_cast< osg::Vec3Array*>(getVertexArray()); + for (osg::Vec3Array::iterator vitr=vertices->begin(); vitr!=vertices->end(); ) + { + if (dco->contains(*vitr)) + { + unsigned int idx=vitr-vertices->begin(); // index of vertex + // remove vertex index from all the primitives + for (unsigned int ipr=0; ipr(prset); + if (dsup) { + for (osg::DrawElementsUShort::iterator usitr=dsup->begin(); usitr!=dsup->end(); ) + { + if ((*usitr)==idx) + { // remove entirely + usitr=dsup->erase(usitr); + } + else + { + if ((*usitr)>idx) (*usitr)--; // move indices down 1 + usitr++; // next index + } + } + } + else + { + osg::notify(osg::WARN) << "Invalid prset " <getType() << " types PrimitiveType,DrawArraysPrimitiveType=1 etc" << std::endl; + } + } + vitr=vertices->erase(vitr); + nrem++; + + } + else + { + vitr++; + } + } +} + +void DelaunayConstraint::merge(DelaunayConstraint *dco) +{ + unsigned int ipr; + osg::Vec3Array *vmerge=dynamic_cast(getVertexArray()); + osg::Vec3Array * varr=dynamic_cast(dco->getVertexArray()); + if (!vmerge) vmerge=new osg::Vec3Array; + setVertexArray(vmerge); + for ( ipr=0; iprgetNumPrimitiveSets(); ipr++) + { + osg::PrimitiveSet* prset=dco->getPrimitiveSet(ipr); + osg::DrawArrays *drarr=dynamic_cast (prset); + if (drarr) + { + // need to add the offset of vmerge->size to each prset indices. + unsigned int noff=vmerge->size(); + unsigned int n1=drarr->getFirst(); // usually 0 + unsigned int numv=drarr->getCount(); // + addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,n1+noff,numv)); + } + } + vmerge->insert(vmerge->end(),varr->begin(),varr->end()); +} + + bool DelaunayTriangulator::triangulate() { // check validity of input array - if (!points_.valid()) { + if (!points_.valid()) + { osg::notify(osg::WARN) << "Warning: DelaunayTriangulator::triangulate(): invalid sample point array" << std::endl; return false; } osg::Vec3Array *points = points_.get(); - if (points->size() < 1) { + if (points->size() < 1) + { osg::notify(osg::WARN) << "Warning: DelaunayTriangulator::triangulate(): too few sample points" << std::endl; return false; } @@ -218,6 +677,31 @@ bool DelaunayTriangulator::triangulate() Triangle_list triangles; Triangle_list discarded_tris; + // GWM July 2005 add constraint vertices to terrain + linelist::iterator linitr; + for (linitr=constraint_lines.begin();linitr!=constraint_lines.end();linitr++) + { + DelaunayConstraint *dc=(*linitr).get(); + const osg::Vec3Array *vercon= dynamic_cast(dc->getVertexArray()); + int nadded=0; + for (unsigned int icon=0;iconsize();icon++) + { + osg::Vec3 p1=(*vercon)[icon]; + int idx=getindex(p1,points_.get()); + if (idx<0) + { // only unique vertices are permitted. + points_->push_back(p1); // add non-unique constraint points to triangulation + nadded++; + } + else + { + osg::notify(osg::WARN) << "DelaunayTriangulator: ignore a duplicate point at "<< p1.x()<< " " << p1.y() << std::endl;; + } + } + // osg::notify(osg::WARN)<< "constraint size "<size()<<" " <begin(), points->end(), Sample_point_compare); @@ -236,15 +720,19 @@ bool DelaunayTriangulator::triangulate() osg::notify(osg::INFO) << "DelaunayTriangulator: finding minimum and maximum Y values\n"; osg::Vec3Array::const_iterator mmi; - for (mmi=points->begin(); mmi!=points->end(); ++mmi) { + for (mmi=points->begin(); mmi!=points->end(); ++mmi) + { if (mmi->y() < miny) miny = mmi->y(); if (mmi->y() > maxy) maxy = mmi->y(); } // add supertriangle vertices to the point list - points->push_back(osg::Vec3(minx - (maxx - minx), miny, 0)); - points->push_back(osg::Vec3(maxx, miny, 0)); - points->push_back(osg::Vec3(maxx, maxy + (maxy - miny), 0)); + // gwm added 1.05* to ensure that supervertices are outside the domain of points. + // original value could make 2 coincident points for regular arrays of x,y,h data. + // this mod allows regular spaced arrays to be used. + points_->push_back(osg::Vec3(minx - 1.05*(maxx - minx), miny - 0.01*(maxy - miny), 0)); + points_->push_back(osg::Vec3(maxx + 0.01*(maxx - minx), miny - 0.01*(maxy - miny), 0)); + points_->push_back(osg::Vec3(maxx + 0.01*(maxx - minx), maxy + 1.05*(maxy - miny), 0)); // add supertriangle to triangle list triangles.push_back(Triangle(last_valid_index+1, last_valid_index+2, last_valid_index+3, points)); @@ -256,7 +744,8 @@ bool DelaunayTriangulator::triangulate() osg::notify(osg::INFO) << "DelaunayTriangulator: triangulating vertex grid (" << (points->size()-3) <<" points)\n"; - for (i=points->begin(); i!=points->end(); ++i, ++pidx) { + for (i=points->begin(); i!=points->end(); ++i, ++pidx) + { // don't process supertriangle vertices if (pidx > last_valid_index) break; @@ -265,7 +754,8 @@ bool DelaunayTriangulator::triangulate() // iterate through triangles Triangle_list::iterator j, next_j; - for (j=triangles.begin(); j!=triangles.end(); j = next_j) { + for (j=triangles.begin(); j!=triangles.end(); j = next_j) + { next_j = j; ++next_j; @@ -276,10 +766,15 @@ bool DelaunayTriangulator::triangulate() // OPTIMIZATION: since points are pre-sorted by the X component, // check whether we can discard this triangle for future operations float xdist = i->x() - cc.x(); - if ((xdist * xdist) > cc.z() && i->x() > cc.x()) { + // this is where the circumcircles radius rather than R^2 is faster. + // original code used r^2 and needed to test xdist*xdist>cc.z && i->x()>cc.x(). + if ((xdist ) > cc.z() ) + { discarded_tris.push_back(*j); triangles.erase(j); - } else { + } + else + { // if the point lies in the triangle's circumcircle then add // its edges to the edge list and remove the triangle @@ -294,7 +789,9 @@ bool DelaunayTriangulator::triangulate() // safe in this case since the set_duplicate is // not used as part of the Less operator. Edge& edge = const_cast(*(result.first)); - edge.set_duplicate(true); + // not clear why this change is needed? But prevents removal of twice referenced edges?? + // edge.set_duplicate(true); + edge.set_duplicate(!edge.get_duplicate()); } } triangles.erase(j); @@ -304,8 +801,10 @@ bool DelaunayTriangulator::triangulate() // remove duplicate edges and add new triangles Edge_set::iterator ci; - for (ci=edges.begin(); ci!=edges.end(); ++ci) { - if (!ci->get_duplicate()) { + for (ci=edges.begin(); ci!=edges.end(); ++ci) + { + if (!ci->get_duplicate()) + { triangles.push_back(Triangle(pidx, ci->ib(), ci->ie(), points)); } } @@ -321,6 +820,157 @@ bool DelaunayTriangulator::triangulate() // rejoin the two triangle lists triangles.insert(triangles.begin(), discarded_tris.begin(), discarded_tris.end()); + // GWM July 2005 eliminate any triangle with an edge crossing a constraint line + // http://www.geom.uiuc.edu/~samuelp/del_project.html + // we could also implement the sourcecode in http://gts.sourceforge.net/reference/gts-delaunay-and-constrained-delaunay-triangulations.html + // this uses the set of lines which are boundaries of the constraints, including points + // added to the contours by tesselation. + for (linelist::iterator dcitr=constraint_lines.begin();dcitr!=constraint_lines.end();dcitr++) + { + //DelaunayConstraint *dc=(*dcitr).get(); + const osg::Vec3Array *vercon= dynamic_cast((*dcitr)->getVertexArray()); + for (unsigned int ipr=0; ipr<(*dcitr)->getNumPrimitiveSets(); ipr++) + { + const osg::PrimitiveSet* prset=(*dcitr)->getPrimitiveSet(ipr); + if (prset->getMode()==osg::PrimitiveSet::LINE_LOOP || + prset->getMode()==osg::PrimitiveSet::LINE_STRIP) + { + // loops or strips + // start with the last point on the loop + unsigned int ip1=getindex((*vercon)[prset->index (prset->getNumIndices()-1)],points_.get()); + for (unsigned int i=0; igetNumIndices(); i++) + { + unsigned int ip2=getindex((*vercon)[prset->index(i)],points_.get()); + if (i>0 || prset->getMode()==osg::PrimitiveSet::LINE_LOOP) + { + // dont check edge from end to start + // for strips + // 2 points on the constraint + bool edgused=false;// first check for exact edge indices are used. + Triangle_list::iterator titr; + const osg::Vec3 curp=(*vercon)[prset->index(i)]; + int it=0; + for (titr=triangles.begin(); titr!=triangles.end() && !edgused; ++titr) + { + //check that the edge ip1-ip2 is not already part of the triangulation. + if (titr->isedge(ip1,ip2)) edgused=true; + if (titr->isedge(ip2,ip1)) edgused=true; + // if (edgused) osg::notify(osg::WARN) << "Edge used in triangle " << it << " " << + // titr->a()<<","<< titr->b()<<","<< titr->c()<< std::endl; + it++; + } + if (!edgused) + { + // then check for intermediate triangles, erase them and replace with constrained triangles. + // find triangle with point ip1 where the 2 edges from ip1 contain the line p1-p2. + osg::Vec2 p1((*points_)[ip1].x(),(*points_)[ip1].y()); // a constraint line joins p1-p2 + osg::Vec2 p2((*points_)[ip2].x(),(*points_)[ip2].y()); + int ntr=0; + std::vector trisToDelete; // array of triangles to delete from terrain. + // form 2 lists of vertices for the edges of the hole created. + // The hole joins vertex ip1 to ip2, and one list of edges lies to the left + // of the line ip1-ip2m the other to the right. + // a list of vertices forming 2 halves of the removed triangles. + // which in turn are filled in with the tesselator. + for (titr=triangles.begin(); titr!=triangles.end(); ) + { + int icut=titr->lineBisects(points_.get(),ip1,p2); + // osg::notify(osg::WARN) << "Testing triangle " << ntr << " "<< ip1 << " ti " << + // titr->a()<< ","<b() <<"," <c() << std::endl; + if (icut>0) + { + // triangle titr starts the constraint edge + std::vector edgeRight, edgeLeft; + edgeRight.push_back(ip1); + edgeLeft.push_back(ip1); + // osg::notify(osg::WARN) << "hole first " << edgeLeft.back()<< " rt " << edgeRight.back()<< std::endl; + trisToDelete.push_back(&(*titr)); + // now find the unique triangle that shares the defined edge + unsigned int e1, e2; // indices of ends of test triangle titr + if (icut==1) + { + // icut=1 implies vertex a is not involved + e1=titr->b(); e2=titr->c(); + } + else if (icut==2) + { + e1=titr->c(); e2=titr->a(); + } + else if (icut==3) + { + e1=titr->a(); e2=titr->b(); + } + edgeRight.push_back(e2); + edgeLeft.push_back(e1); + // osg::notify(osg::WARN) << icut << "hole edges " << edgeLeft.back()<< " rt " << edgeRight.back()<< std::endl; + const Triangle *tradj=getTriangleWithEdge(e2,e1, &triangles); + if (tradj) + { + while (tradj && !tradj->usesVertex(ip2) && trisToDelete.size()<999) + { + trisToDelete.push_back(tradj); + icut=tradj->whichEdge(points_.get(),p1,p2,e1,e2); + // osg::notify(osg::WARN) << ntr << " cur triedge " << icut << " " << ip1 << + // " to " << ip2 << " tadj " << tradj->a()<< ","<b() <<"," + // <c() <b(); e2=tradj->c();} // icut=1 implies vertex a is not involved + else if (icut==2) {e1=tradj->c(); e2=tradj->a();} + else if (icut==3) {e1=tradj->a(); e2=tradj->b();} + if (edgeLeft.back()!=e1 && edgeRight.back()==e2 && e1!=ip2) { + edgeLeft.push_back(e1); + } else if(edgeRight.back()!=e2 && edgeLeft.back()==e1 && e2!=ip2) { + edgeRight.push_back(e2); + } else { + if (!tradj->usesVertex(ip2)) osg::notify(osg::WARN) << "tradj error " << tradj->a()<< " , " << tradj->b()<< " , " << tradj->c()<< std::endl; + } + tradj=getTriangleWithEdge(e2,e1, &triangles); + } + if (trisToDelete.size()>=900) { + osg::notify(osg::WARN) << " found " << trisToDelete.size() << " adjacent tris " <::const_iterator deleteTri=trisToDelete.begin(); + deleteTri!=trisToDelete.end(); deleteTri++) + { + if (&(*tri)==*deleteTri) + { + deleted=true; + tri=triangles.erase(tri); + } + } + if (!deleted) ++tri; + } + } + } // strip test + + ip1=ip2; // next edge of line + } + } + } + } + // GWM Sept 2005 end + + // initialize index storage vector std::vector pt_indices; pt_indices.reserve(triangles.size() * 3); @@ -328,14 +978,17 @@ bool DelaunayTriangulator::triangulate() // build osg primitive osg::notify(osg::INFO) << "DelaunayTriangulator: building primitive(s)\n"; Triangle_list::const_iterator ti; - for (ti=triangles.begin(); ti!=triangles.end(); ++ti) { + for (ti=triangles.begin(); ti!=triangles.end(); ++ti) + { // don't add this triangle to the primitive if it shares any vertex with // the supertriangle // Also don't add degenerate (zero radius) triangles - if (ti->a() <= last_valid_index && ti->b() <= last_valid_index && ti->c() <= last_valid_index && ti->get_circumcircle().z()>0.0) { + if (ti->a() <= last_valid_index && ti->b() <= last_valid_index && ti->c() <= last_valid_index && ti->get_circumcircle().z()>0.0) + { - if (normals_.valid()) { + if (normals_.valid()) + { (normals_.get())->push_back(ti->compute_normal(points)); } @@ -347,10 +1000,211 @@ bool DelaunayTriangulator::triangulate() prim_tris_ = new osg::DrawElementsUInt(GL_TRIANGLES, pt_indices.size(), &(pt_indices.front())); - osg::notify(osg::INFO) << "DelaunayTriangulator: process done, " << (pt_indices.size() / 3) << " triangles created\n"; + osg::notify(osg::WARN) << "DelaunayTriangulator: process done, " << prim_tris_->getNumPrimitives() << " triangles remain\n"; return true; } +void DelaunayTriangulator::removeInternalTriangles(DelaunayConstraint *dc ) +{ + // Triangle_list *triangles + // remove triangle from terrain prim_tris_ internal to each constraintline + // and move to the constraint line to make an alternative geometry, + // possibly with alternative texture, and texture map + int ndel=0; + osg::Vec3Array::iterator normitr = normals_->begin(); + + // osg::notify(osg::WARN) << "DelaunayTriangulator: removeinternals, " << std::endl; + for (osg::DrawElementsUInt::iterator triit=prim_tris_->begin(); triit!=prim_tris_->end(); ) + { + // triangle joins points_[itr, itr+1, itr+2] + Triangle tritest((*triit), *(triit+1), *(triit+2), points_.get()); + if ((*triit==166 && *(triit+1)==162 && *(triit+2)==161) || + (*triit==166 && *(triit+1)==165 && *(triit+2)==164) ) + { + osg::Vec3 ctr=tritest.compute_centroid( points_.get()); + osg::notify(osg::WARN) << "testverts: " << ((*points_)[(*triit)].x()) << "," << ((*points_)[*(triit)].y()) <<","<<((*points_)[*(triit)].z())<windingNumber(ctr))<< std::endl; + } + if ( dc->contains(tritest.compute_centroid( points_.get()) ) ) + { + // centroid is inside the triangle, so IF inside linear, remove + // osg::notify(osg::WARN) << "DelaunayTriangulator: remove, " << (*triit) << "," << *(triit+1) <<","<<*(triit+2)<< std::endl; + dc->addtriangle((*triit), *(triit+1), *(triit+2)); + triit=prim_tris_->erase(triit); + triit=prim_tris_->erase(triit); + triit=prim_tris_->erase(triit); + if (normals_.valid()) + { + // erase the corresponding normal + normitr=normals_->erase(normitr); + } + ndel++; + } + else + { + if (normals_.valid()) + { + normitr++; + } + + triit+=3; + } + } + osg::notify(osg::INFO) << "end of test dc, deleted " << ndel << std::endl; } +//=== DelaunayConstraint functions + +float DelaunayConstraint::windingNumber(const osg::Vec3 testpoint) const +{ + // return winding number of loop around testpoint. Only in 2D, x-y coordinates assumed! + float theta=0; // sum of angles subtended by the line array - the winding number + const osg::Vec3Array *vertices= dynamic_cast(getVertexArray()); + for (unsigned int ipr=0; iprgetMode()==osg::PrimitiveSet::LINE_LOOP) + { + // nothing else loops + // start with the last point on the loop + const osg::Vec3 prev=(*vertices)[prset->index (prset->getNumIndices()-1)]; + osg::Vec3 pi(prev.x()-testpoint.x(),prev.y()-testpoint.y(),0); + pi.normalize(); + for (unsigned int i=0; igetNumIndices(); i++) + { + const osg::Vec3 curp=(*vertices)[prset->index (i)]; + osg::Vec3 edge(curp.x()-testpoint.x(),curp.y()-testpoint.y(),0); + edge.normalize(); + double cth=edge*pi; + if (cth<=-0.99999 ) + { + // testpoint is on edge and between 2 points + return 0; // + } + else + { + if (cth<0.99999) + { + float dang=(cth<1 && cth>-1)?acos(edge*pi):0; // this is unsigned angle + float zsign=edge.x()*pi.y()-pi.x()*edge.y(); // z component of..(edge^pi).z(); + if (zsign>0) theta+=dang; // add the angle subtended appropriately + else if (zsign<0) theta-=dang; + } + } + pi=edge; + } + } + } + + return theta/osg::PI/2.0; // should be 0 or 2 pi. +} +osg::DrawElementsUInt *DelaunayConstraint::makeDrawable() +{ + // initialize index storage vector for internal triangles. + std::vector pt_indices; + pt_indices.reserve(_interiorTris.size() * 3); + trilist::const_iterator ti; + for (ti=_interiorTris.begin(); ti!=_interiorTris.end(); ++ti) + { + + // if (normals_.valid()) { + // (normals_.get())->push_back(ti->compute_normal(points)); + // } + + pt_indices.push_back((*ti)[0]); + pt_indices.push_back((*ti)[1]); + pt_indices.push_back((*ti)[2]); + } + prim_tris_ = new osg::DrawElementsUInt(GL_TRIANGLES, pt_indices.size(), &(pt_indices.front())); + + return prim_tris_.get(); +} +bool DelaunayConstraint::contains(const osg::Vec3 testpoint) const +{ + // true if point is internal to the loop. + float theta=windingNumber(testpoint); // sum of angles subtended by the line array - the winding number + return fabs(theta)>0.9; // should be 0 or 1 (or 2,3,4 for very complex not permitted loops). +} +bool DelaunayConstraint::outside(const osg::Vec3 testpoint) const +{ + // true if point is outside the loop. + float theta=windingNumber(testpoint); // sum of angles subtended by the line array - the winding number + return fabs(theta)<.05; // should be 0 if outside. +} + + +void DelaunayConstraint::addtriangle(const int i1,const int i2, const int i3) +{ + // a triangle joins vertices i1,i2,i3 in the points of the delaunay triangles. + // points is the array of poitns in the triangulator; + // add triangle to the constraint + int *ip=new int[3]; + ip[0]=i1; + ip[1]=i2; + ip[2]=i3; + _interiorTris.push_back(ip); +} +osg::Vec3Array* DelaunayConstraint::getPoints(const osg::Vec3Array *points) +{ + //points_ is the array of points that can be used to render the triangles in this DC. + osg::ref_ptr points_=new osg::Vec3Array; + trilist::iterator ti; + for (ti=_interiorTris.begin(); ti!=_interiorTris.end(); ++ti) { + int idx=0; + int ip[3]={-1,-1,-1}; + // find if points[i1/i2/i3] already in the vertices points_ + for (osg::Vec3Array::iterator ivert=points_->begin(); ivert!=points_->end(); ivert++) + { + if (ip[0]<0 && *ivert==(*points)[(*ti)[0]]) + { + (*ti)[0]=ip[0]=idx; + } + if (ip[1]<0 && *ivert==(*points)[(*ti)[1]]) + { + (*ti)[1]=ip[1]=idx; + } + if (ip[2]<0 && *ivert==(*points)[(*ti)[2]]) + { + (*ti)[2]=ip[2]=idx; + } + idx++; + } + if (ip[0]<0) + { + points_->push_back((*points)[(*ti)[0]]); + (*ti)[0]=ip[0]=points_->size()-1; + } + if (ip[1]<0) + { + points_->push_back((*points)[(*ti)[1]]); + (*ti)[1]=ip[1]=points_->size()-1; + } + if (ip[2]<0) + { + points_->push_back((*points)[(*ti)[2]]); + (*ti)[2]=ip[2]=points_->size()-1; + } + } + makeDrawable(); + return points_.release(); +} + +void DelaunayConstraint::handleOverlaps(void) +{ + // use tesselator to interpolate crossing vertices. + osg::ref_ptr tscx=new osgUtil::Tesselator; // this assembles all the constraints + tscx->setTesselationType(osgUtil::Tesselator::TESS_TYPE_GEOMETRY); + tscx->setBoundaryOnly(true); + tscx->setWindingType( osgUtil::Tesselator::TESS_WINDING_ODD); + // ODD chooses the winding =1, NOT overlapping areas of constraints. + // nb this includes all the edges where delaunay constraints meet + // draw a case to convince yourself!. + + tscx->retesselatePolygons(*this); // find all edges +} + +} // namespace osgutil