Context Navigation

SphereSegment.cpp @ 13041

Revision 13041, 102.1 kB (checked in by robert, 14 months ago)
Ran script to remove trailing spaces and tabs
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`

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

Note: See TracBrowser for help on using the browser.

Context Navigation

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

Download in other formats: