root/OpenSceneGraph/trunk/src/osgUtil/DelaunayTriangulator.cpp @ 9630

Revision 9630, 55.7 kB (checked in by robert, 6 years ago)

Warnings fixes for VS.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
RevLine 
[5328]1/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2006 Robert Osfield
[4611]2 *
3 * This library is open source and may be redistributed and/or modified under 
4 * the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or
5 * (at your option) any later version.  The full license is in LICENSE file
6 * included with this distribution, and on the openscenegraph.org website.
7 *
8 * This library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11 * OpenSceneGraph Public License for more details.
12*/
13
[1889]14#include <osgUtil/DelaunayTriangulator>
[5928]15// NB this algorithm makes heavy use of the osgUtil::Tessellator for constrained triangulation.
[4611]16// truly it is built on the shoulders of giants.
[1889]17
18#include <osg/GL>
19#include <osg/Vec3>
20#include <osg/Array>
21#include <osg/Notify>
22
23#include <algorithm>
24#include <set>
[4611]25#include <map> //GWM July 2005 map is used in constraints.
[5985]26#include <osgUtil/Tessellator> // tessellator triangulates the constrained triangles
[1889]27
28namespace osgUtil
29{
30
31//////////////////////////////////////////////////////////////////////////////////////
32// MISC MATH FUNCTIONS
33
34
35// Compute the circumcircle of a triangle (only x and y coordinates are used),
36// return (Cx, Cy, r^2)
37inline osg::Vec3 compute_circumcircle(
38    const osg::Vec3 &a,
39    const osg::Vec3 &b,
40    const osg::Vec3 &c)
41{
42    float D =
43        (a.x() - c.x()) * (b.y() - c.y()) -
44        (b.x() - c.x()) * (a.y() - c.y());
45
[4358]46    float cx, cy, r2;
47
48    if(D==0.0)
49    {
50        // (Nearly) degenerate condition - either two of the points are equal (which we discount)
51        // or the three points are colinear. In this case we just determine the average of
52        // the three points as the centre for correctness, but squirt out a zero radius.
53        // This method will produce a triangulation with zero area, so we have to check later
54        cx = (a.x()+b.x()+c.x())/3.0;
55        cy = (a.y()+b.y()+c.y())/3.0;
56        r2 = 0.0;
57    }
58    else
59    {
60        cx =
[1889]61        (((a.x() - c.x()) * (a.x() + c.x()) +
62        (a.y() - c.y()) * (a.y() + c.y())) / 2 * (b.y() - c.y()) -
63        ((b.x() - c.x()) * (b.x() + c.x()) +
64        (b.y() - c.y()) * (b.y() + c.y())) / 2 * (a.y() - c.y())) / D;
65
[4358]66        cy =
[1889]67        (((b.x() - c.x()) * (b.x() + c.x()) +
68        (b.y() - c.y()) * (b.y() + c.y())) / 2 * (a.x() - c.x()) -
69        ((a.x() - c.x()) * (a.x() + c.x()) +
70        (a.y() - c.y()) * (a.y() + c.y())) / 2 * (b.x() - c.x())) / D;
71
[4611]72      //  r2 = (c.x() - cx) * (c.x() - cx) + (c.y() - cy) * (c.y() - cy);
73        // the return r square is compared with r*r many times in an inner loop
74        // so for efficiency use the inefficient sqrt once rather than 30* multiplies later.
75        r2 = sqrt((c.x() - cx) * (c.x() - cx) + (c.y() - cy) * (c.y() - cy));
[4358]76    }
[1889]77    return osg::Vec3(cx, cy, r2);
78}
79
80// Test whether a point (only the x and y coordinates are used) lies inside
[4611]81// a circle; the circle is passed as a vector: (Cx, Cy, r).
[1889]82
83inline bool point_in_circle(const osg::Vec3 &point, const osg::Vec3 &circle)
84{
85    float r2 =
86        (point.x() - circle.x()) * (point.x() - circle.x()) +
87        (point.y() - circle.y()) * (point.y() - circle.y());
[4611]88    return r2 <= circle.z()*circle.z();
89//    return r2 <= circle.z();
[1889]90}
91
92
93//
94//////////////////////////////////////////////////////////////////////////////////////
95
96
97// data type for vertex indices
98typedef GLuint Vertex_index;
99
100
101// CLASS: Edge
102// This class describes an edge of a triangle (it stores vertex indices to the two
103// endpoints).
104
105class Edge {
106public:
107
108    // Comparison object (for sorting)
[4611]109    struct Less
110    {
[1889]111        inline bool operator()(const Edge &e1, const Edge &e2) const
112        {
113            if (e1.ibs() < e2.ibs()) return true;
114            if (e1.ibs() > e2.ibs()) return false;
115            if (e1.ies() < e2.ies()) return true;
116            return false;
117        }
118    };
119
[9630]120    Edge(): ib_(0), ie_(0), ibs_(0), ies_(0), duplicate_(false) {}
[1950]121    Edge(Vertex_index ib, Vertex_index ie) : ib_(ib), ie_(ie), ibs_(osg::minimum(ib, ie)), ies_(osg::maximum(ib, ie)), duplicate_(false) {}
[1889]122
123    // first endpoint
124    inline Vertex_index ib() const { return ib_; }
125
126    // second endpoint
127    inline Vertex_index ie() const { return ie_; }
128
129    // first sorted endpoint
130    inline Vertex_index ibs() const { return ibs_; }
131
132    // second sorted endpoint
133    inline Vertex_index ies() const { return ies_; }
134
135    // get the "duplicate" flag
136    inline bool get_duplicate() const { return duplicate_; }
137
138    // set the "duplicate" flag
139    inline void set_duplicate(bool v) { duplicate_ = v; }
140
141private:
142    Vertex_index ib_, ie_;
143    Vertex_index ibs_, ies_;
144    bool duplicate_;
145};
146
147
148// CLASS: Triangle
149
[4611]150class Triangle
151{
[1889]152public:
[9630]153
154    Triangle():
155        a_(0),
156        b_(0),
157        c_(0) {}
158       
159   
[1889]160    Triangle(Vertex_index a, Vertex_index b, Vertex_index c, osg::Vec3Array *points)
161        :    a_(a),
162            b_(b),
163            c_(c),
164            cc_(compute_circumcircle((*points)[a_], (*points)[b_], (*points)[c_]))
165    {
166        edge_[0] = Edge(a_, b_);
167        edge_[1] = Edge(b_, c_);
168        edge_[2] = Edge(c_, a_);
169    }
170
[9630]171    Triangle& operator = (const Triangle& rhs)
172    {
173        if (&rhs==this) return *this;
174       
175        a_ = rhs.a_;
176        b_ = rhs.b_;
177        c_ = rhs.c_;
178        cc_ = rhs.cc_;
179        edge_[0]  = rhs.edge_[0];
180        edge_[1]  = rhs.edge_[1];
181        edge_[2]  = rhs.edge_[2];
182       
183        return *this;
184    }
185
[1889]186    inline Vertex_index a() const { return a_; }
187    inline Vertex_index b() const { return b_; }
188    inline Vertex_index c() const { return c_; }
[5810]189    inline void incrementa(const int delta) { a_+=delta; }
190    inline void incrementb(const int delta) { b_+=delta; }
191    inline void incrementc(const int delta) { c_+=delta; }
[1889]192
193    inline const Edge &get_edge(int idx) const { return edge_[idx];    }
194    inline const osg::Vec3 &get_circumcircle() const { return cc_; }
195
[4611]196    inline osg::Vec3 compute_centroid(const osg::Vec3Array *points) const
197    {
198        return ((*points)[a_] +(*points)[b_] + (*points)[c_])/3;
199    }
200
[1889]201    inline osg::Vec3 compute_normal(osg::Vec3Array *points) const
202    {
203        osg::Vec3 N = ((*points)[b_] - (*points)[a_]) ^ ((*points)[c_] - (*points)[a_]);
204        return N / N.length();
205    }
206
[4611]207    bool isedge(const unsigned int ip1,const unsigned int ip2) const
208    { // is one of the edges of this triangle from ip1-ip2
209        bool isedge=ip1==a() && ip2==b();
210        if (!isedge)
211        {
212            isedge=ip1==b() && ip2==c();
213            if (!isedge)
214            {
215                isedge=ip1==c() && ip2==a();
216            }
217        }
218        return isedge;
219    }
220    // GWM July 2005 add test for triangle intersected by p1-p2.
221    // return true for unused edge
[1889]222
[4611]223    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
224    {
225        // return true if edge iedge of triangle is intersected by ip1,ip2
226        Vertex_index ie1,ie2;
227        if (iedge==0)
228        {
229            ie1=a();
230            ie2=b();
231        }
232        else if (iedge==1)
233        {
234            ie1=b();
235            ie2=c();
236        }
237        else if (iedge==2)
238        {
239            ie1=c();
240            ie2=a();
241        }
242        if (ip1==ie1 || ip2==ie1) return false;
243        if (ip1==ie2 || ip2==ie2) return false;
244   
245        osg::Vec2 tp1((*points)[ie1].x(),(*points)[ie1].y());
246        osg::Vec2 tp2((*points)[ie2].x(),(*points)[ie2].y());
247        return intersect(tp1,tp2,p1,p2);
248    }
249   
250    bool intersectedby(const osg::Vec2 p1,const osg::Vec2 p2,osg::Vec3Array *points) const {
251        // true if line [p1,p2] cuts at least one edge of this triangle
252        osg::Vec2 tp1((*points)[a()].x(),(*points)[a()].y());
253        osg::Vec2 tp2((*points)[b()].x(),(*points)[b()].y());
254        osg::Vec2 tp3((*points)[c()].x(),(*points)[c()].y());
255        bool ip=intersect(tp1,tp2,p1,p2);
256        if (!ip)
257        {
258            ip=intersect(tp2,tp3,p1,p2);
259            if (!ip)
260            {
261                ip=intersect(tp3,tp1,p1,p2);
262            }
263        }
264        return ip;
265    }
266    int whichEdge(osg::Vec3Array *points,const osg::Vec2 p1, const osg::Vec2 p2,
267        const unsigned int e1,const unsigned int e2) const
268    {
269        int icut=0;
270        // find which edge of triangle is cut by line (p1-p2) and is NOT e1-e2 indices.
271        // return 1 - cut is on edge b-c; 2==c-a
272        osg::Vec2 tp1((*points)[a()].x(),(*points)[a()].y()); // triangle vertices
273        osg::Vec2 tp2((*points)[b()].x(),(*points)[b()].y());
274        osg::Vec2 tp3((*points)[c()].x(),(*points)[c()].y());
275        bool ip=intersect(tp2,tp3,p1,p2);
276        if (ip && (a()==e1 || a()==e2)) { return 1;}
277        ip=intersect(tp3,tp1,p1,p2);
278        if (ip && (b()==e1 || b()==e2)) { return 2;}
279        ip=intersect(tp1,tp2,p1,p2);
280        if (ip && (c()==e1 || c()==e2)) { return 3;}
281        return icut;
282    }
283
284    bool usesVertex(const unsigned int ip) const
285    {
286        return ip==a_ || ip==b_ || ip==c_;
287    }
288
[4616]289    int lineBisectTest(const osg::Vec3 apt,const osg::Vec3 bpt,const osg::Vec3 cpt, const osg::Vec2 p2) const
[4611]290    {
291        osg::Vec2 vp2tp=p2-osg::Vec2(apt.x(), apt.y()); // vector from p1 to a.
292        // test is: cross product (z component) with ab,ac is opposite signs
293        // AND dot product with ab,ac has at least one positive value.
294        osg::Vec2 vecba=osg::Vec2(bpt.x(), bpt.y())-osg::Vec2(apt.x(), apt.y());
295        osg::Vec2 vecca=osg::Vec2(cpt.x(), cpt.y())-osg::Vec2(apt.x(), apt.y());
296        float cprodzba=vp2tp.x()*vecba.y() - vp2tp.y()*vecba.x();
297        float cprodzca=vp2tp.x()*vecca.y() - vp2tp.y()*vecca.x();
[4616]298    //    osg::notify(osg::WARN) << "linebisect test " << " tri " << a_<<","<< b_<<","<< c_<<std::endl;
[4611]299        if (cprodzba*cprodzca<0)
300        {
301            // more tests - check dot products are at least partly parallel to test line.
302            osg::Vec2 tp1(bpt.x(),bpt.y()); // triangle vertices
303            osg::Vec2 tp2(cpt.x(),cpt.y());
304            osg::Vec2 tp3(apt.x(),apt.y());
305            bool ip=intersect(tp1,tp2,tp3,p2);
306            if (ip) return 1;
307        }
308        return 0;
309    }
310   
311    int lineBisects(osg::Vec3Array *points, const unsigned int ip1, const osg::Vec2 p2) const
312    {
313        // return true if line starts at vertex <ip1> and lies between the 2 edges which meet at vertex
314        // <vertex> is that which uses index ip1.
315        // line is <vertex> to p2
316        //    return value is 0 - no crossing; 1,2,3 - which edge of the triangle is cut.
317        if (a_==ip1)
318        {
319            // first vertex is the vertex - test that a_ to p2 lies beteen edges a,b and a,c
320            osg::Vec3 apt=(*points)[a_];
321            osg::Vec3 bpt=(*points)[b_];
322            osg::Vec3 cpt=(*points)[c_];
[4616]323            return lineBisectTest(apt,bpt,cpt,p2)?1:0;
[4611]324        }
325        else if (b_==ip1)
326        {
327            // second vertex is the vertex - test that b_ to p2 lies beteen edges a,b and a,c
328            osg::Vec3 apt=(*points)[b_];
329            osg::Vec3 bpt=(*points)[c_];
330            osg::Vec3 cpt=(*points)[a_];
[4616]331            return lineBisectTest(apt,bpt,cpt,p2)?2:0;
[4611]332        }
333        else if (c_==ip1)
334        {
335            // 3rd vertex is the vertex - test that c_ to p2 lies beteen edges a,b and a,c
336            osg::Vec3 apt=(*points)[c_];
337            osg::Vec3 bpt=(*points)[a_];
338            osg::Vec3 cpt=(*points)[b_];
[4616]339            return lineBisectTest(apt,bpt,cpt,p2)?3:0;
[4611]340        }
341        return 0;
342    }
343   
[1889]344private:
[4611]345
[9630]346
[4611]347    bool intersect(const osg::Vec2 p1,const osg::Vec2 p2,const osg::Vec2 p3,const osg::Vec2 p4) const
348    {
349        // intersection point of p1,p2 and p3,p4
350        // test from http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
351        // the intersection must be internal to the lines, not an end point.
352        float det=((p4.y()-p3.y())*(p2.x()-p1.x())-(p4.x()-p3.x())*(p2.y()-p1.y()));
353        if (det!=0)
354        {
355            // point on line is P=p1+ua.(p2-p1) and Q=p3+ub.(p4-p3)
356            float ua=((p4.x()-p3.x())*(p1.y()-p3.y())-(p4.y()-p3.y())*(p1.x()-p3.x()))/det;
357            float ub=((p2.x()-p1.x())*(p1.y()-p3.y())-(p2.y()-p1.y())*(p1.x()-p3.x()))/det;
358            if (ua> 0.00 && ua< 1 && ub> 0.0000  && ub< 1)
359            {
360                return true;
361            }
362        }
363        return false;
364    }
365   
[1889]366    Vertex_index a_;
367    Vertex_index b_;
368    Vertex_index c_;
369    osg::Vec3 cc_;   
370    Edge edge_[3];
371};
372
[4611]373typedef std::list<Triangle> Triangle_list;
[1889]374
375// comparison function for sorting sample points by the X coordinate
376bool Sample_point_compare(const osg::Vec3 &p1, const osg::Vec3 &p2)
377{
[4611]378    // replace pure sort by X coordinate with X then Y.
379    // errors can occur if the delaunay triangulation specifies 2 points at same XY and different Z
380    if (p1.x() != p2.x()) return p1.x() < p2.x();
381    if (p1.y() != p2.y()) return p1.y() < p2.y(); // GWM 30.06.05 - further rule if x coords are same.
382    osg::notify(osg::INFO) << "Two points are coincident at "<<p1.x() <<","<<p1.y() << std::endl;
383    return p1.z() < p2.z(); // never get here unless 2 points coincide
[1889]384}
385
386
387// container types
388typedef std::set<Edge, Edge::Less> Edge_set;
389
390
391DelaunayTriangulator::DelaunayTriangulator():
392    osg::Referenced()
393{
394}
395
396DelaunayTriangulator::DelaunayTriangulator(osg::Vec3Array *points, osg::Vec3Array *normals):
397    osg::Referenced(),
398    points_(points),
399    normals_(normals)
400{
401}
402
403DelaunayTriangulator::DelaunayTriangulator(const DelaunayTriangulator &copy, const osg::CopyOp &copyop):
404    osg::Referenced(copy),
405    points_(static_cast<osg::Vec3Array *>(copyop(copy.points_.get()))),
406    normals_(static_cast<osg::Vec3Array *>(copyop(copy.normals_.get()))),
407    prim_tris_(static_cast<osg::DrawElementsUInt *>(copyop(copy.prim_tris_.get())))
408{
409}
410
411DelaunayTriangulator::~DelaunayTriangulator()
412{
413}
414
[4611]415const Triangle * getTriangleWithEdge(const unsigned int ip1,const unsigned int ip2, const Triangle_list *triangles)
416{ // find triangle in list with edge from ip1 to ip2...
417    Triangle_list::const_iterator trconnitr; // connecting triangle
418    int idx=0;
419    for (trconnitr=triangles->begin(); trconnitr!=triangles->end(); )
420    {
421        if (trconnitr->isedge (ip1,ip2))
422        {
423            // this is the triangle.
424            return &(*trconnitr);
425        }
426        ++trconnitr;
427        idx++;
428    }
429    return NULL; //-1;
430}
431
[6446]432int DelaunayTriangulator::getindex(const osg::Vec3 &pt,const osg::Vec3Array *points)
[4611]433{
434    // return index of pt in points (or -1)
435    for (unsigned int i=0; i<points->size(); i++)
436    {
437        if (pt.x()==(*points)[i].x() &&pt.y()==(*points)[i].y() )
438        {
439            return i;
440        }
441    }
442    return -1;
443}
444
445Triangle_list fillHole(osg::Vec3Array *points,    std::vector<unsigned int> vindexlist)
446{
447    // eg clockwise vertex neighbours around the hole made by the constraint
448    Triangle_list triangles; // returned list
449    osg::ref_ptr<osg::Geometry> gtess=new osg::Geometry; // add all the contours to this for analysis
450    osg::ref_ptr<osg::Vec3Array> constraintverts=new osg::Vec3Array;
[5928]451    osg::ref_ptr<osgUtil::Tessellator> tscx=new osgUtil::Tessellator; // this assembles all the constraints
[4611]452   
453    for (std::vector<unsigned int>::iterator itint=vindexlist.begin(); itint!=vindexlist.end(); itint++)
454    {
455    //    osg::notify(osg::WARN)<< "add point " << (*itint) << " at " << (*points)[*itint].x() << ","<< (*points)[*itint].y() <<std::endl;
456        constraintverts->push_back((*points)[*itint]);
457    }
458   
459    unsigned int npts=vindexlist.size();
460
461    gtess->setVertexArray(constraintverts.get());
462    gtess->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::POLYGON,0,npts));
[5928]463    tscx->setTessellationNormal(osg::Vec3(0.0,0.0,1.0));
464    tscx->setTessellationType(osgUtil::Tessellator::TESS_TYPE_GEOMETRY);
[4611]465    tscx->setBoundaryOnly(false);
[5940]466    tscx->setWindingType( osgUtil::Tessellator::TESS_WINDING_ODD); // the commonest tessellation is default, ODD. GE2 allows intersections of constraints to be found.
[5928]467    tscx->retessellatePolygons(*(gtess.get())); // this should insert extra vertices where constraints overlap
[4611]468
469    // extract triangles from gtess
470   
471    unsigned int ipr;
472    for (ipr=0; ipr<gtess->getNumPrimitiveSets(); ipr++)
473    {
474        unsigned int ic;
475        osg::PrimitiveSet* prset=gtess->getPrimitiveSet(ipr);
476        //                    osg::notify(osg::WARN)<< "gtess set " << ipr << " nprims " << prset->getNumPrimitives() <<
477        //                        " type " << prset->getMode() << std::endl;
478        unsigned int pidx,pidx1,pidx2;
479        switch (prset->getMode()) {
480        case osg::PrimitiveSet::TRIANGLES:
481            for (ic=0; ic<prset->getNumIndices()-2; ic+=3)
482            {
483                if (prset->index(ic)>=npts)
484                {
485                    // this is an added point.
486                    points->push_back((*constraintverts)[prset->index(ic)]);
[5810]487                    pidx=points->size()-1;
[4611]488                }
489                else
490                {
491                    pidx=vindexlist[prset->index(ic)];
492                }
493               
494                if (prset->index(ic+1)>=npts)
495                {
496                    // this is an added point.
497                    points->push_back((*constraintverts)[prset->index(ic+1)]);
[5810]498                    pidx1=points->size()-1;
[4611]499                }
500                else
501                {
502                    pidx1=vindexlist[prset->index(ic+1)];
503                }
504               
505                if (prset->index(ic+2)>=npts)
506                {
507                    // this is an added point.
508                    points->push_back((*constraintverts)[prset->index(ic+2)]);
[5810]509                    pidx2=points->size()-1;
[4611]510                }
511                else
512                {
513                    pidx2=vindexlist[prset->index(ic+2)];
514                }
515                triangles.push_back(Triangle(pidx, pidx1, pidx2, points));
516                //                    osg::notify(osg::WARN)<< "vert " << prset->index(ic) << " in array"<<std::endl;
517            }
518            break;
519        case osg::PrimitiveSet::TRIANGLE_STRIP: // 123, 234, 345...
520
521            for (ic=0; ic<prset->getNumIndices()-2; ic++)
522            {
523                if (prset->index(ic)>=npts)
524                {
525                    // this is an added point.
526                    points->push_back((*constraintverts)[prset->index(ic)]);
[5810]527                    pidx=points->size()-1;
[4611]528                } else {
529                    pidx=vindexlist[prset->index(ic)];
530                }
531                if (prset->index(ic+1)>=npts)
532                {
533                    // this is an added point.
534                    points->push_back((*constraintverts)[prset->index(ic+1)]);
[5810]535                    pidx1=points->size()-1;
[4611]536                }
537                else
538                {
539                    pidx1=vindexlist[prset->index(ic+1)];
540                }
541               
542                if (prset->index(ic+2)>=npts)
543                {
544                    // this is an added point.
545                    points->push_back((*constraintverts)[prset->index(ic+2)]);
[5810]546                    pidx2=points->size()-1;
[4611]547                }
548                else 
549                {
550                    pidx2=vindexlist[prset->index(ic+2)];
551                }
552               
553                if (ic%2==0)
554                {
555                    triangles.push_back(Triangle(pidx, pidx1, pidx2, points));
556                }
557                else
558                {
559                    triangles.push_back(Triangle(pidx1, pidx, pidx2, points));
560                }
561                //                    osg::notify(osg::WARN)<< "vert " << prset->index(ic) << " in array"<<std::endl;
562            }
563            break;
564           
565        case osg::PrimitiveSet::TRIANGLE_FAN:
566            {
567                osg::Vec3 ptest=(*constraintverts)[prset->index(0)];
568                if (prset->index(0)>=npts)
569                {
570                    // this is an added point.
571                    points->push_back((*constraintverts)[prset->index(0)]);
[5810]572                    pidx=points->size()-1;
[4611]573                }
574                else
575                {
576                    pidx=vindexlist[prset->index(0)];
577                }
578                //        osg::notify(osg::WARN)<< "tfan has " << prset->getNumIndices() << " indices"<<std::endl;
579                for (ic=1; ic<prset->getNumIndices()-1; ic++)
580                {
581                    if (prset->index(ic)>=npts)
582                    {
583                        // this is an added point.
584                        points->push_back((*constraintverts)[prset->index(ic)]);
[5810]585                        pidx1=points->size()-1;
[4611]586                    }
587                    else
588                    {
589                        pidx1=vindexlist[prset->index(ic)];
590                    }
591                   
592                    if (prset->index(ic+1)>=npts)
593                    { // this is an added point.
594                        points->push_back((*constraintverts)[prset->index(ic+1)]);
[5810]595                        pidx2=points->size()-1;
[4611]596                    }
597                    else
598                    {
599                        pidx2=vindexlist[prset->index(ic+1)];
600                    }
601                    triangles.push_back(Triangle(pidx, pidx1, pidx2, points));
602                }
603            }
604            break;
605        default:
606            osg::notify(osg::WARN)<< "WARNING set " << ipr << " nprims " << prset->getNumPrimitives() <<
607                " type " << prset->getMode() << " Type not triangle, tfan or strip"<< std::endl;
608            break;
609        }
610    }
611    return triangles;
612}
613
614void DelaunayConstraint::removeVerticesInside(const DelaunayConstraint *dco)
615{    /** remove vertices from this which are internal to dco.
616     * retains potins that are extremely close to edge of dco
617      * defined as edge of dco subtends>acs(0.999999)
618    */
619    int nrem=0;
620    osg::Vec3Array *vertices= dynamic_cast< osg::Vec3Array*>(getVertexArray());
[4728]621    if (vertices)
[4611]622    {
[4728]623        for (osg::Vec3Array::iterator vitr=vertices->begin(); vitr!=vertices->end(); )
[4611]624        {
[4728]625            if (dco->contains(*vitr))
[4611]626            {
[4728]627                unsigned int idx=vitr-vertices->begin(); // index of vertex
628                // remove vertex index from all the primitives
629                for (unsigned int ipr=0; ipr<getNumPrimitiveSets(); ipr++)
630                {
631                    osg::PrimitiveSet* prset=getPrimitiveSet(ipr);
632                    osg::DrawElementsUShort *dsup=dynamic_cast<osg::DrawElementsUShort *>(prset);
633                    if (dsup) {
634                        for (osg::DrawElementsUShort::iterator usitr=dsup->begin(); usitr!=dsup->end(); )
[4611]635                        {
[4728]636                            if ((*usitr)==idx)
637                            { // remove entirely
638                                usitr=dsup->erase(usitr);
639                            }
640                            else
641                            {
642                                if ((*usitr)>idx) (*usitr)--; // move indices down 1
643                                usitr++; // next index
644                            }
[4611]645                        }
646                    }
[4728]647                    else
648                    {
649                        osg::notify(osg::WARN) << "Invalid prset " <<ipr<< " tp " << prset->getType() << " types PrimitiveType,DrawArraysPrimitiveType=1 etc" << std::endl;
650                    }
[4611]651                }
[4728]652                vitr=vertices->erase(vitr);
653                nrem++;
654
[4611]655            }
[4728]656            else
657            {
658                vitr++;
659            }
[4611]660        }
661    }
662}
663
664void DelaunayConstraint::merge(DelaunayConstraint *dco)
665{
666    unsigned int ipr;
[5949]667    if (dco) { // 16 Dec 2006 just in case you pass in a NULL pointer
668        osg::Vec3Array* vmerge=dynamic_cast<osg::Vec3Array*>(getVertexArray());
669        if (!vmerge) vmerge=new osg::Vec3Array;
670        setVertexArray(vmerge);
671        for ( ipr=0; ipr<dco->getNumPrimitiveSets(); ipr++)
[4611]672        {
[5949]673            osg::PrimitiveSet* prset=dco->getPrimitiveSet(ipr);
674            osg::DrawArrays *drarr=dynamic_cast<osg::DrawArrays *> (prset);
675            if (drarr)
676            {
677                // need to add the offset of vmerge->size to each prset indices.
678                unsigned int noff=vmerge->size();
679                unsigned int n1=drarr->getFirst(); // usually 0
680                unsigned int numv=drarr->getCount(); //
681                addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,n1+noff,numv));
682            }
[4611]683        }
[5949]684        osg::Vec3Array* varr=dynamic_cast<osg::Vec3Array*>(dco->getVertexArray());
685        if (varr) vmerge->insert(vmerge->end(),varr->begin(),varr->end());
[4611]686    }
687}
688
[4984]689void DelaunayTriangulator::_uniqueifyPoints()
690{
691    std::sort( points_->begin(), points_->end() );
[4611]692
[4984]693    osg::ref_ptr<osg::Vec3Array> temppts = new osg::Vec3Array;
694    // This won't work... must write our own unique() that compares only the first
695    // two terms of a Vec3 for equivalency
696    //std::insert_iterator< osg::Vec3Array > ti( *(temppts.get()), temppts->begin() );
697    //std::unique_copy( points_->begin(), points_->end(), ti );
698
699    osg::Vec3Array::iterator p = points_->begin();
700    osg::Vec3 v = *p;
[5075]701    // Always push back the first point
702    temppts->push_back( (v = *p));
[4984]703    for( ; p != points_->end(); p++ )
704    {
705        if( v[0] == (*p)[0] && v[1] == (*p)[1] )
706            continue;
707
708        temppts->push_back( (v = *p));
709    }
710
711    points_->clear();
712    std::insert_iterator< osg::Vec3Array > ci(*(points_.get()),points_->begin());
713    std::copy( temppts->begin(), temppts->end(), ci );
714}
715
[5985]716osgUtil::DelaunayConstraint *getconvexhull(osg::Vec3Array *points)
717{ // fits the 'rubberband' around the 2D points for uses as a delaunay constraint.
718    osg::ref_ptr<osgUtil::DelaunayConstraint> dcconvexhull=new osgUtil::DelaunayConstraint; // make
719    // convex hull around all the points
720    // start from first point (at minx); proceed to last x and back
721    osg::Vec3Array *verts=new osg::Vec3Array; // the hull points
722    verts->push_back(*(points->begin()) ); // min x/y point is guaranteed to be on the hull
723    verts->push_back(*(points->begin()+1) ); // second low x/y point is first length to be tested
724    for (osg::Vec3Array::iterator vit=(points->begin()+2); vit!=points->end(); vit++) {
725        // check if point lies outside the current last line segment
726        bool ok=1;
727        while (ok && verts->size()>1) {
728            osg::Vec3 lastseg=*(verts->end()-2)-*(verts->end()-1);
729            osg::Vec3 thisseg=(*vit)-*(verts->end()-1);
730            float cosang=(lastseg^thisseg).z();
731            if (cosang <0.0) { // pop off last point - *vit is further out hull
732                verts->pop_back();
733            } else { ok=0;}
734        }
735        verts->push_back(*vit ); // next low x/y point is next length to be tested
736        // check if the previous external angle is >180 - then remove previous point
737    }
738    for (osg::Vec3Array::reverse_iterator rvit=points->rbegin()+1; rvit!=points->rend(); rvit++) {
739        // check if point lies outside the current last line segment
740        bool ok=1;
741        while (ok && verts->size()>1) {
742            osg::Vec3 lastseg=*(verts->end()-2)-*(verts->end()-1);
743            osg::Vec3 thisseg=(*rvit)-*(verts->end()-1);
744            float cosang=(lastseg^thisseg).z();
745            if (cosang <0.0) { // pop off last point - *rvit is further out hull
746                verts->pop_back();
747            } else { ok=0;}
748        }
749        if ((*rvit)!=*(verts->begin())) verts->push_back(*rvit ); // next low x/y point is next length to be tested
750        // check if the previous external angle is >180 - then remove previous point
751    }
752    dcconvexhull->setVertexArray(verts);
753    dcconvexhull->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP,0,verts->size()) );
754    return dcconvexhull.release();
755}
[4984]756
[1889]757bool DelaunayTriangulator::triangulate()
758{
759    // check validity of input array
[4611]760    if (!points_.valid())
761    {
[1889]762        osg::notify(osg::WARN) << "Warning: DelaunayTriangulator::triangulate(): invalid sample point array" << std::endl;
763        return false;
764    }
765
766    osg::Vec3Array *points = points_.get();
767
[4611]768    if (points->size() < 1)
769    {
[1889]770        osg::notify(osg::WARN) << "Warning: DelaunayTriangulator::triangulate(): too few sample points" << std::endl;
771        return false;
772    }
773
[4984]774    // Eliminate duplicate lat/lon points from input coordinates.
775    _uniqueifyPoints();
776
777
[1889]778    // initialize storage structures
779    Triangle_list triangles;
780    Triangle_list discarded_tris;   
781
[4611]782    // GWM July 2005 add constraint vertices to terrain
783    linelist::iterator linitr;
784    for (linitr=constraint_lines.begin();linitr!=constraint_lines.end();linitr++)
785    {
[4728]786        DelaunayConstraint* dc=(*linitr).get();
787        const osg::Vec3Array* vercon= dynamic_cast<const osg::Vec3Array*>(dc->getVertexArray());
788        if (vercon)
[4611]789        {
[4728]790            int nadded=0;
791            for (unsigned int icon=0;icon<vercon->size();icon++)
[4611]792            {
[4728]793                osg::Vec3 p1=(*vercon)[icon];
794                int idx=getindex(p1,points_.get());
795                if (idx<0)
796                { // only unique vertices are permitted.
797                    points_->push_back(p1); // add non-unique constraint points to triangulation
798                    nadded++;
799                }
800                else
801                {
802                    osg::notify(osg::WARN) << "DelaunayTriangulator: ignore a duplicate point at "<< p1.x()<< " " << p1.y() << std::endl;;
803                }
[4611]804            }
805        }
806    //    osg::notify(osg::WARN)<< "constraint size "<<vercon->size()<<" " <<nadded<< std::endl;
807    }
808        // GWM July 2005 end
809
[1889]810    // pre-sort sample points
811    osg::notify(osg::INFO) << "DelaunayTriangulator: pre-sorting sample points\n";
812    std::sort(points->begin(), points->end(), Sample_point_compare);   
[5985]813    // 24.12.06 add convex hull of points to force sensible outline.
814    osg::ref_ptr<osgUtil::DelaunayConstraint> dcconvexhull=getconvexhull(points);
815    addInputConstraint(dcconvexhull.get());
[1889]816
817    // set the last valid index for the point list
818    GLuint last_valid_index = points->size() - 1;   
819
820    // find the minimum and maximum x values in the point list
821    float minx = (*points)[0].x();
[1895]822    float maxx = (*points)[last_valid_index].x();
[1889]823
[5985]824    // find the minimum and maximum y values in the point list   
[1889]825    float miny = (*points)[0].y();
826    float maxy = miny;
[1895]827   
[1889]828    osg::notify(osg::INFO) << "DelaunayTriangulator: finding minimum and maximum Y values\n";
829    osg::Vec3Array::const_iterator mmi;
[4611]830    for (mmi=points->begin(); mmi!=points->end(); ++mmi)
831    {
[1889]832        if (mmi->y() < miny) miny = mmi->y();
833        if (mmi->y() > maxy) maxy = mmi->y();
834    }
835   
836    // add supertriangle vertices to the point list
[4611]837    // gwm added 1.05* to ensure that supervertices are outside the domain of points.
838    // original value could make 2 coincident points for regular arrays of x,y,h data.
839    // this mod allows regular spaced arrays to be used.
[5949]840    // 16 Dec 2006 this increase in size encourages the supervertex triangles to be long and thin
841    // thus ensuring that the convex hull of the terrain points are edges in the delaunay triangulation
842    // the values do however result in a small loss of numerical resolution.
[5985]843    points_->push_back(osg::Vec3(minx - .10*(maxx - minx), miny - .10*(maxy - miny), 0));
844    points_->push_back(osg::Vec3(maxx + .10*(maxx - minx), miny - .10*(maxy - miny), 0));
845    points_->push_back(osg::Vec3(maxx + .10*(maxx - minx), maxy + .10*(maxy - miny), 0));
846    points_->push_back(osg::Vec3(minx - .10*(maxx - minx), maxy + .10*(maxy - miny), 0));
847                                                                   
848    // add supertriangles to triangle list                           
[1889]849    triangles.push_back(Triangle(last_valid_index+1, last_valid_index+2, last_valid_index+3, points));
[5985]850    triangles.push_back(Triangle(last_valid_index+4, last_valid_index+1, last_valid_index+3, points));
[1889]851
852   
853    // begin triangulation   
854    GLuint pidx = 0;
855    osg::Vec3Array::const_iterator i;   
856   
857    osg::notify(osg::INFO) << "DelaunayTriangulator: triangulating vertex grid (" << (points->size()-3) <<" points)\n";   
858
[4611]859    for (i=points->begin(); i!=points->end(); ++i, ++pidx)
860    {
[1889]861
862        // don't process supertriangle vertices
863        if (pidx > last_valid_index) break;
864
865        Edge_set edges;
866
867        // iterate through triangles
868        Triangle_list::iterator j, next_j;
[4611]869        for (j=triangles.begin(); j!=triangles.end(); j = next_j)
870        {
[1889]871
872            next_j = j;
873            ++next_j;
874
[5985]875            // get the circumcircle (x,y centre & radius)
[1889]876            osg::Vec3 cc = j->get_circumcircle();
877
878            // OPTIMIZATION: since points are pre-sorted by the X component,
879            // check whether we can discard this triangle for future operations
880            float xdist = i->x() - cc.x();
[4611]881            // this is where the circumcircles radius rather than R^2 is faster.
882            // original code used r^2 and needed to test xdist*xdist>cc.z && i->x()>cc.x().
883            if ((xdist ) > cc.z() )
884            {
[5985]885                discarded_tris.push_back(*j); // these are not needed for further tests as no more
886                // points will ever lie inside this triangle.
[1889]887                triangles.erase(j);
[4611]888            }
889            else
890            {
[1889]891
892                // if the point lies in the triangle's circumcircle then add
893                // its edges to the edge list and remove the triangle
894                if (point_in_circle(*i, cc))
895                {
896                    for (int ei=0; ei<3; ++ei)
897                    {
898                        std::pair<Edge_set::iterator, bool> result = edges.insert(j->get_edge(ei));
899                        if (!result.second)
900                        {
901                            // cast away constness of a set element, which is
902                            // safe in this case since the set_duplicate is
903                            // not used as part of the Less operator.
904                            Edge& edge = const_cast<Edge&>(*(result.first));
[4611]905                            // not clear why this change is needed? But prevents removal of twice referenced edges??
906                      //      edge.set_duplicate(true);
907                            edge.set_duplicate(!edge.get_duplicate());
[1889]908                        }
909                    }                   
910                    triangles.erase(j);
911                }
912            }
913        }
914
915        // remove duplicate edges and add new triangles
916        Edge_set::iterator ci;
[4611]917        for (ci=edges.begin(); ci!=edges.end(); ++ci)
918        {
919            if (!ci->get_duplicate())
920            {
[1889]921                triangles.push_back(Triangle(pidx, ci->ib(), ci->ie(), points));
922            }
923        }
924    }
[5810]925    // dec 2006 we used to remove supertriangle vertices here, but then we cant strictly use the supertriangle
926    // vertices to find intersections of constraints with terrain, so moved to later.
[1889]927
928    osg::notify(osg::INFO) << "DelaunayTriangulator: finalizing and cleaning up structures\n";
929
[5810]930     // rejoin the two triangle lists
[1889]931    triangles.insert(triangles.begin(), discarded_tris.begin(), discarded_tris.end());
932
[4611]933        // GWM July 2005 eliminate any triangle with an edge crossing a constraint line
934    // http://www.geom.uiuc.edu/~samuelp/del_project.html
935    // we could also implement the sourcecode in http://gts.sourceforge.net/reference/gts-delaunay-and-constrained-delaunay-triangulations.html
936    // this uses the set of lines which are boundaries of the constraints, including points
[5940]937    // added to the contours by tessellation.
[4611]938    for (linelist::iterator dcitr=constraint_lines.begin();dcitr!=constraint_lines.end();dcitr++)
939    {
940        //DelaunayConstraint *dc=(*dcitr).get();
[4728]941        const osg::Vec3Array* vercon = dynamic_cast<const osg::Vec3Array*>((*dcitr)->getVertexArray());
942        if (vercon)
[4611]943        {
[4728]944            for (unsigned int ipr=0; ipr<(*dcitr)->getNumPrimitiveSets(); ipr++)
[4611]945            {
[4728]946                const osg::PrimitiveSet* prset=(*dcitr)->getPrimitiveSet(ipr);
947                if (prset->getMode()==osg::PrimitiveSet::LINE_LOOP ||
948                    prset->getMode()==osg::PrimitiveSet::LINE_STRIP)
[4611]949                {
[4728]950                    // loops or strips
951                    // start with the last point on the loop
952                    unsigned int ip1=getindex((*vercon)[prset->index (prset->getNumIndices()-1)],points_.get());
953                    for (unsigned int i=0; i<prset->getNumIndices(); i++)
[4611]954                    {
[4728]955                        unsigned int ip2=getindex((*vercon)[prset->index(i)],points_.get());
956                        if (i>0 || prset->getMode()==osg::PrimitiveSet::LINE_LOOP)
[4611]957                        {
[4728]958                            // dont check edge from end to start
959                            // for strips
960                            // 2 points on the constraint
961                            bool edgused=false;// first check for exact edge indices are used.
962                            Triangle_list::iterator titr;
963                            const osg::Vec3 curp=(*vercon)[prset->index(i)];
964                            int it=0;
965                            for (titr=triangles.begin(); titr!=triangles.end() && !edgused; ++titr)
[4611]966                            {
[4728]967                                //check that the edge ip1-ip2 is not already part of the triangulation.
968                                if (titr->isedge(ip1,ip2)) edgused=true;
969                                if (titr->isedge(ip2,ip1)) edgused=true;
970                                //        if (edgused) osg::notify(osg::WARN) << "Edge used in triangle " << it << " " <<
971                                //            titr->a()<<","<< titr->b()<<","<< titr->c()<<  std::endl;
972                                it++;
973                            }
974                            if (!edgused)
975                            {
976                                // then check for intermediate triangles, erase them and replace with constrained triangles.
977                                // find triangle with point ip1 where the 2 edges from ip1 contain the line p1-p2.
978                                osg::Vec2 p1((*points_)[ip1].x(),(*points_)[ip1].y()); // a constraint line joins p1-p2
979                                osg::Vec2 p2((*points_)[ip2].x(),(*points_)[ip2].y());
980                                int ntr=0;
981                                std::vector<const Triangle *> trisToDelete; // array of triangles to delete from terrain.
982                                // form 2 lists of vertices for the edges of the hole created.
983                                // The hole joins vertex ip1 to ip2, and one list of edges lies to the left
984                                // of the line ip1-ip2m the other to the right.
985                                // a list of vertices forming 2 halves of the removed triangles.
[5940]986                                // which in turn are filled in with the tessellator.
[4728]987                                for (titr=triangles.begin(); titr!=triangles.end(); )
[4611]988                                {
[4728]989                                    int icut=titr->lineBisects(points_.get(),ip1,p2);
990                                    //    osg::notify(osg::WARN) << "Testing triangle " << ntr << " "<< ip1 << " ti " <<
991                                    //        titr->a()<< ","<<titr->b() <<"," <<titr->c() << std::endl;
992                                    if (icut>0)
[4611]993                                    {
[4728]994                                        // triangle titr starts the constraint edge
995                                        std::vector<unsigned int> edgeRight, edgeLeft;
996                                        edgeRight.push_back(ip1);
997                                        edgeLeft.push_back(ip1);
998                                        //        osg::notify(osg::WARN) << "hole first " << edgeLeft.back()<<  " rt " << edgeRight.back()<< std::endl;
999                                        trisToDelete.push_back(&(*titr));
1000                                        // now find the unique triangle that shares the defined edge
1001                                        unsigned int e1, e2; // indices of ends of test triangle titr
1002                                        if    (icut==1)
[4611]1003                                        {
[4728]1004                                            // icut=1 implies vertex a is not involved
1005                                            e1=titr->b(); e2=titr->c();
1006                                        }
1007                                        else if (icut==2)
1008                                        {
1009                                            e1=titr->c(); e2=titr->a();
1010                                        }
1011                                        else if (icut==3)
1012                                        {
1013                                            e1=titr->a(); e2=titr->b();
1014                                        }
1015                                        edgeRight.push_back(e2);
1016                                        edgeLeft.push_back(e1);
1017                                        //        osg::notify(osg::WARN) << icut << "hole edges " << edgeLeft.back()<<  " rt " << edgeRight.back()<< std::endl;
1018                                        const Triangle *tradj=getTriangleWithEdge(e2,e1, &triangles);
1019                                        if (tradj)
1020                                        {
1021                                            while (tradj && !tradj->usesVertex(ip2) && trisToDelete.size()<999)
1022                                            {
1023                                                trisToDelete.push_back(tradj);
1024                                                icut=tradj->whichEdge(points_.get(),p1,p2,e1,e2);
1025                                                //    osg::notify(osg::WARN)  << ntr << " cur triedge " << icut << " " << ip1 <<
1026                                                //        " to " << ip2 << " tadj " << tradj->a()<< ","<<tradj->b() <<","
1027                                                //        <<tradj->c() <<std::endl;
1028                                                if        (icut==1) {e1=tradj->b(); e2=tradj->c();} // icut=1 implies vertex a is not involved
1029                                                else if (icut==2) {e1=tradj->c(); e2=tradj->a();}
1030                                                else if (icut==3) {e1=tradj->a(); e2=tradj->b();}
1031                                                if (edgeLeft.back()!=e1 && edgeRight.back()==e2 && e1!=ip2) {
1032                                                    edgeLeft.push_back(e1);
1033                                                } else if(edgeRight.back()!=e2 && edgeLeft.back()==e1 && e2!=ip2) {
1034                                                    edgeRight.push_back(e2);
1035                                                } else {
1036                                                    if (!tradj->usesVertex(ip2)) osg::notify(osg::WARN) << "tradj error " << tradj->a()<<  " , " << tradj->b()<<  " , " << tradj->c()<< std::endl;
1037                                                }
1038                                                tradj=getTriangleWithEdge(e2,e1, &triangles);
[4611]1039                                            }
[4728]1040                                            if (trisToDelete.size()>=900) {
1041                                                osg::notify(osg::WARN) << " found " << trisToDelete.size() << " adjacent tris " <<std::endl;
1042                                            }
[4611]1043                                        }
[4728]1044
1045                                        // both lines end at ip2 point.
1046                                        edgeLeft.push_back(ip2);
1047                                        edgeRight.push_back(ip2);
1048                                        if (tradj) trisToDelete.push_back(tradj);
1049                                        //        osg::notify(osg::WARN) << icut << "hole last " << edgeLeft.back()<<  " rt " << edgeRight.back()<< std::endl;
1050                                        Triangle_list constrainedtris=fillHole(points_.get(),edgeLeft);
1051                                        triangles.insert(triangles.begin(), constrainedtris.begin(), constrainedtris.end());
1052                                        constrainedtris=fillHole(points_.get(),edgeRight);
1053                                        triangles.insert(triangles.begin(), constrainedtris.begin(), constrainedtris.end());
1054
[4611]1055                                    }
[4728]1056                                    ++titr;
1057                                    ntr++;
[4611]1058                                }
[4728]1059                                // remove the triangles list
1060                                Triangle_list::iterator tri; // counts through triangles
1061                                for (tri=triangles.begin(); tri!=triangles.end(); )
[4611]1062                                {
[4728]1063                                    bool deleted=false;
[5985]1064                                    for (std::vector<const Triangle *>::iterator deleteTri=trisToDelete.begin();
1065                                    deleteTri!=trisToDelete.end(); )
[4611]1066                                    {
[4728]1067                                        if (&(*tri)==*deleteTri)
1068                                        {
1069                                            deleted=true;
1070                                            tri=triangles.erase(tri);
[5985]1071                                            deleteTri=trisToDelete.erase(deleteTri); // 24.12.06 remove from delete list.
1072                                        } else {deleteTri++; }
[4611]1073                                    }
[4728]1074                                    if (!deleted) ++tri;
[4611]1075                                }
1076                            }
[4728]1077                        } // strip test
1078
1079                        ip1=ip2; // next edge of line
1080                    }
[4611]1081                }
1082            }
1083        }
1084    }
1085    // GWM Sept 2005 end
[5810]1086    // dec 2006 remove supertriangle vertices - IF we have added some internal vertices (see fillholes)
1087    // then these may not be the last vertices in the list.
1088//    } else { // remove 3 super-triangle vertices more completely, moving any reference indices down.
1089    Triangle_list::iterator tri;
[5985]1090    GLuint supertriend = last_valid_index+4;
[5810]1091    for (tri=triangles.begin(); tri!=triangles.end();)
1092    { // look for triangles using the supertriangle indices or >super indices and move down by 3.
1093        if ((tri->a() > last_valid_index && tri->a() <= supertriend) ||
1094            (tri->b() > last_valid_index && tri->b() <= supertriend ) ||
1095            (tri->c() > last_valid_index && tri->c() <= supertriend )) {
1096            tri=triangles.erase(tri);
1097        } else { // move down by 3 tests
[5985]1098            if (tri->a() > last_valid_index) { // move index down 4
1099                tri->incrementa(-4);
[5810]1100            }
[5985]1101            if (tri->b() > last_valid_index) { // move down 4
1102                tri->incrementb(-4);
[5810]1103            }
[5985]1104            if (tri->c() > last_valid_index) { // move down 4 -- correction 21.12.06- was b() should test index c()
1105                tri->incrementc(-4);
[5810]1106            }
[5985]1107            ++tri; // only increment tri here as the other path increments tri when tri=triangles.erase.
[5810]1108        }
1109    }
1110        // remove 3 supertriangle vertices from points. They may not be the last vertices in points if
1111        // extra points have been inserted by the constraint re-triangulation.
[5985]1112    points->erase(points->begin()+last_valid_index+1,points->begin()+last_valid_index+5);
[5810]1113    //last_valid_index = points->size()-1; // the last valid vertex is last point.
1114    // end of dec 2006 changes.
[4611]1115
[1889]1116    // initialize index storage vector
1117    std::vector<GLuint> pt_indices;
1118    pt_indices.reserve(triangles.size() * 3);
1119
1120    // build osg primitive
1121    osg::notify(osg::INFO) << "DelaunayTriangulator: building primitive(s)\n";
1122    Triangle_list::const_iterator ti;
[4611]1123    for (ti=triangles.begin(); ti!=triangles.end(); ++ti)
1124    {
[1889]1125
1126        // don't add this triangle to the primitive if it shares any vertex with
[5949]1127        // the supertriangle triangles have already been removed.
1128        // Don't need this test: ti->a() <= last_valid_index && ti->b() <= last_valid_index && ti->c() <= last_valid_index &&
1129          // Don't add degenerate (zero radius) triangles
1130        if ( ti->get_circumcircle().z()>0.0)
[4611]1131        {
[1889]1132
[4611]1133            if (normals_.valid())
1134            {
[1889]1135                (normals_.get())->push_back(ti->compute_normal(points));
1136            }
1137
1138            pt_indices.push_back(ti->a());
1139            pt_indices.push_back(ti->b());
1140            pt_indices.push_back(ti->c());
1141        }
1142    }
1143
1144    prim_tris_ = new osg::DrawElementsUInt(GL_TRIANGLES, pt_indices.size(), &(pt_indices.front()));
1145
[6135]1146    osg::notify(osg::INFO) << "DelaunayTriangulator: process done, " << prim_tris_->getNumPrimitives() << " triangles remain\n";
[1889]1147   
1148    return true;
1149}
1150
[4611]1151void DelaunayTriangulator::removeInternalTriangles(DelaunayConstraint *dc )
1152{
[5949]1153    if (dc) { // 16.12.06 just in case....
1154        // Triangle_list *triangles
1155        // remove triangle from terrain prim_tris_ internal to each constraintline
1156        // and move to the constraint line to make an alternative geometry,
1157        // possibly with alternative texture, and texture map
1158        int ndel=0;
1159        osg::Vec3Array::iterator normitr;
1160        if( normals_.valid() )
1161            normitr = normals_->begin();
1162
1163        //        osg::notify(osg::WARN) << "DelaunayTriangulator: removeinternals, " << std::endl;
1164        for (osg::DrawElementsUInt::iterator triit=prim_tris_->begin(); triit!=prim_tris_->end(); )
[4611]1165        {
[5949]1166            // triangle joins points_[itr, itr+1, itr+2]
1167            Triangle tritest((*triit), *(triit+1), *(triit+2), points_.get());
1168            if ( dc->contains(tritest.compute_centroid( points_.get()) ) )
[4611]1169            {
[5949]1170                // centroid is inside the triangle, so IF inside linear, remove
1171                // osg::notify(osg::WARN) << "DelaunayTriangulator: remove, " << (*triit) << "," << *(triit+1) <<","<<*(triit+2)<< std::endl;
1172                dc->addtriangle((*triit), *(triit+1), *(triit+2));
1173                triit=prim_tris_->erase(triit);
1174                triit=prim_tris_->erase(triit);
1175                triit=prim_tris_->erase(triit);
1176                if (normals_.valid())
1177                {
1178                    // erase the corresponding normal
1179                    normitr=normals_->erase(normitr);
1180                }
1181                ndel++;
[4611]1182            }
[5949]1183            else
[4611]1184            {
[5949]1185                if (normals_.valid())
1186                {
1187                    normitr++;
1188                }
1189
1190                triit+=3;
[4611]1191            }
1192        }
[1889]1193
[5949]1194        osg::notify(osg::INFO) << "end of test dc, deleted " << ndel << std::endl;
1195    }
[1889]1196}
[4611]1197//=== DelaunayConstraint functions
1198
[6446]1199float DelaunayConstraint::windingNumber(const osg::Vec3 &testpoint) const
[4611]1200{
1201    // return winding number of loop around testpoint. Only in 2D, x-y coordinates assumed!
1202    float theta=0; // sum of angles subtended by the line array - the winding number
1203    const osg::Vec3Array *vertices= dynamic_cast<const osg::Vec3Array*>(getVertexArray());
[4728]1204    if (vertices)
[4611]1205    {
[4728]1206        for (unsigned int ipr=0; ipr<getNumPrimitiveSets(); ipr++)
[4611]1207        {
[4728]1208            const osg::PrimitiveSet* prset=getPrimitiveSet(ipr);
1209            if (prset->getMode()==osg::PrimitiveSet::LINE_LOOP)
[4611]1210            {
[4728]1211                // nothing else loops
1212                // start with the last point on the loop
1213                const osg::Vec3 prev=(*vertices)[prset->index (prset->getNumIndices()-1)];
1214                osg::Vec3 pi(prev.x()-testpoint.x(),prev.y()-testpoint.y(),0);
1215                pi.normalize();
1216                for (unsigned int i=0; i<prset->getNumIndices(); i++)
[4611]1217                {
[4728]1218                    const osg::Vec3 curp=(*vertices)[prset->index (i)];
1219                    osg::Vec3 edge(curp.x()-testpoint.x(),curp.y()-testpoint.y(),0);
1220                    edge.normalize();
1221                    double cth=edge*pi;
1222                    if (cth<=-0.99999 )
[4611]1223                    {
[4728]1224                        // testpoint is on edge and between 2 points
1225                        return 0; //
[4611]1226                    }
[4728]1227                    else
1228                    {
1229                        if (cth<0.99999)
1230                        {
1231                            float dang=(cth<1 && cth>-1)?acos(edge*pi):0; // this is unsigned angle
1232                            float zsign=edge.x()*pi.y()-pi.x()*edge.y(); // z component of..(edge^pi).z();
1233                            if (zsign>0) theta+=dang; // add the angle subtended appropriately
1234                            else if (zsign<0) theta-=dang;
1235                        }
1236                    }
1237                    pi=edge;
[4611]1238                }
1239            }
1240        }
1241    }
[4728]1242   
[4611]1243    return theta/osg::PI/2.0; // should be 0 or 2 pi.
1244}
1245osg::DrawElementsUInt *DelaunayConstraint::makeDrawable()
1246{
1247    // initialize index storage vector for internal triangles.
1248    std::vector<GLuint> pt_indices;
1249    pt_indices.reserve(_interiorTris.size() * 3);
1250    trilist::const_iterator ti;
1251    for (ti=_interiorTris.begin(); ti!=_interiorTris.end(); ++ti)
1252    {
1253       
1254        //  if (normals_.valid()) {
1255        //        (normals_.get())->push_back(ti->compute_normal(points));
1256        //  }
1257       
1258        pt_indices.push_back((*ti)[0]);
1259        pt_indices.push_back((*ti)[1]);
1260        pt_indices.push_back((*ti)[2]);
1261    }
1262    prim_tris_ = new osg::DrawElementsUInt(GL_TRIANGLES, pt_indices.size(), &(pt_indices.front()));
1263   
1264    return prim_tris_.get();
1265}
[6446]1266bool DelaunayConstraint::contains(const osg::Vec3 &testpoint) const 
[4611]1267{
1268    // true if point is internal to the loop.
1269    float theta=windingNumber(testpoint); // sum of angles subtended by the line array - the winding number
1270    return fabs(theta)>0.9; // should be 0 or 1 (or 2,3,4 for very complex not permitted loops).
1271}
[6446]1272bool DelaunayConstraint::outside(const osg::Vec3 &testpoint) const 
[4611]1273{
1274    // true if point is outside the loop.
1275    float theta=windingNumber(testpoint); // sum of angles subtended by the line array - the winding number
1276    return fabs(theta)<.05; // should be 0 if outside.
1277}
1278
1279
[6446]1280void DelaunayConstraint::addtriangle(int i1, int i2, int i3)
[4611]1281{
1282    // a triangle joins vertices i1,i2,i3 in the points of the delaunay triangles.
1283    // points is the array of poitns in the triangulator;
1284    // add triangle to the constraint
1285    int *ip=new int[3];
1286    ip[0]=i1;
1287    ip[1]=i2;
1288    ip[2]=i3;
1289    _interiorTris.push_back(ip);
1290}
1291osg::Vec3Array* DelaunayConstraint::getPoints(const osg::Vec3Array *points)
1292{
1293    //points_ is the array of points that can be used to render the triangles in this DC.
1294    osg::ref_ptr<osg::Vec3Array> points_=new osg::Vec3Array;
1295    trilist::iterator ti;
1296    for (ti=_interiorTris.begin(); ti!=_interiorTris.end(); ++ti) {
1297        int idx=0;
1298        int ip[3]={-1,-1,-1};
1299        // find if points[i1/i2/i3] already in the vertices points_
1300        for (osg::Vec3Array::iterator ivert=points_->begin(); ivert!=points_->end(); ivert++)
1301        {
1302            if (ip[0]<0 && *ivert==(*points)[(*ti)[0]])
1303            {
1304                (*ti)[0]=ip[0]=idx;
1305            }
1306            if (ip[1]<0 && *ivert==(*points)[(*ti)[1]])
1307            {
1308                (*ti)[1]=ip[1]=idx;
1309            }
1310            if (ip[2]<0 && *ivert==(*points)[(*ti)[2]])
1311            {
1312                (*ti)[2]=ip[2]=idx;
1313            }
1314            idx++;
1315        }
1316        if (ip[0]<0)
1317        {
1318            points_->push_back((*points)[(*ti)[0]]);
1319            (*ti)[0]=ip[0]=points_->size()-1;
1320        }
1321        if (ip[1]<0)
1322        {
1323            points_->push_back((*points)[(*ti)[1]]);
1324            (*ti)[1]=ip[1]=points_->size()-1;
1325        }
1326        if (ip[2]<0)
1327        {
1328            points_->push_back((*points)[(*ti)[2]]);
1329            (*ti)[2]=ip[2]=points_->size()-1;
1330        }
1331    }
1332    makeDrawable();
1333    return points_.release();
1334}
1335
1336void DelaunayConstraint::handleOverlaps(void)
1337{
[5985]1338    // use tessellator to interpolate crossing vertices.
[5928]1339    osg::ref_ptr<osgUtil::Tessellator> tscx=new osgUtil::Tessellator; // this assembles all the constraints
1340    tscx->setTessellationType(osgUtil::Tessellator::TESS_TYPE_GEOMETRY);
[4611]1341    tscx->setBoundaryOnly(true);
[5928]1342    tscx->setWindingType( osgUtil::Tessellator::TESS_WINDING_ODD);
[4611]1343    //  ODD chooses the winding =1, NOT overlapping areas of constraints.
1344    // nb this includes all the edges where delaunay constraints meet
1345    // draw a case to convince yourself!.
1346   
[5928]1347    tscx->retessellatePolygons(*this); // find all edges
[4611]1348}
1349
1350} // namespace osgutil
Note: See TracBrowser for help on using the browser.