TBoundSphere


Croquet-Teapot

Comment:

This is the main bounds test object. It is used for object culling against view planes, ray testing, and for collision detection. It is defined heirarchically for the collision detection and object picking algorithms. 

Things to do:

	Add bidirectional boundsphere references to an octree.
	

Hierarchy:

ProtoObject
Object
TObject
TBoundSphere

Summary:

instance variables:

box children frame globalPosition localPosition normal offset radius radiusSquared side up vertices

Pool:

OpenGLConstants

methods:

instance class
access collision construct mini ball render transform instance creation

Detail:

instance variables:

box
children
frame
globalPosition
localPosition
normal
offset
radius
radiusSquared
side
up
vertices

instance methods:

access
addChild: child


	children ifNil:[children _ OrderedCollection new.].
	children add: child.
box


	^ box.
box: bx


	box _ bx.
children


	^ children.
children: ch


	children _ ch.
frame


	^ frame.
globalPosition


	^ globalPosition.
localPosition


	^ localPosition.
localPosition: pos


	localPosition _ pos.
localPosition: pos radius: rad


	globalPosition _ B3DVector3 new.
	self localPosition: pos.
	self radius: rad.
radius


	^ radius.
radius: rad


	radius _ rad.
	radiusSquared _ rad*rad.
radiusSquared


	^ radiusSquared.
scale: scale


	self radius: radius*scale.
	self localPosition: (self localPosition * scale).
	children do: [:each | (each isMemberOf: TBoundSphere) ifTrue: [each scale: scale]].
sumFaces

	| total |

" Returns the number of faces in the leaves."
	children ifNil:[^ 0].
	vertices ifTrue:[^children size/3].
	total _ 0.
	children do:[:c | total _ total + c sumFaces.].
	^total.
sumNodes

	| total |

" Returns the total number of nodes in the tree."
	children ifNil:[^ 1].
	vertices ifNotNil:[^1].
	total _ 1.
	children do:[:c | total _ total + c sumNodes.].
	^total.
volume

	^4.0 / 3.0 * Float pi * radius * radius * radius

collision
collideFloor: floor transform: trans


	| collideLoc rval vrt |
	
"	(trans localPointToGlobal: localPosition) y - floor > radius ifTrue:[ ^ nil. ]."

	collideLoc _ B3DVector4 x: 0.0 y: Float infinity z: 0.0 w:0.0.

	vertices ifNil:[ 
		children do:[ :c | 
			(trans localPointToGlobal: c localPosition) y - floor > c radius ifFalse:[
				rval _ c collideFloor: floor transform: trans. 
				rval ifNotNil:[ collideLoc y > rval y ifTrue:[collideLoc _ rval.].].].].
			collideLoc y <= floor ifTrue:[^ collideLoc.].
			^ nil.
	] ifNotNil: [
		children ifNil:[
			vertices do:[:v | 
				vrt _ trans localPointToGlobal: v. 
				vrt y < floor ifTrue:[ vrt y < collideLoc y ifTrue:[ 
					collideLoc _ vrt asB3DVector4. 
					collideLoc w: vrt y - floor].].].] 		
		ifNotNil:[
			children do:[:f | 
				vrt _ trans localPointToGlobal: (vertices at: (f+1)).
				vrt y < floor ifTrue:[ vrt y < collideLoc y ifTrue:[ 
					collideLoc _ vrt asB3DVector4. 
					collideLoc w: vrt y - floor.].].].].].
	collideLoc y <= floor ifTrue:[ ^ collideLoc.].
	^ nil.
collidePlane: norm offset: offst
 

	| current |
" The plane is defined in the local coordinate frame - no transforms of vertices are necessary. Compare this to #collideFloor:transform: where the plane is the floor, and the test vertices must be transformed into global coordinate space."

	current _ B3DVector4 new.
	current w: Float infinity.

	current _ self collidePlane: norm offset: offst current: current.
	current w = Float infinity ifTrue:[^ nil. ].
	^ current.
collidePlane: norm offset: offst current: current


	| distance rval vrt |
	
" Presumably, we already know that the current sphere intersects this plane - we just want to test its contents."

	rval _ current.
	vertices ifNil:[ 
		children do:[ :c | 
			(c localPosition dot: norm)-offst < c radius ifTrue:[
				rval _ c collidePlane: norm offset: offst current: rval.].].
			^ rval.
	] ifNotNil: [
		children ifNil:[
			vertices do:[:v | 
				distance _ (v dot: norm) - offst. 
				distance < 0.0 ifTrue:[
					distance < rval w ifTrue:[ rval _ v asB3DVector4. rval w: distance.].].
				].] 		ifNotNil:[
			children do:[:f | 
				vrt _ vertices at: (f+1).
				distance _ (vrt dot: norm) - offst.
				distance < 0.0 ifTrue:[
					distance < rval w ifTrue:[ rval _  vrt asB3DVector4. rval w: distance.].].
				].].].
	^ rval.
