Spatial partitioning structures are often used to accelerate rendering systems by determining what is visible, and what is not currently visible.
I'm sure you've heard of Binary Space Partitioning (BSP) trees, Quadtrees, Octrees, and perhaps you've heard of others, there are many.
In the same way, we can take advantage of spatial partitioning to accelerate the broad phase of our collision detection, where we need to test every entity against every other entity for possible collisions. Spatial partitioning can be used to track smaller groups of proximate objects, which can speed up broadphase testing by several orders, depending on the number of objects being simulated, and the way they are clustered in space.
Today we will be discussing a kind of Binary Tree known generically as a Bounding Volume Hierarchy (BVH), and the specific kind we'll examine is called a Sphere Tree. I will use the term "binary space partitioning", however I am not referring to a BSP Tree.
As you're probably already aware, binary space partitioning schemes work by dividing space into two subspaces, and then recursively subdividing each subspace into smaller subspaces until some HEURISTIC is satisfied. A spheretree works the same way. So how does one go about constructing a spheretree? And what is a heuristic?
Heuristics are a qualitative or quantitative measure which guides an algorithm toward an optimal output.
Which heuristics we choose depends completely on what the phrase "optimal output" means to us - it is purely contextual, it depends what we're doing, and what we're trying to achieve. If we can define what we want our algorithm to produce, then we can determine one or more heuristics which will help us to get there. The most common synonym for the word "heuristic" is known as an "ERROR METRIC"... we measure the amount of error in our calculations, and use that error term to steer subsequent calculations toward perfection. Heuristics are what our algorithm evaluates when it has to make DECISIONS from multiple possible options.
There are basically three methods of constructing a sphere tree... TOP-DOWN, BOTTOM-UP, and INSERTIVE.
I will briefly discuss each of these, but it's worth noting from the outset that only the last one is valuable for a dynamic system, where the tree needs to be rebuilt in realtime.
I will break now, and continue this post a little later today.
I'm sure you've heard of Binary Space Partitioning (BSP) trees, Quadtrees, Octrees, and perhaps you've heard of others, there are many.
In the same way, we can take advantage of spatial partitioning to accelerate the broad phase of our collision detection, where we need to test every entity against every other entity for possible collisions. Spatial partitioning can be used to track smaller groups of proximate objects, which can speed up broadphase testing by several orders, depending on the number of objects being simulated, and the way they are clustered in space.
Today we will be discussing a kind of Binary Tree known generically as a Bounding Volume Hierarchy (BVH), and the specific kind we'll examine is called a Sphere Tree. I will use the term "binary space partitioning", however I am not referring to a BSP Tree.
As you're probably already aware, binary space partitioning schemes work by dividing space into two subspaces, and then recursively subdividing each subspace into smaller subspaces until some HEURISTIC is satisfied. A spheretree works the same way. So how does one go about constructing a spheretree? And what is a heuristic?
Heuristics are a qualitative or quantitative measure which guides an algorithm toward an optimal output.
Which heuristics we choose depends completely on what the phrase "optimal output" means to us - it is purely contextual, it depends what we're doing, and what we're trying to achieve. If we can define what we want our algorithm to produce, then we can determine one or more heuristics which will help us to get there. The most common synonym for the word "heuristic" is known as an "ERROR METRIC"... we measure the amount of error in our calculations, and use that error term to steer subsequent calculations toward perfection. Heuristics are what our algorithm evaluates when it has to make DECISIONS from multiple possible options.
There are basically three methods of constructing a sphere tree... TOP-DOWN, BOTTOM-UP, and INSERTIVE.
I will briefly discuss each of these, but it's worth noting from the outset that only the last one is valuable for a dynamic system, where the tree needs to be rebuilt in realtime.
I will break now, and continue this post a little later today.
Heuristics for tree construction are often "MUTUALLY EXCLUSIVE" - for example, to construct a BSP tree we need to find a Plane which divides the input entities into two subsets that are as equal in number as possible, while cutting through as few entities as possible. We will often need to 'balance' our heuristics to find an outcome which best fits ALL of the heuristics.
Any BVH construction algorithm should observe the following three heuristics:
1. The actual volume enclosed by the bounding volumes should be as small as possible, at each level of the tree.
Note that this is true when a parent node groups together two bounding volumes that are close together.
2. The child bounding volumes of any parent should overlap as little as possible.
Often this will be mutually exclusive to the first heuristic (above), and in general it is better to favor smaller volumes over minimal overlaps. In fact, if you choose a minimal overlap at some low level of the tree, it is likely to cause greater overlaps higher up the tree: so a tree with an overall low overlap is likely to fulfill both requirements.
3. The tree should be as balanced as possible.
You should avoid having some branches of the tree that are very deep while others are very shallow. The worst-case scenario is a tree with one long branch. In this case the advantage of having a hierarchy is minimal. The biggest speed-up is gained when all branches are roughly at the same length.
Let's look at top-down construction first, because it is the most intuitive.
We have a bunch of entities in 3D space, and we want to construct a spheretree according to our three heuristics.
For a top-down approach, we can begin by creating a HUGE bounding sphere around all of the entities in our world.
Now we need to split the entities into two subsets which have the smallest possible boundingspheres while simultaneously holding as equal a number of entities as possible... mutually exclusive goals.
Once we've done that, we repeat this procedure for each subset, until our subsets each only contain a single entity - since we can't continue, we will call this a 'leaf node' - notice that the leaf nodes of our tree contain one entity each, and all the other nodes contain NO entities.
The bottom-up approach is quite different... let's assume once more that we start with a bunch of entities.
We'll begin again by finding the boundingsphere of the entire set, and we'll determine where the origin of this 'world bounding sphere' is located.
Now we'll figure out which entity's origin is closest to this big boundingsphere's origin, and that entity will be the starting point for our algorithm.. we can discard the big boundingsphere, we've chosen the entity which is closest to the midpoint of the distribution.
Now we're looking for the next entity, and our heuristic will be that we want to find an entity which would cause the smallest possible increase in boundingsphere in order to include it. For each 'other entity', we'll determine the midpoint between their boundingspheres, calculate a new boundingsphere using that midpoint as the new origin, and find the sphere which includes these two entities and is the smallest.
This procedure is repeated until the entire input set of entities has been consumed.
The insertive approach, which we will be using, assumes that a tree already exists, although it can be used to generate a tree from scratch.
In the next post, we'll cover the algorithm for maintaining a dynamic spheretree.
Any BVH construction algorithm should observe the following three heuristics:
1. The actual volume enclosed by the bounding volumes should be as small as possible, at each level of the tree.
Note that this is true when a parent node groups together two bounding volumes that are close together.
2. The child bounding volumes of any parent should overlap as little as possible.
Often this will be mutually exclusive to the first heuristic (above), and in general it is better to favor smaller volumes over minimal overlaps. In fact, if you choose a minimal overlap at some low level of the tree, it is likely to cause greater overlaps higher up the tree: so a tree with an overall low overlap is likely to fulfill both requirements.
3. The tree should be as balanced as possible.
You should avoid having some branches of the tree that are very deep while others are very shallow. The worst-case scenario is a tree with one long branch. In this case the advantage of having a hierarchy is minimal. The biggest speed-up is gained when all branches are roughly at the same length.
Let's look at top-down construction first, because it is the most intuitive.
We have a bunch of entities in 3D space, and we want to construct a spheretree according to our three heuristics.
For a top-down approach, we can begin by creating a HUGE bounding sphere around all of the entities in our world.
Now we need to split the entities into two subsets which have the smallest possible boundingspheres while simultaneously holding as equal a number of entities as possible... mutually exclusive goals.
Once we've done that, we repeat this procedure for each subset, until our subsets each only contain a single entity - since we can't continue, we will call this a 'leaf node' - notice that the leaf nodes of our tree contain one entity each, and all the other nodes contain NO entities.
The bottom-up approach is quite different... let's assume once more that we start with a bunch of entities.
We'll begin again by finding the boundingsphere of the entire set, and we'll determine where the origin of this 'world bounding sphere' is located.
Now we'll figure out which entity's origin is closest to this big boundingsphere's origin, and that entity will be the starting point for our algorithm.. we can discard the big boundingsphere, we've chosen the entity which is closest to the midpoint of the distribution.
Now we're looking for the next entity, and our heuristic will be that we want to find an entity which would cause the smallest possible increase in boundingsphere in order to include it. For each 'other entity', we'll determine the midpoint between their boundingspheres, calculate a new boundingsphere using that midpoint as the new origin, and find the sphere which includes these two entities and is the smallest.
This procedure is repeated until the entire input set of entities has been consumed.
The insertive approach, which we will be using, assumes that a tree already exists, although it can be used to generate a tree from scratch.
In the next post, we'll cover the algorithm for maintaining a dynamic spheretree.
The insertion approach is the only one suitable for use in realtime simulations such as 3D games... it can adjust the hierarchy without having to rebuild it completely. Typically, we already have a tree, it can be an empty tree if we are just beginning.
In order to add a new entity to the tree, we need to recursively descend down through the tree: at each node we follow the child (theres only two childs in a binary tree) that would best accommodate the new object. Eventually we'll land in a Leaf, which is then replaced by a new parent sphere which contains both the existing leaf and the new entity.
In order to remove an entity from the tree, we need to replace its parent node with its sibling, then recalculate the boundingspheres in all nodes farther up the hierarchy, until we calculate that there was no change in bounding radius. Since this can (but does not normally) mean recalculating all the bounding spheres all the way back to and including the root node, it is potentially a lot more expensive to delete an entity than it is to add a new one.
Some of the papers you will encounter will go further: instead of the tree stopping with whole entities in the leaf nodes, we could continue to describe spheres which contain FEATURES of the entity (its triangles, its faces, its.. whatever). There is an enormous amount of inter-sphere overlap at this level and so I declare that there is little if any benefit in pursuing this.
Now we have covered spatial partitioning and its relevance to broadphase collision testing... what's next?
In order to add a new entity to the tree, we need to recursively descend down through the tree: at each node we follow the child (theres only two childs in a binary tree) that would best accommodate the new object. Eventually we'll land in a Leaf, which is then replaced by a new parent sphere which contains both the existing leaf and the new entity.
In order to remove an entity from the tree, we need to replace its parent node with its sibling, then recalculate the boundingspheres in all nodes farther up the hierarchy, until we calculate that there was no change in bounding radius. Since this can (but does not normally) mean recalculating all the bounding spheres all the way back to and including the root node, it is potentially a lot more expensive to delete an entity than it is to add a new one.
Some of the papers you will encounter will go further: instead of the tree stopping with whole entities in the leaf nodes, we could continue to describe spheres which contain FEATURES of the entity (its triangles, its faces, its.. whatever). There is an enormous amount of inter-sphere overlap at this level and so I declare that there is little if any benefit in pursuing this.
Now we have covered spatial partitioning and its relevance to broadphase collision testing... what's next?
I don't think I really pointed out the requirements of maintaining a dynamic tree ...
The entities in our tree are represented by the Leaf Nodes.
We've already said that it's ok for SIBLING spheres to overlap.
But if an entity moves in 3D space, eventually it will overlap with its immediate parent sphere, and then with higher parent spheres.
So whenever an entity has moved, we need to check if its sphere is still completely bounded by its immediate parent.
If it isn't, we'll need to delete the entity from the tree and re-insert it, causing at least part of the tree to be rebuilt.
This won't happen every single frame, but it will happen a lot.
Now, let's talk about optimizing the spheretree algorithm(s) which have been outlined so far.
It is important to recognize what we're trying to achieve with this scheme - the spheretree is ONLY a tool for identifying spheres which are proximate in space (close together, and so potentially able to collide). Therefore we're not interested in EXACT bounds, and in fact, our 'hungry sphere tree' is quite a lot of overkill.
Once an entity has been added to the tree, it is acceptable that we consider only the position of its ORIGIN in subsequent checks - we can effectively ignore the radius of the sphere immediately surrounding each entity, and treat them as particles.
Since we've accepted that overlap of sibling spheres is not relevant, then it follows that overlap with the IMMEDIATE PARENT sphere is also not relevant - this can lead to a substantial speedup of the functions used to maintain the tree, without having a severe impact on either the structure of the tree or its value to us as an acceleration tool.
In fact, this same logic CAN be applied right from the very outset - that instead of being interested in spheres which totally enclose our entities, we are only interested in finding spheres which enclose their origins.
I have seen this applied to create acceleration trees for 'point clouds' and particle systems, but never for rapid sphere trees.
However I have implemented this modification in the past, and I did not encounter any problems.
The reason why this can work is best described by using the example of Visibilty Tests based on intersection with the VIEW FRUSTUM - we know that if an entity is PARTLY VISIBLE, we're interested in it - it doesn't need to be TOTALLY VISIBLE, right?
The same logic applies here.
The entities in our tree are represented by the Leaf Nodes.
We've already said that it's ok for SIBLING spheres to overlap.
But if an entity moves in 3D space, eventually it will overlap with its immediate parent sphere, and then with higher parent spheres.
So whenever an entity has moved, we need to check if its sphere is still completely bounded by its immediate parent.
If it isn't, we'll need to delete the entity from the tree and re-insert it, causing at least part of the tree to be rebuilt.
This won't happen every single frame, but it will happen a lot.
Now, let's talk about optimizing the spheretree algorithm(s) which have been outlined so far.
It is important to recognize what we're trying to achieve with this scheme - the spheretree is ONLY a tool for identifying spheres which are proximate in space (close together, and so potentially able to collide). Therefore we're not interested in EXACT bounds, and in fact, our 'hungry sphere tree' is quite a lot of overkill.
Once an entity has been added to the tree, it is acceptable that we consider only the position of its ORIGIN in subsequent checks - we can effectively ignore the radius of the sphere immediately surrounding each entity, and treat them as particles.
Since we've accepted that overlap of sibling spheres is not relevant, then it follows that overlap with the IMMEDIATE PARENT sphere is also not relevant - this can lead to a substantial speedup of the functions used to maintain the tree, without having a severe impact on either the structure of the tree or its value to us as an acceleration tool.
In fact, this same logic CAN be applied right from the very outset - that instead of being interested in spheres which totally enclose our entities, we are only interested in finding spheres which enclose their origins.
I have seen this applied to create acceleration trees for 'point clouds' and particle systems, but never for rapid sphere trees.
However I have implemented this modification in the past, and I did not encounter any problems.
The reason why this can work is best described by using the example of Visibilty Tests based on intersection with the VIEW FRUSTUM - we know that if an entity is PARTLY VISIBLE, we're interested in it - it doesn't need to be TOTALLY VISIBLE, right?
The same logic applies here.
The main problem we have to solve in order to construct our spheretree is called the "Non Weighted Euclidian Single Point Problem", and can be described as follows:
Given two Spheres (or Circles) of known position and radius, find the smallest Sphere (or Circle) which encloses them both.
I've attached a sketch which shows three situations: spheres separated by distance, spheres intersecting, and spheres kissing.
According to my sketch, it appears that we can easily and quickly calculate the Radius of the new boundingsphere: Find the line between the origins of the two input spheres, add both the radii to its length, and we have the diameter of the smallest possible bounding sphere. Its position is a bit more difficult, since we need to extend the line at each end by the appropriate radius - but we've found a fast way to determine the minimal bounding radius (half the extended line's length), and that's all we need MOST of the time.
In order to determine which radius to apply as an extension of which end of our line, we will assume the origin of one input sphere is in fact the origin of our line.
Now we will extend the line in the direction of the other sphere's origin, by the other sphere's radius, and then we know we just have to extend it in the other direction by THIS sphere's radius. This is a naiive solution which works in any direction by taking advantage of our knowledge of the geometries at hand.
Having calculated the full extended line, we now find its midpoint, and that is the position of the origin of our boundingsphere - we can do this by scaling the already-calculated radius in the appropriate direction from either endpoint in order to avoid an unnecessary divide operation.
Armed with this knowledge, we're ready to write our first implementation of the insertion method for a "hungry" spheretree.. let's do that, and then show the much-faster 'particle' version, which I will call "SpeedTree" (no relation to the rendering middleware product).
Today we're going to look at some actual sourcecode.
To begin with, we have this structure:
Every node in our SphereTree will have its own BoundingSphere.
Now let's look at the class which describes a Node in our Tree... please just ignore whatever is not relevant today:
For those of you have trouble understanding that, each node will basically be described by a structure which is roughly equivalent to this:
And finally, we have the code:
It's not perfectly optimized, but it's not ABSOLUTELY horrible, and it will get the job done.
I haven't bothered providing sourcecode for all of the macros that have been used, that is an exercise for OA32 users if they choose to look them up. And please note this code has *NOT* been tested, so there MAY be mistakes - however rest assured that I am indeed implementing this, as this thread documents my ACTUAL coding adventures, not just what's rattling around in my head.
Just looking over that sourcecode, we can see that most of the cost involves calculation of line lengths.
In the next post, I'll show you a slightly different version, which I call "SpeedTree".
The difference is that we will no longer consider the radius of a sphere being inserted into the tree - we'll treat the input sphere as a single particle or point in 3D space. We'll still need to consider the spheres at each node of our tree, but by removing one sphere from the equation, we can reduce the math requirements by approximately half, and so effectively yield a speedup of 200 percent (or thereabouts). In practice it's of course less, but the speedup is nonetheless considerable, and it won't affect the tree's ability to group nearby entities, and so accelerate our broadphase collision testing.
Any speedups here will be incredibly important to our simulation, because we need to do broadphase testing every frame, and because our dynamic tree will be constantly rebuilt as our entities move in worldspace... narrowphase collision testing is VERY expensive, and the whole point of broadphase testing is to try to avoid narrowphase testing !!! So although it may seem like we're optimizing the wrong end of the collision spectrum, if you consider that it is in fact a culling operation, it makes complete sense.
To begin with, we have this structure:
——————————————————————————————————————————————————————————————————————————————————————————————————
; This structure represents a SPHERICAL BOUNDING VOLUME.
; Every Node in our Tree has one of these.
; ——————————————————————————————————————————————————————————————————————————————————————————————————
BoundingSphere struct
vPosition Vec3 <> ;Position of Sphere's Origin, in WorldSpace
fRadius real8 ? ;Sphere's radius
BoundingSphere ends
Every node in our SphereTree will have its own BoundingSphere.
Now let's look at the class which describes a Node in our Tree... please just ignore whatever is not relevant today:
Object BVHNode, 8767843, Primer
; Holds the child nodes of this node.
DefineVariable pChildA, Pointer, NULL
DefineVariable pChildB, Pointer, NULL
; Holds a single bounding volume encompassing all the descendents of this node.
DefineVariable volume, BoundingSphere, {<>}
; Holds the rigid body at this node of the hierarchy.
; Only leaf nodes contain a rigid body.
; Note that it is possible to rewrite the algorithms in this
; class to handle objects at all levels of the hierarchy,
; but the code provided ignores this vector unless ChildA = NULL
DefineVariable body, Pointer , NULL ;-> RigidBody ... only valid in LEAF NODES
RedefineMethod Init, Pointer, ptr RigidBody, ptr BoundingSphere
StaticMethod Insert, ptr RigidBody, ptr BoundingSphere
; Checks the potential contacts from this node downward in
; the hierarchy, writing them to the given array (up to the
; given limit). Returns the number of potential contacts found.
StaticMethod getPotentialContacts, ptr PotentialContact, dword
StaticMethod getPotentialContactsWith, ptr BVHNode,ptr PotentialContact, dword
StaticMethod getBounds, ptr BoundingSphere
StaticMethod isOverlapping, ptr BoundingSphere
StaticMethod makeVolume, ptr BoundingSphere, real8
ObjectEnd
For those of you have trouble understanding that, each node will basically be described by a structure which is roughly equivalent to this:
BVHNode struct
pOwner dd ? ;Pointer to parent node, or NULL=root node
pChildA dd ? ;Pointer to 'left branch' of tree
pChildB dd ? ;Pointer to 'right branch' of tree
volume BoundingSphere <> ;Nested struct describing boundingsphere
pRigidBody dd ? ;Pointer to physics entity, only valid in leaf nodes
BVHNode ends
And finally, we have the code:
;Method: BVHNode.Insert
;Purpose: Insert a new RigidBody into the tree
Method BVHNode.Insert, uses esi, newBody:ptr RigidBody, newVolume:ptr BoundingSphere
LOCAL f1:real8, f2:real8
SetObject esi
; If we are a leaf, then the only option is to spawn two new children and place the new body in one of them (and the current node gets copied into the other, and then this node's sphere is recalculated)
.if .body!=0
; Set Child one = a copy of us.
mov .pChildA,$New(BVHNode,Init,esi,.body,addr .volume)
; Child two holds the new body
mov .pChildB,$New(BVHNode,Init,esi,newBody, newVolume)
; And we now forget the body (we’re no longer a leaf).
mov .body , NULL
; We need to recalculate our bounding volume.
OCall esi.getBounds,newVolume ;calculate boundingsphere of this+newvolume
fsqrt ;we want actual diameter, not dia^2
fstp f1
OCall esi.makeVolume,newVolume,f1 ;
; Otherwise we need to work out which child gets to keep the inserted body.
; We give it to whoever would grow the least to incorporate it.
.else
OCall .pChildA::BVHNode.getBounds, newVolume
fstp f1
OCall .pChildB::BVHNode.getBounds, newVolume
fstp f2
.if $IsLess(f1,f2) ;comparing the sq. diameters of proposed boundingspheres
OCall .pChildA::BVHNode.Insert,newBody,newVolume
.else
OCall .pChildB::BVHNode.Insert,newBody,newVolume
.endif
.endif
MethodEnd
;Method: BVHNode.TwoSphereComparison
;Purpose: Load important data onto FPU for sphere comparison stuff
;Returns: ST(0) = SQUARED Sum of Radii
; ST(1) = SQUARED Distance Between Origins
BVHNode_TwoSphereComparison macro other;ptr BoundingSphere
mov edx,other
lea eax,.volume
fld .BoundingSphere.vPosition.x ;Calculate the sq. length of the line between origins
fsub .BoundingSphere.vPosition.x
fmul st(0),st(0)
fld .BoundingSphere.vPosition.y
fsub .BoundingSphere.vPosition.y
fmul st(0),st(0)
fld .BoundingSphere.vPosition.z
fsub .BoundingSphere.vPosition.z
fmul st(0),st(0)
fadd
fadd
fld .BoundingSphere.fRadius ;Calculate the sq. sum of the radii (ie amount by which to extend our line)
fadd .BoundingSphere.fRadius
fmul st(0),st(0)
endm
;The following method is called while choosing the best (smallest) possible boundingsphere:
;Method: BVHNode.getBounds
;Purpose: Calculate the SQUARED DIAMETER of a proposed new boundingsphere
; which encloses both 'this' sphere, and the given sphere ('other')
; Since we're only looking to find a 'best possible' boundingsphere,
; we can avoid a square root and just compare the sq. diameters.
;Args: other -> BoundingSphere representing the other sphere
;Returns: ST(0) = SQUARED diameter of new sphere
Method BVHNode.getBounds, uses esi, other:ptr BoundingSphere
SetObject esi
BVHNode_TwoSphereComparison other
fadd ;ST0 = squared diameter of minimal bounding sphere
MethodEnd
;Method: BVHNode.makeVolume
;Purpose: Recalculate BoundingSphere of minimal Volume which encloses both 'this' and 'other' sphere
;Args: other -> BVHNode representing 'other sphere'
; fnewDia = Diameter of the new BoundingSphere
;Returns: EAX -> BoundingSphere
Method BVHNode.makeVolume, uses esi, other:ptr BoundingSphere, fnewDia:real8
LOCAL diff:real8
SetObject esi
fld fnewDia ;Precalculate: diff = new diameter - this radius
fsub .volume.fRadius
fstp diff
fld .BoundingSphere.vPosition.x ;Calculate Line between Origins, and its Length
fsub .volume.vPosition.x
fst .volume.vPosition.x
fmul st(0),st(0)
fld .BoundingSphere.vPosition.y
fsub .volume.vPosition.y
fst .volume.vPosition.y
fmul st(0),st(0)
fld .BoundingSphere.vPosition.z
fsub .volume.vPosition.z
fst .volume.vPosition.z
fmul st(0),st(0)
fadd
fadd
fsqrt
fld .volume.vPosition.x ;Use Length to Normalize the Line into a Direction,
fld .volume.vPosition.y ;Scale it by (newDiam - this Radius), and then add our Point A ('this origin')
fld .volume.vPosition.z
fdivr
fmul diff
fadd .volume.vPosition.z
fstp .volume.vPosition.z
fdivr
fmul diff
fadd .volume.vPosition.y
fstp .volume.vPosition.y
fdivr
fmul diff
fadd .volume.vPosition.x
fstp .volume.vPosition.x
fld fnewDia
fmul r4_half
fstp .volume.fRadius
MethodEnd
It's not perfectly optimized, but it's not ABSOLUTELY horrible, and it will get the job done.
I haven't bothered providing sourcecode for all of the macros that have been used, that is an exercise for OA32 users if they choose to look them up. And please note this code has *NOT* been tested, so there MAY be mistakes - however rest assured that I am indeed implementing this, as this thread documents my ACTUAL coding adventures, not just what's rattling around in my head.
Just looking over that sourcecode, we can see that most of the cost involves calculation of line lengths.
In the next post, I'll show you a slightly different version, which I call "SpeedTree".
The difference is that we will no longer consider the radius of a sphere being inserted into the tree - we'll treat the input sphere as a single particle or point in 3D space. We'll still need to consider the spheres at each node of our tree, but by removing one sphere from the equation, we can reduce the math requirements by approximately half, and so effectively yield a speedup of 200 percent (or thereabouts). In practice it's of course less, but the speedup is nonetheless considerable, and it won't affect the tree's ability to group nearby entities, and so accelerate our broadphase collision testing.
Any speedups here will be incredibly important to our simulation, because we need to do broadphase testing every frame, and because our dynamic tree will be constantly rebuilt as our entities move in worldspace... narrowphase collision testing is VERY expensive, and the whole point of broadphase testing is to try to avoid narrowphase testing !!! So although it may seem like we're optimizing the wrong end of the collision spectrum, if you consider that it is in fact a culling operation, it makes complete sense.
Now I'm going to try to describe my proposed SpeedTree, both the theory behind it, and develop the algorithms here, this is NOT implemented but in my mind's eye, and on graph paper, it's looking good.
SpeedTree is a DYNAMIC Binary Space Partitioning tree... it is constructed from Planes, and we can find our entities in the Leaves.
Each plane passes through the origin of an associated entity, with some initial Direction determined by our algorithm. The direction of the plane will then follow that of the linear velocity, if any. And the plane is defined literally by the worldspace position of the entity's center of mass - to define a Plane in 3D, we need a Point (on the Plane), and a Direction. Thus, as the entity moves, the plane moves, and as the entity changes direction of travel, the plane rotates.
Now for the kicker.
Whenever we want to do Point/Plane tests, we ALWAYS give the Plane some "thickness" - normally this is a Numerical Error Tolerance that we simply call "Epsilon".
SpeedTree uses the Radius of the Entity's BoundingSphere as the plane thickness - thus, the plane is actually TWO planes, facing the same direction, and both tangent to the surface of the sphere, on either side of it, defining THREE SUBSPACES , which we usually call "back", "front", and "coplanar".
In fact, when we perform Sphere/Plane tests, it is normal to use this technique of giving the plane the thickness of the radius of the sphere we're testing against - in this case, we're giving it a default thickness of its construction sphere, and we can EXTEND that thickness by "other sphere" in our tests, thus we can effectively then treat the "other sphere" as a simple 3D point in space, so our per-node compare calculation becomes one of an addition followed by a dot product, sounds good!
Our tree construction algorithm will always be "top down insertive" - if we know the positions of all the spheres up front, we can pick one near the middle of the distribution as the starting point.. but I don't think the choice of the root node matters too much, our construction algorithm will assume that NONE of the entities have any initial velocity, and will instead choose the direction from the previously inserted entity to the currently inserted entity. This will cause the space to become most partitioned where clustering is most intense.
In the next post I'll try to develop these ideas further, hopefully providing sketches and pseudocode, and addressing the possible peculiarities of the proposed scheme.
Bounding Volume Hierarchies don't map space - they map the relationship between the positions and sizes of a bunch of spheres.. therefore they are not a true spatial partitioning scheme.
BSP trees map space into binary subspaces (two half-spaces), and normally these spaces are quite static and so the tree is static too.
SPEEDTREE is a hybrid scheme which attempts to fuse elements of both, and is based on planes which can move and rotate - like a BVH the tree is dynamic, but typically requires less updating due to the nature of halfspaces, and uses cheaper tests. Like a BSP, it is based on PLANES, but they are not static.
For the benefit of everyone reading this thread, I would like to go back to the nuts and bolts theory and math behind planar tests in 3D space. So I'll begin with a discussion of planes, point/plane distance test, sphere/plane intersection test, 4D swept test for intersection of a moving sphere and static plane, and finally, we'll derive a swept test for a sphere and a plane which are BOTH moving, and which is the cornerstone of SpeedTree.
PLANE
Imagine that we had an infinitely large, infinitely thin, infinitely flat piece of paper that cuts through the entire universe. It has a front side, and it has a back.. we can rotate it to face into some 3D direction, and we could move it (translate it) so that the plane cuts through some point in WorldSpace (the point is ON the plane, rather than behind it or in front of it).
We can define a Plane using the PLANE EQUATION: Ax + By + Cz + D = 0
where (x,y,z) represents a 3D point resting on our plane, and (A,B,C) represents the Direction which the plane is facing (ie, its NORMAL). So what then, is this unknown D?
D is referred to as the "Plane Distance" (or simply "PlaneD") - literally, it's the distance from the World Origin to our plane, and we can calculate it, given the other values. The full set of values of XYZ and ABCD defines a UNIQUE 3D PLANE.
It is quite common to store a Plane as (N,D) where N=(A,B,C), our Normal, and D was precalculated from some point KNOWN to be on the plane, using the previous equation, which is more clearly rewritten as as: D = - (P . N)
where P = point on plane, and the "dot" represents the DOT PRODUCT EQUATION (Ax + By + Cz)
Given a Plane stored in this form of a Normal vector and a D value, we can test any OTHER point in 3D space to see if it rests apon the plane, in front of the plane, or behind it. To measure the DISTANCE FROM ANY POINT TO A PLANE, we can use the following equation: d = (T - P) . N
where d = distance from testpoint to plane, in the direction of the Normal (shortest distance!!), and T is the point being tested.
We can again rewrite this as: d = N . T + D, allowing us to make use of a plane's precalculated D value.
If we feed values into that equation, let's look at the result: d will either be positive, negative, or zero.
If its zero (or very close - numerical tolerance please!) , then the test point coincides with the plane.
If its positive, then the point is in front of the plane, and otherwise it must be behind the plane, and by the distance of d, when measured along the direction of N.
Now we have an equation that we can use to check whether an arbitrary 3D point is behind, before, or coincides with a given unique 3D plane. This will be the keystone of everything that is to follow, so be sure you can understand how to define a unique plane in 3D, and calculate its D value, and can now write your own DOTPRODUCT function, and your own DistancePointPlane function, and test them.
Next we'll consider a Sphere and a Plane, and modify our existing function slightly to deal with the radius of a sphere.
BSP trees map space into binary subspaces (two half-spaces), and normally these spaces are quite static and so the tree is static too.
SPEEDTREE is a hybrid scheme which attempts to fuse elements of both, and is based on planes which can move and rotate - like a BVH the tree is dynamic, but typically requires less updating due to the nature of halfspaces, and uses cheaper tests. Like a BSP, it is based on PLANES, but they are not static.
For the benefit of everyone reading this thread, I would like to go back to the nuts and bolts theory and math behind planar tests in 3D space. So I'll begin with a discussion of planes, point/plane distance test, sphere/plane intersection test, 4D swept test for intersection of a moving sphere and static plane, and finally, we'll derive a swept test for a sphere and a plane which are BOTH moving, and which is the cornerstone of SpeedTree.
PLANE
Imagine that we had an infinitely large, infinitely thin, infinitely flat piece of paper that cuts through the entire universe. It has a front side, and it has a back.. we can rotate it to face into some 3D direction, and we could move it (translate it) so that the plane cuts through some point in WorldSpace (the point is ON the plane, rather than behind it or in front of it).
We can define a Plane using the PLANE EQUATION: Ax + By + Cz + D = 0
where (x,y,z) represents a 3D point resting on our plane, and (A,B,C) represents the Direction which the plane is facing (ie, its NORMAL). So what then, is this unknown D?
D is referred to as the "Plane Distance" (or simply "PlaneD") - literally, it's the distance from the World Origin to our plane, and we can calculate it, given the other values. The full set of values of XYZ and ABCD defines a UNIQUE 3D PLANE.
It is quite common to store a Plane as (N,D) where N=(A,B,C), our Normal, and D was precalculated from some point KNOWN to be on the plane, using the previous equation, which is more clearly rewritten as as: D = - (P . N)
where P = point on plane, and the "dot" represents the DOT PRODUCT EQUATION (Ax + By + Cz)
Given a Plane stored in this form of a Normal vector and a D value, we can test any OTHER point in 3D space to see if it rests apon the plane, in front of the plane, or behind it. To measure the DISTANCE FROM ANY POINT TO A PLANE, we can use the following equation: d = (T - P) . N
where d = distance from testpoint to plane, in the direction of the Normal (shortest distance!!), and T is the point being tested.
We can again rewrite this as: d = N . T + D, allowing us to make use of a plane's precalculated D value.
If we feed values into that equation, let's look at the result: d will either be positive, negative, or zero.
If its zero (or very close - numerical tolerance please!) , then the test point coincides with the plane.
If its positive, then the point is in front of the plane, and otherwise it must be behind the plane, and by the distance of d, when measured along the direction of N.
Now we have an equation that we can use to check whether an arbitrary 3D point is behind, before, or coincides with a given unique 3D plane. This will be the keystone of everything that is to follow, so be sure you can understand how to define a unique plane in 3D, and calculate its D value, and can now write your own DOTPRODUCT function, and your own DistancePointPlane function, and test them.
Next we'll consider a Sphere and a Plane, and modify our existing function slightly to deal with the radius of a sphere.
OK, what if we wanted to know the distance from a Sphere to a Plane?
If we treated the sphere as a simple Point, we know that: d = (T - P) . N where
d=distance from origin of sphere to plane, measured along plane normal
T=position of sphere origin in world space
P=any Point known to exist on the Plane, say, the one we used when we defined it
. =DotProduct operator
N=Normal of the Plane (Direction it faces)
Given that the Sphere has a Radius (R), ie it is a Point that has some THICKNESS, we can rewrite the equations as: d = ((T - P) . N) - R
In other words, the distance from the surface of the sphere to the plane is that of its origin, minus its radius.... we have taken the radius into account.
Because of the nature of computers, we have to contend with numerical accuracy.
So we will typically give our Plane some Thickness too, by introducing an "Epsilon" value.
The equation for a Point then becomes: d = ((T - P) . N) - Epsilon
Compare this to that for a Sphere: d = ((T - P) . N) - R
Note that we don't need an Epsilon value for the Sphere/Plane distance equation, because we can simply imagine our Sphere as a Point, and the Plane has having +/- the thickness of the Radius - we've simply got a bigger Epsilon :)
The next post will be cool, we'll look at our first 4D intersection test, a MOVING sphere and a Static Plane, over a given TimeStep.
Stare at this image for a while, we are looking at a 3D scene of a moving sphere and a fixed (static) plane, with the plane facing upwards, and looking at it from the side view.. the lower line represents the plane, and the thinner line above it has been shifted by the radius, in the direction of the plane normal. The sphere is shown at three moments in Time, starting at Time T0 and ending at Time T1... could be its current and next position, or its previous and current position, depending on how we think about it. But clearly in the example, it will definitely intersect the plane during its travel, and will first touch the plane at "Unknown Time TU".
We expect that when the sphere first touches the plane, that its origin will coincide with the upper line that represents the plane "shifted along its normal by the radius", and would LAST touch the plane when its origin coincided with a line (not shown) that represented the Plane shifted in the opposite direction by the same amount.
I have to break now, so start staring :D
We expect that when the sphere first touches the plane, that its origin will coincide with the upper line that represents the plane "shifted along its normal by the radius", and would LAST touch the plane when its origin coincided with a line (not shown) that represented the Plane shifted in the opposite direction by the same amount.
I have to break now, so start staring :D
So, we have a moving Sphere, and a fixed Plane.
Let's write down everything we KNOW about this situation:
At Time=T0, the Sphere's Origin is at Position=P0, which is distance=D0 from the Plane.
At Time=T1, the Sphere's Origin is at Position=P1, which is distance=D1 from the Plane.
Then there is U.
At Time=U, the Sphere's Origin is at a distance=R from the Plane.
If we can figure out what Time U equals, then we can interpolate to find the Position :)
The trajectory from P0 to P1 can be parameterized with a variable U, which may be thought of as normalized time, since its value is 0 at T0, P0 and 1 at T1, P1.... that is to say, U is some FRACTION between 0.0 and 1.0 :)
The normalized time at which the sphere first intersects the plane is given by: U=D0-R / (D0-D1)
The position of the sphere's origin at this TIME OF IMPACT (TOI) can then be interpolated with an affine combination of P0 and P1: PU = ((1-U) * P0) + (U * P1)
We are now able to tell exactly what Time the Sphere will touch the Plane, and where it will be.
Note that since the Time is still Normalized, if we want actual time we'll need to interpolate beween the actual times at T0 and T1, with our U coefficient again used as the interpolation factor.
In the next post, we'll look at some example sourcecode which addresses the "singularities" of this formula, ie, the special conditions under which it may fail.
I should have mentioned that our solution is a linear one - it assumes that the trajectory is a simple linear motion for the duration of the timeframe.
It would be more accurate to instead apply this kind of solution to the Acceleration or even the input Forces of the sphere, thus it becomes a second-order or higher (curve) solver. The Runge-Kutta 2nd (midpoint acceleration) or 4th order (RK4) numerical integration techniques could be adopted for this purpose, if such a level of accuracy was demanded. But we're just trying to find the lower bound of the TOI for a BOUNDING sphere, which contains some arbitrary geometry - broadphase, not narrow.
Let's write down everything we KNOW about this situation:
At Time=T0, the Sphere's Origin is at Position=P0, which is distance=D0 from the Plane.
At Time=T1, the Sphere's Origin is at Position=P1, which is distance=D1 from the Plane.
Then there is U.
At Time=U, the Sphere's Origin is at a distance=R from the Plane.
If we can figure out what Time U equals, then we can interpolate to find the Position :)
The trajectory from P0 to P1 can be parameterized with a variable U, which may be thought of as normalized time, since its value is 0 at T0, P0 and 1 at T1, P1.... that is to say, U is some FRACTION between 0.0 and 1.0 :)
The normalized time at which the sphere first intersects the plane is given by: U=D0-R / (D0-D1)
The position of the sphere's origin at this TIME OF IMPACT (TOI) can then be interpolated with an affine combination of P0 and P1: PU = ((1-U) * P0) + (U * P1)
We are now able to tell exactly what Time the Sphere will touch the Plane, and where it will be.
Note that since the Time is still Normalized, if we want actual time we'll need to interpolate beween the actual times at T0 and T1, with our U coefficient again used as the interpolation factor.
In the next post, we'll look at some example sourcecode which addresses the "singularities" of this formula, ie, the special conditions under which it may fail.
I should have mentioned that our solution is a linear one - it assumes that the trajectory is a simple linear motion for the duration of the timeframe.
It would be more accurate to instead apply this kind of solution to the Acceleration or even the input Forces of the sphere, thus it becomes a second-order or higher (curve) solver. The Runge-Kutta 2nd (midpoint acceleration) or 4th order (RK4) numerical integration techniques could be adopted for this purpose, if such a level of accuracy was demanded. But we're just trying to find the lower bound of the TOI for a BOUNDING sphere, which contains some arbitrary geometry - broadphase, not narrow.
The main problem is in our formula for U=D0-R / (D0-D1)
If D0=D1, ie the sphere didn't move at all (relative to the plane), then we will have a DIVIDE BY ZERO ERROR - but we can detect easily when D0=D1, and exit from our test function, since we know that intersection was not possible, because there was no relative velocity, ie displacement in position over time.
Let's see some sourcecode for this function which does eliminate that condition 'for free'..
Note that I've used a mixture of opcodes and macros , your own code can look how you want it to.
In the next post, we'll modify this test for a moving plane as well :)
If D0=D1, ie the sphere didn't move at all (relative to the plane), then we will have a DIVIDE BY ZERO ERROR - but we can detect easily when D0=D1, and exit from our test function, since we know that intersection was not possible, because there was no relative velocity, ie displacement in position over time.
Let's see some sourcecode for this function which does eliminate that condition 'for free'..
Note that I've used a mixture of opcodes and macros , your own code can look how you want it to.
;A macro for distance Point Plane calculation
;P = Point, N = Normal, D = PlaneD
distanceToPoint macro P, N, D
DotProduct P, N
fadd D
endm
;A function for calculating position and time of impact of sphere and plane
;Returns EAX = TRUE, ST(0) = TOI if there is a collision
;Otherwise returns EAX = FALSE
SpherePlaneSweep proc fRadius:real8, pvPos0, pvPos1, pvPlaneNormal, fPlaneD:real8, pvPosU
local fU:real8
local D0:real8
local D1:real8
mov eax,pvPlaneNormal
mov edx, pvPos0
distanceToPoint .Vec3, .Vec3, fPlaneD
fstp D0
mov edx, pvPos1
distanceToPoint .Vec3, .Vec3, fPlaneD
fstp D1
;check if it was touching on previous frame
fld D0
fabs
.if $IsLessOrEqual(ST(0), fRadius)
mov eax,pvPos0
mov edx,pvPosU
Vec3Copy .Vec3, .Vec3
fldz
return TRUE
.endif
;check if the sphere penetrated during this frame (and if indeed it moved at all)
.if ($IsGreater(D0,fRadius)) && ($IsLess(D1,fRadius))
; Pu = (1-u)*P0 + (u*P1) = point of first contact
fld D0
fsub fRadius
fld D0
fsub D1
fdiv ;leave U on fpu stack
fst fU
mov eax,pvPosU
mov edx,pvPos0
fld .Vec3.x
fld .Vec3.y
fld .Vec3.z
fld1
fsub fU
fmul st(3),st(0)
fmul st(2),st(0)
fmul
fstp .Vec3.z
fstp .Vec3.y
fstp .Vec3.x
mov edx,pvPos1
fld .Vec3.x
fld .Vec3.y
fld .Vec3.z
fld fU
fmul st(3),st(0)
fmul st(2),st(0)
fmul
fadd .Vec3.z
fstp .Vec3.z
fadd .Vec3.y
fstp .Vec3.y
fadd .Vec3.x
fstp .Vec3.x
return TRUE
.endif
return FALSE
SpherePlaneSweep endp
In the next post, we'll modify this test for a moving plane as well :)
In this post, we will consider a moving sphere and a moving plane. For the purpose of this post, we will consider only that the Plane is moving, and moving along the axis defined by its Normal (either direction)... rotation is not considered, we'll get to that though!
In order to calculate our U coefficient for a moving sphere and a moving plane, we need to, mathematically speaking, treat the plane as TWO planes, which we'll call Plane0 and Plane1.
As mentioned, they will both have the same Direction, ie, the Normal does not change.
But the position of the plane does, and this directly affects the PlaneD value.
Previously, we stated that "PlaneD is the distance from the World Origin to the Plane" (I forgot to mention that this is the distance measured in the direction of N, but that should be obvious by now).
Given that we know the Velocity of our moving plane, or that we can calculate it from P0, P1, and the timestep (T1-T0), we can determine DIRECTLY how much PlaneD has to be adjusted for Plane1.
And whether we add or subtract this amount can be determined by the Sign of Plane0.PlaneD, which tells us which side of the Plane we can find the world origin, and thus whether the Plane is moving toward or away from the world origin :)
So, all we need to do in order to calculate U for a moving sphere and a moving plane is to determine the "Plane1.PlaneD" value, pass this as an extra parameterm and use that when we calculate D1:
We're in business!
Now we can calculate the TOI and Position of first Contact of a moving sphere and a moving plane.
This can be very useful in narrowphase testing, but we're not yet interested in that.
What we care about is that this will allow us to construct and maintain a DYNAMIC BSP TREE.
But only as long as the Planes are not able to rotate...
In the next post, we'll see how to introduce Rotation of the Plane in our Swept Sphere test.
In order to calculate our U coefficient for a moving sphere and a moving plane, we need to, mathematically speaking, treat the plane as TWO planes, which we'll call Plane0 and Plane1.
As mentioned, they will both have the same Direction, ie, the Normal does not change.
But the position of the plane does, and this directly affects the PlaneD value.
Previously, we stated that "PlaneD is the distance from the World Origin to the Plane" (I forgot to mention that this is the distance measured in the direction of N, but that should be obvious by now).
Given that we know the Velocity of our moving plane, or that we can calculate it from P0, P1, and the timestep (T1-T0), we can determine DIRECTLY how much PlaneD has to be adjusted for Plane1.
And whether we add or subtract this amount can be determined by the Sign of Plane0.PlaneD, which tells us which side of the Plane we can find the world origin, and thus whether the Plane is moving toward or away from the world origin :)
So, all we need to do in order to calculate U for a moving sphere and a moving plane is to determine the "Plane1.PlaneD" value, pass this as an extra parameterm and use that when we calculate D1:
mov eax,pvPlaneNormal
mov edx, pvPos0
distanceToPoint .Vec3, .Vec3, fPlaneD0
fstp D0
mov edx, pvPos1
distanceToPoint .Vec3, .Vec3, fPlaneD1
fstp D1
We're in business!
Now we can calculate the TOI and Position of first Contact of a moving sphere and a moving plane.
This can be very useful in narrowphase testing, but we're not yet interested in that.
What we care about is that this will allow us to construct and maintain a DYNAMIC BSP TREE.
But only as long as the Planes are not able to rotate...
In the next post, we'll see how to introduce Rotation of the Plane in our Swept Sphere test.
Now we wish to introduce angular velocity into the equation.
Angular velocity is a 3D vector which can be decomposed into an axis and a rate of change of orientation about that axis, by convention in radians per second.
If our plane has a known angular velocity about its 'construction point', then we can easily rotate the Direction vector to find out where it is facing at the end of the timestep.
Simply put, we now have a situation where the two Planes (Plane0 and Plane1) change completely between T0 and T1. We can no longer calculate the value of Plane1.D as previously, because the Direction is no longer constant... we'll need to completely calculate the value of PlaneD at T1, as we did at T0.
Plane1.PlaneD is now dynamic, we can track the PlaneD value for use in T0, but we'll be calculating a new value for T1, and also for TU if there is intersection, at each timestep. So we drop the extra PlaneD param that we added to our function in the previous post, since we'll only know Plane0.D on entry to the function.
Once we've calculated U, interpolating the rotation for TU is extremely simple, and there's two ways to do it.The first is simply to interpolate between Plane0.N and Plane1.N, by the U coefficient.
The second way is to interpolate the magnitude of the angular velocity.. since we're going to assume that the angular velocity remains constant for the duration of the timestep, we simply need to extract the rate of angular change from the angular velocity vector (ie its magnitude), and scale it by U, then rotate the Plane0.N vector by that amount about the rotation axis. I'm not sure which is cheaper, I have a feeling they work out about the same. If I remember right, there's 12 multiplications and 9 additions in a trigonometric rotation formula, noting that we don't have a suitable transform matrix here. I guess interpolating the vector N is cheaper.
Finally, having calculated a "PlaneU.N", we can calculate a "PlaneU.D" to match.
So very little needs to change in our function.
But again we have a new singularity to deal with.
Our formula can fail if the angular velocity is allowed to exceed "PI radians per timestep" (180 degrees per timestep) ... this is , of course, incredibly fast rotation if our timesteps are sanely small, so it's fine that we place this upper limit on the angular velocity of all entities in the simulation.
We now have a function that can determine the TOI and position of first contact for a moving sphere and a moving, rotating plane :)
This is a VERY useful test, and it's not really costing a lot more than the previous one, especially since we'll remember the new PlaneD value at the end of each timestep for use in the next timestep (COHERENCE, remember?)
Angular velocity is a 3D vector which can be decomposed into an axis and a rate of change of orientation about that axis, by convention in radians per second.
If our plane has a known angular velocity about its 'construction point', then we can easily rotate the Direction vector to find out where it is facing at the end of the timestep.
Simply put, we now have a situation where the two Planes (Plane0 and Plane1) change completely between T0 and T1. We can no longer calculate the value of Plane1.D as previously, because the Direction is no longer constant... we'll need to completely calculate the value of PlaneD at T1, as we did at T0.
Plane1.PlaneD is now dynamic, we can track the PlaneD value for use in T0, but we'll be calculating a new value for T1, and also for TU if there is intersection, at each timestep. So we drop the extra PlaneD param that we added to our function in the previous post, since we'll only know Plane0.D on entry to the function.
Once we've calculated U, interpolating the rotation for TU is extremely simple, and there's two ways to do it.The first is simply to interpolate between Plane0.N and Plane1.N, by the U coefficient.
The second way is to interpolate the magnitude of the angular velocity.. since we're going to assume that the angular velocity remains constant for the duration of the timestep, we simply need to extract the rate of angular change from the angular velocity vector (ie its magnitude), and scale it by U, then rotate the Plane0.N vector by that amount about the rotation axis. I'm not sure which is cheaper, I have a feeling they work out about the same. If I remember right, there's 12 multiplications and 9 additions in a trigonometric rotation formula, noting that we don't have a suitable transform matrix here. I guess interpolating the vector N is cheaper.
Finally, having calculated a "PlaneU.N", we can calculate a "PlaneU.D" to match.
So very little needs to change in our function.
But again we have a new singularity to deal with.
Our formula can fail if the angular velocity is allowed to exceed "PI radians per timestep" (180 degrees per timestep) ... this is , of course, incredibly fast rotation if our timesteps are sanely small, so it's fine that we place this upper limit on the angular velocity of all entities in the simulation.
We now have a function that can determine the TOI and position of first contact for a moving sphere and a moving, rotating plane :)
This is a VERY useful test, and it's not really costing a lot more than the previous one, especially since we'll remember the new PlaneD value at the end of each timestep for use in the next timestep (COHERENCE, remember?)
Just wanted to add this:
Given a typical physics timestep of 0.1 seconds, the maximum rate of rotation is then 10*PI radians per second, which is 5 complete rotations of the plane per second, as the upper limit of the angular velocity before our formula can begin to fail.
Due to constraints apon my time, I am going to shelf the SpeedTree temporarily and instead implement a rather specialized Dynamic AABB Tree.
It's specialized because we will work with the assumption that all of the entities that it contains are Spheres (every RigidBody will have, at least, a BoundingSphere). We will use our knowledge of these spheres when building and maintaining the tree.
Also, it will depart from the traditional AABB tree in several other ways - I'll describe the system in depth, but you're going to have to wait a few hours, I need sleep :)
The traditional AABB tree is a straightforward Leafy Binary Tree where each node in the tree represents an axis-aligned bounding box (AABB), with nonleaf nodes having two child nodes, and leaf nodes having at least one entity.
Each non-leaf node has the property that its AABB encloses all of its child nodes’ AABBs, and
the tree is built by successively splitting the 'root' AABB in half (cutting through its PRINCIPAL AXIS), and then inserting the entities into one of the two child nodes based on which side of the split it is on (more on this later).
For a dynamic tree, we'll relax some of those constraints.. As previously stated, we don’t limit the number of child nodes that a node can have, and we allow entities to be stored in nonleaf nodes. And rather than build the entire tree based on static data, we only build subtrees when they are actually being queried (again, more on this later).
And we're going to have to choose a nesting heuristic for the entire tree, which dictates how nodes are split and how objects are stored in the child nodes. We'll take a look at a few candidates.
Finally, our implementation will only ever deal with Sphere entities, so we can make a few optimizations based on that fact.
So let's get started by covering the most important three functions in a Dynamic Partitioning Tree: INSERTING, DELETING AND UPDATING OUR SPHERE ENTITIES :)
Since this is about Dynamic trees, I am going to assume (as you should) that we already have a tree - if we don't, the first Insertion of an entity should generate a Root Node storing our entity.
I'm just going to cover the theory first, this probably won't make complete sense until we see some of the structures and code examples, but you should be able to see a few similarities to the SphereTree - this is in fact just another kind of BVH.
Inserting an entity into the tree is really easy... Starting at the root node, and using whatever nesting heuristic we have decided to use, choose the child which is the most appropriate container for the entity until either we land in a leaf node, or the nesting heuristic demands that the entity should be stored in the current nonleaf node. Associate the final tree node with the object, then expand the AABBs of all nodes from the current node, all the way up to and including the root node, as necessary.
Note that this can be done during the return from this recursion, with some care... that's a much more elegant way than walking back up through the tree via "parent pointers", or worse, RECURSING BACK UP the tree, I've seen both of those crap ideas implemented by people who are smarter than me.
Deleting is also quite simple :)
First we will remove the entity from the node which contains it.
Next, starting with the current node, check to see if the object’s bounding volume was 'kissing' against the edge of the AABB (ie, the AABB could potentially shrink as a result of this deletion).. if so, mark the AABB for recalculation, walk to the parent, and repeat this process until the object is no longer flush against the AABB.
Finally, we perform a simple garbage collection: if the node is empty of entities and has no children, we can delete the node, and delete its reference from its parent node. If the parent node now has no
children, then it should be UNSPLIT, and if it is ALSO devoid of entities, it can ALSO be deleted, and we simply repeat this all the way back to root, or until we hit a node that contains entities.
Updating is again, pretty straightward.. Since our spheres are rotationally invariant, the only time we need to UPDATE is when an entity moves, causing its BoundingBox to change.
First, locate the Node containing the entity which has Moved, and walk back up the tree until we find a node whose AABB will completely contain the sphere's boundingbox, or Root node, whichever we reach first. Next, use our nesting heuristic to determine the child node which should now contain the sphere. Finally, if this destination node is different than the sphere’s current node, remove the sphere from the current node and insert it into the destination node. Either way, recalculate the container node’s AABB.
This technique exploits both spatial and temporal coherence: if an object does not move between queries, then we have no work to do... and if an object does move, we've limited the search for its new container to the space within which it has actually moved. This will rarely require adjusting all the way back to root, that only happens when an entity steps outside the root node's box.
In the next post, we'll get into some of the dirty details of implementation specifics and we'll talk about a few nesting heuristics that we might choose from. We'll also talk about Node Splitting - when and why, and how to avoid the possible "Frgamentation" that can result from doing so.
Each non-leaf node has the property that its AABB encloses all of its child nodes’ AABBs, and
the tree is built by successively splitting the 'root' AABB in half (cutting through its PRINCIPAL AXIS), and then inserting the entities into one of the two child nodes based on which side of the split it is on (more on this later).
For a dynamic tree, we'll relax some of those constraints.. As previously stated, we don’t limit the number of child nodes that a node can have, and we allow entities to be stored in nonleaf nodes. And rather than build the entire tree based on static data, we only build subtrees when they are actually being queried (again, more on this later).
And we're going to have to choose a nesting heuristic for the entire tree, which dictates how nodes are split and how objects are stored in the child nodes. We'll take a look at a few candidates.
Finally, our implementation will only ever deal with Sphere entities, so we can make a few optimizations based on that fact.
So let's get started by covering the most important three functions in a Dynamic Partitioning Tree: INSERTING, DELETING AND UPDATING OUR SPHERE ENTITIES :)
Since this is about Dynamic trees, I am going to assume (as you should) that we already have a tree - if we don't, the first Insertion of an entity should generate a Root Node storing our entity.
I'm just going to cover the theory first, this probably won't make complete sense until we see some of the structures and code examples, but you should be able to see a few similarities to the SphereTree - this is in fact just another kind of BVH.
Inserting an entity into the tree is really easy... Starting at the root node, and using whatever nesting heuristic we have decided to use, choose the child which is the most appropriate container for the entity until either we land in a leaf node, or the nesting heuristic demands that the entity should be stored in the current nonleaf node. Associate the final tree node with the object, then expand the AABBs of all nodes from the current node, all the way up to and including the root node, as necessary.
Note that this can be done during the return from this recursion, with some care... that's a much more elegant way than walking back up through the tree via "parent pointers", or worse, RECURSING BACK UP the tree, I've seen both of those crap ideas implemented by people who are smarter than me.
Deleting is also quite simple :)
First we will remove the entity from the node which contains it.
Next, starting with the current node, check to see if the object’s bounding volume was 'kissing' against the edge of the AABB (ie, the AABB could potentially shrink as a result of this deletion).. if so, mark the AABB for recalculation, walk to the parent, and repeat this process until the object is no longer flush against the AABB.
Finally, we perform a simple garbage collection: if the node is empty of entities and has no children, we can delete the node, and delete its reference from its parent node. If the parent node now has no
children, then it should be UNSPLIT, and if it is ALSO devoid of entities, it can ALSO be deleted, and we simply repeat this all the way back to root, or until we hit a node that contains entities.
Updating is again, pretty straightward.. Since our spheres are rotationally invariant, the only time we need to UPDATE is when an entity moves, causing its BoundingBox to change.
First, locate the Node containing the entity which has Moved, and walk back up the tree until we find a node whose AABB will completely contain the sphere's boundingbox, or Root node, whichever we reach first. Next, use our nesting heuristic to determine the child node which should now contain the sphere. Finally, if this destination node is different than the sphere’s current node, remove the sphere from the current node and insert it into the destination node. Either way, recalculate the container node’s AABB.
This technique exploits both spatial and temporal coherence: if an object does not move between queries, then we have no work to do... and if an object does move, we've limited the search for its new container to the space within which it has actually moved. This will rarely require adjusting all the way back to root, that only happens when an entity steps outside the root node's box.
In the next post, we'll get into some of the dirty details of implementation specifics and we'll talk about a few nesting heuristics that we might choose from. We'll also talk about Node Splitting - when and why, and how to avoid the possible "Frgamentation" that can result from doing so.
Let's press forward.
Since we're now dealing with Boxes, let's describe a naiive approach to (re)calculating the AABB of a given node. For this, we will typically begin at Leaf nodes.
First, determine the AABB which encloses all the Spheres in the current node.
Next, get each child node’s AABB, and grow the current node’s AABB to accommodate it.
Finally, mark the AABB as current, and signal an update to the parent node, which may then require an update.
Note that once we've calculated a node's AABB, we only have to recalculate it when a sphere is deleted whose AABB was 'kissing' against the node's AABB... if it wasn't, then deleting it has no effect on the AABB. This will become more clear with some descriptive pictures and sourcecode, be patient :)
If we're prepared to add some more pointer variables to each node, we can greatly optimize this operation.
We'll do this by only considering the Origins of our spheres and child boxes (ie, treat everything as a Point in 3D space). For any given node directly containing one or more entities, we'll keep track which sphere's origin is furthest to the north, south, east, west, top and bottom of the node's box, taking into account the radius. This will allow us to more quickly recalculate a node's box based on knowledge of which contained elements have changed.
I'll definitely be posting a picture of that, because it's a fundamental optimization that I'll be implementing.
Assuming that we're able to efficiently manage the growing and shrinking tree of boxes, let's introduce our first heuristic: we will place a cap on the number of entities that can be stored in a node.
If we detect that a node has reached capacity, it will need to be Split, which involves distributing the node's contents into child nodes. This brings us to our next topic, FRAGMENTATION.
It's not only possible but LIKELY for entities to be added in such a way that it will “fragment” the tree, causing them to be inserted at much deeper levels of the tree than they would if the tree were built
containing them to begin with. This also causes the tree as a whole to become imbalanced, since we end up with one branch that is much deeper than the average case. To stop this from occurring, we will allow the AABBs of SIBLINGS (ie the immediate child nodes of a given node) to overlap in 3D space.
A good example of this is a situation where we have two spheres that are close to one another and offset diagonally, such that no Major Axis can be used to separate them without intersecting one or both.
In this case, we will place each child into the chld box which would grow least to accommodate it, and then grow the box accordingly, and finally alert the Parent node that the child box(es) grew (since the parent may also need to grow as it must CONTAIN all children).
I'll cover the finer points of tree management as the sourcecode is rolled out.
Let's take a look at some of the possible Nesting Heuristics. Although these are named after various kinds of Spatial Partitioning algorithm, we are ONLY using them to determine how and when to split up the nodes of our DAABB tree, so don't be confused by the names.
KDTree
This is the one used for traditional AABB trees, it reminds me a lot of a BSPTree heuristic: Given a set of Points (ie our sphere origins), we'll determine where the midpoint of the Distribution is located, and which of the BoundingBox's major axes is the longest one... then we'll determine a theoretical plane whose Normal is that longest axis, and which passes through that midpoint, thus cutting the box into two child boxes.
We will then distribute our Spheres into whichever child box contains their origin. Points on the splitting plane can be placed into either child, and most implementations will just use a convention, for example the Front halfspace... we can and should place them with a view to keeping our tree well-balanced, so we'll make a choice at runtime on that basis.
TernaryTree
This is just an extension of the previous heuristic: if a Sphere is intersected by (ie straddles) our axial cutting plane (and so rightfully belongs to BOTH halfspaces), we will place it into a THIRD subtree, or alternatively, we can just leave it in the current node. The latter is the technique I chose for previous work with BSP trees. I've never personally implemented the 'third subtree techique' before now, and it's still quite unintuitive to me as it makes spatial queries a lot more complex (requiring dual tree walk), whereas simply leaving those objects where they are makes perfect sense and can be walked linearly.
Octree
This is also an extended KD heuristic: we find the KD midpoint, then we cut the box using all three axial planes passing thru that point, thus we don't care which axis is 'longest' anymore.
This cuts the box into eight subspaces, we then distribute our Spheres as per KD.
Icoseptree:
Wow, what a name, huh? It's the Octree heuristic, extended with the same logic as a Ternary heuristic.
There are now three subtrees for each major axis, and so each (nonleaf) node has a total of 27 subspaces available for distributing our spheres.
It is my hope that I have been eloquent rather than simply brief in these descriptions.
If you really want to hurt yourself, go study each of these Static tree structures and tell me if I missed anything important.
In the next post, I hope to provide some pictures to back up these ideas, and other than that we're ready to start putting some code together :)
Since we're now dealing with Boxes, let's describe a naiive approach to (re)calculating the AABB of a given node. For this, we will typically begin at Leaf nodes.
First, determine the AABB which encloses all the Spheres in the current node.
Next, get each child node’s AABB, and grow the current node’s AABB to accommodate it.
Finally, mark the AABB as current, and signal an update to the parent node, which may then require an update.
Note that once we've calculated a node's AABB, we only have to recalculate it when a sphere is deleted whose AABB was 'kissing' against the node's AABB... if it wasn't, then deleting it has no effect on the AABB. This will become more clear with some descriptive pictures and sourcecode, be patient :)
If we're prepared to add some more pointer variables to each node, we can greatly optimize this operation.
We'll do this by only considering the Origins of our spheres and child boxes (ie, treat everything as a Point in 3D space). For any given node directly containing one or more entities, we'll keep track which sphere's origin is furthest to the north, south, east, west, top and bottom of the node's box, taking into account the radius. This will allow us to more quickly recalculate a node's box based on knowledge of which contained elements have changed.
I'll definitely be posting a picture of that, because it's a fundamental optimization that I'll be implementing.
Assuming that we're able to efficiently manage the growing and shrinking tree of boxes, let's introduce our first heuristic: we will place a cap on the number of entities that can be stored in a node.
If we detect that a node has reached capacity, it will need to be Split, which involves distributing the node's contents into child nodes. This brings us to our next topic, FRAGMENTATION.
It's not only possible but LIKELY for entities to be added in such a way that it will “fragment” the tree, causing them to be inserted at much deeper levels of the tree than they would if the tree were built
containing them to begin with. This also causes the tree as a whole to become imbalanced, since we end up with one branch that is much deeper than the average case. To stop this from occurring, we will allow the AABBs of SIBLINGS (ie the immediate child nodes of a given node) to overlap in 3D space.
A good example of this is a situation where we have two spheres that are close to one another and offset diagonally, such that no Major Axis can be used to separate them without intersecting one or both.
In this case, we will place each child into the chld box which would grow least to accommodate it, and then grow the box accordingly, and finally alert the Parent node that the child box(es) grew (since the parent may also need to grow as it must CONTAIN all children).
I'll cover the finer points of tree management as the sourcecode is rolled out.
Let's take a look at some of the possible Nesting Heuristics. Although these are named after various kinds of Spatial Partitioning algorithm, we are ONLY using them to determine how and when to split up the nodes of our DAABB tree, so don't be confused by the names.
KDTree
This is the one used for traditional AABB trees, it reminds me a lot of a BSPTree heuristic: Given a set of Points (ie our sphere origins), we'll determine where the midpoint of the Distribution is located, and which of the BoundingBox's major axes is the longest one... then we'll determine a theoretical plane whose Normal is that longest axis, and which passes through that midpoint, thus cutting the box into two child boxes.
We will then distribute our Spheres into whichever child box contains their origin. Points on the splitting plane can be placed into either child, and most implementations will just use a convention, for example the Front halfspace... we can and should place them with a view to keeping our tree well-balanced, so we'll make a choice at runtime on that basis.
TernaryTree
This is just an extension of the previous heuristic: if a Sphere is intersected by (ie straddles) our axial cutting plane (and so rightfully belongs to BOTH halfspaces), we will place it into a THIRD subtree, or alternatively, we can just leave it in the current node. The latter is the technique I chose for previous work with BSP trees. I've never personally implemented the 'third subtree techique' before now, and it's still quite unintuitive to me as it makes spatial queries a lot more complex (requiring dual tree walk), whereas simply leaving those objects where they are makes perfect sense and can be walked linearly.
Octree
This is also an extended KD heuristic: we find the KD midpoint, then we cut the box using all three axial planes passing thru that point, thus we don't care which axis is 'longest' anymore.
This cuts the box into eight subspaces, we then distribute our Spheres as per KD.
Icoseptree:
Wow, what a name, huh? It's the Octree heuristic, extended with the same logic as a Ternary heuristic.
There are now three subtrees for each major axis, and so each (nonleaf) node has a total of 27 subspaces available for distributing our spheres.
It is my hope that I have been eloquent rather than simply brief in these descriptions.
If you really want to hurt yourself, go study each of these Static tree structures and tell me if I missed anything important.
In the next post, I hope to provide some pictures to back up these ideas, and other than that we're ready to start putting some code together :)
Well I tried, and I'm just a terrible artist so I will just keep moving and if you're really lucky some images will magically appear in these posts at a later date.
I am a really BAD artist, but I can code.
So I will show you code instead, since that's why we're all here, isn't it?
We need to start somewhere, let's start with the spine and get this thread moving :)
I am a really BAD artist, but I can code.
So I will show you code instead, since that's why we're all here, isn't it?
We need to start somewhere, let's start with the spine and get this thread moving :)
; ——————————————————————————————————————————————————————————————————————————————————————————————————
;Some constants we'll reference in our 'nesting heuristic':
UNDEFINED equ -1
X_AXIS equ 0
Y_AXIS equ 1
Z_AXIS equ 2
; ——————————————————————————————————————————————————————————————————————————————————————————————————
; ——————————————————————————————————————————————————————————————————————————————————————————————————
;The BBNode class implements the nodes of a Dynamic Axially-Aligned BoundingBox (DAABB) Tree
;which has been optimized specifically to store Sphere entities for BroadPhase collision detection.
Object BBNode, 345342, Primer
RedefineMethod Init, Pointer
RedefineMethod Done
DefineVariable vMin, Vec3, {<>} ;Current extents of this node's AABB
DefineVariable vMax, Vec3, {<>}
DefineVariable vSplit, Vec3, {<>} ;Point at which the Node was "split" into subnodes
;
Definevariable dAxis, dword, UNDEFINED
;We use these pointers to keep track of which spheres are pressing on the sides of the node's AABB
;This in turn lets us keep track of the dynamic size of an AABB without having to evaluate all of its child spheres.
;Since this optimization is not implemented yet, I will leave these commented out for the moment.
; DefineVariable pFurthestEast, Pointer
; DefineVariable pFurthestWest, Pointer
; DefineVariable pFurthestNorth, Pointer
; DefineVariable pFurthestSouth, Pointer
; DefineVariable pFurthestDown, Pointer
; DefineVariable pFurthestUp, Pointer
Embed children, Collection ;holds array of child nodes (ie subtrees, branches, as you prefer)
Embed bodies, Collection ;holds array of RigidBody objects
ObjectEnd
LPBBNODE typedef ptr BBNode
; ——————————————————————————————————————————————————————————————————————————————————————————————————
;The BBTree class is mostly just a container for the Root Node of our DAABB tree.
;It implements some "entry-level" methods for recursive functions.
Object BBTree, 234566, Primer
DefineVariable pRootNode, LPBBNODE, NULL
ObjectEnd
; ——————————————————————————————————————————————————————————————————————————————————————————————————
;Method: BBNode.Init
;Purpose: Initialize BBNode object
;Args: pOwnerParent -> BBNode which is the immediate parent of this node
;Returns: nothing
Method BBNode.Init, uses esi, pOwnerParent
SetObject esi
ACall Init, pOwnerParent
OCall .children::Collection.Init, esi, 16, 256, UNDEFINED
OCall .bodies ::Collection.Init, esi, 16, 256, UNDEFINED
MethodEnd
;Method: BBNode.Done
;Purpose: Destroy this node
Method BBNode.Done, uses esi
SetObject esi
;Release this node's resource objects
OCall .children::Collection.Done ;subtrees will be destroyed recursively, this is automated by ObjAsm
OCall .bodies ::Collection.Done
MethodEnd
I must apologize for the lack of recent posts, it's not because I've lost interest, rather the opposite.
There's no point in blindly surging forward with this, I've done mountains of research and studied everything I could lay my hands on related to this work... it's important stuff !!
As mentioned previously, there are quite a number of possible heuristics that we might employ for the purpose of determing If, When and How a Node will be split up into child nodes (subtrees).
And there has been quite a LOT of research done into which ones are best under various conditions.
The goal of our BVH is to cull pairs of entities which cannot possibly collide during the BroadPhase collision testing, and the benefits of having a broadphase are proportional to the number of entities.
IE, with low numbers of entities, there's almost no point having a BVH at all, this is PURELY about HIGH NUMBERS OF ENTITIES !!! I am determined to beat the Bullet Engine, which already is better and faster than the PhysX engine and others. Bullet implements only a Leafy, Binary Tree using the KD-Tree heuristic, which means that ALL entities end up in the Leaf Nodes, and that any non-leaf node has only TWO subspaces (child nodes) - this can create very deep trees so we end up wasting a LOT of time recursing and adjusting the tree.. think we can do better.
Bearing all that in mind, the heuristic which scales best for LOTS of entities is ICOSEPTREE.
And so, this is the one that I will be implementing and talking about, I will simply ignore all the others.
The Icoseptree divides each Non-Leaf node into 27 subspaces.
What follows is perhaps the most important function in our BVH code, the implementation of that heuristic: given a Sphere represented by a Point and Radius, this method must decide which of the 27 subspaces contains our Sphere. Note that this is a Recursive function, it will "Walk the Tree" - at each Node it will make the decision which Child Node (subspace) contains our Sphere entity, and it will then Recurse that node, until it reaches a Leaf Node, or until it reaches a Node whose BoundingBox is quite a lot smaller than the entity... we've already said that we are willing to have entities in Non-Leaf nodes, that decision is not arbitrary - we don't want large entities to be shoved too far down the tree because there's a high chance that they will straddle multiple subspaces and we would rather they fall neatly into exactly one subspace - we're trying to avoid 'growing' the child nodes to encompass them, because that would force the Parent nodes to expand such that their children don't violate their bounds.
The code provided mentions some Macros which are not provided here since they would only serve to confuse - the comments explain what's going on - any questions?
There's no point in blindly surging forward with this, I've done mountains of research and studied everything I could lay my hands on related to this work... it's important stuff !!
As mentioned previously, there are quite a number of possible heuristics that we might employ for the purpose of determing If, When and How a Node will be split up into child nodes (subtrees).
And there has been quite a LOT of research done into which ones are best under various conditions.
The goal of our BVH is to cull pairs of entities which cannot possibly collide during the BroadPhase collision testing, and the benefits of having a broadphase are proportional to the number of entities.
IE, with low numbers of entities, there's almost no point having a BVH at all, this is PURELY about HIGH NUMBERS OF ENTITIES !!! I am determined to beat the Bullet Engine, which already is better and faster than the PhysX engine and others. Bullet implements only a Leafy, Binary Tree using the KD-Tree heuristic, which means that ALL entities end up in the Leaf Nodes, and that any non-leaf node has only TWO subspaces (child nodes) - this can create very deep trees so we end up wasting a LOT of time recursing and adjusting the tree.. think we can do better.
Bearing all that in mind, the heuristic which scales best for LOTS of entities is ICOSEPTREE.
And so, this is the one that I will be implementing and talking about, I will simply ignore all the others.
The Icoseptree divides each Non-Leaf node into 27 subspaces.
What follows is perhaps the most important function in our BVH code, the implementation of that heuristic: given a Sphere represented by a Point and Radius, this method must decide which of the 27 subspaces contains our Sphere. Note that this is a Recursive function, it will "Walk the Tree" - at each Node it will make the decision which Child Node (subspace) contains our Sphere entity, and it will then Recurse that node, until it reaches a Leaf Node, or until it reaches a Node whose BoundingBox is quite a lot smaller than the entity... we've already said that we are willing to have entities in Non-Leaf nodes, that decision is not arbitrary - we don't want large entities to be shoved too far down the tree because there's a high chance that they will straddle multiple subspaces and we would rather they fall neatly into exactly one subspace - we're trying to avoid 'growing' the child nodes to encompass them, because that would force the Parent nodes to expand such that their children don't violate their bounds.
The code provided mentions some Macros which are not provided here since they would only serve to confuse - the comments explain what's going on - any questions?
;Method: BBNode.WhichChild
;Purpose: Determine which Icoseptree Child rightfully owns the given SPHERE (origin, radius)
Method BBNode.WhichChild,uses esi, pos, radius:real8
LOCAL _c,sv:real8,bv:real8
SetObject esi
;If this is a Leaf node, we HAVE to leave the sphere here in this node, theres no choice
;We might decide to SPLIT this node, but we sure won't do it NOW - we'll defer that decision.
.if !.bIsSplit
return -1
.endif
and _c,0
; If the sphere’s volume is at least 1/8 of this node's AABB volume
; then we'll just leave the sphere here in this node.
; I think this should probably be 1/27 , what do you think?
.if .bExtentsAreCorrect ;better only do this when the extents are good
fldpi ;volume of sphere
fmul r4_4
fdiv r4_3
fmul radius
fstp sv ;=4/3 pi r
fld .vMax.x ;volume of box / 8
fsub .vMin.x
fld .vMax.y
fsub .vMin.y
fld .vMax.z
fsub .vMin.z
fmul
fmul
fdiv r4_8
fstp bv ;= (length*width*height) / 8
.if $IsGreater(sv, bv)
return -1 ;"leave it here"
.endif
.endif
;ICOSEPTREE CHILD SELECTION:
;The following three tests will yield a value from 0 to 26 inclusive (ie 27 possible values)
;This indicates which Icoseptree Child owns the input sphere
mov edx,pos
;We'll perform tests for each Major Axis...
;Is the Sphere Origin plus the Radius less than this node's 'split point'
.if $IsAddLess(.Vec3.x, radius, .vSplit.x)
add _c,1
;Otherwise, is Origin minus Radius greater than split point?
.elseif $IsSubGreater(.Vec3.x, radius, .vSplit.x)
add _c,2
.endif
.if $IsAddLess(.Vec3.y, radius, .vSplit.y)
add _c,3
.elseif $IsSubGreater(.Vec3.y, radius, .vSplit.y)
add _c,6
.endif
.if $IsAddLess(.Vec3.z, radius, .vSplit.z)
add _c,9
.elseif $IsSubGreater(.Vec3.z, radius, .vSplit.z)
add _c,18
.endif
;Assert that the indicated child index is valid in this node
mov eax,_c
.if eax >= .children.dCount
DbgWarning "Error in BBNode.WhichChild - the indicated child node does not exist"
DbgDec _c,"Icoseptree child index"
;We COULD take this opportunity to create the necessary Child Node,
;but since its NOT a leaf, we expected it to already have been done?
int 3
.endif
MethodEnd