Microsoft Patent | User-guided surface reconstruction

Patent: User-guided surface reconstruction

Drawings: Click to check drawins

Publication Number: 20100085353

Publication Date: 20100408

Assignee: Microsoft Corporation

Abstract

Described is a technology by which a user interacts with a surface representative of a point cloud data to correct for imperfect scan data. The surface is reconstructed based on the interaction. Real time viewing of the image is facilitated by parallel surface reconstruction. For example, the user may draw strokes to reduce topological ambiguities in poorly-sampled areas. An algorithm automatically adds new oriented sample points to the original point cloud based on the user interaction. Then a new isosurface is generated for the augmented point cloud. The user also may specify the geometry of missing areas of the surface. The user copies a set of points from another point cloud, and places the points around the target area. A new isosurface is then generated.

Claims

1. In a computing environment, a method comprising, receiving user interaction with a surface representative of point cloud data, modifying the point cloud data based on the user interaction, and generating a new surface from the modified point cloud data by parallel operations on a graphics processing unit.

2. The method of claim 1 wherein receiving the user interaction comprises detecting a stroke.

3. The method of claim 1 wherein receiving the user interaction with the surface comprises detecting interaction at a location corresponding to one part of the surface, and wherein modifying the point cloud data comprises adding points corresponding to the location so as to be visible in the new surface.

4. The method of claim 1 wherein receiving the user interaction with the surface comprises detecting interaction at a location corresponding to one part of the surface, and wherein modifying the point cloud data comprises removing points corresponding to the location so as to not be visible in the new surface.

5. The method of claim 1 wherein modifying the point cloud data based on the user interaction comprises computing a depth range of a bounding box of a current view of a point set of point cloud data, and processing the interaction along the current view direction in the depth range to determine where to modify the point cloud data.

6. The method of claim 5 wherein the user interaction is a stroke, and wherein processing the interaction along the current view direction comprises determining a set of points distributed along the stroke, and for each point, computing a function value.

7. The method of claim 6 further comprising, removing a point based on its function value relative to a threshold value.

8. The method of claim 1 wherein generating the new surface comprises using Poisson equation for surface reconstruction.

9. In a computing environment having a graphics processing unit, a system comprising, a surface construction mechanism coupled to the graphics processing unit that communicates with the graphics processing unit to construct a surface based on processing a set of points, an interactive program that modifies the set of points at a location corresponding to interaction with the surface into a set of modified points, the surface construction mechanism communicating with the graphics processing unit to reconstruct a surface based on processing the set of modified points.

10. The system of claim 9 wherein the user interaction comprises a stroke, or a shape, or any combination of a stroke or a shape.

11. The system of claim 9 wherein the set of points is contained with nodes of an octree, in which the octree is built by computing nodes in parallel via the graphics processing unit.

12. The system of claim 11 further comprising, means for determining which node or set of nodes correspond to the modified points based on the location and a current depth of the surface.

13. The system of claim 9 wherein the interactive program modifies the set of points by adding points with so as to become visible in the constructed surface.

14. The system of claim 9 wherein the interactive program modifies the set of points by adding points copied from another surface.

15. The system of claim 9 wherein the interactive program modifies the set of points by removing points so as to not be visible in the constructed surface.

16. The system of claim 9 wherein the user interaction is a stroke, and further comprising, means for determining a set of points distributed along the stroke.

17. One or more computer-readable media having computer-executable instructions, which when executed perform steps, comprising, building an octree of nodes for a given set of oriented points, each node containing a set of points enclosed by the node and pointers to related nodes, constructing a surface based on the nodes, receiving user interaction that changes the points of at least one of the nodes, and reconstructing the surface based on the changes to the points.

18. The one or more computer-readable media of claim 17 having further computer-executable instructions comprising, determining which node or nodes in which to change point data based on a current view of the surface.

19. The one or more computer-readable media of claim 17 having further computer-executable instructions comprising, modifying the points to add visible information to the surface.

20. The one or more computer-readable media of claim 17 wherein receiving the user interaction that changes the points comprises, adding points copied from another surface, or removing points, or both adding points copied from another surface and removing points.

Description

BACKGROUND

[0001] Surface reconstruction from point clouds is actively researched in computer graphics. This reconstruction approach is widely used for fitting three-dimensional scanned data, filling holes on surfaces, and re-meshing existing surfaces.

[0002] Although a number of algorithms exist that are capable of producing high-quality surfaces, because of the computational complexity, each such algorithm operates as an offline process, e.g., cannot be used for real time graphics processing. In other words, none of the surface reconstruction algorithms can achieve interactive performance for editing and the like.

[0003] Moreover, such algorithms generally can only handle noise-free and uniformly-sampled point clouds. For noisy data, this method may fail to produce an acceptable surface. For example, with real-world scanned data, some areas of the surface may be under-sampled or completely missing. Automatic techniques fail to faithfully reconstruct the topology of the surface around these areas.

SUMMARY

[0004] This Summary is provided to introduce a selection of representative concepts in a simplified form that are further described below in the Detailed Description. This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used in any way that would limit the scope of the claimed subject matter.

[0005] Briefly, various aspects of the subject matter described herein are directed towards a technology by which a user interacts with a surface representative of a point cloud image to correct for imperfect scan data. The surface is reconstructed based on the interaction. Real time viewing of the image is facilitated by parallel surface reconstruction.

[0006] In one aspect, a mechanism allows the user to draw strokes to reduce topological ambiguities in poorly-sampled areas. An algorithm of the mechanism automatically adds new oriented sample points to the original point cloud based on the user interaction. Then a new isosurface is generated for the augmented point cloud.

[0007] In one aspect, the user may specify the geometry of missing areas of the surface. The user copies a set of points from another point cloud, and places the points around the target area. The new isosurface is then generated.

[0008] Other advantages may become apparent from the following detailed description when taken in conjunction with the drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

[0009] The present invention is illustrated by way of example and not limited in the accompanying figures in which like reference numerals indicate similar elements and in which:

[0010] FIG. 1 is a block diagram showing example components for octree building and parallel surface reconstruction.