collideSphere: aSphere

	"Answer whether I collilde into aSphere."
	| myPos itsPos dist pt |
	"Note: For some reason the global position of some spheres is broken. E.g.,
		(frame globalTransform localPointToGlobal: localPosition) distanceTo: globalPosition.
	is NOT zero."
	myPos := self frame globalTransform localPointToGlobal: localPosition.
	itsPos := aSphere frame globalTransform localPointToGlobal: localPosition.
	dist := myPos - itsPos.
	dist := dist dot: dist.
	dist > (self radiusSquared + aSphere radiusSquared) ifTrue:[^nil]. "no intersection"
	"The global spheres intersect. If I have any local spheres, then recurse into them."
	normal ifNil:[ children ifNotNil:[ children size > 0 ifTrue:[
		children do:[:each|
			pt := each collideSphere: aSphere.
			pt ifNotNil:[^pt].
		].
		^nil.
	]]].
	"I am a leaf sphere"
	aSphere isLeafSphere ifTrue:[
		"It's a leaf sphere as well. Answer an approximate collision point (not accurate!)"
		^(self globalPosition + aSphere globalPosition) * 0.5.
	].
	"Otherwise have the other sphere collide"
	^aSphere collideSphere: self.
isLeafSphere

	"Answer whether I am a leaf sphere"
	^normal notNil
ray: framePtr tri: p1 tri: p2 tri: p3



" This method works by using the orthogonal picking matrix as orthogonal planes going through the origin at 'origin'. The points of the triangle are tested to the two planes orthoganal to the direction of picking. If they are on either side of both planes, then there is a very high probability that our picking ray interesects the triangle, and we do the actual calculation of the intersection point. Note the comments below. I use a vertical bar to indicate which side of the plane a point is on.  ***** p1 p2 | ***** means that p1 and p2 are on the same side, 
***** p1 | p2 ***** indicates that they are on opposite sides. p1 is always on the left side, because it is the first one I test. This really just means - same side as p1 - no meaning other than that. 

This can be optimized for tri-strips, tri-fans, and even full face based meshes. When I have the time... "

	| d1 d2 d3 pp1 pp2 dd1 dd2 xPlane yPlane  |
		
	xPlane _ framePtr row1.
	d1 _ xPlane dot: p1.
	d2 _ xPlane dot: p2.
	d3 _ xPlane dot: p3.
	dd1 _ d1*d2.
	dd2 _ d1*d3.

" test if all points are on the same side of the x plane."
	dd1 < 0 ifFalse:[ 				"***** p1 p2 |  *****"
		dd2 < 0  ifFalse:[  ^ nil		"***** p1 p2 p3 | *****" ] 
				ifTrue:[ 			"***** p1 p2 | p3 *****"
						pp1 _ p1 + ((p3-p1) * (d1/(d1-d3))).
						pp2 _ p2 + ((p3-p2) * (d2/(d2-d3))).]]
		ifTrue:[ 				"***** p1 | p2 *****"
			dd2 < 0  ifTrue:[ 		"***** p1 | p2 p3 *****" 
						pp1 _ p1 + ((p2-p1) * (d1/(d1-d2))).
						pp2 _ p1 + ((p3-p1) * (d1/(d1-d3))).] 
					ifFalse:[ 		"***** p1 p3 | p2 *****"
						pp1 _ p1 + ((p2-p1) * (d1/(d1-d2))).
						pp2 _ p3 + ((p2-p3) * (d3/(d3-d2))).].].
" At this point, the triangle is intersected by the x plane and we know where. If both points are on the same side of the y plane, we are done. "
	yPlane _ framePtr row2.

	d1 _ yPlane dot: pp1.
	d2 _ yPlane dot: pp2.
" test if both points are on the same side of the y plane."
	d1 * d2 < 0 ifFalse:[ 
	 ^ nil ].
	pp1 _ pp1 + ((pp2-pp1) * ( d1/ (d1-d2))).
	^ pp1.
	
	


	

construct
calcBoundSphere: verts


	|  iP iN xmx xmn ymx ymn zmx zmn xspan yspan zspan maxSpan d1 d2 center radSq rad |
	
" This routine is from the Graphics Gems 1, C code on page 723. verts is a B3DVector3Array. It constructs a bounding sphere from the vertices."

	iP _ Float infinity.
	iN _ Float infinity negated.

	xmx _ B3DVector3 x: iN y: iN z: iN.
	ymx _ B3DVector3 x: iN y: iN z: iN.
	zmx _ B3DVector3 x: iN y: iN z: iN.
	xmn _ B3DVector3 x: iP y: iP z: iP.
	ymn _ B3DVector3 x: iP y: iP z: iP.
	zmn _ B3DVector3 x: iP y: iP z: iP.

" First pass, find a pair of maximally distant points between each axis."
	verts do:[ :v |
		v x > xmx x ifTrue:[ xmx _ v].
		v x < xmn x ifTrue:[ xmn _ v].
		v y > ymx y ifTrue:[ ymx _ v].
		v y < ymn y ifTrue:[ ymn _ v].
		v z > zmx z ifTrue:[ zmx _ v].
		v z < zmn z ifTrue:[ zmn _ v].].

	xspan _ (xmx-xmn) squaredLength.
	yspan _ (ymx-ymn) squaredLength.
	zspan _ (zmx-zmn) squaredLength.
		
	maxSpan _ xspan.
	d1 _ xmn.
	d2 _ xmx.
	yspan > maxSpan ifTrue:[maxSpan _ yspan. d1 _ ymn. d2 _ ymx.].
	zspan > maxSpan ifTrue:[maxSpan _ zspan. d1 _ zmn. d2 _ zmx.].
	
	center _ (d1 + d2)/2.0.
	radSq _ (d1 - center) squaredLength.
	rad _ radSq sqrt.
	
"Second pass, expand sphere to include outlying points."

	verts do:[ :v |
		d1 _ (v-center)squaredLength.
		d1 > radSq ifTrue:[
			d1 _ d1 sqrt.
			rad _ (d1 + rad)/2.0.
			radSq _ rad*rad.
			d2 _ d1 - rad.
			center _ ((rad * center) + (d2*v))/d1.
			].
		].

	localPosition _ center.
	radius _ rad.
	radiusSquared _ radSq.
"	vertices _ verts. " "don't need this and it screws up other things later."
	^  self.
calcBoundSphere: verts faces: faces


	|  v iP iN xmx xmn ymx ymn zmx zmn xspan yspan zspan maxSpan d1 d2 center radSq rad |
	
" This routine is from the Graphics Gems 1, C code on page 723. verts is a B3DVector3Array. It constructs a bounding sphere based upon the faces and vertices."

	iP _ Float infinity.
	iN _ Float infinity negated.

	xmx _ B3DVector3 x: iN y: iN z: iN.
	ymx _ B3DVector3 x: iN y: iN z: iN.
	zmx _ B3DVector3 x: iN y: iN z: iN.
	xmn _ B3DVector3 x: iP y: iP z: iP.
	ymn _ B3DVector3 x: iP y: iP z: iP.
	zmn _ B3DVector3 x: iP y: iP z: iP.

" First pass, find a pair of maximally distant points between each axis."
	faces do:[ :f |
		v _ verts at:(f+1).
		v x > xmx x ifTrue:[ xmx _ v].
		v x < xmn x ifTrue:[ xmn _ v].
		v y > ymx y ifTrue:[ ymx _ v].
		v y < ymn y ifTrue:[ ymn _ v].
		v z > zmx z ifTrue:[ zmx _ v].
		v z < zmn z ifTrue:[ zmn _ v].].

	xspan _ (xmx-xmn) squaredLength.
	yspan _ (ymx-ymn) squaredLength.
	zspan _ (zmx-zmn) squaredLength.
		
	maxSpan _ xspan.
	d1 _ xmn.
	d2 _ xmx.
	yspan > maxSpan ifTrue:[maxSpan _ yspan. d1 _ ymn. d2 _ ymx.].
	zspan > maxSpan ifTrue:[maxSpan _ zspan. d1 _ zmn. d2 _ zmx.].
	
	center _ (d1 + d2)/2.0.
	radSq _ (d1 - center) squaredLength.
	rad _ radSq sqrt.
	
"Second pass, expand sphere to include outlying points."

	faces do:[ :f |
		v _ verts at: (f+1).
		d1 _ (v-center)squaredLength.
		d1 > radSq ifTrue:[
			d1 _ d1 sqrt.
			rad _ (d1 + rad)/2.0.
			radSq _ rad*rad.
			d2 _ d1 - rad.
			center _ ((rad * center) + (d2*v))/d1.
			].
		].

	localPosition _ center.
	radius _ rad.
	radiusSquared _ radSq.
	^  self.
calcSurface: verts faces: faces


"This is for a first approximation - I may need a better solution later."

	| n nabs m nlist ray t |

	normal _ B3DVector3 new.
	nlist _ OrderedCollection new.

	1 to: faces size by: 3 do:[ :i |
		
		n _ 
		(((verts at:1+(faces at: i)) - (verts at:1+(faces at: i+1))) cross:
		(((verts at:1+(faces at: i)) - (verts at:1+(faces at: i+2))))).
		n squaredLength = 0.0 ifFalse:[
			nlist add: n.
			normal _ normal + n. ].].

	normal _ normal/( faces size/3).
	normal squaredLength = 0.0 ifTrue:[normal _ B3DVector3 x: 0.0 y: 1.0 z: 0.0].
	normal normalize.
	nabs _ normal abs.
	up _ B3DVector3 new.

	nabs x < nabs y ifTrue:[
		nabs x < nabs z ifTrue:[up x:1.0] 
		ifFalse:[up z:1.0].]ifFalse:[
		nabs y < nabs z ifTrue:[up y: 1.0] ifFalse:[up z: 1.0]].
	side _ (up cross: normal) normalized.
	up _ (normal cross: side) normalized.

	m _ B3DMatrix4x4 identity.
	
	m at: 1 at: 1 put: side x.
	m at: 1 at: 2 put: side y.
	m at: 1 at: 3 put: side z.
	m at: 2 at: 1 put: up x.
	m at: 2 at: 2 put: up y.
	m at: 2 at: 3 put: up z.
	m at: 3 at: 1 put: normal x.
	m at: 3 at: 2 put: normal y.
	m at: 3 at: 3 put: normal z.

	offset _ 0.
	1 to: faces size by: 3 do:[ :i |
		ray _ self ray: m tri: localPosition - (verts at: 1+(faces at:i))
					tri: localPosition - (verts at: 1+(faces at: i+1))
					tri: localPosition - (verts at: 1+(faces at: i+2)).
		
		ray ifNotNil:[ 
			nabs x > nabs y ifTrue:[
				nabs x > nabs z ifTrue:[t _ ray x/normal x] 
				ifFalse:[ t _ ray z/normal z].]ifFalse:[
				nabs y > nabs z ifTrue:[t _ ray y/normal y] ifFalse:[t _ ray z/normal z]].


		t abs > offset abs ifTrue:[ offset _ t negated.]]].
calcTree: verts faces: faces box: bx depth: depth


	| child faceBox boxes flist |

" See TBoundSphere >>#calcTree:faces:depth for more info."

	vertices _ nil.
	(depth = 0 or:[ faces size <= 12])ifFalse:[
	"depth = 0 ifFalse:["
		boxes _ self splitBoxes: bx.
" Calculate the list of faces that overlap each box, and recurse."
		boxes do:[:b |
			flist _ OrderedCollection new.
			1 to: faces size by:3 do:[:i |
				(b intersectT1: (verts at:1+(faces at: i)) 
					T2:(verts at: 1+(faces at: i+1)) 
					T3:(verts at: 1+(faces at: i+2))) ifTrue:[
						flist add: (faces at: i).
						flist add: (faces at: i+1).
						flist add: (faces at: i+2).].].
			flist size > 1 ifTrue:[
				faceBox _ TBox new.
				flist do:[:f | faceBox growVertex: (verts at: 1+f).].
				faceBox _ faceBox intersectBox: b.
				child _ (TBoundSphere calcTree: verts faces: flist box: faceBox depth: depth-1).
				child children ifNotNil:[self addChild: child].					
				].].]
	ifTrue:[
		vertices _ verts. 
		children _ faces asIntegerArray. 
		].

" Now calculate the minimal bounding sphere. We do it this way because the polygons actually overlap between boxes. If we were to use the mtfBall approach, the resulting radius could well be larger than this one."
	self box: bx.
	self localPosition: bx center radius: bx diagonal/2.0.
	vertices ifNotNil:[self calcSurface: verts faces: faces.].

				
		

	
calcTree: verts faces: faces depth: depth


	| bx epsilon |

" 
CONSTRUCTION FOR SPHERE HIERARCHY

