Context Navigation

SphereSegment.cpp @ 10658

Revision 10658, 103.0 kB (checked in by robert, 4 years ago)
Converted all glBegin/glEnd code across to using osg::GLBeginEndAdapter
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
28	using namespace osgSim;
29
30	// Define the collection of nested classes, all Drawables, which make
31	// up the parts of the sphere segment.
32
33	/**
34	SphereSegment::Surface is the Drawable which represents the specified area of the
35	sphere's surface.
36	*/
37	class SphereSegment::Surface: public osg::Drawable
38	{
39	public:
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
60	protected:
61
62	private:
63
64	SphereSegment* _ss;
65	};
66
67	void SphereSegment::Surface::drawImplementation(osg::RenderInfo& renderInfo) const
68	{
69	_ss->Surface_drawImplementation(*renderInfo.getState());
70	}
71
72	osg:: BoundingBox SphereSegment::Surface::computeBound() const
73	{
74	osg:: BoundingBox bbox;
75	_ss->Surface_computeBound(bbox);
76	return bbox;
77	}
78
79
80	/**
81	SphereSegment::EdgeLine is the Drawable which represents the line around the edge
82	of the specified area of the sphere's EdgeLine.
83	*/
84	class SphereSegment::EdgeLine: public osg::Drawable
85	{
86	public:
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
106	protected:
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
119	private:
120
121	SphereSegment* _ss;
122	};
123
124	void SphereSegment::EdgeLine::drawImplementation(osg::RenderInfo& renderInfo) const
125	{
126	_ss->EdgeLine_drawImplementation(*renderInfo.getState());
127	}
128
129	osg::BoundingBox SphereSegment::EdgeLine::computeBound() const
130	{
131	osg::BoundingBox bbox;
132	_ss->EdgeLine_computeBound(bbox);
133	return bbox;
134	}
135
136
137
138	/**
139	SphereSegment::Side is a Drawable which represents one of the
140	planar areas, at either the minimum or maximum azimuth.
141	*/
142	class SphereSegment::Side: public osg::Drawable
143	{
144	public:
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
164	protected:
165
166	virtual osg::BoundingBox computeBound() const;
167
168	private:
169	SphereSegment* _ss;
170	SphereSegment::SideOrientation _planeOrientation;
171	SphereSegment::BoundaryAngle _BoundaryAngle;
172	};
173
174
175	void SphereSegment::Side::drawImplementation(osg::RenderInfo& renderInfo) const
176	{
177	_ss->Side_drawImplementation(*renderInfo.getState(), _planeOrientation, _BoundaryAngle);
178	}
179
180	osg::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	/**
190	SphereSegment::Spoke is a Drawable which represents a spoke.
191	*/
192	class SphereSegment::Spoke: public osg::Drawable
193	{
194	public:
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
215	protected:
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
227	private:
228	SphereSegment* _ss;
229	SphereSegment::BoundaryAngle _azAngle, _elevAngle;
230	};
231
232	void SphereSegment::Spoke::drawImplementation(osg::RenderInfo& renderInfo) const
233	{
234	_ss->Spoke_drawImplementation(*renderInfo.getState(), _azAngle, _elevAngle);
235	}
236
237	osg::BoundingBox SphereSegment::Spoke::computeBound() const
238	{
239	osg::BoundingBox bbox;
240	_ss->Spoke_computeBound(bbox, _azAngle, _elevAngle);
241	return bbox;
242	}
243
244	SphereSegment::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
258	void SphereSegment::setCentre(const osg::Vec3& c)
259	{
260	_centre = c;
261	dirtyAllDrawableDisplayLists();
262	dirtyAllDrawableBounds();
263	dirtyBound();
264	}
265
266	const osg::Vec3& SphereSegment::getCentre() const
267	{
268	return _centre;
269	}
270
271	void SphereSegment::setRadius(float r)
272	{
273	_radius = r;
274	dirtyAllDrawableDisplayLists();
275	dirtyAllDrawableBounds();
276	dirtyBound();
277	}
278
279	float SphereSegment::getRadius() const
280	{
281	return _radius;
282	}
283
284
285	void 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
312	void 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
322	void 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
335	void 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
344	void SphereSegment::setDensity(int density)
345	{
346	_density = density;
347	dirtyAllDrawableDisplayLists();
348	}
349
350	int SphereSegment::getDensity() const
351	{
352	return _density;
353	}
354
355	void 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
372	void SphereSegment::dirtyAllDrawableDisplayLists()
373	{
374	for(DrawableList::iterator itr = _drawables.begin();
375	itr != _drawables.end();
376	++itr)
377	{
378	(*itr)->dirtyDisplayList();
379	}
380	}
381
382	void SphereSegment::dirtyAllDrawableBounds()
383	{
384	for(DrawableList::iterator itr = _drawables.begin();
385	itr != _drawables.end();
386	++itr)
387	{
388	(*itr)->dirtyBound();
389	}
390	}
391
392	void 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
503	bool 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() + _radiuscos(elev)sin(az),
520	_centre.y() + _radiuscos(elev)cos(az),
521	_centre.z() + _radius*sin(elev))
522	);
523	}
524	}
525	return true;
526	}
527
528	void 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() + _radiuscos(_elevMax)sin(az),
549	_centre.y() + _radiuscos(_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() + _radiuscos(_elevMin)sin(az),
561	_centre.y() + _radiuscos(_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() + _radiuscos(elev)sin(_azMin),
574	_centre.y() + _radiuscos(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() + _radiuscos(elev)sin(_azMax),
586	_centre.y() + _radiuscos(elev)cos(_azMax),
587	_centre.z() + _radius*sin(elev));
588	}
589	gl.End();
590	}
591	}
592
593	bool 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() + _radiuscos(_elevMax)sin(az),
607	_centre.y() + _radiuscos(_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() + _radiuscos(_elevMin)sin(az),
617	_centre.y() + _radiuscos(_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() + _radiuscos(elev)sin(_azMin),
628	_centre.y() + _radiuscos(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() + _radiuscos(elev)sin(_azMax),
638	_centre.y() + _radiuscos(elev)cos(_azMax),
639	_centre.z() + _radius*sin(elev));
640	}
641
642	return true;
643	}
644
645	void 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() + _radiuscos(elev)sin(az),
697	_centre.y() + _radiuscos(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() + _radiuscos(elev)sin(az),
724	_centre.y() + _radiuscos(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() + _radiuscos(elev)sin(az),
765	_centre.y() + _radiuscos(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() + _radiuscos(elev)sin(az),
794	_centre.y() + _radiuscos(elev)cos(az),
795	_centre.z() + _radius*sin(elev));
796	}
797	gl.End();
798	}
799
800	}
801	}
802	}
803	}
804
805	bool 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() + _radiuscos(elev)sin(az),
822	_centre.y() + _radiuscos(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() + _radiuscos(elev)sin(az),
836	_centre.y() + _radiuscos(elev)cos(az),
837	_centre.z() + _radius*sin(elev));
838	}
839	}
840
841	return true;
842	}
843
844	void 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() + _radiuscos(elev)sin(az),
858	_centre.y() + _radiuscos(elev)cos(az),
859	_centre.z() + _radius*sin(elev));
860	gl.End();
861	}
862	}
863
864	bool 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() + _radiuscos(elev)sin(az),
871	_centre.y() + _radiuscos(elev)cos(az),
872	_centre.z() + _radius*sin(elev));
873
874	return true;
875	}
876
877	void SphereSegment::setDrawMask(DrawMask dm)
878	{
879	_drawMask=dm;
880	dirtyAllDrawableDisplayLists();
881	dirtyAllDrawableBounds();
882	dirtyBound();
883	}
884
885	struct 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
905	protected:
906
907	ActivateTransparencyOnType& operator = (const ActivateTransparencyOnType&) { return *this; }
908	};
909
910	struct 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
927	protected:
928
929	DeactivateTransparencyOnType& operator = (const DeactivateTransparencyOnType&) { return *this; }
930	};
931
932	void 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
940	void 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
948	void 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
956	void 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
964	void 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
976	class 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
1068	SphereSegment::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
1113	osg::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
1157	namespace 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 = bb - 4.0a*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 = v1one_minus_r + v2r;
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 v1one_minus_r + v2r;
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 = dzdz-tt(dxdx + dydy);
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 = v1one_minus_r + v2r;
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 = dzdz-tt(dxdx + dydy);
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 v1one_minus_r + v2r;
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 = dxdx + dydy + 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 = v1one_minus_r + v2r;
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 = dxdx + dydy + 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 v1one_minus_r + v2r;
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
2875	using namespace SphereSegmentIntersector;
2876
2877	SphereSegment::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
3004	osg::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.

Context Navigation

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

Download in other formats: