-------------------------------------------------------------------- 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.