This method is used to compute a heirarchy of bounding spheres for a complex object. This allows for extremely fast picking and collision detection. The object is broken into octrees if the ratio of sides permits. Long thin objects tend to be split into two parts, flat objects into four.

	1.	This is a top down approach. I use aligned boxes (TBox) as the basis of the
		construction, given the assumption that when the boxes are small enough
		their shape and relative orientation doesn't matter as long as they
		completely cover the object. We will construct spheres from these boxes
		(or cubes) that only contain the contents of the cubes, hence, overlap
		is not much of an issue.
	2.	Once the top level containing box is calculated from the underlying objects
		vertices (TBox>>#max is the maximum value corner and TBox>>#min is the
		minimal value corner), we slice the object into smaller boxes and recurse.
		This slicing is only done along the 'longer' dimensions, hence a long
		thin object may only be sliced once perpendicular to its long axis, 
		a flat object may only be sliced twice each perpendicular to each other.
		We determine whether a dimension is 'sliceable' based upon the ratio of
		its length relative to the other dimensions.
	3.	The recursion first determines which surfaces of the containing box 
		overlap the resulting bounding box. These resulting values are stored as 
		an array for use in the same way by the children of this box.
	4.	We test the depth of the recursion at this stage. If the depth is zero or
		the number of surfaces that we overlap is zero we stop the recursion. If
		this cube does not contain any surfaces its parent will delete it.
	5.	Inside the recursion, we determine which of the set of faces overlap
		this box. Once we know this, we calculate the bounding box of these
		surfaces. We then take the intersection between this box and the 
		previously constructed box from the parent. This can only be the same
		size or smaller, which improves the accuracy of the bounding box 
		representation of the object and further improves its value in 
		collision detection.
	6.	Once we construct the child cubes, we calculate the minimal bounding boxes
		of these in attempt to make the current cube size even smaller.
	7.	The last thing to do here is if this cube covers four faces or less, we
		shouldn't bother testing the children, and go directly to these surfaces.
		We set this cube to be a leaf node and destroy all of the children.
	8.	And now, the things we have been waiting for. We construct the sphere
		tree from the cube tree recursively (of course). The boolean leaf value indicates
		whether or not this is a leaf sphere, hence its children are faces.


We use the TBox class to help us out here. The TBox class defines a min and max corner and has a number of useful methods to simplify the code here."

	bx _ TBox new.
	faces do:[:f| bx growVertex: (verts at: 1+f).].
" Grow by epsilon to make sure the box has some depth - otherwise, it will fail to construct good boxes."
	epsilon _ 0.000001.
	bx growVertex: bx min - (B3DVector3 x: epsilon y: epsilon z: epsilon).
	bx growVertex: bx max + (B3DVector3 x: epsilon y: epsilon z: epsilon).
	self calcTree: verts faces: faces box: bx depth: depth.				
	globalPosition _ B3DVector3 new.

	
frame: frm


	frame _ frm.
	children ifNotNil:[ vertices ifNil:[ 
		children do:[:c | c frame: frm.].].].
splitBoxes: bx


	| epsilon ratio boxes splitBoxes extent |
" Split along the axis where the ratio with the other axes is greater than 0.6. I am sure there is a more efficient number. Just don't know what it is."

"	ratio _ 0.737."
	ratio _ 0.6.
	extent _ bx extent.

" Make sure we don't divide by zero. "
	epsilon _ 0.000001.
	extent _ extent max: (B3DVector3 x: epsilon y: epsilon z: epsilon).
	boxes _ OrderedCollection new.
	boxes add: bx.

	(((extent x/extent y) > ratio) and:[(extent x/extent z) > ratio]) ifTrue:[ 
		splitBoxes _ OrderedCollection new.
		boxes do:[:b | 
			splitBoxes add: b splitXMin.
			splitBoxes add: b splitXMax.
			].
		boxes _ splitBoxes.
		].

	(((extent y/extent x) > ratio) and:[(extent y/extent z) > ratio]) ifTrue:[ 
		splitBoxes _ OrderedCollection new.
		boxes do:[:b | 
			splitBoxes add: b splitYMin.
			splitBoxes add: b splitYMax.
			].
		boxes _ splitBoxes.
		].

	(((extent z/extent x) > ratio) and:[(extent z/extent y) > ratio]) ifTrue:[ 
		splitBoxes _ OrderedCollection new.
		boxes do:[:b | 
			splitBoxes add: b splitZMin.
			splitBoxes add: b splitZMax.
			].
		boxes _ splitBoxes.
		].

	^ boxes.
union: sphere


" We calculate the minimal sphere that contains both this sphere and the argument sphere."

" This needs to be performed in global coordinates."

	| vector length norm p1 p2 |
	sphere ifNil:[ ^ TBoundSphere localPosition: (self globalPosition) radius: (self radius).].

	vector _ sphere globalPosition - self globalPosition.
	vector length = 0.0 ifTrue:[
		^ TBoundSphere localPosition: (self globalPosition) radius: (self radius max: sphere radius).
		].
	length _ vector length.
	norm _ vector/length.
	length + sphere radius > self radius ifTrue:[ p2 _ norm * (length + sphere radius).] 
		ifFalse:[ ^ TBoundSphere localPosition: (self globalPosition) radius: (self radius).].
	length + self radius > sphere radius ifTrue:[ p1 _ norm * (self radius) negated ] 
		ifFalse:[ ^ TBoundSphere localPosition: (sphere globalPosition) radius: (sphere radius).].

	length _ ((p1-p2) length)/2.0.
	p1 _ ((p1+p2)/2.0) + self globalPosition.
	^ TBoundSphere localPosition: p1 radius: length.

	

mini ball
mtfBall: aCollection

	"Create a minimum enclosing sphere from aCollection of points.
	See Emo Welzl, 'Smallest enclosing disks (balls and ellipsoids)' (1991),
	H. Maurer (Ed.), New Results and New Trends in Computer Science, LNCS 555
	http://citeseer.nj.nec.com/welzl91smallest.html"
	| border idxList |
	idxList := WordArray new: aCollection size.
	1 to: idxList size do:[:i| idxList at: i put: i].
	self privateMtfRandomize: idxList.
	border := Array new: 4.
	self privateMtfBall0: aCollection with: idxList with: 1 with: border.
	radius := radiusSquared sqrt.

	box _ TBox new.
	aCollection do:[ :v | box growVertex:(v).].
mtfBall: aCollection faces: faceCollection

	"Create a minimum enclosing sphere from aCollection of points.
	See Emo Welzl, 'Smallest enclosing disks (balls and ellipsoids)' (1991),
	H. Maurer (Ed.), New Results and New Trends in Computer Science, LNCS 555
	http://citeseer.nj.nec.com/welzl91smallest.html"
	| border idxList |
	idxList := WordArray new: faceCollection size.
	1 to: idxList size do:[:i| idxList at: i put: (faceCollection at: i) + 1].
	self privateMtfRandomize: idxList.
	border := Array new: 4.
	self privateMtfBall0: aCollection with: idxList with: 1 with: border.
	radius := radiusSquared sqrt.
privateMtfBall0: pointList with: indexList with: minIndex with: borderList

	"Note: In this variant the default recursion path is inlined. The recursion exclusively affects changes of the border vertices so its depth is bounded by the max. number of vertices on the border (which is four)."
	| idx p dist2 index |
	localPosition := B3DVector3 new.
	radiusSquared := -1.0.
	index := indexList size.
	[index >= minIndex] whileTrue:[
		idx := indexList at: index.
		p := pointList at: idx.
		p -= localPosition.
		dist2 := (p dot: p) * 0.999. "numerical accuracy"
		dist2 <= radiusSquared ifFalse:[
			p := pointList at: idx.
			borderList at: 1 put: p.
			self privateMtfBall1: pointList with: indexList with: index+1 with: borderList.
			"Move p to front"
			indexList replaceFrom: index to: indexList size-1 with: indexList startingAt: index+1.
			indexList at: indexList size put: idx.
		].
		index := index - 1.
	].
privateMtfBall1: pointList with: indexList with: minIndex with: borderList

	"Note: In this variant the default recursion path is inlined. The recursion exclusively affects changes of the border vertices so its depth is bounded by the max. number of vertices on the border (which is four)."
	| idx p dist2 index |
	localPosition := (borderList at: 1).
	radiusSquared := 0.0.
	index := indexList size.
	[index >= minIndex] whileTrue:[
		idx := indexList at: index.
		p := pointList at: idx.
		p -= localPosition.
		dist2 := (p dot: p) * 0.999. "numerical accuracy"
		dist2 <= radiusSquared ifFalse:[
			p := pointList at: idx.
			borderList at: 2 put: p.
			self privateMtfBall2: pointList with: indexList with: index+1 with: borderList.
			"Move p to front"
			indexList replaceFrom: index to: indexList size-1 with: indexList startingAt: index+1.
			indexList at: indexList size put: idx.
		].
		index := index - 1.
	].
privateMtfBall2: pointList with: indexList with: minIndex with: borderList

	"Note: In this variant the default recursion path is inlined. The recursion exclusively affects changes of the border vertices so its depth is bounded by the max. number of vertices on the border (which is four)."
	| idx p dist2 index |
	self privateMtfSphereFromP1: (borderList at: 1) p2: (borderList at: 2).
	index := indexList size.
	[index >= minIndex] whileTrue:[
		idx := indexList at: index.
		p := pointList at: idx.
		p -= localPosition.
		dist2 := (p dot: p) * 0.999. "numerical accuracy"
		dist2 <= radiusSquared ifFalse:[
			p := pointList at: idx.
			borderList at: 3 put: p.
			self privateMtfBall3: pointList with: indexList with: index+1 with: borderList.
			"Move p to front"
			indexList replaceFrom: index to: indexList size-1 with: indexList startingAt: index+1.
			indexList at: indexList size put: idx.
		].
		index := index - 1.
	].
privateMtfBall3: pointList with: indexList with: minIndex with: borderList

	"Note: In this variant the default recursion path is inlined. The recursion exclusively affects changes of the border vertices so its depth is bounded by the max. number of vertices on the border (which is four)."
	| idx p dist2 index |
	self privateMtfSphereFromP1: (borderList at: 1) p2: (borderList at: 2) p3: (borderList at: 3).
	index := indexList size.
	[index >= minIndex] whileTrue:[
		idx := indexList at: index.
		p := pointList at: idx.
		p -= localPosition.
		dist2 := (p dot: p) * 0.999. "numerical accuracy"
		dist2 <= radiusSquared ifFalse:[
			p := pointList at: idx.
			borderList at: 4 put: p.
			self privateMtfSphereFromP1: (borderList at: 1) p2: (borderList at: 2) 
								p3: (borderList at: 3) p4: (borderList at: 4).
			"Move p to front"
			indexList replaceFrom: index to: indexList size-1 with: indexList startingAt: index+1.
			indexList at: indexList size put: idx.
		].
		index := index - 1.
	].
privateMtfRandomize: idxList

	"Randomize input"
	| max rnd j tmp |
	max := idxList size.
	rnd := Random new" seed: 16r2134567".
	1 to: max do:[:i|
		j := (rnd next * i) truncated +1.
		tmp := idxList at: i.
		idxList at: i put: (idxList at: j).
		idxList at: j put: tmp.
	].
privateMtfSphereFromP1: p1 p2: p2

	"Two points"
	| delta center rad2 |
	center := (p1 + p2).
	center *= 0.5.
	delta := center - p1.
	rad2 := delta dot: delta.
	localPosition := center.
	radiusSquared := rad2.
privateMtfSphereFromP1: p1 p2: p2 p3: p3

	| d1 d2 d3 v1 v2 v3 m b rad2 center |
	v1 := (p1 + p2). v1 *= 0.5.
	v2 := (p1 + p3). v2 *= 0.5.
	v3 := p1.

	d1 := (p2 - p1).
	d2 := (p3 - p1).
	d3 := d1 cross: d2.
	m := B3DMatrix4x4 new.
	m
		a11: d1 x; a12: d1 y; a13: d1 z;
		a21: d2 x; a22: d2 y; a23: d2 z;
		a31: d3 x; a32: d3 y; a33: d3 z.
	b := B3DVector3
			x: (v1 dot: d1)
			y: (v2 dot: d2)
			z: (v3 dot: d3).
	center := m solve3x3: b.
	center ifNil:[^self privateMtfSphereFromP1: p1 p2: p2].
	d1 := center - p1.
	rad2 := d1 dot: d1.
	localPosition := center.
	radiusSquared := rad2.
privateMtfSphereFromP1: p1 p2: p2 p3: p3 p4: p4

	| d1 d2 d3 v1 v2 v3 m b rad2 center |
	v1 := (p1 + p2). v1 *= 0.5.
	v2 := (p1 + p3). v2 *= 0.5.
	v3 := (p1 + p4). v3 *= 0.5.

	d1 := (p2 - p1).
	d2 := (p3 - p1).
	d3 := (p4 - p1).
	m := B3DMatrix4x4 new.
	m
		a11: d1 x; a12: d1 y; a13: d1 z;
		a21: d2 x; a22: d2 y; a23: d2 z;
		a31: d3 x; a32: d3 y; a33: d3 z.
	b := B3DVector3
			x: (v1 dot: d1)
			y: (v2 dot: d2)
			z: (v3 dot: d3).
	center := m solve3x3: b.
	center ifNil:[^self privateMtfSphereFromP1: p1 p2: p2 p3: p3].
	d1 := center - p1.
	rad2 := d1 dot: d1.
	localPosition := center.
	radiusSquared := rad2.

render
pick: pointer


	| rval |

" There is a problem here. I do not have the norm for these faces, which means that this will
 allow us to pick a triangle on the wrong side of the object. This must be fixed.

