--------------------------------------------------------------------
On the KD-tree construction algorithms with surface area heuristics
--------------------------------------------------------------------
The original text written by Vlastimil Havran, 19th November 2005.
Version of the text: 1.0, 19th November 2005.
The text is available at
http://www.cgg.cvut.cz/~havran/kdtreecontstruction.html
-----------------------------------------------------------------------
This text is an addendum of the doctoral thesis:
"Heuristics Ray Shooting Algorithms",
written by Vlastimil Havran,
submitted on 30th November 2000 and defended on 20th of April 2001 at
the Czech Technical University, Prague.
The thesis is available at
http://www.cgg.cvut.cz/~havran/phdthesis.html
The text below is an addendum A1 to the thesis, explaining the kd-tree
construction in O(N log N) time based on the surface area heuristics,
part 4.9 - Time complexity Analysis.
---------------------------------------------
Introduction:
David MacDonald and Kellogg S. Booth (thesis supervisor, web page at
http://www.cs.ubc.ca/~ksbooth/ ), did quite a good work at that time
concerning the efficiency of the data structures for ray tracing. The
construction time for their kd-tree was probably not O(N*log(N)), as
they mention in their paper, but much worse since this was not the
problem they have studied, but still they had constructed fairly good
kd-trees. The kd-tree construction was also improved later, see
chapter 4 of the thesis.
Since the issue of the kd-tree construction time complexity was raised
by a few people during four years, including
Colin Withers in 2002,
Ingo Wald (http://graphics.cs.uni-sb.de/~wald) in 2003,
Daniel Neilson in 2004,
Cyril Soler in 2004, and
Eric Haines during SIGGRAPH 2005,
I have decided to put more explanation to the web page, while the
efficient solution is sketched briefly in part 4.9 - Time complexity
Analysis, in the thesis revision 1.1 on page 80, the third paragraph.
So let us describe three algorithms for kd-tree construction with
decreasing complexity, the efficient algorithm in O(N * log N) time
referred as Algorithm~3.
-------------------------------------------------------
Algorithm~1:
A naive approach for a single axis is to take every possible object
boundary as a splitting plane in the current box, there are 2*N
boundaries, find the objects are on the left (N_L), which objects are
on the right (N_R), and the number of objects straddling the tested
splitting plane. Then you need for a single axis time O(2*N*N), which
is terrible. If the objects are uniformly distributed in the scene
(the easiest case for this construction), you have then O(N-1)
interior nodes and O(N) leaves, each leaf just one object, then the
total time for kd-tree construction is O(3*N*N). In this text I put
the constants to the big-Oh notation, since it says something.
--------------------------------------------------------
Algorithm~2:
So, better approach is this: Let us improve by sorting. Before
evaluation of the cost for every axis you do sorting in
O(2*N*log(2*N)) time, but let us again assume a uniform distribution
of objects for simplicity of the analysis. So after splitting we have
N/2 objects on the left and N/2 objects on the right. Assuming that we
have again (N-1) interior nodes and N leaves, then the total sum is as
you can derive O(N*2 * log(N*2) * (log(N*2) - 1/2)). So in fact less
than quadratic but more than O(N * log (N)), without constants it is
O(N * log(N) * log(N)). All the logarithms are base of 2.
--------------------------------------------------------
Algorithm~3:
The Algorithm 2 seems to be what most people do when doing their first
implementation of kd-tree construction, most likely they reveal that
Algorithm 1 is really too inefficient during early desgin. This includes
the problems and complaints of the people that showed their interest
in the problem I have communicated to, explaining the things either
by an e-mail or by telling that.
--------
As it is written in my thesis, part 4.9 - Time complexity Analysis, in
the thesis revision 1.1 on page 80, the third paragraph, the time
complexity of such a construction is O(N * log N). Since I have never
designed the Algorithm~1 or Algorithm~2, but directly algorithm~3 in
a tight collaboration with Martin Devera
(http://luxik.cdi.cz/~devik/qos/htb/) around the year 1996 and 1997, I
have never considered to make it a publication, since we have assumed
that people has to do it in the same way as we did. As I have found later,
that assumption was true, quite a similar algorithm was published in
computational geometry for sorting in K-dimensional space:
"An O(n * log n) Algorithm for the All-Nearest-Neighbors Problem",
written by Pravin M. Vaidya, Discrete Computational Geometry 4:
101-115, 1989.
The paper of Vaidya is cited in various contexts and this is not the
only algorithm of O(N log N) time complexity for the problem. Our
method is a direct application of Vaidya's algorithm in fact on the
construction of kd-trees with surface area heuristics. The time
complexity analysis holds for arbitrary distribution of objects.
----
Let us now describe the algorithm in detail, from a simple version in
1D for points, then for line segments in 1D, and finally for 3D boxes
in 3D space.
-------
Version for points in 1D space
-------
Let us describe first for a single axis in informal way when you
construct a tree over numbers only in one-dimensional space and you
evaluate some cost function for every possible position. So you sort
the numbers in an array in O(N * log N) time. Then you have an array
sorted, you go through the array linearly, evaluate the cost function
and decide, where you put your "splitting plane" in 1D. You have now
two subarrays, which are still sorted, then you recurse, until you
finish with leaves. You are finished when you got a single number to
the leaves of a tree. This is fairly simple to understand, design, and
implement for such a situation.
-----
Version for line segments in 1D space
-----
Now we can make it slightly more complicated. Instead of numbers (1D
points) we will have 1D line segments (=1D boxes or 1D intervals). The
line segments are represented by points, one point for the left
boundary of the line segment, one point for the right boundary of the
line segment. At every point we keep the bidirectional pointers; from
the left boundary to the right boundary and from the right boundary to
the left boundary. We keep an information if the point represents left
or right boundary of the line segment. Also we keep a pointer from the
left boundary point to some data connected to that line segment. Then
the rest of the algorithm is the same as already described, but
instead of an array we keep the point data in a bidirectionally linked
list. So for every point we have the link to the predecessor and to
the successor on 1D axis. The structure is like that:
struct Data1D {
bool leftOrRight; // either left or right
void *data; // the object data
// a pointer from left boundary to right boundary or a pointer from
// the right boundary to the left boundary
Data1D *otherBoundary;
Data1D *previous; // on the x-axis the previous number when sorted
Data1D *next; // on the x-axis the next number when sorted
};
We sort all the points in 1D using array, during sorting we update the
pointers for boundaries (Data1D.otherBoundary). After sorting we set
the links Data1D.previous and Data1D.next. We find a candidate for a
splitting plane evaluating the cost function for all the possible
entries in the list and make two sublists, the left sublist and the
right sublist. However, this is not finished, the answer to the
question comes: why we need the bidirectional links Data1D.previous
and Data1D.next ?
If after decision on splitting plane position we create two sublists,
it can happen that any line segment straddles the splitting
plane. This is easy to detect, but we have to traverse all the entries
in the left sublist. If going from the left side to the right side of
the left sublist using Data1D.next pointer, we check for every point
representing the left boundary using a pointer Data1D.otherBoundary,
if its right boundary is on the left or right side of the splitting
plane. If the right boundary is on the left side, it means that also
the whole object is on the left side of the splitting plane. In the
opposite case, the line segment straddles the splitting plane. In the
case we duplicate the right boundary of this line segment creating a
new entry Data1D in the list, we set the position of this right
boundary to the 1D splitting plane value. The duplicated boundary is
appended at the end of the left sublist. The pointers from the left to
the right boundary have to be reinitialized in the left boundary entry.
Similarly for this line segment straddling the splitting plane we also
create a new left boundary entry and prepend it at the beginning of
the right sublist using Data1D.previous and Data1D.next pointers. To
summarise, we have to create for a line-segment two new boundaries at
the position of the splitting plane and insert the right boundary at
the end of the left sublist and the left boundary at the beginning of
the right sublist.
-----
Version for 3D boxes in 3D space
(also holds to N-D boxes in N-D space)
-----
So we have solved it in linear time for that 1D axis in one pass, then
the whole construction remains O(N * log N). The construction for 3D
boxes, so for kd-tree over 3D objects simplified to 3D boxes is in
fact the same as for 1D line segments in 1D. After splitting we are
only forced to keep updated also the links for other two axes, which
means to go through the lists of two remaining axes and deciding if
a current entry belongs to the left sublist in that axis or the right
sublist. So in addition to 1D, our data structures are extended by new
pointers for two other axes:
struct Data3D {
bool leftOrRight; // either left or right
void *data; // the object data
// a pointer from left boundary to right boundary or a pointer from
// the right boundary to the left boundary
Data3D *otherBoundary;
Data3D *previous[3]; // on the x/y/z-axis the previous boundary when sorted
Data3D *next[3]; // on the x/y/z-axis the next boundary when sorted
};
In other words, during the first traversal through axis to be
subdivided we have created the left sublist and the right sublist for
this axis. During the second traversal using 'nextAxis' (the other
than splitting axis) we just updated links to data.previous[nextAxis],
data.next[nextAxis] creating two other sublists, for the left subarray
and for the right subarray. Similarly, during the third traversal for
"nextnextAxis" we update the data data.previous[nextnextAxis],
data.previous[nextnextAxis], creating next two sublists.
Then we have resorted data for both children of the currently splitted
node for all three axes in O(N) time, getting in total six sublists,
three sublists for the left subtree to be build, and the next three
sublists for the right subtree to be build. Then we can apply any
search for minima of the cost function to any axis and recurse.
In case you do not use split clipping described in the section 4.6.2
of the thesis and you want to really optimize the construction time,
then there is futher way to increase the efficiency of the whole
kd-tree construction by a constant factor only. So assuming you do not
need some additional precomputation for your cost function, you can in
fact evaluate the cost function already during creating two child
lists of Data3D entries, when creating the sublists for the left child
and also during resorting phase for other two axes.
A special note for the cost function evaluation is the fact that for
the boxes of objects sharing the same values for the left and right
boundaries at the same time, it is advisable to put to the lists first
the right boundaries and then left boundaries. This decreases the
value of cost function locally to minimum.
------------------------------------------------------------
Question: Why there are bidirectional lists in Algorithm~3 ?
The bounding boxes that are straddling the splitting plane, they can
be reduced usually also in other two axes. This complicates the
algorithm only a bit, and this is why it is preferable to use a list
instead of array. In principle, you have to resort new positions of
clipped bounding boxes by going through the the list towards a new
position of the reinserted boundary, the left boundary along axis is
reinserted backwards in the list using Data3D.prev[nextAxis], the left
boundary is reinserted forwards in the list, using
Data3D.next[nextAxis]. The same is carried out for "nextnextAxis".
For resorting we need searching. Searching in the list takes
complexity in O(N) time at most for a single boundary entry, however,
in an array it could have been carried out in principle in O(log N)
time or even in O(log log N) time. In practice, the split clipped
boundary lie closely to the splitting plane, so only O(1) steps are
needed to find a new position. If the size of the boxes around the
objects is limited in size with respect to the bounding box of the
whole scene containing all objects and the objects are fat, it can be
proved in formal way.
However, there is another problem than searching for a new boundary
than just searching. For resorting also the deletion and insertion of
the boundary at a new position has to be done. This is impossible for
an array in an efficient way, since then the resorting is implemented
by copying in O(N) time. Obviously, you cannot make an insertion to
the array in constant time and copying is quite costly, so performance
wise at the end the use of list is better and under assumption above
the complexity remains O(N * log N) at the end.
If you do not make split clip operation for the objects straddling the
boundary, you can stay with unidirectional links using Data3D.next
only, or you have to copy the whole array, if you work with arrays, in
O(N) time. Both solutions have their advantages and disadvantages
concerning memory access pattern and overall efficiency.
---------------------------------------------------------------
For further questions on the algorithm if not clear, contact
Vlastimil Havran
, the e-mail address can be found on the web page:
http://www.google.de/search?hl=de&q=Vlastimil+Havran&btnG=Google-Suche&meta=
---------------------------------------------------------------
End of the document.