root/OpenSceneGraph/trunk/src/osgSim/SphereSegment.cpp @ 13041

Revision 13041, 102.1 kB (checked in by robert, 2 years ago)

Ran script to remove trailing spaces and tabs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2006 Robert Osfield
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
14#include <osgSim/SphereSegment>
15
16#include <osg/Notify>
17#include <osg/CullFace>
18#include <osg/LineWidth>
19#include <osg/Transform>
20#include <osg/Geometry>
21#include <osg/TriangleIndexFunctor>
22#include <osg/ShapeDrawable>
23#include <osg/io_utils>
24
25#include <algorithm>
26#include <list>
27
28using namespace osgSim;
29
30// Define the collection of nested classes, all Drawables, which make
31// up the parts of the sphere segment.
32
33/**
34SphereSegment::Surface is the Drawable which represents the specified area of the
35sphere's surface.
36 */
37class SphereSegment::Surface: public osg::Drawable
38{
39public:
40    Surface(SphereSegment* ss): osg::Drawable(), _ss(ss) {}
41
42    Surface():_ss(0)
43    {
44        OSG_WARN<<
45            "Warning: unexpected call to osgSim::SphereSegment::Surface() default constructor"<<std::endl;
46    }
47
48    Surface(const Surface& rhs, const osg::CopyOp& co=osg::CopyOp::SHALLOW_COPY):osg::Drawable(rhs,co), _ss(0)
49    {
50        OSG_WARN<<
51            "Warning: unexpected call to osgSim::SphereSegment::Surface() copy constructor"<<std::endl;
52    }
53
54    META_Object(osgSim,Surface)
55
56    void drawImplementation(osg::RenderInfo& renderInfo) const;
57
58    virtual osg::BoundingBox computeBound() const;
59
60protected:
61
62private:
63
64    SphereSegment* _ss;
65};
66
67void SphereSegment::Surface::drawImplementation(osg::RenderInfo& renderInfo) const
68{
69    _ss->Surface_drawImplementation(*renderInfo.getState());
70}
71
72osg:: BoundingBox SphereSegment::Surface::computeBound() const
73{
74    osg:: BoundingBox bbox;
75    _ss->Surface_computeBound(bbox);
76    return bbox;
77}
78
79
80/**
81SphereSegment::EdgeLine is the Drawable which represents the line around the edge
82of the specified area of the sphere's EdgeLine.
83 */
84class SphereSegment::EdgeLine: public osg::Drawable
85{
86public:
87    EdgeLine(SphereSegment* ss): osg::Drawable(), _ss(ss) { init(); }
88
89    EdgeLine():_ss(0)
90    {
91        init();
92        OSG_WARN<<
93            "Warning: unexpected call to osgSim::SphereSegment::EdgeLine() default constructor"<<std::endl;
94    }
95
96    EdgeLine(const EdgeLine& rhs, const osg::CopyOp& co=osg::CopyOp::SHALLOW_COPY):osg::Drawable(rhs,co), _ss(0)
97    {
98        OSG_WARN<<
99            "Warning: unexpected call to osgSim::SphereSegment::EdgeLine() copy constructor"<<std::endl;
100    }
101
102    META_Object(osgSim,EdgeLine)
103
104    void drawImplementation(osg::RenderInfo& renderInfo) const;
105
106protected:
107
108    void init()
109    {
110        // switch off lighting.
111        getOrCreateStateSet()->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
112
113        //getOrCreateStateSet()->setAttributeAndModes(new osg::LineWidth(2.0),osg::StateAttribute::OFF);
114    }
115
116
117    virtual osg::BoundingBox computeBound() const;
118
119private:
120
121    SphereSegment* _ss;
122};
123
124void SphereSegment::EdgeLine::drawImplementation(osg::RenderInfo& renderInfo) const
125{
126    _ss->EdgeLine_drawImplementation(*renderInfo.getState());
127}
128
129osg::BoundingBox SphereSegment::EdgeLine::computeBound() const
130{
131    osg::BoundingBox bbox;
132    _ss->EdgeLine_computeBound(bbox);
133    return bbox;
134}
135
136
137
138/**
139SphereSegment::Side is a Drawable which represents one of the
140planar areas, at either the minimum or maximum azimuth.
141 */
142class SphereSegment::Side: public osg::Drawable
143{
144public:
145    Side(SphereSegment* ss, SphereSegment::SideOrientation po, SphereSegment::BoundaryAngle pa):
146            osg::Drawable(), _ss(ss), _planeOrientation(po), _BoundaryAngle(pa) {}
147
148    META_Object(osgSim,Side)
149
150    void drawImplementation(osg::RenderInfo& renderInfo) const;
151
152protected:
153
154    Side():_ss(0), _planeOrientation(SphereSegment::AZIM), _BoundaryAngle(SphereSegment::MIN)
155    {
156        OSG_WARN<<
157            "Warning: unexpected call to osgSim::SphereSegment::Side() default constructor"<<std::endl;
158    }
159
160    Side(const Side& rhs, const osg::CopyOp& co=osg:: CopyOp::SHALLOW_COPY):
161        osg::Drawable(rhs,co),
162        _ss(0),
163        _planeOrientation(rhs._planeOrientation),
164        _BoundaryAngle(rhs._BoundaryAngle)
165    {
166        OSG_WARN<<
167            "Warning: unexpected call to osgSim::SphereSegment::Side() copy constructor"<<std::endl;
168    }
169
170    virtual osg::BoundingBox computeBound() const;
171
172private:
173    SphereSegment* _ss;
174    SphereSegment::SideOrientation _planeOrientation;
175    SphereSegment::BoundaryAngle _BoundaryAngle;
176};
177
178
179void SphereSegment::Side::drawImplementation(osg::RenderInfo& renderInfo) const
180{
181    _ss->Side_drawImplementation(*renderInfo.getState(), _planeOrientation, _BoundaryAngle);
182}
183
184osg::BoundingBox SphereSegment::Side::computeBound() const
185{
186    osg::BoundingBox bbox;
187    _ss->Side_computeBound(bbox, _planeOrientation, _BoundaryAngle);
188    return bbox;
189}
190
191
192
193/**
194SphereSegment::Spoke is a Drawable which represents a spoke.
195 */
196class SphereSegment::Spoke: public osg::Drawable
197{
198public:
199    Spoke(SphereSegment* ss, SphereSegment::BoundaryAngle azAngle, SphereSegment::BoundaryAngle elevAngle):
200            osg::Drawable(), _ss(ss), _azAngle(azAngle), _elevAngle(elevAngle) { init(); }
201
202    META_Object(osgSim,Spoke)
203
204    void drawImplementation(osg::RenderInfo& renderInfo) const;
205
206protected:
207
208    Spoke():_ss(0)
209    {
210        init();
211        OSG_WARN<<
212            "Warning: unexpected call to osgSim::SphereSegment::Spoke() default constructor"<<std::endl;
213    }
214
215    Spoke(const Spoke& rhs, const osg::CopyOp& co=osg:: CopyOp::SHALLOW_COPY):
216        osg::Drawable(rhs,co),
217        _ss(0),
218        _azAngle(rhs._azAngle), _elevAngle(rhs._elevAngle)
219    {
220        OSG_WARN<<
221            "Warning: unexpected call to osgSim::SphereSegment::Spoke() copy constructor"<<std::endl;
222    }
223
224
225    void init()
226    {
227        // switch off lighting.
228        getOrCreateStateSet()->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
229
230        //getOrCreateStateSet()->setAttributeAndModes(new osg::LineWidth(2.0),osg::StateAttribute::OFF);
231    }
232
233    virtual osg::BoundingBox computeBound() const;
234
235private:
236    SphereSegment* _ss;
237    SphereSegment::BoundaryAngle _azAngle, _elevAngle;
238};
239
240void SphereSegment::Spoke::drawImplementation(osg::RenderInfo& renderInfo) const
241{
242    _ss->Spoke_drawImplementation(*renderInfo.getState(), _azAngle, _elevAngle);
243}
244
245osg::BoundingBox SphereSegment::Spoke::computeBound() const
246{
247    osg::BoundingBox bbox;
248    _ss->Spoke_computeBound(bbox, _azAngle, _elevAngle);
249    return bbox;
250}
251
252SphereSegment::SphereSegment(const osg::Vec3& centre, float radius, const osg::Vec3& vec, float azRange,
253                float elevRange, int density):
254    osg::Geode(),
255    _centre(centre), _radius(radius),
256    _density(density),
257    _drawMask(DrawMask(ALL))
258{
259    // Rather than store the vector, we'll work out the azimuth boundaries and elev
260    // boundaries now, rather than at draw time.
261    setArea(vec, azRange, elevRange);
262
263    init();
264}
265
266void SphereSegment::setCentre(const osg::Vec3& c)
267{
268    _centre = c;
269    dirtyAllDrawableDisplayLists();
270    dirtyAllDrawableBounds();
271    dirtyBound();
272}
273
274const osg::Vec3& SphereSegment::getCentre() const
275{
276    return _centre;
277}
278
279void SphereSegment::setRadius(float r)
280{
281    _radius = r;
282    dirtyAllDrawableDisplayLists();
283    dirtyAllDrawableBounds();
284    dirtyBound();
285}
286
287float SphereSegment::getRadius() const
288{
289    return _radius;
290}
291
292
293void SphereSegment::setArea(const osg::Vec3& v, float azRange, float elevRange)
294{
295    osg::Vec3 vec(v);
296
297    vec.normalize();    // Make sure we're unit length
298
299    // Calculate the elevation range
300    float xyLen = sqrtf(vec.x()*vec.x() + vec.y()*vec.y());
301    float elev = atan2(vec.z(), xyLen);   // Elevation angle
302
303    elevRange /= 2.0f;
304    _elevMin = elev - elevRange;
305    _elevMax = elev + elevRange;
306
307    // Calculate the azimuth range, cater for trig ambiguities
308    float az = atan2(vec.x(), vec.y());
309
310    azRange /= 2.0f;
311    _azMin = az - azRange;
312    _azMax = az + azRange;
313
314    dirtyAllDrawableDisplayLists();
315    dirtyAllDrawableBounds();
316    dirtyBound();
317}
318
319void SphereSegment::getArea(osg::Vec3& vec, float& azRange, float& elevRange) const
320{
321    azRange = _azMax - _azMin;
322    elevRange = _elevMax - _elevMin;
323
324    float az = azRange/2.0f;
325    float elev = elevRange/2.0f;
326    vec.set(cos(elev)*sin(az), cos(elev)*cos(az), sin(elev));
327}
328
329void SphereSegment::setArea(float azMin, float azMax,
330    float elevMin, float elevMax)
331{
332    _azMin=azMin;
333    _azMax=azMax;
334    _elevMin=elevMin;
335    _elevMax=elevMax;
336
337    dirtyAllDrawableDisplayLists();
338    dirtyAllDrawableBounds();
339    dirtyBound();
340}
341
342void SphereSegment::getArea(float &azMin, float &azMax,
343    float &elevMin, float &elevMax) const
344{
345    azMin=_azMin;
346    azMax=_azMax;
347    elevMin=_elevMin;
348    elevMax=_elevMax;
349}
350
351void SphereSegment::setDensity(int density)
352{
353    _density = density;
354    dirtyAllDrawableDisplayLists();
355}
356
357int SphereSegment::getDensity() const
358{
359    return _density;
360}
361
362void SphereSegment::init()
363{
364    addDrawable(new Surface(this));
365
366    addDrawable(new EdgeLine(this));
367
368    addDrawable(new Side(this,AZIM,MIN));
369    addDrawable(new Side(this,AZIM,MAX));
370    addDrawable(new Side(this,ELEV,MIN));
371    addDrawable(new Side(this,ELEV,MAX));
372
373    addDrawable(new Spoke(this,MIN,MIN));
374    addDrawable(new Spoke(this,MIN,MAX));
375    addDrawable(new Spoke(this,MAX,MIN));
376    addDrawable(new Spoke(this,MAX,MAX));
377}
378
379void SphereSegment::dirtyAllDrawableDisplayLists()
380{
381    for(DrawableList::iterator itr = _drawables.begin();
382        itr != _drawables.end();
383        ++itr)
384    {
385        (*itr)->dirtyDisplayList();
386    }
387}
388
389void SphereSegment::dirtyAllDrawableBounds()
390{
391    for(DrawableList::iterator itr = _drawables.begin();
392        itr != _drawables.end();
393        ++itr)
394    {
395        (*itr)->dirtyBound();
396    }
397}
398
399void SphereSegment::Surface_drawImplementation(osg::State& state) const
400{
401    const float azIncr = (_azMax - _azMin)/_density;
402    const float elevIncr = (_elevMax - _elevMin)/_density;
403
404    // Draw the area on the sphere surface if needed
405    // ---------------------------------------------
406    if(_drawMask & SURFACE)
407    {
408        osg::GLBeginEndAdapter& gl = (state.getGLBeginEndAdapter());
409
410        gl.Color4fv(_surfaceColor.ptr());
411
412        bool drawBackSide = true;
413        bool drawFrontSide = true;
414
415        // draw back side.
416        if (drawBackSide)
417        {
418            for(int i=0; i+1<=_density; i++)
419            {
420                // Because we're drawing quad strips, we need to work out
421                // two azimuth values, to form each edge of the (z-vertical)
422                // strips
423                float az1 = _azMin + (i*azIncr);
424                float az2 = _azMin + ((i+1)*azIncr);
425
426                gl.Begin(GL_QUAD_STRIP);
427                for (int j=0; j<=_density; j++)
428                {
429                    float elev = _elevMin + (j*elevIncr);
430
431                    // QuadStrip Edge formed at az1
432                    // ----------------------------
433
434                    // Work out the sphere normal
435                    float x = cos(elev)*sin(az1);
436                    float y = cos(elev)*cos(az1);
437                    float z = sin(elev);
438
439                    gl.Normal3f(-x, -y, -z);
440                    gl.Vertex3f(_centre.x() + _radius*x,
441                            _centre.y() + _radius*y,
442                            _centre.z() + _radius*z);
443
444                    // QuadStrip Edge formed at az2
445                    // ----------------------------
446
447                    // Work out the sphere normal
448                    x = cos(elev)*sin(az2);
449                    y = cos(elev)*cos(az2);
450                    // z = sin(elev);   z doesn't change
451
452                    gl.Normal3f(-x, -y, -z);
453                    gl.Vertex3f(_centre.x() + _radius*x,
454                            _centre.y() + _radius*y,
455                            _centre.z() + _radius*z);
456                }
457                gl.End();
458            }
459        }
460
461        // draw front side
462        if (drawFrontSide)
463        {
464            for(int i=0; i+1<=_density; i++)
465            {
466                // Because we're drawing quad strips, we need to work out
467                // two azimuth values, to form each edge of the (z-vertical)
468                // strips
469                float az1 = _azMin + (i*azIncr);
470                float az2 = _azMin + ((i+1)*azIncr);
471
472                gl.Begin(GL_QUAD_STRIP);
473                for (int j=0; j<=_density; j++)
474                {
475                    float elev = _elevMin + (j*elevIncr);
476
477                    // QuadStrip Edge formed at az1
478                    // ----------------------------
479
480                    // Work out the sphere normal
481                    float x = cos(elev)*sin(az2);
482                    float y = cos(elev)*cos(az2);
483                    float z = sin(elev);
484
485                    gl.Normal3f(x, y, z);
486                    gl.Vertex3f(_centre.x() + _radius*x,
487                            _centre.y() + _radius*y,
488                            _centre.z() + _radius*z);
489
490                    // QuadStrip Edge formed at az2
491                    // ----------------------------
492
493                    // Work out the sphere normal
494                    // z = sin(elev);   z doesn't change
495
496                    x = cos(elev)*sin(az1);
497                    y = cos(elev)*cos(az1);
498
499                    gl.Normal3f(x, y, z);
500                    gl.Vertex3f(_centre.x() + _radius*x,
501                            _centre.y() + _radius*y,
502                            _centre.z() + _radius*z);
503                }
504                gl.End();
505            }
506        }
507    }
508}
509
510bool SphereSegment::Surface_computeBound(osg::BoundingBox& bbox) const
511{
512    bbox.init();
513
514    float azIncr = (_azMax - _azMin)/_density;
515    float elevIncr = (_elevMax - _elevMin)/_density;
516
517    for(int i=0; i<=_density; i++){
518
519        float az = _azMin + (i*azIncr);
520
521        for(int j=0; j<=_density; j++){
522
523            float elev = _elevMin + (j*elevIncr);
524
525            bbox.expandBy(
526                osg::Vec3(_centre.x() + _radius*cos(elev)*sin(az),
527                        _centre.y() + _radius*cos(elev)*cos(az),
528                        _centre.z() + _radius*sin(elev))
529            );
530        }
531    }
532    return true;
533}
534
535void SphereSegment::EdgeLine_drawImplementation(osg::State& state) const
536{
537    const float azIncr = (_azMax - _azMin)/_density;
538    const float elevIncr = (_elevMax - _elevMin)/_density;
539
540    // Draw the edgeline if necessary
541    // ------------------------------
542    if(_drawMask & EDGELINE)
543    {
544        osg::GLBeginEndAdapter& gl = (state.getGLBeginEndAdapter());
545
546        gl.Color4fv(_edgeLineColor.ptr());
547
548        // Top edge
549        gl.Begin(GL_LINE_STRIP);
550        int i;
551        for(i=0; i<=_density; i++)
552        {
553            float az = _azMin + (i*azIncr);
554            gl.Vertex3f(
555                _centre.x() + _radius*cos(_elevMax)*sin(az),
556                _centre.y() + _radius*cos(_elevMax)*cos(az),
557                _centre.z() + _radius*sin(_elevMax));
558        }
559        gl.End();
560
561        // Bottom edge
562        gl.Begin(GL_LINE_STRIP);
563        for(i=0; i<=_density; i++)
564        {
565            float az = _azMin + (i*azIncr);
566            gl.Vertex3f(
567                _centre.x() + _radius*cos(_elevMin)*sin(az),
568                _centre.y() + _radius*cos(_elevMin)*cos(az),
569                _centre.z() + _radius*sin(_elevMin));
570        }
571        gl.End();
572
573        // Left edge
574        gl.Begin(GL_LINE_STRIP);
575        int j;
576        for(j=0; j<=_density; j++)
577        {
578            float elev = _elevMin + (j*elevIncr);
579            gl.Vertex3f(
580                _centre.x() + _radius*cos(elev)*sin(_azMin),
581                _centre.y() + _radius*cos(elev)*cos(_azMin),
582                _centre.z() + _radius*sin(elev));
583        }
584        gl.End();
585
586        // Right edge
587        gl.Begin(GL_LINE_STRIP);
588        for(j=0; j<=_density; j++)
589        {
590            float elev = _elevMin + (j*elevIncr);
591            gl.Vertex3f(
592                _centre.x() + _radius*cos(elev)*sin(_azMax),
593                _centre.y() + _radius*cos(elev)*cos(_azMax),
594                _centre.z() + _radius*sin(elev));
595        }
596        gl.End();
597    }
598}
599
600bool SphereSegment::EdgeLine_computeBound(osg::BoundingBox& bbox) const
601{
602    bbox.init();
603
604    float azIncr = (_azMax - _azMin)/_density;
605    float elevIncr = (_elevMax - _elevMin)/_density;
606
607    // Top edge
608    int i;
609    for(i=0; i<=_density; i++)
610    {
611        float az = _azMin + (i*azIncr);
612        bbox.expandBy(
613            _centre.x() + _radius*cos(_elevMax)*sin(az),
614            _centre.y() + _radius*cos(_elevMax)*cos(az),
615            _centre.z() + _radius*sin(_elevMax));
616    }
617
618    // Bottom edge
619    for(i=0; i<=_density; i++)
620    {
621        float az = _azMin + (i*azIncr);
622        bbox.expandBy(
623            _centre.x() + _radius*cos(_elevMin)*sin(az),
624            _centre.y() + _radius*cos(_elevMin)*cos(az),
625            _centre.z() + _radius*sin(_elevMin));
626    }
627
628    // Left edge
629    int j;
630    for(j=0; j<=_density; j++)
631    {
632        float elev = _elevMin + (j*elevIncr);
633        bbox.expandBy(
634            _centre.x() + _radius*cos(elev)*sin(_azMin),
635            _centre.y() + _radius*cos(elev)*cos(_azMin),
636            _centre.z() + _radius*sin(elev));
637    }
638
639    // Right edge
640    for(j=0; j<=_density; j++)
641    {
642        float elev = _elevMin + (j*elevIncr);
643        bbox.expandBy(
644            _centre.x() + _radius*cos(elev)*sin(_azMax),
645            _centre.y() + _radius*cos(elev)*cos(_azMax),
646            _centre.z() + _radius*sin(elev));
647    }
648
649    return true;
650}
651
652void SphereSegment::Side_drawImplementation(osg::State& state,
653                                SphereSegment::SideOrientation orientation,
654                                SphereSegment::BoundaryAngle boundaryAngle) const
655{
656    // Draw the planes if necessary
657    // ----------------------------
658    if(_drawMask & SIDES)
659    {
660        osg::GLBeginEndAdapter& gl = (state.getGLBeginEndAdapter());
661
662        bool drawBackSide = true;
663        bool drawFrontSide = true;
664        int start, end, delta;
665
666        gl.Color4fv(_planeColor.ptr());
667
668        // draw back side.
669        if (drawBackSide)
670        {
671
672            if(orientation == AZIM)      // This is a plane at a given azimuth
673            {
674                const float az = (boundaryAngle==MIN?_azMin:_azMax);
675                const float elevIncr = (_elevMax - _elevMin)/_density;
676
677                // Normal
678                osg::Vec3 normal = osg::Vec3(cos(_elevMin)*sin(az), cos(_elevMin)*cos(az), sin(_elevMin))
679                                    ^ osg::Vec3(cos(_elevMax)*sin(az), cos(_elevMax)*cos(az), sin(_elevMax));
680
681                if (boundaryAngle==MIN)
682                {
683                    start = _density;
684                    end = 0;
685                }
686                else
687                {
688                    start = 0;
689                    end = _density;
690                    normal = -normal;   // Make sure normals oriented 'outwards'
691                }
692                delta = end>start?1:-1;
693
694                if (drawBackSide)
695                {
696                    // Tri fan
697                    gl.Normal3f(-normal.x(),-normal.y(),-normal.z());
698                    gl.Begin(GL_TRIANGLE_FAN);
699                    gl.Vertex3fv(_centre.ptr());
700                    for (int j=start; j!=end+delta; j+=delta)
701                    {
702                        float elev = _elevMin + (j*elevIncr);
703                        gl.Vertex3f( _centre.x() + _radius*cos(elev)*sin(az),
704                                    _centre.y() + _radius*cos(elev)*cos(az),
705                                    _centre.z() + _radius*sin(elev));
706                    }
707                    gl.End();
708                }
709
710                if (boundaryAngle==MIN)
711                {
712                    start = 0;
713                    end = _density;
714                }
715                else
716                {
717                    start = _density;
718                    end = 0;
719                }
720                delta = end>start?1:-1;
721
722                if (drawFrontSide)
723                {
724                    gl.Normal3fv(normal.ptr());
725                    gl.Begin(GL_TRIANGLE_FAN);
726                    gl.Vertex3fv(_centre.ptr());
727                    for (int j=start; j!=end+delta; j+=delta)
728                    {
729                        float elev = _elevMin + (j*elevIncr);
730                        gl.Vertex3f( _centre.x() + _radius*cos(elev)*sin(az),
731                                    _centre.y() + _radius*cos(elev)*cos(az),
732                                    _centre.z() + _radius*sin(elev));
733                    }
734                    gl.End();
735                }
736
737            }
738            else if(orientation == ELEV) // This is a plane at a given elevation
739            {
740                const float elev = (boundaryAngle==MIN?_elevMin:_elevMax);
741                const float azIncr = (_azMax - _azMin)/_density;
742
743                // Normal
744                osg::Vec3 normal = osg::Vec3(cos(elev)*sin(_azMax), cos(elev)*cos(_azMax), sin(elev))
745                                    ^ osg::Vec3(cos(elev)*sin(_azMin), cos(elev)*cos(_azMin), sin(elev));
746
747
748                if (boundaryAngle==MIN)
749                {
750                    start = _density;
751                    end = 0;
752                    normal = -normal;   // Make sure normals orientated 'outwards'
753                }
754                else
755                {
756                    start = 0;
757                    end = _density;
758                }
759                delta = end>start?1:-1;
760
761                if (drawBackSide)
762                {
763                    gl.Normal3f(-normal.x(),-normal.y(),-normal.z());
764
765                    // Tri fan
766                    gl.Begin(GL_TRIANGLE_FAN);
767                    gl.Vertex3fv(_centre.ptr());
768                    for (int j=start; j!=end+delta; j+=delta)
769                    {
770                        float az = _azMin + (j*azIncr);
771                        gl.Vertex3f( _centre.x() + _radius*cos(elev)*sin(az),
772                                    _centre.y() + _radius*cos(elev)*cos(az),
773                                    _centre.z() + _radius*sin(elev));
774                    }
775                    gl.End();
776                }
777
778                if (boundaryAngle==MIN)
779                {
780                    start = 0;
781                    end = _density;
782                }
783                else
784                {
785                    start = _density;
786                    end = 0;
787                }
788                delta = end>start?1:-1;
789
790                if (drawFrontSide)
791                {
792                    gl.Normal3fv(normal.ptr());
793
794                    // Tri fan
795                    gl.Begin(GL_TRIANGLE_FAN);
796                    gl.Vertex3fv(_centre.ptr());
797                    for (int j=start; j!=end+delta; j+=delta)
798                    {
799                        float az = _azMin + (j*azIncr);
800                        gl.Vertex3f( _centre.x() + _radius*cos(elev)*sin(az),
801                                    _centre.y() + _radius*cos(elev)*cos(az),
802                                    _centre.z() + _radius*sin(elev));
803                    }
804                    gl.End();
805                }
806
807            }
808        }
809    }
810}
811
812bool SphereSegment::Side_computeBound(osg::BoundingBox& bbox,
813                                SphereSegment::SideOrientation orientation,
814                                SphereSegment::BoundaryAngle boundaryAngle) const
815{
816    bbox.init();
817    bbox.expandBy(_centre);
818
819    if(orientation == AZIM)      // This is a plane at a given azimuth
820    {
821        const float az = (boundaryAngle==MIN?_azMin:_azMax);
822        const float elevIncr = (_elevMax - _elevMin)/_density;
823
824        for (int j=0; j<=_density; j++)
825        {
826            float elev = _elevMin + (j*elevIncr);
827            bbox.expandBy(
828                _centre.x() + _radius*cos(elev)*sin(az),
829                _centre.y() + _radius*cos(elev)*cos(az),
830                _centre.z() + _radius*sin(elev));
831        }
832    }
833    else if(orientation == ELEV) // This is a plane at a given elevation
834    {
835        const float elev = (boundaryAngle==MIN?_elevMin:_elevMax);
836        const float azIncr = (_azMax - _azMin)/_density;
837
838        for(int i=0; i<=_density; i++)
839        {
840            float az = _azMin + (i*azIncr);
841            bbox.expandBy(
842                _centre.x() + _radius*cos(elev)*sin(az),
843                _centre.y() + _radius*cos(elev)*cos(az),
844                _centre.z() + _radius*sin(elev));
845        }
846    }
847
848    return true;
849}
850
851void SphereSegment::Spoke_drawImplementation(osg::State& state, BoundaryAngle azAngle, BoundaryAngle elevAngle) const
852{
853    if(_drawMask & SPOKES)
854    {
855        osg::GLBeginEndAdapter& gl = (state.getGLBeginEndAdapter());
856
857        gl.Color4fv(_spokeColor.ptr());
858
859        const float az = (azAngle==MIN?_azMin:_azMax);
860        const float elev = (elevAngle==MIN?_elevMin:_elevMax);
861
862        gl.Begin(GL_LINES);
863            gl.Vertex3fv(_centre.ptr());
864            gl.Vertex3f( _centre.x() + _radius*cos(elev)*sin(az),
865                        _centre.y() + _radius*cos(elev)*cos(az),
866                        _centre.z() + _radius*sin(elev));
867        gl.End();
868    }
869}
870
871bool SphereSegment::Spoke_computeBound(osg::BoundingBox& bbox, BoundaryAngle azAngle, BoundaryAngle elevAngle) const
872{
873    const float az = (azAngle==MIN?_azMin:_azMax);
874    const float elev = (elevAngle==MIN?_elevMin:_elevMax);
875
876    bbox.expandBy(_centre);
877    bbox.expandBy(  _centre.x() + _radius*cos(elev)*sin(az),
878                    _centre.y() + _radius*cos(elev)*cos(az),
879                    _centre.z() + _radius*sin(elev));
880
881    return true;
882}
883
884void SphereSegment::setDrawMask(int dm)
885{
886    _drawMask=dm;
887    dirtyAllDrawableDisplayLists();
888    dirtyAllDrawableBounds();
889    dirtyBound();
890}
891
892struct ActivateTransparencyOnType
893{
894    ActivateTransparencyOnType(const std::type_info& t): _t(t) {}
895
896    void operator()(osg::ref_ptr<osg::Drawable>& dptr) const
897    {
898        if(typeid(*dptr)==_t)
899        {
900            osg::StateSet* ss = dptr->getOrCreateStateSet();
901            ss->setRenderingHint(osg::StateSet::TRANSPARENT_BIN);
902
903            ss->setAttributeAndModes(new osg::CullFace(osg::CullFace::BACK),osg::StateAttribute::ON);
904            ss->setMode(GL_BLEND,osg::StateAttribute::ON);
905
906            dptr->dirtyDisplayList();
907        }
908    }
909
910    const std::type_info&  _t;
911
912protected:
913
914    ActivateTransparencyOnType& operator = (const ActivateTransparencyOnType&) { return *this; }
915};
916
917struct DeactivateTransparencyOnType
918{
919    DeactivateTransparencyOnType(const std::type_info& t): _t(t) {}
920
921    void operator()(osg::ref_ptr<osg::Drawable>& dptr) const
922    {
923        if(typeid(*dptr)==_t)
924        {
925            osg::StateSet* ss = dptr->getStateSet();
926            if(ss) ss->setRenderingHint(osg::StateSet::OPAQUE_BIN);
927
928            dptr->dirtyDisplayList();
929        }
930    }
931
932    const std::type_info&  _t;
933
934protected:
935
936    DeactivateTransparencyOnType& operator = (const DeactivateTransparencyOnType&) { return *this; }
937};
938
939void SphereSegment::setSurfaceColor(const osg::Vec4& c)
940{
941    _surfaceColor=c;
942
943    if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(Surface)));
944    else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(Surface)));
945}
946
947void SphereSegment::setSpokeColor(const osg::Vec4& c)
948{
949    _spokeColor=c;
950
951    if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(Spoke)));
952    else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(Spoke)));
953}
954
955void SphereSegment::setEdgeLineColor(const osg::Vec4& c)
956{
957    _edgeLineColor=c;
958
959    if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(EdgeLine)));
960    else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(EdgeLine)));
961}
962
963void SphereSegment::setSideColor(const osg::Vec4& c)
964{
965    _planeColor=c;
966
967    if(c.w() != 1.0) std::for_each(_drawables.begin(), _drawables.end(), ActivateTransparencyOnType(typeid(Side)));
968    else std::for_each(_drawables.begin(), _drawables.end(), DeactivateTransparencyOnType(typeid(Side)));
969}
970
971void SphereSegment::setAllColors(const osg::Vec4& c)
972{
973    setSurfaceColor(c);
974    setSpokeColor(c);
975    setEdgeLineColor(c);
976    setSideColor(c);
977}
978
979////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
980//
981// SphereSegment intersection code.
982
983class PolytopeVisitor : public osg::NodeVisitor
984{
985    public:
986
987        typedef std::pair<osg::Matrix, osg::Polytope> MatrixPolytopePair;
988        typedef std::vector<MatrixPolytopePair> PolytopeStack;
989
990        struct Hit
991        {
992            Hit(const osg::Matrix& matrix, osg::NodePath& nodePath, osg::Drawable* drawable):
993                _matrix(matrix),
994                _nodePath(nodePath),
995                _drawable(drawable) {}
996
997            osg::Matrix                 _matrix;
998            osg::NodePath               _nodePath;
999            osg::ref_ptr<osg::Drawable> _drawable;
1000        };
1001
1002        typedef std::vector<Hit> HitList;
1003
1004
1005        PolytopeVisitor(const osg::Matrix& matrix, const osg::Polytope& polytope):
1006            osg::NodeVisitor(osg::NodeVisitor::TRAVERSE_ACTIVE_CHILDREN)
1007        {
1008            _polytopeStack.push_back(MatrixPolytopePair());
1009            _polytopeStack.back().first = matrix;
1010            _polytopeStack.back().second.setAndTransformProvidingInverse(polytope, _polytopeStack.back().first);
1011        }
1012
1013        META_NodeVisitor("osgSim","PolytopeVisitor")
1014
1015        void reset()
1016        {
1017            _polytopeStack.clear();
1018            _hits.clear();
1019        }
1020
1021        void apply(osg::Node& node)
1022        {
1023            if (_polytopeStack.back().second.contains(node.getBound()))
1024            {
1025                traverse(node);
1026            }
1027        }
1028
1029        void apply(osg::Transform& transform)
1030        {
1031            if (_polytopeStack.back().second.contains(transform.getBound()))
1032            {
1033                const osg::Polytope& polytope = _polytopeStack.front().second;
1034
1035                osg::Matrix matrix = _polytopeStack.back().first;
1036                transform.computeLocalToWorldMatrix(matrix, this);
1037
1038                _polytopeStack.push_back(MatrixPolytopePair());
1039                _polytopeStack.back().first = matrix;
1040                _polytopeStack.back().second.setAndTransformProvidingInverse(polytope, matrix);
1041
1042                traverse(transform);
1043
1044                _polytopeStack.back();
1045            }
1046        }
1047
1048        void apply(osg::Geode& node)
1049        {
1050            if (_polytopeStack.back().second.contains(node.getBound()))
1051            {
1052                for(unsigned int i=0; i<node.getNumDrawables(); ++i)
1053                {
1054                    if (_polytopeStack.back().second.contains(node.getDrawable(i)->getBound()))
1055                    {
1056                        _hits.push_back(Hit(_polytopeStack.back().first,getNodePath(),node.getDrawable(i)));
1057                    }
1058
1059                }
1060
1061                traverse(node);
1062            }
1063        }
1064
1065        HitList& getHits() { return _hits; }
1066
1067    protected:
1068
1069        PolytopeStack   _polytopeStack;
1070        HitList         _hits;
1071
1072};
1073
1074
1075SphereSegment::LineList SphereSegment::computeIntersection(const osg::Matrixd& transform, osg::Node* subgraph)
1076{
1077    OSG_INFO<<"Creating line intersection between sphere segment and subgraph."<<std::endl;
1078
1079    osg::BoundingBox bb = getBoundingBox();
1080
1081    osg::Polytope polytope;
1082    polytope.add(osg::Plane(1.0,0.0,0.0,-bb.xMin()));
1083    polytope.add(osg::Plane(-1.0,0.0,0.0,bb.xMax()));
1084    polytope.add(osg::Plane(0.0,1.0,0.0,-bb.yMin()));
1085    polytope.add(osg::Plane(0.0,-1.0,0.0,bb.yMax()));
1086    polytope.add(osg::Plane(0.0,0.0,1.0,-bb.zMin()));
1087    polytope.add(osg::Plane(0.0,0.0,-1.0,bb.zMax()));
1088
1089    osg::Plane pl;
1090    pl.set(osg::Vec3(-1.0,0.0,0.0), bb.corner(0));
1091    PolytopeVisitor polytopeVisitor(transform, polytope);
1092
1093    subgraph->accept(polytopeVisitor);
1094
1095    if (polytopeVisitor.getHits().empty())
1096    {
1097        OSG_INFO<<"No hits found."<<std::endl;
1098        return LineList();
1099    }
1100
1101    // create a LineList to store all the compute line segments
1102    LineList all_lines;
1103
1104    // compute the line intersections with each of the hit drawables
1105    OSG_INFO<<"Hits found. "<<polytopeVisitor.getHits().size()<<std::endl;
1106    PolytopeVisitor::HitList& hits = polytopeVisitor.getHits();
1107    for(PolytopeVisitor::HitList::iterator itr = hits.begin();
1108        itr != hits.end();
1109        ++itr)
1110    {
1111        SphereSegment::LineList lines = computeIntersection(itr->_matrix, itr->_drawable.get());
1112        all_lines.insert(all_lines.end(), lines.begin(), lines.end());
1113    }
1114
1115    // join all the lines that have ends that are close together..
1116
1117    return all_lines;
1118}
1119
1120osg::Node* SphereSegment::computeIntersectionSubgraph(const osg::Matrixd& transform, osg::Node* subgraph)
1121{
1122    OSG_INFO<<"Creating line intersection between sphere segment and subgraph."<<std::endl;
1123
1124    osg::BoundingBox bb = getBoundingBox();
1125
1126    osg::Polytope polytope;
1127    polytope.add(osg::Plane(1.0,0.0,0.0,-bb.xMin()));
1128    polytope.add(osg::Plane(-1.0,0.0,0.0,bb.xMax()));
1129    polytope.add(osg::Plane(0.0,1.0,0.0,-bb.yMin()));
1130    polytope.add(osg::Plane(0.0,-1.0,0.0,bb.yMax()));
1131    polytope.add(osg::Plane(0.0,0.0,1.0,-bb.zMin()));
1132    polytope.add(osg::Plane(0.0,0.0,-1.0,bb.zMax()));
1133
1134    osg::Plane pl;
1135    pl.set(osg::Vec3(-1.0,0.0,0.0), bb.corner(0));
1136    PolytopeVisitor polytopeVisitor(transform, polytope);
1137
1138    subgraph->accept(polytopeVisitor);
1139
1140    if (polytopeVisitor.getHits().empty())
1141    {
1142        OSG_INFO<<"No hits found."<<std::endl;
1143        return 0;
1144    }
1145
1146    // create a LineList to store all the compute line segments
1147    osg::Group* group = new osg::Group;
1148
1149    // compute the line intersections with each of the hit drawables
1150    OSG_INFO<<"Hits found. "<<polytopeVisitor.getHits().size()<<std::endl;
1151    PolytopeVisitor::HitList& hits = polytopeVisitor.getHits();
1152    for(PolytopeVisitor::HitList::iterator itr = hits.begin();
1153        itr != hits.end();
1154        ++itr)
1155    {
1156        group->addChild(computeIntersectionSubgraph(itr->_matrix, itr->_drawable.get()));
1157    }
1158
1159    // join all the lines that have ends that are close together..
1160
1161    return group;
1162}
1163
1164namespace SphereSegmentIntersector
1165{
1166
1167    struct dereference_less
1168    {
1169        template<class T, class U>
1170        inline bool operator() (const T& lhs,const U& rhs) const
1171        {
1172            return *lhs < *rhs;
1173        }
1174    };
1175
1176    struct SortFunctor
1177    {
1178        typedef std::vector< osg::Vec3 > VertexArray;
1179
1180        SortFunctor(VertexArray& vertices):
1181            _vertices(vertices) {}
1182
1183        bool operator() (unsigned int p1, unsigned int p2) const
1184        {
1185            return _vertices[p1]<_vertices[p2];
1186        }
1187
1188        VertexArray& _vertices;
1189
1190    protected:
1191
1192        SortFunctor& operator = (const SortFunctor&) { return *this; }
1193    };
1194
1195
1196    struct TriangleIntersectOperator
1197    {
1198
1199        TriangleIntersectOperator():
1200            _radius(-1.0f),
1201            _azMin(0.0f),
1202            _azMax(0.0f),
1203            _elevMin(0.0f),
1204            _elevMax(0.0f),
1205            _numOutside(0),
1206            _numInside(0),
1207            _numIntersecting(0) {}
1208
1209        enum SurfaceType
1210        {
1211            NO_SURFACE = 0,
1212            RADIUS_SURFACE,
1213            LEFT_SURFACE,
1214            RIGHT_SURFACE,
1215            BOTTOM_SURFACE,
1216            TOP_SURFACE
1217        };
1218
1219        struct Triangle;
1220
1221        struct Edge : public osg::Referenced
1222        {
1223            typedef std::vector<Triangle*> TriangleList;
1224
1225            enum IntersectionType
1226            {
1227                NO_INTERSECTION,
1228                POINT_1,
1229                POINT_2,
1230                MID_POINT,
1231                BOTH_ENDS
1232            };
1233
1234            Edge(unsigned int p1, unsigned int p2, SurfaceType intersectionEdge=NO_SURFACE):
1235                _intersectionType(NO_INTERSECTION),
1236                _intersectionVertexIndex(0),
1237                _p1Outside(false),
1238                _p2Outside(false),
1239                _intersectionEdge(intersectionEdge)
1240            {
1241                if (p1>p2)
1242                {
1243                    _p1 = p2;
1244                    _p2 = p1;
1245                }
1246                else
1247                {
1248                    _p1 = p1;
1249                    _p2 = p2;
1250                }
1251            }
1252
1253            bool operator < (const Edge& edge) const
1254            {
1255                if (_p1<edge._p1) return true;
1256                else if (_p1>edge._p1) return false;
1257                else return _p2<edge._p2;
1258            }
1259
1260            inline void addTriangle(Triangle* tri)
1261            {
1262                TriangleList::iterator itr = std::find(_triangles.begin(), _triangles.end(), tri);
1263                if (itr==_triangles.end()) _triangles.push_back(tri);
1264            }
1265
1266            void removeFromToTraverseList(Triangle* tri)
1267            {
1268                TriangleList::iterator itr = std::find(_toTraverse.begin(), _toTraverse.end(), tri);
1269                if (itr!=_toTraverse.end()) _toTraverse.erase(itr);
1270            }
1271
1272
1273            bool completlyOutside() const { return _p1Outside && _p2Outside; }
1274
1275            unsigned int _p1;
1276            unsigned int _p2;
1277
1278            TriangleList _triangles;
1279
1280            // intersection information
1281            IntersectionType    _intersectionType;
1282            osg::Vec3           _intersectionVertex;
1283            unsigned int        _intersectionVertexIndex;
1284            bool                _p1Outside;
1285            bool                _p2Outside;
1286            TriangleList        _toTraverse;
1287            SurfaceType         _intersectionEdge;
1288        };
1289
1290        typedef std::list< osg::ref_ptr<Edge> >                     EdgeList;
1291
1292        struct Triangle : public osg::Referenced
1293        {
1294
1295            Triangle(unsigned int p1, unsigned int p2, unsigned int p3):
1296                _p1(p1), _p2(p2), _p3(p3),
1297                _e1(0), _e2(0), _e3(0)
1298            {
1299                sort();
1300            }
1301
1302            bool operator < (const Triangle& rhs) const
1303            {
1304                if (_p1 < rhs._p1) return true;
1305                else if (_p1 > rhs._p1) return false;
1306                else if (_p2 < rhs._p2) return true;
1307                else if (_p2 > rhs._p2) return false;
1308                else return (_p3 < rhs._p3);
1309            }
1310
1311            bool operator == (const Triangle& rhs) const
1312            {
1313                return (_p1 == rhs._p1) && (_p2 != rhs._p2) && (_p3 != rhs._p3);
1314            }
1315
1316            bool operator != (const Triangle& rhs) const
1317            {
1318                return (_p1 != rhs._p1) || (_p2 != rhs._p2) || (_p3 != rhs._p3);
1319            }
1320
1321            void sort()
1322            {
1323                if (_p1>_p2) std::swap(_p1,_p2);
1324                if (_p1>_p3) std::swap(_p1,_p3);
1325                if (_p2>_p3) std::swap(_p2,_p3);
1326            }
1327
1328            Edge* oppositeActiveEdge(Edge* edge)
1329            {
1330                if (edge!=_e1 &&  edge!=_e2 &&  edge!=_e3)
1331                {
1332                    OSG_INFO<<"Edge problem"<<std::endl;
1333                    return 0;
1334                }
1335
1336                if (edge!=_e1 && _e1 && _e1->_intersectionType!=Edge::NO_INTERSECTION) return _e1;
1337                if (edge!=_e2 && _e2 && _e2->_intersectionType!=Edge::NO_INTERSECTION) return _e2;
1338                if (edge!=_e3 && _e3 && _e3->_intersectionType!=Edge::NO_INTERSECTION) return _e3;
1339                return 0;
1340            }
1341
1342
1343            unsigned int _p1;
1344            unsigned int _p2;
1345            unsigned int _p3;
1346
1347            Edge* _e1;
1348            Edge* _e2;
1349            Edge* _e3;
1350        };
1351
1352
1353        struct Region
1354        {
1355            enum Classification
1356            {
1357                INSIDE = -1,
1358                INTERSECTS = 0,
1359                OUTSIDE = 1
1360            };
1361
1362            Region():
1363                _radiusSurface(OUTSIDE),
1364                _leftRightSurfaces(OUTSIDE),
1365                _leftSurface(OUTSIDE),
1366                _rightSurface(OUTSIDE),
1367                _bottomSurface(OUTSIDE),
1368                _topSurface(OUTSIDE) {}
1369
1370
1371
1372            void classify(const osg::Vec3& vertex, double radius2, double azimMin, double azimMax, double elevMin, double elevMax)
1373            {
1374                double azimCenter = (azimMax+azimMin)*0.5;
1375                double azimRange = (azimMax-azimMin)*0.5;
1376
1377                double rad2 = vertex.length2();
1378                double length_xy = sqrtf(vertex.x()*vertex.x() + vertex.y()*vertex.y());
1379                double elevation = atan2((double)vertex.z(),length_xy);
1380
1381                // radius surface
1382                if      (rad2 > radius2) _radiusSurface = OUTSIDE;
1383                else if (rad2 < radius2) _radiusSurface = INSIDE;
1384                else                     _radiusSurface = INTERSECTS;
1385
1386                // bottom surface
1387                if      (elevation<elevMin) _bottomSurface = OUTSIDE;
1388                else if (elevation>elevMin) _bottomSurface = INSIDE;
1389                else                        _bottomSurface = INTERSECTS;
1390
1391                // top surface
1392                if      (elevation>elevMax) _topSurface = OUTSIDE;
1393                else if (elevation<elevMax) _topSurface = INSIDE;
1394                else                        _topSurface = INTERSECTS;
1395
1396
1397                double dot_alphaMin = cos(azimMin)*vertex.x() - sin(azimMin)*vertex.y();
1398                if (dot_alphaMin<0.0) _leftSurface = OUTSIDE;
1399                else if (dot_alphaMin>0.0) _leftSurface = INSIDE;
1400                else _leftSurface = INTERSECTS;
1401
1402                double dot_alphaMax = cos(azimMax)*vertex.x() - sin(azimMax)*vertex.y();
1403                if (dot_alphaMax>0.0) _rightSurface = OUTSIDE;
1404                else if (dot_alphaMax<0.0) _rightSurface = INSIDE;
1405                else _rightSurface = INTERSECTS;
1406
1407                double azim = atan2(vertex.x(),vertex.y());
1408                double azimDelta1 = azim-azimCenter;
1409                double azimDelta2 = 2.0*osg::PI + azim-azimCenter;
1410                double azimDelta = std::min(fabs(azimDelta1), fabs(azimDelta2));
1411                if (azimDelta>azimRange)
1412                {
1413                    _leftRightSurfaces = OUTSIDE;
1414                }
1415                else if (azimDelta<azimRange)
1416                {
1417                    _leftRightSurfaces = INSIDE;
1418                }
1419                else if (azimDelta==azimRange)
1420                {
1421                    _leftRightSurfaces = INTERSECTS;
1422                }
1423            }
1424
1425            Classification _radiusSurface;
1426            Classification _leftRightSurfaces;
1427            Classification _leftSurface;
1428            Classification _rightSurface;
1429            Classification _bottomSurface;
1430            Classification _topSurface;
1431        };
1432
1433        struct RegionCounter
1434        {
1435            RegionCounter():
1436                _numVertices(0),
1437                _outside_radiusSurface(0),
1438                _inside_radiusSurface(0),
1439                _intersects_radiusSurface(0),
1440                _outside_leftRightSurfaces(0),
1441                _inside_leftRightSurfaces(0),
1442                _intersects_leftRightSurfaces(0),
1443                _outside_leftSurface(0),
1444                _inside_leftSurface(0),
1445                _intersects_leftSurface(0),
1446                _outside_rightSurface(0),
1447                _inside_rightSurface(0),
1448                _intersects_rightSurface(0),
1449                _outside_bottomSurface(0),
1450                _inside_bottomSurface(0),
1451                _intersects_bottomSurface(0),
1452                _outside_topSurface(0),
1453                _inside_topSurface(0),
1454                _intersects_topSurface(0) {}
1455
1456
1457            void add(const Region& region)
1458            {
1459                ++_numVertices;
1460
1461                if (region._radiusSurface == Region::OUTSIDE)       ++_outside_radiusSurface;
1462                if (region._radiusSurface == Region::INSIDE)        ++_inside_radiusSurface;
1463                if (region._radiusSurface == Region::INTERSECTS)    ++_intersects_radiusSurface;
1464
1465                if (region._leftRightSurfaces == Region::OUTSIDE)    ++_outside_leftRightSurfaces;
1466                if (region._leftRightSurfaces == Region::INSIDE)     ++_inside_leftRightSurfaces;
1467                if (region._leftRightSurfaces == Region::INTERSECTS) ++_intersects_leftRightSurfaces;
1468
1469                if (region._leftSurface == Region::OUTSIDE)         ++_outside_leftSurface;
1470                if (region._leftSurface == Region::INSIDE)          ++_inside_leftSurface;
1471                if (region._leftSurface == Region::INTERSECTS)      ++_intersects_leftSurface;
1472
1473                if (region._rightSurface == Region::OUTSIDE)        ++_outside_rightSurface;
1474                if (region._rightSurface == Region::INSIDE)         ++_inside_rightSurface;
1475                if (region._rightSurface == Region::INTERSECTS)     ++_intersects_rightSurface;
1476
1477                if (region._bottomSurface == Region::OUTSIDE)       ++_outside_bottomSurface;
1478                if (region._bottomSurface == Region::INSIDE)        ++_inside_bottomSurface;
1479                if (region._bottomSurface == Region::INTERSECTS)    ++_intersects_bottomSurface;
1480
1481                if (region._topSurface == Region::OUTSIDE)          ++_outside_topSurface;
1482                if (region._topSurface == Region::INSIDE)           ++_inside_topSurface;
1483                if (region._topSurface == Region::INTERSECTS)       ++_intersects_topSurface;
1484            }
1485
1486            Region::Classification overallClassification() const
1487            {
1488                // if all vertices are outside any of the surfaces then we are completely outside
1489                if (_outside_radiusSurface==_numVertices ||
1490                    _outside_leftRightSurfaces==_numVertices ||
1491                    _outside_topSurface==_numVertices ||
1492                    _outside_bottomSurface==_numVertices) return Region::OUTSIDE;
1493
1494                // if all the vertices on all the sides and inside then we are completely inside
1495                if (_inside_radiusSurface==_numVertices &&
1496                    _inside_leftRightSurfaces==_numVertices &&
1497                    _inside_topSurface==_numVertices &&
1498                    _inside_bottomSurface==_numVertices) return Region::INSIDE;
1499
1500                return Region::INTERSECTS;
1501            }
1502
1503            bool intersecting(SurfaceType surfaceType) const
1504            {
1505                // if all vertices are outside any of the surfaces then we are completely outside
1506                if ((surfaceType!=RADIUS_SURFACE && _outside_radiusSurface!=0) ||
1507                    (surfaceType!=LEFT_SURFACE && _outside_leftSurface!=0) ||
1508                    (surfaceType!=RIGHT_SURFACE && _outside_rightSurface!=0) ||
1509                    (surfaceType!=TOP_SURFACE && _outside_topSurface!=0) ||
1510                    (surfaceType!=BOTTOM_SURFACE && _outside_bottomSurface!=0)) return false;
1511
1512                // if all the vertices on all the sides and inside then we are completely inside
1513                if ((surfaceType!=RADIUS_SURFACE && _inside_radiusSurface!=0) &&
1514                    (surfaceType!=LEFT_SURFACE && _inside_leftSurface!=0) &&
1515                    (surfaceType!=RIGHT_SURFACE && _inside_rightSurface!=0) &&
1516                    (surfaceType!=TOP_SURFACE && _inside_topSurface!=0) &&
1517                    (surfaceType!=BOTTOM_SURFACE && _inside_bottomSurface!=0)) return false;
1518
1519                return true;
1520            }
1521
1522            int numberOfIntersectingSurfaces() const
1523            {
1524                int sidesThatIntersect = 0;
1525                if (_outside_radiusSurface!=_numVertices && _inside_radiusSurface!=_numVertices) ++sidesThatIntersect;
1526                if (_outside_leftSurface!=_numVertices && _inside_leftSurface!=_numVertices) ++sidesThatIntersect;
1527                if (_outside_rightSurface!=_numVertices && _inside_rightSurface!=_numVertices) ++sidesThatIntersect;
1528                if (_outside_topSurface!=_numVertices && _inside_topSurface!=_numVertices) ++sidesThatIntersect;
1529                if (_outside_bottomSurface!=_numVertices && _inside_bottomSurface!=_numVertices) ++sidesThatIntersect;
1530                return sidesThatIntersect;
1531            }
1532
1533            unsigned int _numVertices;
1534            unsigned int _outside_radiusSurface;
1535            unsigned int _inside_radiusSurface;
1536            unsigned int _intersects_radiusSurface;
1537
1538            unsigned int _outside_leftRightSurfaces;
1539            unsigned int _inside_leftRightSurfaces;
1540            unsigned int _intersects_leftRightSurfaces;
1541
1542            unsigned int _outside_leftSurface;
1543            unsigned int _inside_leftSurface;
1544            unsigned int _intersects_leftSurface;
1545
1546            unsigned int _outside_rightSurface;
1547            unsigned int _inside_rightSurface;
1548            unsigned int _intersects_rightSurface;
1549
1550            unsigned int _outside_bottomSurface;
1551            unsigned int _inside_bottomSurface;
1552            unsigned int _intersects_bottomSurface;
1553
1554            unsigned int _outside_topSurface;
1555            unsigned int _inside_topSurface;
1556            unsigned int _intersects_topSurface;
1557
1558        };
1559
1560
1561        typedef std::vector< osg::Vec3 >                            VertexArray;
1562        typedef std::vector< Region >                               RegionArray;
1563        typedef std::vector< bool >                                 BoolArray;
1564        typedef std::vector< unsigned int >                         IndexArray;
1565        typedef std::vector< osg::ref_ptr<Triangle> >               TriangleArray;
1566        typedef std::set< osg::ref_ptr<Edge>, dereference_less >    EdgeSet;
1567
1568        VertexArray     _originalVertices;
1569        RegionArray     _regions;
1570        BoolArray       _vertexInIntersectionSet;
1571        IndexArray      _candidateVertexIndices;
1572        IndexArray      _remapIndices;
1573        TriangleArray   _triangles;
1574        EdgeSet         _edges;
1575
1576        osg::Vec3       _centre;
1577        double          _radius;
1578        double          _azMin, _azMax, _elevMin, _elevMax;
1579
1580        unsigned int    _numOutside;
1581        unsigned int    _numInside;
1582        unsigned int    _numIntersecting;
1583
1584        SphereSegment::LineList _generatedLines;
1585
1586        void computePositionAndRegions(const osg::Matrixd& matrix, osg::Vec3Array& array)
1587        {
1588            _originalVertices.resize(array.size());
1589            _regions.resize(array.size());
1590            _vertexInIntersectionSet.resize(array.size(), false);
1591            _candidateVertexIndices.clear();
1592
1593            double radius2 = _radius*_radius;
1594
1595            for(unsigned int i=0; i<array.size(); ++i)
1596            {
1597                osg::Vec3 vertex = array[i]*matrix - _centre;
1598                _originalVertices[i] = vertex;
1599                _regions[i].classify(vertex, radius2, _azMin, _azMax, _elevMin, _elevMax);
1600            }
1601
1602        }
1603
1604        inline void operator()(unsigned int p1, unsigned int p2, unsigned int p3)
1605        {
1606            RegionCounter rc;
1607            rc.add(_regions[p1]);
1608            rc.add(_regions[p2]);
1609            rc.add(_regions[p3]);
1610
1611            Region::Classification classification = rc.overallClassification();
1612
1613            // reject if outside.
1614            if (classification==Region::OUTSIDE)
1615            {
1616                ++_numOutside;
1617                return;
1618            }
1619
1620            if (rc.numberOfIntersectingSurfaces()==0)
1621            {
1622                ++_numInside;
1623                return;
1624            }
1625
1626            ++_numIntersecting;
1627
1628            _triangles.push_back(new Triangle(p1,p2,p3));
1629
1630            if (!_vertexInIntersectionSet[p1])
1631            {
1632                _vertexInIntersectionSet[p1] = true;
1633                _candidateVertexIndices.push_back(p1);
1634            }
1635
1636            if (!_vertexInIntersectionSet[p2])
1637            {
1638                _vertexInIntersectionSet[p2] = true;
1639                _candidateVertexIndices.push_back(p2);
1640            }
1641
1642            if (!_vertexInIntersectionSet[p3])
1643            {
1644                _vertexInIntersectionSet[p3] = true;
1645                _candidateVertexIndices.push_back(p3);
1646            }
1647
1648        }
1649
1650        void removeDuplicateVertices()
1651        {
1652            OSG_INFO<<"Removing duplicates : num vertices in "<<_candidateVertexIndices.size()<<std::endl;
1653
1654            if (_candidateVertexIndices.size()<2) return;
1655
1656            std::sort(_candidateVertexIndices.begin(), _candidateVertexIndices.end(), SortFunctor(_originalVertices));
1657
1658            _remapIndices.resize(_originalVertices.size());
1659            for(unsigned int i=0; i< _originalVertices.size(); ++i)
1660            {
1661                _remapIndices[i] = i;
1662            }
1663
1664            bool verticesRemapped = false;
1665            IndexArray::iterator itr = _candidateVertexIndices.begin();
1666            unsigned int lastUniqueIndex = *(itr++);
1667            for(; itr != _candidateVertexIndices.end(); ++itr)
1668            {
1669                //unsigned int i = *itr;
1670                // OSG_INFO<<"  i="<<i<<" lastUniqueIndex="<<lastUniqueIndex<<std::endl;
1671                if (_originalVertices[*itr]==_originalVertices[lastUniqueIndex])
1672                {
1673                    OSG_INFO<<"Combining vertex "<<*itr<<" with "<<lastUniqueIndex<<std::endl;
1674                    _remapIndices[*itr] = lastUniqueIndex;
1675                    verticesRemapped = true;
1676                }
1677                else
1678                {
1679                    lastUniqueIndex = *itr;
1680                }
1681            }
1682
1683            if (verticesRemapped)
1684            {
1685                OSG_INFO<<"Remapping triangle vertices "<<std::endl;
1686                for(TriangleArray::iterator titr = _triangles.begin();
1687                    titr != _triangles.end();
1688                    ++titr)
1689                {
1690                    (*titr)->_p1 = _remapIndices[(*titr)->_p1];
1691                    (*titr)->_p2 = _remapIndices[(*titr)->_p2];
1692                    (*titr)->_p3 = _remapIndices[(*titr)->_p3];
1693                    (*titr)->sort();
1694                }
1695            }
1696        }
1697
1698        void removeDuplicateTriangles()
1699        {
1700            OSG_INFO<<"Removing duplicate triangles : num triangles in "<<_triangles.size()<<std::endl;
1701
1702            if (_triangles.size()<2) return;
1703
1704            std::sort(_triangles.begin(), _triangles.end(), dereference_less());
1705
1706            unsigned int lastUniqueTriangle = 0;
1707            unsigned int numDuplicates = 0;
1708            for(unsigned int i=1; i<_triangles.size(); ++i)
1709            {
1710                if ( *(_triangles[lastUniqueTriangle]) != *(_triangles[i]) )
1711                {
1712                    ++lastUniqueTriangle;
1713                    if (lastUniqueTriangle!=i)
1714                    {
1715                        _triangles[lastUniqueTriangle] = _triangles[i];
1716                    }
1717                }
1718                else
1719                {
1720                    ++numDuplicates;
1721                }
1722            }
1723            if (lastUniqueTriangle<_triangles.size()-1)
1724            {
1725                _triangles.erase(_triangles.begin()+lastUniqueTriangle+1, _triangles.end());
1726            }
1727
1728            OSG_INFO<<"Removed duplicate triangles : num duplicates found "<<numDuplicates<<std::endl;
1729            OSG_INFO<<"Removed duplicate triangles : num triangles out "<<_triangles.size()<<std::endl;
1730        }
1731
1732        void buildEdges(Triangle* tri)
1733        {
1734            tri->_e1 = addEdge(tri->_p1, tri->_p2, tri);
1735            tri->_e2 = addEdge(tri->_p2, tri->_p3, tri);
1736            tri->_e3 = addEdge(tri->_p1, tri->_p3, tri);
1737        }
1738
1739        void buildEdges()
1740        {
1741            _edges.clear();
1742            for(TriangleArray::iterator itr = _triangles.begin();
1743                itr != _triangles.end();
1744                ++itr)
1745            {
1746                Triangle* tri = itr->get();
1747
1748                RegionCounter rc;
1749                rc.add(_regions[tri->_p1]);
1750                rc.add(_regions[tri->_p2]);
1751                rc.add(_regions[tri->_p3]);
1752                int numIntersections = rc.numberOfIntersectingSurfaces();
1753
1754                if (numIntersections>=1)
1755                {
1756                    buildEdges(tri);
1757                }
1758
1759            }
1760            OSG_INFO<<"Number of edges "<<_edges.size()<<std::endl;
1761
1762            unsigned int numZeroConnections = 0;
1763            unsigned int numSingleConnections = 0;
1764            unsigned int numDoubleConnections = 0;
1765            unsigned int numMultiConnections = 0;
1766            OSG_INFO<<"Number of edges "<<_edges.size()<<std::endl;
1767            for(EdgeSet::iterator eitr = _edges.begin();
1768                eitr != _edges.end();
1769                ++eitr)
1770            {
1771                const Edge* edge = eitr->get();
1772                unsigned int numConnections = edge->_triangles.size();
1773                if (numConnections==0) ++numZeroConnections;
1774                else if (numConnections==1) ++numSingleConnections;
1775                else if (numConnections==2) ++numDoubleConnections;
1776                else ++numMultiConnections;
1777            }
1778
1779            OSG_INFO<<"Number of numZeroConnections "<<numZeroConnections<<std::endl;
1780            OSG_INFO<<"Number of numSingleConnections "<<numSingleConnections<<std::endl;
1781            OSG_INFO<<"Number of numDoubleConnections "<<numDoubleConnections<<std::endl;
1782            OSG_INFO<<"Number of numMultiConnections "<<numMultiConnections<<std::endl;
1783        }
1784
1785        Edge* addEdge(unsigned int p1, unsigned int p2, Triangle* tri)
1786        {
1787
1788            osg::ref_ptr<Edge> edge = new Edge(p1, p2);
1789            EdgeSet::iterator itr = _edges.find(edge);
1790            if (itr==_edges.end())
1791            {
1792                edge->addTriangle(tri);
1793                _edges.insert(edge);
1794                return edge.get();
1795            }
1796            else
1797            {
1798                Edge* edge = const_cast<Edge*>(itr->get());
1799                edge->addTriangle(tri);
1800                return edge;
1801            }
1802        }
1803
1804        SphereSegment::LineList connectIntersections(EdgeList& hitEdges)
1805        {
1806            SphereSegment::LineList lineList;
1807
1808            OSG_INFO<<"Number of edge intersections "<<hitEdges.size()<<std::endl;
1809
1810            if (hitEdges.empty()) return lineList;
1811
1812            // now need to build the toTraverse list for each hit edge,
1813            // but should only contain triangles that actually hit the intersection surface
1814            EdgeList::iterator hitr;
1815            for(hitr = hitEdges.begin();
1816                hitr != hitEdges.end();
1817                ++hitr)
1818            {
1819                Edge* edge = hitr->get();
1820                edge->_toTraverse.clear();
1821                //OSG_INFO<<"edge= "<<edge<<std::endl;
1822                for(Edge::TriangleList::iterator titr = edge->_triangles.begin();
1823                    titr != edge->_triangles.end();
1824                    ++titr)
1825                {
1826                    Triangle* tri = *titr;
1827
1828                    // count how many active edges there are on this triangle
1829                    unsigned int numActiveEdges = 0;
1830                    unsigned int numEdges = 0;
1831                    if (tri->_e1 && tri->_e1->_intersectionType!=Edge::NO_INTERSECTION) ++numActiveEdges;
1832                    if (tri->_e2 && tri->_e2->_intersectionType!=Edge::NO_INTERSECTION) ++numActiveEdges;
1833                    if (tri->_e3 && tri->_e3->_intersectionType!=Edge::NO_INTERSECTION) ++numActiveEdges;
1834
1835                    if (tri->_e1) ++numEdges;
1836                    if (tri->_e2) ++numEdges;
1837                    if (tri->_e3) ++numEdges;
1838
1839                    // if we have one or more then add it into the edges to traverse list
1840                    if (numActiveEdges>1)
1841                    {
1842                        //OSG_INFO<<"   adding tri="<<tri<<std::endl;
1843                        edge->_toTraverse.push_back(tri);
1844                    }
1845
1846                    // OSG_INFO<<"Number active edges "<<numActiveEdges<<" num original edges "<<numEdges<<std::endl;
1847                }
1848            }
1849
1850            while(!hitEdges.empty())
1851            {
1852                // find the an open edge
1853                for(hitr = hitEdges.begin();
1854                    hitr != hitEdges.end();
1855                    ++hitr)
1856                {
1857                    Edge* edge = hitr->get();
1858                    if (edge->_toTraverse.size()==1) break;
1859                }
1860
1861                if (hitr == hitEdges.end())
1862                {
1863                    hitr = hitEdges.begin();
1864                }
1865
1866                // OSG_INFO<<"New line "<<std::endl;
1867
1868
1869                osg::Vec3Array* newLine = new osg::Vec3Array;
1870                lineList.push_back(newLine);
1871
1872                Edge* edge = hitr->get();
1873                while (edge)
1874                {
1875                    // OSG_INFO<<"   vertex "<<edge->_intersectionVertex<<std::endl;
1876                    newLine->push_back(edge->_intersectionVertex+_centre/*+osg::Vec3(0.0f,0.0f,200.0f)*/);
1877
1878                    Edge* newEdge = 0;
1879
1880                    Triangle* tri = !(edge->_toTraverse.empty()) ? edge->_toTraverse.back() : 0;
1881                    if (tri)
1882                    {
1883
1884                        newEdge = tri->oppositeActiveEdge(edge);
1885
1886                        edge->removeFromToTraverseList(tri);
1887                        newEdge->removeFromToTraverseList(tri);
1888
1889                        // OSG_INFO<<"   tri="<<tri<<" edge="<<edge<<" newEdge="<<newEdge<<std::endl;
1890
1891                        if (edge==newEdge)
1892                        {
1893                            OSG_INFO<<"   edge returned to itself problem "<<std::endl;
1894                        }
1895                    }
1896                    else
1897                    {
1898                        newEdge = 0;
1899                    }
1900
1901                    if (edge->_toTraverse.empty())
1902                    {
1903                        edge->_intersectionType = Edge::NO_INTERSECTION;
1904
1905                        // remove edge for the hitEdges.
1906                        hitr = std::find(hitEdges.begin(), hitEdges.end(), edge);
1907                        if (hitr!=hitEdges.end()) hitEdges.erase(hitr);
1908                    }
1909
1910                    // move on to next edge in line.
1911                    edge = newEdge;
1912
1913                }
1914            }
1915            return lineList;
1916        }
1917
1918        template<class I>
1919        SphereSegment::LineList computeIntersections(I intersector)
1920        {
1921            // collect all the intersecting edges
1922            EdgeList hitEdges;
1923            for(EdgeSet::iterator itr = _edges.begin();
1924                itr != _edges.end();
1925                ++itr)
1926            {
1927                Edge* edge = const_cast<Edge*>(itr->get());
1928                if (intersector(edge))
1929                {
1930                    hitEdges.push_back(edge);
1931                }
1932            }
1933
1934            return connectIntersections(hitEdges);
1935        }
1936
1937        template<class I>
1938        void trim(SphereSegment::LineList& lineList, osg::Vec3Array* sourceLine, I intersector)
1939        {
1940            if (sourceLine->empty()) return;
1941
1942            // OSG_INFO<<"Testing line of "<<sourceLine->size()<<std::endl;
1943
1944            unsigned int first=0;
1945            while (first<sourceLine->size())
1946            {
1947                // find first valid vertex.
1948                for(; first<sourceLine->size(); ++first)
1949                {
1950                    if (intersector.distance((*sourceLine)[first]-_centre)>=0.0) break;
1951                }
1952
1953                if (first==sourceLine->size())
1954                {
1955                    // OSG_INFO<<"No valid points found"<<std::endl;
1956                    return;
1957                }
1958
1959                // find last valid vertex.
1960                unsigned int last = first+1;
1961                for(; last<sourceLine->size(); ++last)
1962                {
1963                    if (intersector.distance((*sourceLine)[last]-_centre)<0.0) break;
1964                }
1965
1966                if (first==0 && last==sourceLine->size())
1967                {
1968                    // OSG_INFO<<"Copying complete line"<<std::endl;
1969                    lineList.push_back(sourceLine);
1970                }
1971                else
1972                {
1973                    // OSG_INFO<<"Copying partial line line"<<first<<" "<<last<<std::endl;
1974
1975                    osg::Vec3Array* newLine = new osg::Vec3Array;
1976
1977                    if (first>0)
1978                    {
1979                        newLine->push_back(intersector.intersectionPoint((*sourceLine)[first-1]-_centre, (*sourceLine)[first]-_centre)+_centre);
1980                    }
1981
1982                    for(unsigned int i=first; i<last; ++i)
1983                    {
1984                        newLine->push_back((*sourceLine)[i]);
1985                    }
1986                    if (last<sourceLine->size())
1987                    {
1988                        newLine->push_back(intersector.intersectionPoint((*sourceLine)[last-1]-_centre, (*sourceLine)[last]-_centre)+_centre);
1989                    }
1990
1991                    lineList.push_back(newLine);
1992                }
1993
1994                first = last;
1995            }
1996        }
1997
1998
1999        // handle a paired of surfaces that work to enclose a convex region, which means that
2000        // points can be inside either surface to be valid, and be outside both surfaces to be invalid.
2001        template<class I>
2002        void trim(SphereSegment::LineList& lineList, osg::Vec3Array* sourceLine, I intersector1, I intersector2)
2003        {
2004            if (sourceLine->empty()) return;
2005
2006            // OSG_INFO<<"Testing line of "<<sourceLine->size()<<std::endl;
2007
2008            unsigned int first=0;
2009            while (first<sourceLine->size())
2010            {
2011                // find first valid vertex.
2012                for(; first<sourceLine->size(); ++first)
2013                {
2014                    osg::Vec3 v = (*sourceLine)[first]-_centre;
2015                    if (intersector1.distance(v)>=0.0 || intersector2.distance(v)>=0.0 ) break;
2016                }
2017
2018                if (first==sourceLine->size())
2019                {
2020                    // OSG_INFO<<"No valid points found"<<std::endl;
2021                    return;
2022                }
2023
2024                // find last valid vertex.
2025                unsigned int last = first+1;
2026                for(; last<sourceLine->size(); ++last)
2027                {
2028                    osg::Vec3 v = (*sourceLine)[last]-_centre;
2029                    if (intersector1.distance(v)<0.0 && intersector2.distance(v)<0.0 ) break;
2030                }
2031
2032                if (first==0 && last==sourceLine->size())
2033                {
2034                    // OSG_INFO<<"Copying complete line"<<std::endl;
2035                    lineList.push_back(sourceLine);
2036                }
2037                else
2038                {
2039                    OSG_INFO<<"Copying partial line line"<<first<<" "<<last<<std::endl;
2040
2041                    osg::Vec3Array* newLine = new osg::Vec3Array;
2042
2043                    if (first>0)
2044                    {
2045                        osg::Vec3 start = (*sourceLine)[first-1]-_centre;
2046                        osg::Vec3 end = (*sourceLine)[first]-_centre;
2047
2048                        float end1 = intersector1.distance(end);
2049                        float end2 = intersector2.distance(end);
2050
2051
2052                        // work out which intersector to use by discounting the one that
2053                        // isn't a plausible candidate.
2054                        bool possible1 = end1>=0.0;
2055                        bool possible2 = end2>=0.0;
2056                        if (possible1 && possible2)
2057                        {
2058
2059                            double start1 = intersector1.distance(start);
2060                            double start2 = intersector2.distance(start);
2061
2062                            // need to check which intersection is latest.
2063                            double end1 = intersector1.distance(end);
2064                            double delta1 = (end1-start1);
2065
2066                            double end2 = intersector2.distance(end);
2067                            double delta2 = (end2-start2);
2068
2069                            double r1 = fabs(delta1)>0.0 ? start1/delta1 : 0.0;
2070                            double r2 = fabs(delta2)>0.0 ? start2/delta2 : 0.0;
2071
2072                            // choose intersection which is nearest the end point.
2073                            if (r1<r2)
2074                            {
2075                                OSG_INFO<<"start point, 1 near to end than 2"<<r1<<" "<<r2<<std::endl;
2076                                possible1 = true;
2077                                possible2 = false;
2078                            }
2079                            else
2080                            {
2081                                OSG_INFO<<"start point, 2 near to end than 1"<<std::endl;
2082                                possible1 = false;
2083                                possible2 = true;
2084                            }
2085                        }
2086
2087                        if (possible1)
2088                        {
2089                            newLine->push_back(intersector1.intersectionPoint(start, end)+_centre);
2090                        }
2091                        else
2092                        {
2093                            newLine->push_back(intersector2.intersectionPoint(start, end)+_centre);
2094                        }
2095
2096                    }
2097
2098                    for(unsigned int i=first; i<last; ++i)
2099                    {
2100                        newLine->push_back((*sourceLine)[i]);
2101                    }
2102
2103                    if (last<sourceLine->size())
2104                    {
2105                        osg::Vec3 start = (*sourceLine)[last-1]-_centre;
2106                        osg::Vec3 end = (*sourceLine)[last]-_centre;
2107
2108                        double start1 = intersector1.distance(start);
2109                        double start2 = intersector2.distance(start);
2110
2111                        // work out which intersector to use by discounting the one that
2112                        // isn't a plausible candidate.
2113                        bool possible1 = start1>=0.0;
2114                        bool possible2 = start2>=0.0;
2115                        if (possible1 && possible2)
2116                        {
2117                            double end1 = intersector1.distance(end);
2118                            double end2 = intersector2.distance(end);
2119
2120                            possible1 = end1<0.0;
2121                            possible2 = end2<0.0;
2122
2123                            if (possible1 && possible2)
2124                            {
2125                                // need to check which intersection is latest.
2126                                double end1 = intersector1.distance(end);
2127                                double delta1 = (end1-start1);
2128
2129                                double end2 = intersector2.distance(end);
2130                                double delta2 = (end2-start2);
2131
2132                                double r1 = fabs(delta1)>0.0 ? start1/delta1 : 0.0;
2133                                double r2 = fabs(delta2)>0.0 ? start2/delta2 : 0.0;
2134
2135                                // choose intersection which is nearest the end point.
2136                                if (r1>r2)
2137                                {
2138                                    OSG_INFO<<"end point, 1 near to end than 2"<<r1<<" "<<r2<<std::endl;
2139                                    possible1 = true;
2140                                    possible2 = false;
2141                                }
2142                                else
2143                                {
2144                                    OSG_INFO<<"end point, 2 near to end than 1"<<std::endl;
2145                                    possible1 = false;
2146                                    possible2 = true;
2147                                }
2148                            }
2149                        }
2150
2151                        if (possible1)
2152                        {
2153                            newLine->push_back(intersector1.intersectionPoint(start, end)+_centre);
2154                        }
2155                        else
2156                        {
2157                            newLine->push_back(intersector2.intersectionPoint(start, end)+_centre);
2158                        }
2159
2160                    }
2161
2162                    lineList.push_back(newLine);
2163                }
2164
2165                first = last;
2166            }
2167        }
2168
2169        template<class I>
2170        void trim(SphereSegment::LineList& lineList, I intersector)
2171        {
2172            SphereSegment::LineList newLines;
2173
2174            // collect all the intersecting edges
2175            for(SphereSegment::LineList::iterator itr = lineList.begin();
2176                itr != lineList.end();
2177                ++itr)
2178            {
2179                osg::Vec3Array* line = itr->get();
2180                trim(newLines, line, intersector);
2181            }
2182            lineList.swap(newLines);
2183        }
2184
2185
2186        template<class I>
2187        void trim(SphereSegment::LineList& lineList, I intersector1, I intersector2)
2188        {
2189            SphereSegment::LineList newLines;
2190
2191            // collect all the intersecting edges
2192            for(SphereSegment::LineList::iterator itr = lineList.begin();
2193                itr != lineList.end();
2194                ++itr)
2195            {
2196                osg::Vec3Array* line = itr->get();
2197                trim(newLines, line, intersector1, intersector2);
2198            }
2199            lineList.swap(newLines);
2200        }
2201
2202        struct LinePair
2203        {
2204            LinePair(osg::Vec3Array* line):
2205                _line(line),
2206                _lineEnd(0),
2207                _neighbourLine(0),
2208                _neighbourLineEnd(0),
2209                _distance(FLT_MAX) {}
2210
2211            bool operator < (const LinePair& linePair) const
2212            {
2213                return _distance < linePair._distance;
2214            }
2215
2216            void consider(osg::Vec3Array* testline)
2217            {
2218                if (_neighbourLine.valid())
2219                {
2220                    float distance = ((*_line)[0]-(*testline)[0]).length();
2221                    if (distance<_distance)
2222                    {
2223                        _lineEnd = 0;
2224                        _neighbourLine = testline;
2225                        _neighbourLineEnd = 0;
2226                        _distance = distance;
2227                    }
2228
2229                    distance = ((*_line)[0]-(*testline)[testline->size()-1]).length();
2230                    if (distance<_distance)
2231                    {
2232                        _lineEnd = 0;
2233                        _neighbourLine = testline;
2234                        _neighbourLineEnd = testline->size()-1;
2235                        _distance = distance;
2236                    }
2237
2238                    distance = ((*_line)[_line->size()-1]-(*testline)[0]).length();
2239                    if (distance<_distance)
2240                    {
2241                        _lineEnd = _line->size()-1;
2242                        _neighbourLine = testline;
2243                        _neighbourLineEnd = 0;
2244                        _distance = distance;
2245                    }
2246
2247                    distance = ((*_line)[_line->size()-1]-(*testline)[testline->size()-1]).length();
2248                    if (distance<_distance)
2249                    {
2250                        _lineEnd = _line->size()-1;
2251                        _neighbourLine = testline;
2252                        _neighbourLineEnd = testline->size()-1;
2253                        _distance = distance;
2254                    }
2255
2256                }
2257                else
2258                {
2259                    _neighbourLine = testline;
2260                    if (_neighbourLine==_line)
2261                    {
2262                        _lineEnd = 0;
2263                        _neighbourLineEnd = _neighbourLine->size()-1;
2264                        _distance = ((*_line)[_lineEnd]-(*_neighbourLine)[_neighbourLineEnd]).length();
2265                    }
2266                    else
2267                    {
2268                        _distance = ((*_line)[0]-(*_neighbourLine)[0]).length();
2269                        _lineEnd = 0;
2270                        _neighbourLineEnd = 0;
2271
2272                        float distance = ((*_line)[0]-(*_neighbourLine)[_neighbourLine->size()-1]).length();
2273                        if (distance<_distance)
2274                        {
2275                            _lineEnd = 0;
2276                            _neighbourLineEnd = _neighbourLine->size()-1;
2277                            _distance = distance;
2278                        }
2279
2280                        distance = ((*_line)[_line->size()-1]-(*_neighbourLine)[0]).length();
2281                        if (distance<_distance)
2282                        {
2283                            _lineEnd = _line->size()-1;
2284                            _neighbourLineEnd = 0;
2285                            _distance = distance;
2286                        }
2287
2288                        distance = ((*_line)[_line->size()-1]-(*_neighbourLine)[_neighbourLine->size()-1]).length();
2289                        if (distance<_distance)
2290                        {
2291                            _lineEnd = _line->size()-1;
2292                            _neighbourLineEnd = _neighbourLine->size()-1;
2293                            _distance = distance;
2294                        }
2295                    }
2296                }
2297            };
2298
2299            bool contains(osg::Vec3Array* line) const
2300            {
2301                return _line==line || _neighbourLine==line;
2302            }
2303
2304
2305            osg::ref_ptr<osg::Vec3Array>    _line;
2306            unsigned int                    _lineEnd;
2307            osg::ref_ptr<osg::Vec3Array>    _neighbourLine;
2308            unsigned int                    _neighbourLineEnd;
2309            float                           _distance;
2310        };
2311
2312        void joinEnds(float fuseDistance, bool doFuse, bool allowJoinToSelf)
2313        {
2314            SphereSegment::LineList fusedLines;
2315            SphereSegment::LineList unfusedLines;
2316
2317            // first separate the already fused lines from the unfused ones.
2318            for(SphereSegment::LineList::iterator itr = _generatedLines.begin();
2319                itr != _generatedLines.end();
2320                ++itr)
2321            {
2322                osg::Vec3Array* line = itr->get();
2323                if (line->size()>=2)
2324                {
2325                    if ((*line)[0]==(*line)[line->size()-1])
2326                    {
2327                        fusedLines.push_back(line);
2328                    }
2329                    else
2330                    {
2331                        unfusedLines.push_back(line);
2332                    }
2333                }
2334            }
2335
2336            while (unfusedLines.size()>=1)
2337            {
2338                // generate a set of line pairs to establish which
2339                // line pair has the minimum distance.
2340                typedef std::multiset<LinePair> LinePairSet;
2341                LinePairSet linePairs;
2342                for(unsigned int i=0; i<unfusedLines.size(); ++i)
2343                {
2344                    unsigned int j = allowJoinToSelf ? i : i+1;
2345                    if (j<unfusedLines.size())
2346                    {
2347                        LinePair linePair(unfusedLines[i].get());
2348                        for(; j<unfusedLines.size(); ++j)
2349                        {
2350                            linePair.consider(unfusedLines[j].get());
2351                        }
2352                        linePairs.insert(linePair);
2353                    }
2354                }
2355
2356                if (linePairs.empty())
2357                {
2358                    OSG_INFO<<"Line Pairs empty"<<std::endl;
2359                    break;
2360                }
2361
2362                for(LinePairSet::iterator itr = linePairs.begin();
2363                    itr != linePairs.end();
2364                    ++itr)
2365                {
2366                    OSG_INFO<<"Line "<<itr->_line.get()<<" "<<itr->_lineEnd<<"  neighbour "<<itr->_neighbourLine.get()<<" "<<itr->_neighbourLineEnd<<" distance="<<itr->_distance<<std::endl;
2367                }
2368
2369                LinePair linePair = *linePairs.begin();
2370                if (linePair._distance > fuseDistance)
2371                {
2372                    OSG_INFO<<"Completed work, shortest distance left is "<<linePair._distance<<std::endl;
2373                    break;
2374                }
2375
2376                if (linePair._line == linePair._neighbourLine)
2377                {
2378                    OSG_INFO<<"Fusing line to itself"<<std::endl;
2379                    osg::Vec3Array* line = linePair._line.get();
2380                    osg::Vec3 average = ((*line)[0]+(*line)[line->size()-1])*0.5f;
2381
2382                    if (doFuse)
2383                    {
2384                        (*line)[0] = average;
2385                        (*line)[line->size()-1] = average;
2386                    }
2387                    else
2388                    {
2389                        // add start of line to end.
2390                        line->push_back((*line)[0]);
2391                    }
2392
2393                    fusedLines.push_back(line);
2394
2395                    SphereSegment::LineList::iterator fitr = std::find(unfusedLines.begin(), unfusedLines.end(), line);
2396                    if (fitr != unfusedLines.end())
2397                    {
2398                        unfusedLines.erase(fitr);
2399                    }
2400                    else
2401                    {
2402                        OSG_INFO<<"Error couldn't find line in unfused list, exiting fusing loop."<<std::endl;
2403                        break;
2404                    }
2405                }
2406                else
2407                {
2408
2409                    osg::Vec3Array* line1 = linePair._line.get();
2410                    int fuseEnd1 =  linePair._lineEnd;
2411                    int openEnd1 = fuseEnd1==0 ? line1->size()-1 : 0;
2412                    int direction1 = openEnd1<fuseEnd1 ? 1 : -1;
2413
2414                    osg::Vec3Array* line2 = linePair._neighbourLine.get();
2415                    int fuseEnd2 =  linePair._neighbourLineEnd;
2416                    int openEnd2 = fuseEnd2==0 ? line2->size()-1 : 0;
2417                    int direction2 = fuseEnd2<openEnd2 ? 1 : -1;
2418
2419                    osg::Vec3Array* newline = new osg::Vec3Array;
2420
2421                    // copy across all but fuse end of line1
2422                    for(int i=openEnd1;
2423                        i != fuseEnd1;
2424                        i += direction1)
2425                    {
2426                        newline->push_back((*line1)[i]);
2427                    }
2428
2429                    // add the average of the two fused ends
2430
2431                    if (doFuse)
2432                    {
2433                        osg::Vec3 average = ((*line1)[fuseEnd1] + (*line2)[fuseEnd2])*0.5f;
2434                        newline->push_back(average);
2435                    }
2436                    else
2437                    {
2438                        newline->push_back((*line1)[fuseEnd1]);
2439                        newline->push_back((*line2)[fuseEnd2]);
2440                    }
2441
2442                    // copy across from the next point in from fuseEnd2 to the openEnd2.
2443                    for(int j=fuseEnd2 + direction2;
2444                        j != openEnd2 + direction2;
2445                        j += direction2)
2446                    {
2447                        newline->push_back((*line2)[j]);
2448                    }
2449
2450                    // remove line1 from unfused list.
2451                    SphereSegment::LineList::iterator fitr = std::find(unfusedLines.begin(), unfusedLines.end(), line1);
2452                    if (fitr != unfusedLines.end())
2453                    {
2454                        unfusedLines.erase(fitr);
2455                    }
2456
2457                    // remove line2 from unfused list.
2458                    fitr = std::find(unfusedLines.begin(), unfusedLines.end(), line2);
2459                    if (fitr != unfusedLines.end())
2460                    {
2461                        unfusedLines.erase(fitr);
2462                    }
2463
2464                    // add the newline into the unfused for further processing.
2465                    unfusedLines.push_back(newline);
2466
2467                    OSG_INFO<<"Fusing two separate lines "<<newline<<std::endl;
2468                }
2469
2470                _generatedLines = fusedLines;
2471                _generatedLines.insert(_generatedLines.end(), unfusedLines.begin(), unfusedLines.end());
2472
2473            }
2474        }
2475
2476    };
2477
2478    bool computeQuadraticSolution(double a, double b, double c, double& s1, double& s2)
2479    {
2480        // avoid division by zero.
2481        if (a==0.0)
2482        {
2483            s1 = 0.0;
2484            s2 = 0.0;
2485            return false;
2486        }
2487
2488        double inside_sqrt = b*b - 4.0*a*c;
2489
2490        // avoid sqrt of negative number
2491        if (inside_sqrt<0.0)
2492        {
2493            s1 = 0.0;
2494            s2 = 0.0;
2495            return false;
2496        }
2497
2498        double rhs = sqrt(inside_sqrt);
2499        s1 = (-b + rhs)/(2.0*a);
2500        s2 = (-b - rhs)/(2.0*a);
2501
2502        return true;
2503    }
2504
2505    struct AzimPlaneIntersector
2506    {
2507        AzimPlaneIntersector(TriangleIntersectOperator& tif, double azim, bool lowerOutside):
2508            _tif(tif),
2509            _lowerOutside(lowerOutside)
2510        {
2511            _plane.set(cos(azim),-sin(azim),0.0,0.0);
2512            _endPlane.set(sin(azim),cos(azim),0.0,0.0);
2513        }
2514
2515        TriangleIntersectOperator& _tif;
2516        osg::Plane _plane;
2517        osg::Plane _endPlane;
2518        bool _lowerOutside;
2519
2520        inline bool operator() (TriangleIntersectOperator::Edge* edge)
2521        {
2522            edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2523
2524            osg::Vec3& v1 = _tif._originalVertices[edge->_p1];
2525            osg::Vec3& v2 = _tif._originalVertices[edge->_p2];
2526
2527            double d1 = _plane.distance(v1);
2528            double d2 = _plane.distance(v2);
2529
2530            edge->_p1Outside = _lowerOutside ? (d1<0.0) : (d1>0.0);
2531            edge->_p2Outside = _lowerOutside ? (d2<0.0) : (d2>0.0);
2532
2533            // if both points inside then discard
2534            if (d1<0.0 && d2<0.0) return false;
2535
2536            // if both points outside then discard
2537            if (d1>0.0 && d2>0.0) return false;
2538
2539            if (d1==0.0)
2540            {
2541                if (d2==0.0)
2542                {
2543                    edge->_intersectionType = TriangleIntersectOperator::Edge::BOTH_ENDS;
2544                }
2545                else
2546                {
2547                    edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_1;
2548                }
2549            }
2550            else if (d2==0.0)
2551            {
2552                edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_2;
2553            }
2554            else
2555            {
2556
2557                double div = d2-d1;
2558                if (div==0.0)
2559                {
2560                    edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2561                    return false;
2562                }
2563
2564                double r =  -d1 / div;
2565                if (r<0.0 || r>1.0)
2566                {
2567                    edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2568                    return false;
2569                }
2570
2571                // OSG_INFO<<"r = "<<r<<std::endl;
2572
2573                double one_minus_r = 1.0-r;
2574
2575                edge->_intersectionType = TriangleIntersectOperator::Edge::MID_POINT;
2576                edge->_intersectionVertex = v1*one_minus_r + v2*r;
2577            }
2578
2579            return true;
2580        }
2581
2582        // compute the intersection between line segment and surface
2583        osg::Vec3 intersectionPoint(const osg::Vec3& v1, const osg::Vec3& v2)
2584        {
2585            double d1 = _plane.distance(v1);
2586            double d2 = _plane.distance(v2);
2587
2588            double div = d2-d1;
2589            if (div==0.0)
2590            {
2591                return v1;
2592            }
2593
2594            double r =  -d1 / div;
2595            double one_minus_r = 1.0-r;
2596
2597            return v1*one_minus_r + v2*r;
2598        }
2599
2600        // positive distance to the inside.
2601        double distance(const osg::Vec3& v)
2602        {
2603            return _lowerOutside ? _plane.distance(v) : -_plane.distance(v) ;
2604        }
2605
2606    protected:
2607
2608        AzimPlaneIntersector& operator = (const AzimPlaneIntersector&) { return *this; }
2609    };
2610
2611    struct ElevationIntersector
2612    {
2613        ElevationIntersector(TriangleIntersectOperator& tif, double elev, bool lowerOutside):
2614            _tif(tif),
2615            _elev(elev),
2616            _lowerOutside(lowerOutside) {}
2617
2618        TriangleIntersectOperator& _tif;
2619        double _elev;
2620        bool _lowerOutside;
2621
2622        inline bool operator() (TriangleIntersectOperator::Edge* edge)
2623        {
2624            edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2625
2626            osg::Vec3& v1 = _tif._originalVertices[edge->_p1];
2627            osg::Vec3& v2 = _tif._originalVertices[edge->_p2];
2628
2629            double length_xy1 = sqrt(v1.x()*v1.x() + v1.y()*v1.y());
2630            double elev1 = atan2((double)v1.z(),length_xy1);
2631
2632            double length_xy2 = sqrt(v2.x()*v2.x() + v2.y()*v2.y());
2633            double elev2 = atan2((double)v2.z(),length_xy2);
2634
2635            edge->_p1Outside = _lowerOutside ? (elev1<_elev) : (elev1>_elev);
2636            edge->_p2Outside = _lowerOutside ? (elev2<_elev) : (elev2>_elev);
2637
2638            // if both points inside then discard
2639            if (elev1<_elev && elev2<_elev) return false;
2640
2641            // if both points outside then discard
2642            if (elev1>_elev && elev2>_elev) return false;
2643
2644            if (elev1==_elev)
2645            {
2646                if (elev2==_elev)
2647                {
2648                    edge->_intersectionType = TriangleIntersectOperator::Edge::BOTH_ENDS;
2649                }
2650                else
2651                {
2652                    edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_1;
2653                }
2654            }
2655            else if (elev2==_elev)
2656            {
2657                edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_2;
2658            }
2659            else
2660            {
2661                double dx = v2.x()-v1.x();
2662                double dy = v2.y()-v1.y();
2663                double dz = v2.z()-v1.z();
2664                double t = tan(_elev);
2665                double tt = t*t;
2666                double a = dz*dz-tt*(dx*dx + dy*dy);
2667                double b = 2.0*(v1.z()*dz - tt*(v1.x()*dx + v1.y()*dy));
2668                double c = v1.z()*v1.z() - tt*(v1.x()*v1.x() + v1.y()*v1.y());
2669
2670                double s1, s2;
2671                if (!computeQuadraticSolution(a,b,c,s1,s2))
2672                {
2673                    edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2674                    return false;
2675                }
2676                double r = 0.0;
2677                if (s1>=0.0 && s1<=1.0)
2678                {
2679                    r = s1;
2680                }
2681                else if (s2>=0.0 && s2<=1.0)
2682                {
2683                    r = s2;
2684                }
2685                else
2686                {
2687                    OSG_INFO<<"neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
2688
2689                    edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2690                    return false;
2691                }
2692
2693                double one_minus_r = 1.0-r;
2694                edge->_intersectionType = TriangleIntersectOperator::Edge::MID_POINT;
2695                edge->_intersectionVertex = v1*one_minus_r + v2*r;
2696            }
2697
2698            return true;
2699        }
2700
2701        // compute the intersection between line segment and surface
2702        osg::Vec3 intersectionPoint(const osg::Vec3& v1, const osg::Vec3& v2)
2703        {
2704            double dx = v2.x()-v1.x();
2705            double dy = v2.y()-v1.y();
2706            double dz = v2.z()-v1.z();
2707            double t = tan(_elev);
2708            double tt = t*t;
2709            double a = dz*dz-tt*(dx*dx + dy*dy);
2710            double b = 2.0*(v1.z()*dz - tt*(v1.x()*dx + v1.y()*dy));
2711            double c = v1.z()*v1.z() - tt*(v1.x()*v1.x() + v1.y()*v1.y());
2712
2713            double s1, s2;
2714            if (!computeQuadraticSolution(a,b,c,s1,s2))
2715            {
2716                OSG_INFO<<"Warning::neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
2717                return v1;
2718            }
2719            double r = 0.0;
2720            if (s1>=0.0 && s1<=1.0)
2721            {
2722                r = s1;
2723            }
2724            else if (s2>=0.0 && s2<=1.0)
2725            {
2726                r = s2;
2727            }
2728            else
2729            {
2730                OSG_INFO<<"Warning::neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
2731                return v1;
2732            }
2733
2734            double one_minus_r = 1.0-r;
2735            return v1*one_minus_r + v2*r;
2736        }
2737
2738        // positive distance to the inside.
2739        double distance(const osg::Vec3& v)
2740        {
2741            double length_xy = sqrt(v.x()*v.x() + v.y()*v.y());
2742            double computedElev = atan2((double)v.z(),length_xy);
2743
2744            return _lowerOutside ? computedElev-_elev : _elev-computedElev ;
2745        }
2746
2747    protected:
2748
2749        ElevationIntersector& operator = (const ElevationIntersector&) { return *this; }
2750
2751    };
2752
2753    struct RadiusIntersector
2754    {
2755        RadiusIntersector(TriangleIntersectOperator& tif):
2756            _tif(tif) {}
2757
2758        TriangleIntersectOperator& _tif;
2759
2760        inline bool operator() (TriangleIntersectOperator::Edge* edge)
2761        {
2762            edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2763
2764            osg::Vec3& v1 = _tif._originalVertices[edge->_p1];
2765            osg::Vec3& v2 = _tif._originalVertices[edge->_p2];
2766            double radius1 = v1.length();
2767            double radius2 = v2.length();
2768
2769            edge->_p1Outside = radius1>_tif._radius;
2770            edge->_p2Outside = radius2>_tif._radius;
2771
2772            // if both points inside then discard
2773            if (radius1<_tif._radius && radius2<_tif._radius) return false;
2774
2775            // if both points outside then discard
2776            if (radius1>_tif._radius && radius2>_tif._radius) return false;
2777
2778            if (radius1==_tif._radius)
2779            {
2780                if (radius2==_tif._radius)
2781                {
2782                    edge->_intersectionType = TriangleIntersectOperator::Edge::BOTH_ENDS;
2783                }
2784                else
2785                {
2786                    edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_1;
2787                }
2788            }
2789            else if (radius2==_tif._radius)
2790            {
2791                edge->_intersectionType = TriangleIntersectOperator::Edge::POINT_2;
2792            }
2793            else
2794            {
2795                double dx = v2.x()-v1.x();
2796                double dy = v2.y()-v1.y();
2797                double dz = v2.z()-v1.z();
2798                double a = dx*dx + dy*dy + dz*dz;
2799                double b = 2.0*(v1.x()*dx + v1.y()*dy + v1.z()*dz);
2800                double c = v1.x()*v1.x() +  v1.y()*v1.y() + v1.z()*v1.z() - _tif._radius*_tif._radius;
2801
2802                double s1, s2;
2803                if (!computeQuadraticSolution(a,b,c,s1,s2))
2804                {
2805                    edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2806                    return false;
2807                }
2808                double r = 0.0;
2809                if (s1>=0.0 && s1<=1.0)
2810                {
2811                    r = s1;
2812                }
2813                else if (s2>=0.0 && s2<=1.0)
2814                {
2815                    r = s2;
2816                }
2817                else
2818                {
2819                    OSG_INFO<<"neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
2820
2821                    edge->_intersectionType = TriangleIntersectOperator::Edge::NO_INTERSECTION;
2822                    return false;
2823                }
2824
2825                double one_minus_r = 1.0-r;
2826                edge->_intersectionType = TriangleIntersectOperator::Edge::MID_POINT;
2827                edge->_intersectionVertex = v1*one_minus_r + v2*r;
2828
2829            }
2830
2831            return true;
2832        }
2833
2834        // compute the intersection between line segment and surface
2835        osg::Vec3 intersectionPoint(const osg::Vec3& v1, const osg::Vec3& v2)
2836        {
2837            double dx = v2.x()-v1.x();
2838            double dy = v2.y()-v1.y();
2839            double dz = v2.z()-v1.z();
2840            double a = dx*dx + dy*dy + dz*dz;
2841            double b = 2.0*(v1.x()*dx + v1.y()*dy + v1.z()*dz);
2842            double c = v1.x()*v1.x() +  v1.y()*v1.y() + v1.z()*v1.z() - _tif._radius*_tif._radius;
2843
2844            double s1, s2;
2845            if (!computeQuadraticSolution(a,b,c,s1,s2))
2846            {
2847                OSG_INFO<<"Warning: neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
2848                return v1;
2849            }
2850            double r = 0.0;
2851            if (s1>=0.0 && s1<=1.0)
2852            {
2853                r = s1;
2854            }
2855            else if (s2>=0.0 && s2<=1.0)
2856            {
2857                r = s2;
2858            }
2859            else
2860            {
2861                OSG_INFO<<"Warning: neither segment intersects s1="<<s1<<" s2="<<s2<<std::endl;
2862                return v1;
2863            }
2864
2865            double one_minus_r = 1.0-r;
2866            return v1*one_minus_r + v2*r;
2867
2868        }
2869
2870        // positive distance to the inside.
2871        double distance(const osg::Vec3& v)
2872        {
2873            return _tif._radius-v.length();
2874        }
2875
2876
2877    protected:
2878
2879        RadiusIntersector& operator = (const RadiusIntersector&) { return *this; }
2880
2881    };
2882}
2883
2884using namespace SphereSegmentIntersector;
2885
2886SphereSegment::LineList SphereSegment::computeIntersection(const osg::Matrixd& matrix, osg::Drawable* drawable)
2887{
2888    // cast to Geometry, return empty handed if Drawable not a Geometry.
2889    osg::Geometry* geometry = dynamic_cast<osg::Geometry*>(drawable);
2890    if (!geometry) return SphereSegment::LineList();
2891
2892    // get vertices from geometry, return empty handed if a Vec3Array not present.
2893    osg::Vec3Array* vertices = dynamic_cast<osg::Vec3Array*>(geometry->getVertexArray());
2894    if (!vertices) return SphereSegment::LineList();
2895
2896    typedef osg::TriangleIndexFunctor<TriangleIntersectOperator> TriangleIntersectFunctor;
2897    TriangleIntersectFunctor tif;
2898
2899    tif._centre = _centre;
2900    tif._radius = _radius;
2901    tif._azMin = _azMin;
2902    tif._azMax = _azMax;
2903    tif._elevMin = _elevMin;
2904    tif._elevMax = _elevMax;
2905
2906    tif.computePositionAndRegions(matrix, *vertices);
2907
2908    // traverse the triangles in the Geometry dedicating intersections
2909    geometry->accept(tif);
2910
2911    OSG_INFO<<"_numOutside = "<<tif._numOutside<<std::endl;
2912    OSG_INFO<<"_numInside = "<<tif._numInside<<std::endl;
2913    OSG_INFO<<"_numIntersecting = "<<tif._numIntersecting<<std::endl;
2914
2915    tif.removeDuplicateVertices();
2916    tif.removeDuplicateTriangles();
2917    tif.buildEdges();
2918
2919
2920    RadiusIntersector radiusIntersector(tif);
2921    AzimPlaneIntersector azMinIntersector(tif,_azMin, true);
2922    AzimPlaneIntersector azMinEndIntersector(tif,_azMin-osg::PI*0.5, true);
2923    AzimPlaneIntersector azMaxIntersector(tif,_azMax, false);
2924    AzimPlaneIntersector azMaxEndIntersector(tif,_azMax-osg::PI*0.5, true);
2925    ElevationIntersector elevMinIntersector(tif,_elevMin, true);
2926    ElevationIntersector elevMaxIntersector(tif,_elevMax, false);
2927
2928    // create the line intersections with the terrain
2929    SphereSegment::LineList radiusLines = tif.computeIntersections(radiusIntersector);
2930    SphereSegment::LineList elevMinLines = tif.computeIntersections(elevMinIntersector);
2931    SphereSegment::LineList elevMaxLines = tif.computeIntersections(elevMaxIntersector);
2932    SphereSegment::LineList azMinLines;
2933    SphereSegment::LineList azMaxLines;
2934
2935    double azimRange = _azMax-_azMin;
2936    if (azimRange<2.0*osg::PI)
2937    {
2938        azMinLines = tif.computeIntersections(azMinIntersector);
2939        azMaxLines = tif.computeIntersections(azMaxIntersector);
2940
2941        // trim the azimuth intersection lines by the radius
2942        tif.trim(azMinLines,radiusIntersector);
2943        tif.trim(azMaxLines,radiusIntersector);
2944
2945        // trim the azim intersection lines by the elevation
2946        tif.trim(azMinLines, elevMinIntersector);
2947        tif.trim(azMaxLines, elevMinIntersector);
2948
2949        // trim the azim intersection lines by the elevation
2950        tif.trim(azMinLines, elevMaxIntersector);
2951        tif.trim(azMaxLines, elevMaxIntersector);
2952
2953        // trim the centeral ends of the azim lines
2954        tif.trim(azMinLines,azMinEndIntersector);
2955        tif.trim(azMaxLines,azMaxEndIntersector);
2956
2957        if (azimRange<=osg::PI)
2958        {
2959            // trim the radius and elevation intersection lines by the azimMin
2960            tif.trim(radiusLines, azMinIntersector);
2961            tif.trim(elevMinLines, azMinIntersector);
2962            tif.trim(elevMaxLines, azMinIntersector);
2963
2964            // trim the radius and elevation intersection lines by the azimMax
2965            tif.trim(radiusLines, azMaxIntersector);
2966            tif.trim(elevMinLines, azMaxIntersector);
2967            tif.trim(elevMaxLines, azMaxIntersector);
2968        }
2969        else
2970        {
2971            // need to have new intersector which handles convex azim planes
2972            tif.trim(radiusLines, azMinIntersector, azMaxIntersector);
2973            tif.trim(elevMinLines, azMinIntersector, azMaxIntersector);
2974            tif.trim(elevMaxLines, azMinIntersector, azMaxIntersector);
2975        }
2976
2977    }
2978
2979    // trim elevation intersection lines by radius
2980    tif.trim(elevMinLines,radiusIntersector);
2981    tif.trim(elevMaxLines,radiusIntersector);
2982
2983    // trim the radius and elevation intersection lines by the elevMin
2984    tif.trim(radiusLines, elevMinIntersector);
2985
2986    // trim the radius and elevation intersection lines by the elevMax
2987    tif.trim(radiusLines, elevMaxIntersector);
2988
2989    // collect all lines together.
2990    tif._generatedLines.insert(tif._generatedLines.end(), radiusLines.begin(), radiusLines.end());
2991    tif._generatedLines.insert(tif._generatedLines.end(), azMinLines.begin(), azMinLines.end());
2992    tif._generatedLines.insert(tif._generatedLines.end(), azMaxLines.begin(), azMaxLines.end());
2993    tif._generatedLines.insert(tif._generatedLines.end(), elevMinLines.begin(), elevMinLines.end());
2994    tif._generatedLines.insert(tif._generatedLines.end(), elevMaxLines.begin(), elevMaxLines.end());
2995
2996    OSG_INFO<<"number of separate lines = "<<tif._generatedLines.size()<<std::endl;
2997
2998    float fuseDistance = 1.0;
2999    tif.joinEnds(fuseDistance, true, true);
3000
3001    OSG_INFO<<"number of separate lines after fuse = "<<tif._generatedLines.size()<<std::endl;
3002
3003    float joinDistance = 1e8;
3004    tif.joinEnds(joinDistance, false, false);
3005    OSG_INFO<<"number of separate lines after join = "<<tif._generatedLines.size()<<std::endl;
3006
3007    tif.joinEnds(joinDistance, false, true);
3008    OSG_INFO<<"number of separate lines after second join = "<<tif._generatedLines.size()<<std::endl;
3009
3010    return tif._generatedLines;
3011}
3012
3013osg::Node* SphereSegment::computeIntersectionSubgraph(const osg::Matrixd& matrix, osg::Drawable* drawable)
3014{
3015    SphereSegment::LineList generatedLines = computeIntersection(matrix, drawable);
3016
3017    osg::Geode* geode = new osg::Geode;
3018
3019    geode->getOrCreateStateSet()->setMode(GL_LIGHTING,osg::StateAttribute::OFF);
3020
3021    for(osgSim::SphereSegment::LineList::iterator itr = generatedLines.begin();
3022       itr != generatedLines.end();
3023       ++itr)
3024    {
3025        osg::Geometry* geom = new osg::Geometry;
3026        geode->addDrawable(geom);
3027
3028        osg::Vec3Array* vertices = itr->get();
3029        geom->setVertexArray(vertices);
3030        geom->addPrimitiveSet(new osg::DrawArrays(GL_LINE_STRIP, 0, vertices->getNumElements()));
3031    }
3032
3033#if 0
3034    float radius = 0.1f;
3035    for(unsigned int i=0; i<tif._originalVertices.size(); ++i)
3036    {
3037        osg::ShapeDrawable* sd = new osg::ShapeDrawable(new osg::Sphere(tif._originalVertices[i]+tif._centre,radius));
3038
3039        TriangleIntersectFunctor::RegionCounter rc;
3040        rc.add(tif._regions[i]);
3041
3042        TriangleIntersectFunctor::Region::Classification region = rc.overallClassification();
3043        if (region==TriangleIntersectFunctor::Region::OUTSIDE)
3044        {
3045            sd->setColor(osg::Vec4(1.0,0.0,0.0,1.0));
3046        }
3047        else if (region==TriangleIntersectFunctor::Region::INSIDE)
3048        {
3049            sd->setColor(osg::Vec4(1.0,1.0,0.0,1.0));
3050        }
3051        else if (region==TriangleIntersectFunctor::Region::INTERSECTS)
3052        {
3053            sd->setColor(osg::Vec4(1.0,1.0,1.0,1.0));
3054        }
3055
3056        geode->addDrawable(sd);
3057    }
3058#endif
3059
3060    return geode;
3061}
Note: See TracBrowser for help on using the browser.