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

Revision 11164, 103.0 kB (checked in by robert, 5 years ago)

Replaced use of unsigned int/enum mask combinations with int/enum mask combinations to avoid the need for casting enums to unsigned ints,
and to avoid associated warnings.

Update wrappers to reflect these changes.

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