[0011] FIGS. 2A-2C are representations of elements of data structures including the ordering of vertices and edges in a node (FIG. 2A), the ordering of a node's children and the ordering of nodes sharing a vertex, (FIG. 2B), and ordering of a node's neighboring nodes (FIG. 2C).

[0012] FIG. 3 is a flow diagram showing example steps for building a node array.

[0013] FIGS. 4A and 4B are representations of lookup tables (FIG. 4A) and the use of the lookup tables to compute neighboring nodes (FIG. 4B).

[0014] FIG. 5 is a flow diagram showing example steps for parallel surface reconstruction based upon user interaction.

[0015] FIGS. 6-9 are representations of images showing how a scanned image is reconstructed and modified based upon user interaction.

[0016] FIG. 10 shows an illustrative example of a computing environment into which various aspects of the present invention may be incorporated.

DETAILED DESCRIPTION

[0017] Various aspects of the technology described herein are generally directed towards a user-guided reconstruction mechanism (algorithm) to correct imperfect scan data, e.g., where areas of the surface are under-sampled or completely missing. The algorithm allows the user to draw strokes to reduce topological ambiguities in poorly-sampled areas. The algorithm automatically adds new oriented sample points to the original point cloud based on user interaction. Then a new isosurface is generated for the augmented point cloud.

[0018] The algorithm also allows the user to specify the geometry of missing areas of the surface. The user first copies a set of points from another point cloud and places the points around the target area. The new isosurface is then generated.

[0019] While some of the examples described herein are directed towards a parallel surface reconstruction algorithm implemented using NVIDIA Corporation's CUDA programming framework, e.g., to leverage a set of convenient parallel primitives such as scan, compact and sort, it is understood that these are only examples. Other frameworks may be used, along with graphics processors from other vendors, and other code and functions may be substituted for such primitives. As such, the present invention is not limited to any particular embodiments, aspects, concepts, structures, functionalities or examples described herein. Rather, any of the embodiments, aspects, concepts, structures, functionalities or examples described herein are non-limiting, and the present invention may be used various ways that provide benefits and advantages in computing and surfacing in general.

[0020] Turning to FIG. 1, there is shown a general block diagram exemplifying how a tree building mechanism 102 of a parallel surface construction mechanism (algorithm) 104 may work with an interactive application 106 (which may be an operating system component) and its point data 108 to dynamically render a scene 110. As described herein, one such application is directed to user guided surface reconstruction.

[0021] In general, the application 106 calls an API or the like (e.g., with supplied parameters) to work with the parallel surface construction mechanism (algorithm) 104, which uses the tree construction mechanism 102 to build data structures 112 (e.g., octrees and/or other suitable data structures) as needed to produce the data for displaying the surface. The tree construction mechanism 102 provides code and data for execution in the GPU 114, which executes the code to build the data structures 112.

[0022] In one aspect, there is provided fast octree construction of the GPU, including by building octrees in real-time via the GPU's parallelism, and then using level-order traversals. In contrast to conventional CPU octree builders which often construct trees by depth-first traversals, the GPU-based technique processes the octree nodes at the same tree level in parallel, one level at a time. With level-order traversals, the technique leverages the parallelism by spawning a new thread for every node at the same tree level. In one current implementation, 256 threads are used for each block, wherein the block number is computed by dividing the total number of parallel processes by the thread number per block. For example, in Step 2 (line 5) of Listing 1 (described below), the block number is N/256.

[0023] In another aspect, octrees are constructed that supply the information necessary for GPU surface reconstruction. In particular, the octree data structure provides fast access to tree nodes as well as the neighborhood information of each node (e.g., links to the neighbors of the node). While information of individual nodes is relatively easy to collect, computing the neighborhood information requires a large number of searches for every single node. Collecting neighborhood information for the nodes of the tree is thus extremely expensive, even on the GPU.

[0024] To address this problem, the algorithm takes advantage of the observation that a node's neighbors are determined by the relative position of the node with respect to its parent and its parent's neighbors. Based on this observation, two look up tables (LUT) are built, which record the indirect pointers to a node's relatives. Unlike direct pointers, indirect pointers are independent of specific instances of octrees and thus can be precomputed. At runtime, the actual pointers are quickly generated by querying the LUTs.

[0025] Further, based on octrees built in this manner, there is provided a GPU algorithm for the known Poisson surface reconstruction method; the Poisson method can reconstruct high quality surfaces through a global optimization. As part of the GPU algorithm, the algorithm derives an efficient procedure for evaluating the divergence vector in the Poisson equation and an adaptive marching cubes procedure for extracting isosurfaces from an implicit function defined over the volume spanned by an octree. As will be seen, these procedures are designed to exploit a modern GPU's fine-grained parallel architecture and use the octree neighborhood information. Note that GPU algorithms can also be readily designed for other implicit reconstruction methods by using the described octree construction technique and the adaptive marching cubes procedure for extracting isosurfaces on the GPU.

[0026] To build an octree O with maximum depth D from a given set of sample points Q={q.sub.i|i=1, . . . N}, the design of the octree data structure comprises four arrays, namely a vertex array, edge array, face array, and node array. The vertex, edge, and face arrays record the vertices, edges, and faces of the octree nodes, respectively. These arrays are relatively straightforward. In the vertex array, each vertex v records v.nodes, the pointers to the octree nodes that share vertex v. Following v.nodes, related elements may be reached, such as any edges sharing v. In the edge array, each edge records the pointers to its two vertices. Similarly in the face array each face records the pointers to its four edges.

[0027] FIG. 2A-2C show element ordering for quadtrees; the case with octrees is analogous. FIG. 2A shows the ordering of vertices (outside the box) and edges (inside the box) in a node. FIG. 2B shows the ordering of a node's children as well as the ordering of nodes sharing a vertex. FIG. 2C shows the ordering of a node's neighboring nodes.

[0028] The node array, which records the octree nodes, is more complex. Each node t in the node array NodeArray contains three pieces of information, namely (1) the shuffled xyz key named t.key, (2) the sample points contained in t, and (3) pointers to related data including its parent, children, neighbors, and other information as described below.

[0029] The shuffled xyz key is used because each octree node has eight children, whereby it is convenient to number a child node using a 3-bit code ranging from zero to seven. This 3-bit code encodes the sub-region covered by each child. Using the xyz convention, if the x bit is one, the child covers an octant that is "right in x"; otherwise the child covers an octant that is "left in x". The y and z bits are similarly set. The shuffled xyz key of a node at tree depth D is defined as the bit string x.sub.1y.sub.1z.sub.1x.sub.2y.sub.2z.sub.2 . . . x.sub.Dy.sub.Dz.sub.D, indicating the path from the root to this node in the octree. Therefore a shuffled xyz key at depth D has three D bits. One implementation uses 32 bits to represent the key, allowing a maximum tree depth of 10. The unused bits are set to zero.

[0030] Each octree node records the sample points enclosed by the node. The sample points are stored in a point array and sorted such that all points in the same node are contiguous. Therefore, for each node t, the point array only needs to store the number of points enclosed, t.pnum, and the index of the first point, t.pidx.

[0031] Connectivity pointers are stored for each node, that is, the pointers to the parent node, 8 child nodes, 27 neighboring nodes including itself, 8 vertices, 12 edges, and 6 faces are recorded. The pointers are represented as indices to the corresponding arrays. For example, t's parent node is NodeArray[t.parent] and t's first neighboring node is NodeArray[t.neighs[0]]. If the pointed element does not exist, the corresponding pointer is set to -1. Since each node has 27 neighbors at the same depth, the array t.neighs is of size 27. For consistent ordering of the related elements, these elements are ordered according to their shuffled xyz keys. For example, t's first child node t.children[0] has the smallest key among t's eight children and the last child t.children[7] has the largest key. For a vertex, its key value is defined as the sum of the keys of all nodes sharing the vertex. In this way, vertices can also be sorted; similarly, edges and faces can be sorted.

[0032] One example node array is built using a reverse level-order traversal of the octree, starting from the finest depth D and moving towards the root, one depth at a time. Listing 1 provides the pseudo code for the construction of NodeArrayD, the node array at depthD:

TABLE-US-00001 Listing 1: 1: // Step 1: compute bounding box 2: Compute Q`s the bounding box using Reduce primitive 3: // Step 2: compute shuffled xyz key and sorting code 4: code .rarw.new array 5: for each i = 0 to N - 1 in parallel 6: Compute key, q.sub.i`s shuffled xyz key at depth D 7: code[i] = key << 32 + i 8: // Step 3: sort all sample points 9: sortCode .rarw. new array 10: Sort(sortCode, code) 11: Generate the new point array according to sortCode 12: // Step 4: find the unique nodes 13: mark .rarw. new array 14: uniqueCode .rarw.new array 15: for each element i in sortcode in parallel 16: if sortCode[i].key .noteq. sortCode[i - 1].key then 17: mark[i] = true 18: else 19: mark[i] = false 20: Compact(uniqueCode, mark, sortCode) 21: Create uniqueNode according to uniqueCode 22: // Step 5: augment uniqueNode 23: nodeNums .rarw. new array 24: nodeAddress .rarw. new array 25: for each element i in uniqueNode in parallel 26: if element i - 1 and i share the same parent then 27: nodeNums[i] = 0 28: else 29: nodeNums[i] = 8 30: Scan(nodeAddress, nodeNums, +) 31: // Step 6: create NodeArray.sub.D 32: Create NodeArray.sub.D 33: for each element i in uniqueNode in parallel 34: t = uniqueNode[i] 35: address = nodeAddress[i] + t.x.sub.Dy.sub.Dz.sub.D 36: NodeArray.sub.D[address] = t

[0033] This construction process comprises six steps generally represented in FIG. 3. In a first step 302, the bounding box of the point set Q is computed. This is done by carrying out parallel reduction operations on the coordinates of all sample points. The Reduce primitive performs a scan on an input array and outputs the result of a binary associative operator, such as min or max, applied to all elements of the input array.

[0034] In step 304, the process computes the 32-bit shuffled xyz keys at depth D for all sample points, in parallel. Given a point p, its shuffled xyz key is computed in a top-down manner. The x bit at depth d, 1.ltoreq.d.ltoreq.D, is computed as:

x d = { 0 , if p x < C d x , 1 , otherwise , ##EQU00001##

where Cd is the centroid of the node that contains p at depth d-1. The y and z bits y.sub.d and z.sub.d are similarly computed. Unused bits are set to zero. The shuffled xyz key and the 32-bit point index are concatenated to a 64-bit code for the subsequent sorting operation.

[0035] At step 306, the sample points are sorted using a known sort primitive, which first performs a split-based radix sort per block and then a parallel merge sort of blocks. After sorting, points having the same key are contiguous in the sorted array. Then the index of each sample point in the original point array is computed by extracting the lower 32 bits of the point's code. The new point array is then constructed by copying the positions and normals from the original point array using the extracted indices.

[0036] At step 308, a node array is generated by removing duplicate keys in the sorted array, as follows. First, for each element of the sorted array, the element is marked as invalid if its key value equals that of its preceding element in the array. Then, a compact primitive is used to generate the node array which does not contain invalid elements. During this process, the relationship between the point array and the node array can be easily built. More particularly, for each element of the node array, the number of points contained by this node and the index of the first point in the point array are recorded.

[0037] At step 310, the node array obtained via step 308 is augmented to ensure that each node's seven siblings are also included, since each octree node has either eight or zero children. Note that in lines 25.about.29 of the pseudo code, each element in the node array is checked to see if it shares the same parent with the preceding element. This is done by comparing their keys. If the result is yes, nodeNums[i] is set to zero; otherwise it is set to eight. Then a parallel prefix sum/scan primitive is performed on the array nodeNums, and the result is stored in the array nodeAddress. Each element of nodeAddress thus holds the sum of all its preceding elements in nodeNums. In other words, nodeAddress contains the starting address of every unique node in the final node array.

[0038] Via step 312, the node array NodeArrayD is created. For each node that is added at step 310, only the key value is computed and the number of points contained is set to zero. For each node in uniqueNode, its corresponding element in NodeArrayD is located through nodeAddress and its 3-bit x.sub.Dy.sub.Dz.sub.D key, the node's data copied to this element. For each sample point in the point array, the pointer to the octree node that contains it is saved.

[0039] At other depths, the node array (e.g., at depth D-1) can be easily built from NodeArrayD. Recall that the eight siblings having the same parent are contiguous in NodeArrayD. For every eight sibling nodes in NodeArrayD, a parent node is generated by setting the last three bits of the keys of these nodes to zero. Again, the resulting parent nodes are augmented (as in step 310) to generate the final array NodeArrayD-1. At this point, each node in NodeArrayD can get the index of its parent node in NodeArrayD-1. For each node t in NodeArrayD-1, the pointers to its children are saved. The number of points t.pnum is computed as the sum of that of t's children, and the index of the first point t.pidx is set to be that of t's first child. The node arrays at other depths can be built the same way. The node arrays of all depths are then concatenated to form a single node array, NodeArray. Another array BaseAddressArray is also created, with each element of the array recording the index of the first node at each depth in NodeArray.

[0040] With respect to computing neighborhood Information, for each octree node in NodeArray, the process finds its neighboring octree nodes at the same depth. This neighborhood information used for computing the implicit function and running the marching cubes procedure (described below), and is also used for building the vertex, edge, and face arrays. Listing 2 shows the computing of neighboring nodes:

TABLE-US-00002 Listing 2: 1: for each node t at depth d in parallel 2: for each j = 0 to 26 in parallel 3: i .rarw. t`s 3-bit xyz key 4: p .rarw. NodeArray[t.parent] 5: if p.neighs[LUTparent[i][j]] .noteq. -1 then 6: h .rarw. NodeArray[p.neighs[LUTparent[i][j]]] 7: t.neighs[j] = h.children[LUTchild[i][j]] 8: else 9: t.neighs[j] = -1

[0041] Further provided is a lookup table (LUT)-based technique for computing the neighborhood information of every octree node efficiently, in parallel. Each node has up to 26 neighbors at the same depth, distributed among its sibling nodes and the child nodes of its parent's neighbors. One way to compute the neighbors is to enumerate all these candidate nodes, which requires 26.times.27.times.8=5616 searches for each node (26 neighbors, its parent and 26 neighbors of its parent, each neighbor having 8 children). An alternative approach provided herein is based on the observation is that a node's neighbors are determined by the relative position of the node with respect to its parent and its parent's neighbors; based on this observation, two lookup tables (LUTs) may be pre-computed to significantly speed up this neighborhood computation. Two such tables are shown in FIG. 4A, and are defined as a parent table 440, LUTparent, and a child table 442 LUTchild.

[0042] The parent table comprises a two-dimensional array providing the following information: for an octree node t whose parent is p, if t's index (or xyz key) in p.children is i, then the index of t.neighs[j]'s parent in p.neighs is LUTparent[i][j].

[0043] The child table LUTchild is a two-dimensional array with the following information: For the node t with parent p and index i in p.children as above, if node t's j-th neighbor t.neighs[j], whose parent node is h, the index of t.neigh[j] in h.children is LUTchild[i][j]. The size of both tables is 8.times.27.

[0044] For convenience a node is regarded as a neighbor of itself with index 13 in neighs. Note that two kinds of pointers are distinguished. Direct pointers are those represented as indices into one of the "global" arrays: the node, vertex, edge, and face arrays. For example, t.parent is a direct pointer. Indirect pointers are those represented as indices into one of the "local" arrays of a node: t.neighs, t.children, t.vertices, t.edges, and t.faces. The above-described tables record only indirect pointers, which are independent of specific instances of octrees and thus can be precomputed.

[0045] Listing 2 provides the pseudo-code for computing the neighboring nodes for each node t at depth d in parallel. First, it fetches t's parent p and its xyz key, which is t's index in p.children. To compute t's j-th neighbor t.neighs[j], this neighbor's parent node h is found by querying LUTparent, getting the neighbor using a second query to LUTchild; (this technique only needs 27 searches and is over two orders of magnitude faster that conventional approaches).

[0046] For clarity quadtrees are used in FIGS. 4A and 4B to illustrate Listing 2. The two tables for quadtrees, LUTparent and LUTchild, are of size 4.times.9 as shown in FIG. 4A. As shown in FIG. 4B, the quadtree node t's parent is p, and t's index in p.children is 0 (i.e., i=0). To compute t's 2-th neighbor (i.e., j=2), we first get p's 1-th neighbor, which is b, according to LUTparent[0][2].ident.1. Because LUTchild[0][2].ident.3, b's 3-th child, which is r, is the neighboring node desired. Therefore, t.neighs[2]=b.children[3]=r.

[0047] To compute t's 7-th neighbor (i.e., j=7), the process gets first get p's 4-th neighbor, which is p itself, according to LUTparent[0][7].ident.4. Because LUTchild[0][7].ident.1, p's 1-th child, which is s, is the node that is wanted. Therefore, t.neighs[7]=p.children[1]=s.

[0048] When computing a node's neighbors, its parent's neighbors are required. For this reason code corresponding to Listing 2 is preformed for all depths using a (forward) level-order traversal of the octree. If node t's j-th neighbor does not exist, t.neighs[j] is set as -1. For the root node, all its neighbors are -1, except its 13-th neighbor which is the root itself.

[0049] With respect to collecting information of vertices, edges, and faces of octree nodes by computing vertex, edge, and face arrays, for a vertex array, each octree node has eight corner vertices. Adding the eight vertices of every node into the vertex array introduces many duplications because a corner may be shared by up to eight nodes. One way to create a duplication-free vertex array is to sort all the candidate vertices by their keys and then remove duplicate keys, however, thus is inefficient due to the large number of nodes, e.g., on the order of 670 thousand nodes at depth 8, with the number of candidate vertices over five million.

[0050] A more efficient way to create the vertex array makes use of node neighbors computed as described above. Building the vertex array at octree depth d first finds, in parallel, a unique owner node for every corner vertex. The owner node of a corner is defined as the node that has the smallest shuffled xyz key among all nodes sharing the corner. Observing that all nodes that share corners with node t have to be t's neighbors allows quickly locating the owner of each corner from t's neighbors. Second, for each node t in parallel, all corner vertices whose owner is t itself are collected. The unique vertex array is then created. During this process, the vertex pointers t.vertices are saved. For each vertex v in the vertex array, the node pointers v.nodes are also appropriately set.

[0051] To build the vertex array of all octree nodes, the above process is performed at each depth independently, and the resulting vertex arrays are concatenated to form a single vertex array. Unlike the node array, this vertex array still has duplicate vertices between different depths. However, since this does not affect subsequent surface reconstruction, these duplicate vertices may be left as is.

[0052] With respect to other arrays, the edge and face arrays can be built in a similar way. For each edge/face of each node, its owner node is found. Then the unique edge/face array is created by collecting edges/faces from the owner nodes.

[0053] With respect to GPU surface reconstruction, to reconstruct surfaces from sample points using the octree constructed as above, the reconstruction generally comprises two steps. First, an implicit function .phi. over the volume spanned by the octree nodes is computed using Poisson surface reconstruction. Then, an adaptive marching cubes procedure extracts a watertight mesh as an isosurface of the implicit function. Note that, instead of Poisson surface reconstruction, other known methods for GPU surface reconstruction may be used. As mentioned above, however, the Poisson approach reconstructs high quality surfaces through a global optimization, and moreover, only requires solving a well-conditioned sparse linear system, which can be efficiently done on the GPU.

[0054] To this end, the following steps are performed on the GPU: [0055] 1. Build a linear system Lx=b, where L is the Laplacian matrix and b is the divergence vector. [0056] 2. Solve the above linear system using a multigrid solver. [0057] 3. Compute the isovalue as an average of the implicit function values at sample points [0058] 4. Extract the isosurface using marching cubes.

[0059] With respect to step (1) and computing the Laplacian Matrix L, as is known, the implicit function .phi. is a weighted linear combination of a set of blending functions {F.sub.o} with each function F.sub.o corresponding to a node of the octree. An entry of the Laplacian matrix L.sub.o,o'=(F.sub.o,.DELTA.F.sub.o') is the inner product of blending function F.sub.o and the Laplacian of F.sub.o. The blending function F.sub.o is given by a fixed basis function F:

F o ( q ) = F ( q - o c o w ) 1 o w 3 , ##EQU00002##

where o,c and o,w are the center and width of the octree node o. F is non-zero only inside the cube [-1, 1].sup.3. As is known, F is a separable function of x, y and z. As a result, the blending function F.sub.o is separable as well and can be expressed as:

F.sub.o(x, y, z)=f.sub.o.x,o.w(x)f.sub.o.y,o.w(y)f.sub.o.z,o.w(z)

[0060] Given the definition of Laplacian

.DELTA. F o ' = .differential. 2 F o ' .differential. x 2 + .differential. 2 F o ' .differential. y 2 + .differential. 2 F o ' .differential. z 2 , ##EQU00003##

the Laplacian matrix entry L.sub.o,o' can be computed as:

L o , o ' = F o , .differential. 2 F o ' .differential. x 2 + F o , .differential. 2 F o ' .differential. y 2 + F o , .differential. 2 F o ' .differential. z 2 = f o x , o w , f o ' x , o ' w '' f o y , o w , f o ' y , o ' w f o z , o w , f o ' z , o ' w + f o x , o w , f o ' x , o ' w f o y , o w , f o ' y , o ' w '' f o z , o w , f o ' z , o ' w + f o x , o w , f o ' x , o ' w f o y , o w , f o ' y , o ' w f o z , o w , f o ' z , o ' w '' . ##EQU00004##

[0061] The above inner products can be efficiently computed by looking up two pre-computed two-dimensional tables: one for f.sub.o, f.sub.o'and the other for f.sub.o, f''.sub.o'. These two tables are queried using the x-bits, y-bits, or z-bits of the shuffled xyz keys of node o and o'. This reduces the table size significantly. For a maximal octree depth 9, the table size is (2.sup.10-1).times.(2.sup.10-1). The table size may be further reduced because the entries of the tables are symmetric.

[0062] For evaluating the divergence vector b, the divergence coefficients b.sub.o can be computed as:

b.sub.o=.SIGMA..sub.o'.epsilon..sub.D{right arrow over (v)}.sub.o'{right arrow over (u)}.sub.o,o',

where {right arrow over (u)}.sub.o,o'=F.sub.o(q), .gradient.F.sub.o' .sup.D is the set of all octree nodes at depth D. The inner productF.sub.o(q), .gradient.F.sub.o'can be efficiently computed using a precomputed lookup table for f.sub.o, f.sub.o'as in the computation of L.sub.o,o'. As for {right arrow over (v)}.sub.o', it is computed as:

{right arrow over (v)}.sub.o'=.SIGMA..sub.q.sub.i.sub..epsilon.Q.alpha..sub.o',q.sub.i{righ- t arrow over (n)}.sub.i, (1)

where .alpha..sub.o,qi is the weight by which each sampling point q.sub.i distributes the normal {right arrow over (n)}.sub.i to its eight closest octree nodes at depth-D. Listing 3 provides the pseudocode for computing the divergence vector b.

TABLE-US-00003 Listing 3 Compute Divergence Vector b: 1: // Step 1: compute vector field 2: for each node o at depth D in parallel 3: = 0 4: for j = 0 to 26 5: t .rarw. NodeArray[o.neighs[j]] 6: for k = 0 to t.pnum 7: i = t.pidx + k 8: + = 9: // Step 2: compute divergence for finer depth nodes 10: for d = D to 5 11: for each node o at depth d in parallel 12: b.sub.o = 0 13: for j = 0 to 26 14: t .rarw. NodeArray[o.neighs[j]] 15: for k = 0 to t.dnum 16: idx = t.didx + k 17: o` .rarw. NodeArray[idx] 18: b.sub.o + = 19: // Step 3: compute divergence for coarser depth nodes 20: for d = 4 to 0 21: divg .rarw.new array 22: for node o at depth d 23: for each depth-D node o` covered by all nodes in o.neighs in parallel 24: divg[i] = 25: b.sub.o = Reduce(divg, +)

[0063] This computation takes three steps. In a first step, the vector field {right arrow over (v)}.sub.o' is computed for each octree node o' according to Equation (1). Since Equation (1) essentially distributes sample point q.sub.i's normal {right arrow over (n)}.sub.i to its eight nearest octree nodes at depth D, vector {right arrow over (v)}.sub.o' is only affected by the sample points that are contained in either node o' or its twenty-six neighbors. The pointers to the node neighbors (recorded as described above) are used to locate these neighbors.

[0064] In a second step, the divergence at every finer depth, which is defined as any depth greater than four, is computed in parallel for all nodes, as shown in Step 2 of Listing 3. One way to accumulate b.sub.o for each octree node o is to iterate through all nodes o' at depth D. However, this costly full iteration is actually not necessary` because the basis function F's domain of support is the cube [-1, 1].sup.3, {right arrow over (u)}.sub.o,o' equals zero for a large number node pairs (o, o'). Specifically, for node o, only the depth-D nodes whose ancestors are either o or o's neighbors have nonzero {right arrow over (u)}.sub.o,o'. These nodes can be located by iterating over o's neighbors. Note that t.dnum and t.didx are the number of depth-D nodes covered by t and the pointer to t's first depth-D node respectively. These data can be easily obtained and recorded during tree construction.

[0065] In a third step, the divergence at every coarser depth, which in one implementation is defined as any depth not greater than four, is computed. For nodes at a coarser depth, the approach taken in the second step is not appropriate because it cannot exploit the fine-grained parallelism of GPUs.

[0066] The node number at coarser depths is much smaller than that at finer depths, and the divergence of a node at a coarser depth may be affected by many depth-D nodes. For example, at depth zero, there is only one root node and all depth-D nodes contribute to its divergence. To maximize parallelism, the computation is parallelized over all covered depth-D nodes for nodes at coarser depths. As shown in step 3 of Listing 3, the divergence contribution is first computed for each depth-D node in parallel, followed by a reduction operation to sum up the contributions.

[0067] The GPU multigrid solver is rather straightforward. For each depth d from coarse to fine, the linear system L.sup.dx.sup.d=b.sup.d is solved using a conjugate gradient solver for sparse matrices as is known. L.sup.d contains as many as 27 nonzero entries in a row. For each row, the values and column indices of nonzero entries are stored in a fixed-sized array. The number of the nonzero entries is also recorded.

[0068] Note that the divergence coefficients at depth d need to be updated using solutions at coarser depths. For the blending function F.sub.o of an arbitrary octree node o, only the blending functions of o's ancestors and their 26 neighbors may overlap with F.sub.o. Therefore, these nodes may be visited through the pointers stored in parent and neighs fields of node o.

[0069] To evaluate the implicit function value at an arbitrary point q in the volume, the octree is traversed. Listing 4 shows the pseudo code of a depth-first traversal for this purpose:

TABLE-US-00004 Listing 4: compute implicit function value .phi..sub.q for point q: 1: .phi..sub.q = 0 2: nodestack .rarw.new stack 3: nodestack.push(proot) 4: while nodestack is not empty 5: o .rarw. NodeArray[nodestack.pop( )] 6: .phi..sub.q + = F.sub.o(q).phi..sub.o 7: for i = 0 to 7 8: t .rarw. NodeArray[o.children[i]] 9: if q.x-t.x < t.w and q.y-t.y < t.w and q.z-t.z < t.w then 10: nodestack.push(o.children[i])

[0070] A stack is used to store the pointers to the nodes to be traversed. For this traversal, a stack size of 8D is enough for octrees with a maximal depth D.

[0071] Note that the implicit function value of a sample point q.sub.i can be evaluated in a more efficient way, because the process already knows the depth-D node o where q.sub.i is located. In other words, only octree nodes whose blending function may overlap with that of o need to be traversed. These nodes include o itself, o's neighbors, o's ancestors, and the neighbors of o's ancestors. With the implicit function values at all sample points, the isovalue is computed as an average:

.phi.=.SIGMA..sub.i.phi.(q.sub.i)/N

[0072] A marching cubes technique on the leaf nodes of the octree is used to extract the isosurface. The output is a vertex array and a triangle array which can be rendered directly.

[0073] The depth-D nodes are processed in five steps. First, the implicit function values are computed for all octree vertices in parallel. As in the case with the sample points, each vertex v's implicit function value can be efficiently computed by traversing only the related nodes, which can be located through the pointers stored in v.nodes. Second, the number of output vertices is computed with a single pass over the octree edges and the output address is computed by performing a scan operation. Third, each node's cube category is calculated and the number and addresses of output triangles are computed.

[0074] In Step 4 and 5 the vertices and triangles are generated and saved. During this process, for each face of each node, if one of its four edges has a surface-edge intersection, the face is deemed to contain surface-edge intersections, and the face is marked. This information is propagated to the node's ancestors. For leaf nodes at other depths, nodes that do not produce triangles in parallel are filtered out. For each node, if the implicit function values at its eight corners have the same sign and none of its six faces contain surface-edge intersections, the node does not need any further processing. Otherwise, the node is subdivided to depth D. The depth-D nodes generated by this subdivision are collected to build the new node, vertex and edge arrays.

[0075] Then, code corresponding to Listing 5 is performed to generate vertices and triangles:

TABLE-US-00005 Listing 5 Marching Cubes 1: // Step 1: compute implicit function values for octree vertices 2: vvalue .rarw.new array 3: for each octree vertex i at depth-D in parallel 4: Compute the implicit function value vvalue[i] 5: vvalue[i] - = .phi. 6: // Step 2: compute vertex number and address 7: vexNums .rarw.new array 8: vexAddress .rarw.new array 9: for each edge i at depth-D in parallel 10: if the values of i`s two vertices have different sign then 11: vexNums[i] = 1 12: else 13: vexNums[i] = 0 14: Scan(vexAddress, vexNums, +) 15: // Step 3: compute triangle number and address 16: triNums .rarw.new array 17: triAddress .rarw.new array 18: for each node i at depth-D in parallel 19: Compute the cube category based the values of i`s vertices 20: Compute triNums[i] according to the cube category 21: Scan(triAddress, triNums, +) 22: // Step 4: generate vertices 23: Create VertexBuffer according to vexAddress 24: for each edge i at depth-D in parallel 25: if vexNums[i] == 1 then 26: Compute the surface-edge intersection point q 27: VertexBuffer[vexAddress[i]] = q 28: // Step 5: generate triangles 29: Create TriangleBuffer according to triAddress 30: for each node i at depth-D in parallel 31: Generate triangles based on the cube category 32: Save triangles to TriangleBuffer[triAddress[i]]

