Efficient linear octree generation from voxels

Efficient linear octree generation from voxels Renben Shu and Mohan S Kankanhalli

is an efficient data structure for compact representation of voxel data. There are several algorithms to generate a linear octree from voxel data with time complexity O(2) for an input of size 2 voxels. We present a new algorithm which first extracts the surface of the object. Based on this surface data, the object is partitioned into a s:t of parallelepipeds where each parallelepiped is a contiguous run of voxels along one axis. Starting from the lowest level of the octree, the algorithm proceeds iteratively to the highest level. computing maximal overlaps of the parallelepipeds at each level. For any level, the voxels which are not in the overlap are octree nodes and are output at that level. The maximal overlapped parallelepipeds form the input to the next higher level in the algorithm. For a connected object having n3 voxels, the algorithm has a time complexity of O(S) where S is the size of the surface of the object. The algorithm has been implemented and tested for a variety of medical data.

The linear

Keywords:

octree

quadtreeloctree

construction,

linear

quadtrees/

octrees

Medical visualization involves a huge amount of data. The interactive needs of medical visualization’.’ mandate fast processing of data. There is a need for compact storage and efficient handling of the voxel data. The octree structure is an efficient representation technique, which leads to less storage3 and is amenable to different kinds of geometric operations. The octree representation of a structure is a hierarchical 2” x 2” x 2” binary voxel array. The octree is obtained by recursively dividing the space into eight octants. Each octant is fully white, fully black or partially full. Fully white or black octants are not further subdivided, but the partially full octants are recursively divided until maximal fully white or black octantsm are obtained.

Institute

of Systems Science, National University of Singapore. Kent Ridge. Singapore 051 I (email: mohankc iss.nus.sg) Prpr rrceid: 10 Aqust 1991; re~~i.sedpup receiwtl 25 .Nowmhrr 1993

structure is represented by a tree in which the root corresponds to the entire array. Every non-terminal node has eight children, representing the octants. The leaves of the tree are either fully white or fully black. Converting the voxel representation of data into an octree representation has traditionally been a slow process4. In this paper, we present a new algorithm to convert the voxel representation of a connected object into an octree representation which takes O(S) time, where S is the size of the surface of the object. Octree methods have been successfully used in a variety of applications’ ‘. Samet’s workx,’ is a comprehensive reference in the tield.There have been several pixel-to-quadtree and voxel-to-octree conversion algorithms”. Chen and Huang”’ presented a survey of algorithms for construction and manipulation of octrees. Samet” presented an O(n’) pixel-to-quadtree conversion algorithm which does a raster-order traversal of the image and maintains a valid quadtree throughout the construction procedure. Shaffer and Samet I2 then developed a better pixel-to-quadtree conversion algorithm which has a dominant time complexity proportional to the number of nodes in the quadtree (the total time complexity being O(n*), where r? is the number of pixels). This algorithm makes at most one insertion for each node in the quadtree at the extra expense of keeping track of the active nodes during the raster-order traversal.The Shaffer-Samet algorithm has recently been further improved by Holroyd and Mason4 by doing a Morton-order traversal.Clearly, all these algorithms can be extended to three dimensions, in which case the time taken will be O(n’). Franklin and Akman” have an algorithm to convert a set of parallelepipeds to an octree by stacking and splitting them. The time complexity of their algorithm is O(nlogn) for 71 parallelepipeds in a n x n x n universe. Yin et. ~1.‘~ presented a parallel algorithm for generating linear octree from voxel data. This algorithm does a straightforward bottom-up merging of sub-cubes into cubes using the Mortonorder traversal. The algorithm has a sequential time complexity of 0(n3) for a voxel set of size n3 and a parallel time complexity of O(n3/p) using p processors. This

0262-8856/94/05/0297-07 (_: 1994 Butterworth-Heinemann Image and Vision Computing Volume 12 Number 5 June

Ltd 1994

297

Efficient linear octree

Our algorithm has a time complexity of O(S), where S is the surface area of the object, for n3 voxels of an connected object. The time complexity is also O(S) for a general object with connected components if we know at least one voxel of each connected component. PRELIMINARIES To present concepts clearly, we first present the 2-dimensional pixel-to-quadtree conversion algorithm in the next section. We subsequently present the generalization to octrees. We assume that the input is a pixel array of size n x n where n = 2’, where each pixel is either black (full) or white (empty). The input is assumed to be connected. This condition can be relaxed to allow for disconnected components as long as we know the location of at least one pixel in each connected component. The seed-pixel for a connected component could be selected interactively for a visualization system, or can be found by exhaustive search. The output is a set of nodes which represents the linear quadtree of the input object. The linear quadtreet5 is a pointerless quadtree which only stores the black leaf nodes of a quadtree. The location of a node and all its ancestors are encoded into a single label. So the location, size and the path from the root to the node are implicitly encoded in the label. We now present the linear quadtree representation shown in Tin et a1.14,which we use in our al orithm. Assume that the input universe is a set of 2/i x 29 . The decomposition of the object into quadrants is shown in Figure 1. Each quadrant which belongs to the object is labelled (L, X, y) where L is the level of the node and x, y are the relative coordinates of the node. The root node is at level 0 and the label of root node is (O,O,O).Suppose the label of a quadrant is (L, x, y), and suppose this quadrant is not fully white or black. In this case, it is further subY

4 I I

(2.2.3)

(2.3.3)

(2.2.2)

(2.3.2)

(l,O,l)

f

(l,l,O)

(l,O,O)

--+ X Figure 1

298

Decomposition

method

Image and Vision

Computing

Volume

12 Number

divided into four sub-quadrants. Then the label of the sub-quadrant in the lower-left position is (L’, x’, y’), where L’ = L + 1; x’ = 2x; y’ = 2~. The labels of the are (L’, x’ + 1, y’), three sub-quadrants other (L’, x’, y’ + 1) and (L’, x’ + 1, y’ + 1). In general, for a quadrant with label (L, x, y), the label of the lower-left sub-quadrant at level L’ = L + 1, is (L’, x’, y’), where x’ = x ‘2’; y’ = y. 2’. The labels for its three sibling sub-quadrants can be derived as shown earlier. Note that in the linear quadtree, only the labels of quadrants which are black (full) are stored in (x, y) sorted order. Labels of empty and grey nodes are not stored. Pointers are not used to represent the tree, but that information is encoded in the label. An example of a simple object and its linear quadtree is shown in Figure 2. THE ALGORITHM Assume that the input is a 2-dimensional scene consisting of an n x n array of pixels, where n = 2k. The coordinates range 0,. . . , n - 1 in the X and Y directions. Each pixel is either black or white. Black pixels represent the object in the scene. The algorithm consists of two stages - surface extraction and merging: Input: 2D pixel data (n x n array) Output: Set {(L, x, y)} representing the linear quadtreee

For each y, determine the starts and ends of the runs of black pixels. A run is a contiguous set of black pixels in a row of the image. We will be using X-runs which are runs parallel to the X axis. At the end of this step, we will have a sorted set (possibly empty) of {(s;, ei)} runs of black pixels for every row of the image, where s;, ei = start, end coordinate of the ith run. We present the boundary tracking method for obtaining the runs: In this method, we can find the contour of the objects represented by the pixel data. The contour is a set of edges of the pixels. Only the edges parallel to the Y-axis need to be considered to determine the X-runs. The start and end of runs can be determined by sorting these edges first in Y and then in X. A bucket-sort can be used for this purpose. If the object is connected or if we know at-least one pixel of every connected component, then we can use the boundary-trafkiyg alg,orithm presented by Shu and Krueger 6 which 1s an improvement of the algorithm proposed by Gordon and Udupa”. This algorithm exploits the surface coherence of a connected object to compute the boundary in time O(S), where S is the size of the boundary. Loop for every level of the tree to perform merging: 2.1 For each consecutive pair of rows of the input image, find the overlaps of the runs of the two rows. The overlaps can be easily found by implementing a finite-state machine with four states - no runs, first row run, second row run and overlap (both row runs). Note that the pixels which are not in the overlap are the quadtree nodes at the current level. So the pixels which are not in the overlap are output as the linear quadtree nodes (L, x, y), where 5 June 1994

Efficient

linear octree

generation:

R Shu and M S Kankanhalli

Decomposition

Object

(l,l,l)

Linear Quadtree = { (1,0,0),(1,1,0),(2,0,2), (2,0,3),f2,1,2)]

2.2

L = level, .Yand _r are coordinates of the pixel. Now the overlapping runs can potentially be merged to form nodes at a higher level. However, not all pixels in the overlap can be passed to the higher level. Only the overlap starting on an even X-coordinate and ending on an even X-coordinate can be considered. So the overlaps have to be truncated and the overlapped pixels which are not in the truncated overlaps must be output as quadtree nodes at the current level of the tree. Divide the starts and ends of the truncated overlaps by two to be passed to the higher level. We also divide the Y-coordinate by two. So, conceptually, the resolution of the input is haived at the end of each loop iteration. At the end of this step, we have another set of {(s;, e,)} starts & ends, which are the input contiguous runs for the next iteration.

The looping routine below in pseudo-code:

described

above is presented

image

Figure 2

Linear quadtree example

size-2n; for (level=k;level>=O;level--1 I

size=size/ lim=

size-

2; 1;

for (j=O;j
merge(row~j),row(j+l)); /*merge

runsoftwo

rows

*/

I This merging process can be better understood by following through the entire construction shown in the example of Figure 3. We now summarize the merging algorithm below:

(if (ii) (iii)

and Vision

find the overlaps of the runs of the two rows for all rows output non-overlapping pixels as quadtree nodes at that level truncate overlaps such that the ends are even numbers Computing

Volume

72 Number

5 June

1994

299

Efficient linear octree generation:

R Shu and M S Kankanhalli

Runs

(1)

8

Nodes Output

_‘.‘.‘,‘.‘.‘_‘_‘_‘.‘.‘.~,‘,‘.‘.’ t

-

-

r

I

I

-

I

F,T.7

-

i:::::::]

‘.~.~.‘.1.‘.‘.‘.~.‘.~,~.‘,~.‘.‘.~ _..,,.____. .,.,.~.~.~.~.~.,.,.,.~.~.,.,.,., .~.~i:~.‘;.‘P.‘.I:.r:.i’. ,~,~.‘.~.‘.~,~,~.~.‘.‘.~.~.~. .‘.‘.-.‘l’.‘.‘.‘1’.‘,~.‘.t.‘.‘.‘.

.~.~.‘.‘.‘.‘.‘.~.~,~,~.‘.‘.~.‘.~. - * r.._.i_:.~:.~,,~,.~.~~.~~.‘~.‘~.’~.~.~..~’., ~~~~~~~~~‘~‘~~~~~‘~‘~‘~~~~~~~~~~, .~.‘.‘.‘.~.~.‘.‘.‘.‘.‘.‘. 1.‘~.~‘~,).,‘~.~.,.,‘,‘.‘.’.,.,. _~,‘.‘.‘.‘,‘.‘,‘.‘.,.‘.~.~. .~.~.‘.‘,‘,‘.‘.‘.‘.‘.‘.‘. .‘.‘.‘.‘,‘,‘,‘.‘.‘.‘.~.~. ;;;. . _ r’ ‘d ‘~,:~.:r t .:::j.:L...*.:r;.:.s.:.c:.~:. , _.__.,..___ll. ,‘,‘.‘.‘,‘_‘.‘_‘,‘,‘,~.~.~.~.‘.~. ‘_‘_... . I ,::::I .~.~.~.~.‘.‘_~.~_~.‘_‘.~.‘.~,~.~. _,_,... .,.,.,. _,_..~...,.,.,., ‘.‘.‘.‘.; L L < ;::: I I I ::::I I I .I.,‘_’ - r - r _ . _ r _ r - r .,.r,.+

6

4

4

2

0

6

8

(2) 4 .

.‘.‘_‘.‘_‘.‘.‘_

1

_‘,‘,‘,‘_‘_‘_‘_’

(y.:.: :.:.:.::.:::::_:.:.:.,~,,.,.,.~.~.,...,. . ... . . .. . i.:.:.:.~.:.:.:.‘.‘.‘.:.:.:.:.:.:~:.:.:.:.:.:.:.~. ,.I._._.,._.,.I. .~.~.,.,.,._._.~~,~,.,.,._._._.,_ 0

2

1

3

Cl*41

4

(3)

W1)

Figure 3

(iv) (v) (vi) (vii)

Example

of merging

process

o

1

pixels which are not in the truncated overlaps as nodes divide the bounds of the truncated overlaps by 2 to form new runs output this set of overlaps as input runs for the next level iteration go to step (i) if all levels have not been processed. output

The generalization of the 2-dimensional algorithm is presented now. The input for the 3-D case would be an

300

Image and Vision

Computing

Volume

12 Number

n x n x n array of voxels, where n = 2k. The output would be a set {(L, X, y, z)}, which represents the linear octree. The representation scheme is the same as described in the 2-D case with an extra ‘z’ coordinate. Instead of considering runs of pixels, we will consider runs of voxels which are parallelepipeds in the 2 direction with unit thickness in X and Y directions. Of course, the algorithm would also work with parallelepipeds in the A’ or Y directions. These runs can be determined by 3-dimensional surface tracking 5 June 1994

Efficient

algorithm16. Alternatively, ray-casting could also be used to determine these runsi3. The merging procedure for the 3D case is exactly like that for the two-dimensional case. Instead of merging two rows at a time, we now have to merge four Z parallelepipeds at a time. Otherwise, this is a direct extension of the two-dimensional procedure. Like in the 2-dimensional case, the only extra data structure required will be that needed for holding this sorted set of runs of black voxels {(si, ei)}, where si, ei = start, end coordinate of the ith run.We present the pseudo-code for the 3D voxel-to-linear octree conversion algorithm.The looping routine described above is presented below in pseudo-code: size for

= 2n; (level=k;

level--)

level>=O;

{ size=size/ lim=size-1; for {

(i=O; for

2; i
(j=O;

i=i+2)

j
j-j+21

{

merge(p(i,j),p(i,j+l), p(i+l,j),p(i+l,j+l)); /*merge4parallelepipedsp(,) */

linear octree generation:

R Shu and M S Kankanhalli

where S is the size of the boundary of the object. Now, in the second stage of the algorithm, these X-runs are merged and nodes are output. The number of times the merge procedure is executed is: WP,

= n/2 + n/4 +

. + n/210pzn

= O(n) The number of merge operations required to obtain the linear quadtree is O(S). This is because we merge the parallelepipeds at the boundary of the object. So the dominant time complexity for the second step will be O(S). On average, the size of the boundary of the object is O(n)‘. We can argue similarly for the threedimensional case and show that the average-case time complexity for voxel-to-octree conversion would be O(S), where S is the size of the surface of the object. Obviously, an adversary can create a worst-case input which would take 0(n2) for 2D and O(n”) for the 3D algorithm. A possible worst-case input to the algorithm is the checker-board pattern. But our algorithm performs very well for the average case.

REMARK

1 1 1 The 3D merging

(9

(ii) (iii) (iv) (v) (vi) (vii)

procedure

is now described

briefly:

find the overlaps of the runs of four parallelepipeds. Again, this step can be easily implemented using a finite-state machine with sixteen states - 1 state indicating no runs, 4 states indicating one parallelepiped run, 6 states indicating two parallelepiped overlap, 4 states indicating three parallelepiped overlap, and 1 state indicating all four parallelepipeds overlapping. Only the last case is output as runs to the next level after truncation output non-overlapping voxels as octree nodes at that level truncate overlaps of the runs such that the ends are even numbers output voxels which are not in the truncated overlaps as linear octree nodes at the current level divide the bounds of the truncated overlaps by 2 to form new runs output this set of runs as input parallelepiped runs for the next level of iteration go to step (i) if all levels have not been processed.

ANALYSIS

OF THE

ALGORITHM

The first step of the algorithm determines the contiguous X-runs of black pixels for each row. The two-dimensional version of surface tracking16, ” can be used to find the boundary of the object, and the X-runs could be easily determined from this information by sorting. If the object is simply connected or if there is a seed pixel for every connected component, then the time taken will be

Image

Since our algorithm is similar in form to the Franklin Akman13,we will indicate the differences between the two approaches. The reader is referred elsewhereI for details of the FranklinAkman algorithm.Both algorithms need parallelepipeds as input data, so both need n x n lists to store the original parallelepipeds. The final goal of both approaches is the same, i.e. to combine all voxels in the parallelepipeds into cubesas big as possible. But the basic ideas of the algorithms are different. The space and time complexities are quite different, too. The FranklinAkman algorithm is based on a number of combine/split operations on rows or squares to form octants. First, every parallelepiped is decomposed into a series of rows ‘maximum called components’ such that every component satisfies a divisibility property. This procedure for all n x n parallelepipeds will take a time & space complexity of O(n* logn). Then the rows in every list are either combined into squares or split by half into two sub-rows. Before doing this, a sorting must be performed. Assuming the use of a radix sort, the time complexity of this operation is O(n2 logn + tspiit/combine). After this step, the squares are similarly combined/split into cubes. Finally, the elements in all lists are output as octree nodes. The overall dominant time complexity is 0(n2 logn) or O(n log n) for rc parallelepipeds. In contrast, our algorithm uses only one extra data structure to hold all the parallelepipeds, including the input and combined parallelepipeds. There is no necessity of decomposing the input parallelepipeds into rows, and we do not use combining/splitting of rows (squares) into squares (cubes). The octant generation is a one-step procedure which outputs the octants based on the comparison of the start/end coordinates of every four neighbouring parallelepipeds. As we show in the paper, the time complexity is thus O(S), where S is the size of the surface of the object and the space complexity is

and Vision

Computing

Volume

12 Number

5 June

1994

301

Efficient linear octree generation:

Figure 4 Example human torso

R Shu and M S Kankanhalli

data - CT scan of

O(n’), which is less than that of the FranklinAkman algorithm. Another difference is that the output of our algorithm is a linear octree unlike a full octree as in the case of the Franklin-Akman algorithm.

in Table 1. The times shown in the last two columns indicate the total time required for the octree conversion without including the disk-read time. The data sets represent diverse objects, and we can see from the experiments that our algorithm takes roughly 60% of the time taken by the Yin algorithm.

IMPLEMENTATION AND RESULTS We implemented the three-dimensional algorithm on an IBM RS-6000/320 workstation and timed the algorithm for several voxel data sets. An example of the data set used is shown in Figure 4. We also implemented the Yin algorithmI and tested it for the data sets. The input voxel-data had an 8-bit density value for each voxel, but we thresholded the input and considered it to have only two density values - 0 and 1. The results are summarized Table 1 Timings algorithms

comparison

for

two

voxel-to-octree

conversion

Data set

Resolution

Yin algorithnt timing (s)

Our algorithm timing (s)

Cylhead CT Torso Molecule NerveCell 3dhead CTknee MRbrain CThead

256x256~ 113 256x256~ 113 128 x 128 x 113 256x256~ 113 256x256~ 113 512 x 512 x 113 256x256~ 113 256x256~ 113

37.69 14.23 34.18 13.62 122.30 173.94 123.46 154.74

26.22 9.65 20.64 8.07 87.03 127.90 89.51 101.62

302

Image and Vision

Computing

Volume

12 Number

CONCLUSIONS We have presented a new algorithm for constructing a linear octree from voxel data. This algorithm efficiently combines parallelepipeds to obtain the octree nodes. For connected objects, this algorithm performs significantly better than the existing algorithms. This algorithm has been used in an octree-based data compression utility for the medical visualization system developed for the National University of Singapore Hospital at ISS. The algorithm can also be easily parallelized so that it can be used efficiently on a parallel machine”.

ACKNOWLEDGEMENTS We thank Richard Krueger for some of the code used in the implementation. We are grateful to Kuan-Tsae Huang and Zuo-Ming Yin for many helpful discussions. We would also like to thank the anonymous

5 June 1994

Efficient linear octree generation:

referee whose insightful suggestions have significantly improved the presentation of the material in the paper.

8 9 IO

REFERENCES

I1

Westover, L ‘Interactive volume rendering’, Proc. Chapel Hill Volume Visualization Workshop (May 1989) pp9-16 Wilhelms, J and Van Gelder, A ‘Octrees for faster isosurface generation’, Compur. Graph., Vol 24 No 5 (November 1990) pp 57-6 I Tamminen, M and Samet, H ‘Efficient octree conversion by connectivity labeling’, Comput. Graph., Vol I8 No 3 (July 1984) pp 43351 Holroyd. F C and Mason, D C ‘Efficient linear quadtree construction algorithm’, Image & Vision Cornput.. Vol 8 No 3 (August 1990) pp 2 188224 Meagher, D ‘Geometric modeling using octree encoding’, Compu/. Graph. & Image Process., Vol 19 No 2 (June 1982) pp 1299147 Samet, H and Webber, R E ‘Hierarchical data structures for computer graphics. Part 1. Fundamentals’, IEEE Comput. Graph. & Applic., Vol 8 No 3 (May 1988) ~~48-68 Samet, H and Webber, R E ‘Hierarchical data structures for computer graphics. Part II. Applications’, IEEE Comput. Graph. & App/ic.. Vol 8 No 4 (July 1988) pp 59-75

12

13

14

15 16

17

18

Image and Vision

R Shu and M S Kankanhalli

Samet, H The Design and Analysis of Spatial Data Structures, Addison-Wesley, Reading, MA (1990) Samet, H Applications qf Spatial Data Struc~turcs, AddisonWesley, Reading, MA (1990) Chen, H H and Huang, T S ‘A survey of construction and manipulation of octrees’, Compur. Vision, Graph. & Image Process., Vol43 No 3 (September 1988) pp 40943 I Samet. H ‘An algorithm for converting rasters to quadtrees’, IEEE Trans. PAMI, Vol 3 No 1 (January 1981) ~~93-95 Shaffer, C A and Samet, H ‘Optimal quadtree construction algorithms’, Comput. Vision. Graph. & Imugc Process.. Vol 37 No 3 (March 1987) ~~402-419 Franklin. W R and Akman. V ‘Buildine. an octree from a set of parallelepipeds’, IEEE Comput. Graph.-& Applic.. Vol 5 No IO (October 1985) pp 58-64 Yin, 2 M, Huang, K T, Shu, Rand Tan. P H ‘A new optimal parallel algorithm for generating linear level octree from voxels’. Proc. Singupore Supercompu/ing Cottf:, Singapore (December 1990) Gargantim. I ‘An effective way to store quadtrees’, Commun. ACM, Vol 25 No 12 (December 1982) ~~905~.910 Shu, R and Krueger, R ‘A 3D surface construction algorithm for volume data’, Proc. Computer Gruphics In!. 1991 Cortfl, Boston, MA (June 1991) Gordon, D and Udupa, J ‘Fast surface tracking in threedimensional binary images’, Comput. Vision, Graph. & Image Process.. Vol 45 No 2 (February 1989) pp 196214 Shu, R and Kankanhalli, M S ‘Generating a linear octree from voxel data for a connected object’, Proc. SPIE VIth Int. Conf: on Medical fmuging. Newport Beach, CA (February 1992) pp 17-25

Computing

Volume

12 Number

5 June 1994

303