Another thing that needs to be looked at is just because we find a triangle to pick does NOT mean that there isn't a closer one. This code exits on the very first triangle - but I am leaving it -as is- for performance reasons. To fix it, look at the next 'clause', where I had the same problem with groups of polys."

	(pointer pickLocalBoundSphere: self) ifFalse:[^false].
	vertices ifNotNil:[^pointer pickTriangles: vertices list: children].
	rval _ false.
	children do:[:c | rval _ (c pick: pointer) or:[ rval]. ].
	^rval.		
pickChildren: pointer


	| rval |

" This method exists to avoid testing the top-most sphere twice. It is just a special case of the #pick: method. Notice the commented out code."

"	(pointer pickLocalBoundSphere: self)ifTrue:["
	vertices ifNotNil:[^pointer pickTriangles: vertices list: children].
	rval _ false.
	children do:[:c | rval _ (c pick: pointer) or:[ rval]. ].
	^ rval
render: ogl box: frm


	children ifNotNil:[
		vertices ifNil: [ children do:[ :c | c render: ogl box: frm].].].
	self box ifNotNil:[frm render: ogl box: self box].
render: ogl box: frm depth: depth


	children ifNotNil:[
		vertices ifNil: [ children do:[ :c | c render: ogl box: frm depth: depth-1. ].].].
	depth = 0 ifTrue:[ box ifNotNil:[
		frm render: ogl box: self box.].].
render: ogl frame: frm segments: segs


	children ifNotNil:[
		vertices ifNil: [ children do:[ :c | c render: ogl frame: frm segments: segs. ].].].
	frm render: ogl sphere: self segments: segs.
renderBoundSphere: frm

	| mid ogl |
	children ifNotNil:[ normal ifNil: [ 
		children do:[ :c | c renderBoundSphere: frm. ].
		^self.
	].].
	ogl _ Croquet world ogl.
	normal ifNotNil:[
		ogl glColor3fv: #(1.0 0 0)asFloatArray.
		ogl glPolygonMode: GLFrontAndBack with: GLLine.
		1 to: 2 do:[:i|
		ogl glBegin: GLTriangleFan.
		mid := i == 1 ifTrue:[normal] ifFalse:[normal negated].
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := side negated.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := (side + up) negated normalized.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := up negated.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := (side - up) normalized.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := side.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := (side+up) normalized.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := up.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		mid := (up - side) normalized.
			ogl glNormal3fv: mid.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).
		ogl glEnd.
		].
		ogl glPolygonMode: GLFrontAndBack with: GLFill.
	].
renderSurface: ogl


	| mid |
	children ifNotNil:[
		normal ifNil: [ children do:[ :c | c renderSurface: ogl. ].].].
	normal ifNotNil:[
	 		ogl glColor3fv: #(1.0 0 0)asFloatArray.
			ogl glBegin: GLLineStrip.
    
			ogl glVertex3fv:localPosition + (normal*offset) + normal.
			ogl glVertex3fv:localPosition + (normal * offset).
			ogl glVertex3fv:localPosition + (normal*offset) + (up*radius).
		mid _ (up+side) normalized.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).		
			ogl glVertex3fv:localPosition + (normal*offset) + (side*radius).
		mid _ (side - up) normalized.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).		
			ogl glVertex3fv:localPosition + (normal*offset) + (up*radius negated).
		mid _( side+up)normalized negated.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).		
			ogl glVertex3fv:localPosition + (normal*offset) + (side*radius negated).
		mid _ (up-side)normalized.
			ogl glVertex3fv:localPosition + (normal*offset) + (mid*radius).		
			ogl glVertex3fv:localPosition + (normal*offset) + (up*radius).
		ogl glEnd.].

transform
transform: trans


	trans ifNotNil:[
		globalPosition _ trans localPointToGlobal: localPosition.] 
	ifNil: [globalPosition _ localPosition.].

" As far as I can tell, we never need to transform the children boundSpheres.

	children ifNotNil:[
		leaf ifFalse:[
		children do:[:cs | cs transform: trans.].].]."

class methods:

instance creation
calcBoundSphere: verts


	 ^ TBoundSphere new calcBoundSphere: verts.
calcBoundSphere: verts faces: faces


	 ^ TBoundSphere new calcBoundSphere: verts faces: faces.
calcTree: vertices faces: faces box: bx depth: depth


	^ TBoundSphere new calcTree: vertices faces: faces box: bx depth: depth.
calcTree: vertices faces: faces depth: depth


	^ TBoundSphere new calcTree: vertices faces: faces depth: depth.
localPosition: pos radius: rad


	^ TBoundSphere new localPosition: pos radius: rad.
mtfBall: verts


	 ^ self new mtfBall: verts.
mtfBall: verts faces: faces


	 ^ self new mtfBall: verts faces: faces.

^top


- made by Dandelion -