[0076] This procedure is carried out iteratively until no new triangles are produced. Note that in each iteration, there is no need to handle the nodes subdivided in previous iterations. To remove duplicate surface vertices and merge vertices located closely to each other, the shuffled xyz key is computed for each vertex, with the keys used to sort all vertices. Vertices having the same key values are merged by performing a parallel compact operation. The elements in the triangle array are updated accordingly and all degenerated triangles are removed. Each triangle's normal is also computed.

[0077] The GPU-based surface reconstruction can be employed in new or existing applications. As one example, a user-guided reconstruction algorithm is provided to compensate for imperfect scan data where many areas of the surface are either under-sampled or completely missing. In this algorithm, a user draws strokes (or otherwise interacts) around poorly-sampled areas to reduce topological ambiguities. Benefiting from the high performance of GPU reconstruction, the user can view the reconstructed mesh in fractions of a second after drawing a stroke, unlike other applications that take on the order of several minutes to update the reconstructed mesh.

[0078] To this end, new oriented sample points are added to the original point cloud based on user interaction. Then a new isosurface is generated for the augmented point cloud. For example, suppose Q is the original point set and Q' is the current point set after each user interaction; after the user draws a stroke, the following steps are taken to generate the new surface, as generally represented in FIG. 5:

[0079] Step 502 computes the depth range of Q's bounding box under the current view. Step 504 iteratively extrudes the stroke along the current view direction in the depth range, with a user-specified interval w. For each extruded stroke, a set of points are uniformly distributed along the stroke, also with interval w. This point set is denoted as S.

[0080] For the points in S, step 506 computes their implicit function values in parallel using the procedure in Listing 4. Step 508 remove points from S whose implicit function values are not less than the current isovalue threshold .phi.. Step 510 compute normals for all points in S, after which step 512 adds S to the current point set, forming Q'.

[0081] Step 514 perform GPU reconstruction with Q' as its input, and generates the new isosurface. In Step 504, the interval w is set to be the width of an octree node at depth D by default. Note that step 508 removes points outside of the current reconstructed surface because the process only wants to add new points in inner regions, where topological ambiguity is found.

[0082] Further, note that unwanted points may be accidentally introduced in some inner regions. When this happens, the user can remove those points interactively. In Step 514, the new isovalue is computed as the average of the implicit function values of points in the original point set Q to restrict the influence of newly-added points to local areas. The new points are only used to change the local vector field.

[0083] One current system provides two ways to compute the normals for points in S in Step 510. One is based on normal interpolation; for each point s.sub.i.epsilon. S, the octree of Q' is traversed to find all points of Q' that are enclosed by a box centered at s.sub.i. Then s.sub.i's normal is computed as an interpolation of the normals of these points.

[0084] The interpolation weight of a point q' is proportional to the reciprocal of the squared distance between q' and s.sub.i. The box size is a userspecified parameter. If no point is found given the current box size, the algorithm automatically increases the box size and traverses the octree again.

[0085] The other scheme for computing the normals is relatively straightforward. The normals are restricted to be orthogonal to both the current viewing direction and the tangents of the stroke. On implementation lets the normals point to the right side of the stroke. Note that for the first normal computation scheme, the user's interaction is not limited to drawing strokes. Users may draw a rectangle or any closed shape to define an area where they want to insert new points. This shape is then extruded and a set of points is uniformly distributed inside the extruded shape.

[0086] FIGS. 6-9 show how user-guided reconstructions works with a scanned elephant model. FIG. 6 shows the input scan, while FIG. 7 shows a result from automatic reconstruction. As can be seen, the elephant's head and trunk are mistakenly connected.

[0087] To correct this, the user interactively draws a stroke 770 to indicate to the interactive application that there is a mistake at that location. The surface is then redrawn based on the stroke along the current view direction in the depth range, with the user-specified interval w.

[0088] FIG. 8 shows the improved surface after the user draws the stroke 770 shown in FIG. 7. As can be seen by the area 880, the head and trunk are now separated, reducing the topological ambiguity from the scan.

[0089] In another example, a tail 990 copied from another surface is added to the elephant via user interaction. A new elephant surface with the new tail is immediately reconstructed for the user to view.

Exemplary Operating Environment

[0090] FIG. 10 illustrates an example of a suitable computing and networking environment 1000 on which the examples of FIGS. 1-9 may be implemented. The computing system environment 1000 is only one example of a suitable computing environment and is not intended to suggest any limitation as to the scope of use or functionality of the invention. Neither should the computing environment 1000 be interpreted as having any dependency or requirement relating to any one or combination of components illustrated in the exemplary operating environment 1000.

[0091] The invention is operational with numerous other general purpose or special purpose computing system environments or configurations. Examples of well known computing systems, environments, and/or configurations that may be suitable for use with the invention include, but are not limited to: personal computers, server computers, hand-held or laptop devices, tablet devices, multiprocessor systems, microprocessor-based systems, set top boxes, programmable consumer electronics, network PCs, minicomputers, mainframe computers, distributed computing environments that include any of the above systems or devices, and the like.

[0092] The invention may be described in the general context of computer-executable instructions, such as program modules, being executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, and so forth, which perform particular tasks or implement particular abstract data types. The invention may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in local and/or remote computer storage media including memory storage devices.

[0093] With reference to FIG. 10, an exemplary system for implementing various aspects of the invention may include a general purpose computing device in the form of a computer 1010. Components of the computer 1010 may include, but are not limited to, a processing unit 1020, a system memory 1030, and a system bus 1021 that couples various system components including the system memory to the processing unit 1020. The system bus 1021 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures. By way of example, and not limitation, such architectures include Industry Standard Architecture (ISA) bus, Micro Channel Architecture (MCA) bus, Enhanced ISA (EISA) bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnect (PCI) bus also known as Mezzanine bus.

[0094] The computer 1010 typically includes a variety of computer-readable media. Computer-readable media can be any available media that can be accessed by the computer 1010 and includes both volatile and nonvolatile media, and removable and non-removable media. By way of example, and not limitation, computer-readable media may comprise computer storage media and communication media. Computer storage media includes volatile and nonvolatile, removable and non-removable media implemented in any method or technology for storage of information such as computer-readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can accessed by the computer 1010. Communication media typically embodies computer-readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism and includes any information delivery media. The term "modulated data signal" means a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media includes wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, RF, infrared and other wireless media. Combinations of the any of the above may also be included within the scope of computer-readable media.

[0095] The system memory 1030 includes computer storage media in the form of volatile and/or nonvolatile memory such as read only memory (ROM) 1031 and random access memory (RAM) 1032. A basic input/output system 1033 (BIOS), containing the basic routines that help to transfer information between elements within computer 1010, such as during start-up, is typically stored in ROM 1031. RAM 1032 typically contains data and/or program modules that are immediately accessible to and/or presently being operated on by processing unit 1020. By way of example, and not limitation, FIG. 10 illustrates operating system 1034, application programs 1035, other program modules 1036 and program data 1037.

[0096] The computer 1010 may also include other removable/non-removable, volatile/nonvolatile computer storage media. By way of example only, FIG. 10 illustrates a hard disk drive 1041 that reads from or writes to non-removable, nonvolatile magnetic media, a magnetic disk drive 1051 that reads from or writes to a removable, nonvolatile magnetic disk 1052, and an optical disk drive 1055 that reads from or writes to a removable, nonvolatile optical disk 1056 such as a CD ROM or other optical media. Other removable/non-removable, volatile/nonvolatile computer storage media that can be used in the exemplary operating environment include, but are not limited to, magnetic tape cassettes, flash memory cards, digital versatile disks, digital video tape, solid state RAM, solid state ROM, and the like. The hard disk drive 1041 is typically connected to the system bus 1021 through a non-removable memory interface such as interface 1040, and magnetic disk drive 1051 and optical disk drive 1055 are typically connected to the system bus 1021 by a removable memory interface, such as interface 1050.

[0097] The drives and their associated computer storage media, described above and illustrated in FIG. 10, provide storage of computer-readable instructions, data structures, program modules and other data for the computer 1010. In FIG. 10, for example, hard disk drive 1041 is illustrated as storing operating system 1044, application programs 1045, other program modules 1046 and program data 1047. Note that these components can either be the same as or different from operating system 1034, application programs 1035, other program modules 1036, and program data 1037. Operating system 1044, application programs 1045, other program modules 1046, and program data 1047 are given different numbers herein to illustrate that, at a minimum, they are different copies. A user may enter commands and information into the computer 1010 through input devices such as a tablet, or electronic digitizer, 1064, a microphone 1063, a keyboard 1062 and pointing device 1061, commonly referred to as mouse, trackball or touch pad. Other input devices not shown in FIG. 10 may include a joystick, game pad, satellite dish, scanner, or the like. These and other input devices are often connected to the processing unit 1020 through a user input interface 1060 that is coupled to the system bus, but may be connected by other interface and bus structures, such as a parallel port, game port or a universal serial bus (USB). A monitor 1091 or other type of display device is also connected to the system bus 1021 via an interface, such as a video interface 1090. The monitor 1091 may also be integrated with a touch-screen panel or the like. Note that the monitor and/or touch screen panel can be physically coupled to a housing in which the computing device 1010 is incorporated, such as in a tablet-type personal computer. In addition, computers such as the computing device 1010 may also include other peripheral output devices such as speakers 1095 and printer 1096, which may be connected through an output peripheral interface 1094 or the like.

[0098] The computer 1010 may operate in a networked environment using logical connections to one or more remote computers, such as a remote computer 1080. The remote computer 1080 may be a personal computer, a server, a router, a network PC, a peer device or other common network node, and typically includes many or all of the elements described above relative to the computer 1010, although only a memory storage device 1081 has been illustrated in FIG. 10. The logical connections depicted in FIG. 10 include one or more local area networks (LAN) 1071 and one or more wide area networks (WAN) 1073, but may also include other networks. Such networking environments are commonplace in offices, enterprise-wide computer networks, intranets and the Internet.

[0099] When used in a LAN networking environment, the computer 1010 is connected to the LAN 1071 through a network interface or adapter 1070. When used in a WAN networking environment, the computer 1010 typically includes a modem 1072 or other means for establishing communications over the WAN 1073, such as the Internet. The modem 1072, which may be internal or external, may be connected to the system bus 1021 via the user input interface 1060 or other appropriate mechanism. A wireless networking component 1074 such as comprising an interface and antenna may be coupled through a suitable device such as an access point or peer computer to a WAN or LAN. In a networked environment, program modules depicted relative to the computer 1010, or portions thereof, may be stored in the remote memory storage device. By way of example, and not limitation, FIG. 10 illustrates remote application programs 1085 as residing on memory device 1081. It may be appreciated that the network connections shown are exemplary and other means of establishing a communications link between the computers may be used.

[0100] An auxiliary subsystem 1099 (e.g., for auxiliary display of content) may be connected via the user interface 1060 to allow data such as program content, system status and event notifications to be provided to the user, even if the main portions of the computer system are in a low power state. The auxiliary subsystem 1099 may be connected to the modem 1072 and/or network interface 1070 to allow communication between these systems while the main processing unit 1020 is in a low power state.

CONCLUSION

[0101] While the invention is susceptible to various modifications and alternative constructions, certain illustrated embodiments thereof are shown in the drawings and have been described above in detail. It should be understood, however, that there is no intention to limit the invention to the specific forms disclosed, but on the contrary, the intention is to cover all modifications, alternative constructions, and equivalents falling within the spirit and scope of the invention.