LINEBURG

If the data contain noise, we can use a simple postprocessing for the

implicit surface. There are many ways to view this process, derived by

Whitaker [172], but perhaps the most relevant here is based on TV denois-

13.4 Numerical Implementation

ing of images described in Chapter 11. Consider φ0 , the level set function

whose zero isocontour is the surface we wish to smooth. Then we let

There are three key numerical ingredients in our implicit surface recon-

u0 = H(φ0 ) (H is the Heaviside function) be the noisy image, which we

struction. First, we need a fast algorithm to compute the distance function

input into the TV denoising algorithm. Then we minimize

to an arbitrary data set on rectangular grids. Second, we need to ¬nd a good

1 initial surface for our gradient ¬‚ow. Third, we have to solve time-dependent

(H(φ) ’ H(φ0 ))2 dx,

µT V (H(φ)) +

2 partial di¬erential equations for the level set function.

We can use an arbitrary initial surface that contains the data set such as

where µ > 0 is the regularization parameter that balances between ¬delity

a rectangular bounding box, since we do not have to assume any a priori

and regularization. The variational level set method of Zhao et al. [175]

knowledge for the topology of the reconstructed surface. However, a good

gives

initial surface is important for the e¬ciency of our method. We start from

φt = [µκ ’ (H(φ) ’ H(φ0 ))]| φ|, (13.5) any initial exterior region that is a subset of the true exterior region. All

grid points that are not in the initial exterior region are labeled as interior

and we take φ0 as the initial guess.

points. Those interior grid points that have at least one exterior neighbor

13.4. Numerical Implementation 143 144 13. Reconstruction of Surfaces from Unorganized Data Points

are labeled as temporary boundary points. Then we use the following pro- Figure 13.1 shows the reconstruction of a torus with missing data. The

cedure to march the temporary boundary inward toward the data set. We hole is ¬lled nicely with a patch of minimal surface. Figure 13.2 shows the

put all the temporary boundary points in a heap-sort binary tree structure, reconstruction of a rat brain from MRI data, which is both noisy and highly

sorting according to distance values. Take the temporary boundary point nonuniform between slices. Next we show the reconstruction of a dragon on

that has the largest distance (on the heap top) and check to see whether a 300—212—136 grid using high-resolution data in Figure 13.3(a) and much

it has an interior neighbor that has a larger or equal distance value. If it lower resolution data in ¬gure 13.3(b). Figure 13.4 shows the reconstruction

does not have such an interior neighbor, turn this temporary boundary of a statuette of the Buddha on two di¬erent grids using the same data set

point into an exterior point, take this point out of the heap, add all this composed of 543,652 points.

point™s interior neighbors into the heap, and re-sort according to distance Other extensions are possible. Suppose we are given values of the normal

values. If it does have such an interior neighbor, we turn this temporary to the surface at the same or di¬erent set S of points. The ¬rst step,

boundary point into a ¬nal boundary point, take it out of the heap, and analogous to the fast computation of unsigned distance, is to construct

re-sort the heap. None of its neighbors are added to the heap. We repeat a unit vector de¬ned throughout the grid that interpolates this set. One

this procedure on the temporary boundary points until the the maximum possibility involves the construction of a harmonic map, which is easier

distance of the temporary boundary points is smaller than some tolerance, than in sounds using any of the techniques developed by Vese and Osher

e.g., the size of a grid cell, which means that all the temporary boundary [170], Alouges [4], E and Wang [57], or Tang et al. [162]. Given this unit

points in the heap are close enough to the data set. Finally, we turn these vector N (x) we add to our energy E(“) another quantity cE (“), where

temporary boundary points into the ¬nal set of boundary points, and our c > 0 is a constant whose dimension is length,

tagging procedure is ¬nished. Since we visit each interior grid point at most 1

p

φ p

once, the procedure will be completed in no more than O(N log N ) oper- 1’N ·

E (“) = ds . (13.7)

| φ|

ations, where log N comes from the heap-sort algorithm. Moreover, since “

the maximum distance for the boundary heap is strictly decreasing, we can Again using our variational level set calculus, we see that the gradient

prove that those interior points that have a distance no smaller than the descent evolution associated with equation (13.7) is

maximum distance of the temporary boundary heap at any time will re-

φ

main as interior points; i.e., there is a nonempty interior region when the φt = | φ| · ’ ·n (13.8)

| φ|

tagging algorithm is ¬nished. We can also show that at least one of the

¬nal boundary points is within the tolerance distance to the data set. for p = 1. See Burchard et al. [21].

Starting from an arbitrary exterior region that is a subset of the ¬nal

exterior region, the furthest point on the temporary boundary is tangent to

a distance contour and does not have an interior point that is farther away.

The furthest point will be tagged as an exterior point, and the boundary

will move inward at that point. Now another point on the temporary bound-

ary becomes the furthest point, and hence the whole temporary boundary

moves inward. After a while the temporary boundary is close to a distance

contour and moves closer and closer to the data set, following the distance

contours until the distance contours begin to break into spheres around

data points. The temporary boundary point at the breaking point of the

distance contour, which is equally distant from distinct data points, will

have neighboring interior points that have a larger distance. So this tem-

porary boundary point will be tagged as a ¬nal boundary point by our

procedure, and the temporary boundary will stop moving inward at this

breaking point. The temporary boundary starts deviating from the distance

contours and continues moving closer to the data set until all temporary

boundary points either have been tagged as ¬nal boundary points or are

close to the data points. The ¬nal boundary is approximately a polyhedron

with vertices belonging to the data set.

13.4. Numerical Implementation 145 146 13. Reconstruction of Surfaces from Unorganized Data Points

(a) data points (b) reconstruction (c) reconstruction

Figure 13.1. Hole-¬lling for a torus.

(a) data points (b) initial guess (c) ¬nal reconstruction

Figure 13.2. Reconstruction of a rat brain from data of MRI slices.

(a) 146x350x146 grid (b) 63x150x64 grid

Figure 13.4. Reconstruction of a “Happy Buddha” from 543,652 data points on

di¬erent grid resolutions.

(a) 437,645 points (b) 100,250 points

Figure 13.3. Reconstruction of a dragon using data sets of di¬erent resolution on

a 300 — 212 — 136 grid.

95 95

90

90

85

85

80

80

75

75

y

y

70

70

65

65

60

60

55 55

30 35 40 45 50 55 60 65 70

30 35 40 45 50 55 60 65 70

x x

Plate 1. (Figure 9.5). Initial Plate 2. (Figure 9.5).

placement of both types of Particle positions after the

particles on both sides of the initial attraction step is used

interface. to place them on the appro-

priate side of the interface.

Plate 4 (Figure 12.8). Recovered objects without well-

defined boundaries, using the multi-channel version of the

active contour model without edges from [32]. The three

Plate 3 (Figure 12.7). Color image, its gray-level version, objects could not be recovered using only one channel or

and the three RGB channels. the intensity image.

log(den) vel

0

3

500

2.5

2

1000

1.5

1

1500

0.5

0

2000

“ 0.5

0 2 4 6 8 10 0 2 4 6 8 10

entropy press

5 7

x 10 x 10

3

10

2.5

8

2

6

1.5

4 1

2 0.5

0

0

0 2 4 6 8 10 0 2 4 6 8 10

Plate 6 (Figure 15.9). The gamma-law gas is depicted in

red, while the stiff Tait equation of state water is depicted

in green. Note that the log of the density is shown, since

the density ratio is approximately 1000 to 1. This calcula-

tion uses only 100 grid cells.

Plate 5 (Figure 12.12). Color picture with junctions. Three

level set functions representing up to eight regions. Six seg-

ments are detected. We show the final zero level sets [36].

log(den) vel log(den) vel

500

3

3 0

2.5

400

2.5

500

2

300

1.5 2

1000

1 200

1.5

1500

0.5

100

0 1

2000

0

“ 0.5

0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 8 10

entropy press

entropy press

5 7 4 6

x 10 x 10 x 10 x 10

12 6

10

10

10 5 9

8 8

8 4

7

6

6 6

3

5

4

4 2 4

3

2

2 1

2

0 1

0 0

0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 8 10

Plate 7 (Figure 15.10). This is the same calculation as in Plate 8 (Figure 15.11). In this calculation two interfaces are

Figure 15.9, except that 500 grid cells are used. present, since the air surrounds the water on both sides.

This calculation uses only 100 grid cells.

velocity field

log(den) vel

500

3 0.7

400 0.6

2.5

300 0.5

2

200 0.4

1.5

100 0.3

1

0 0.2

0 2 4 6 8 10 0 2 4 6 8 10

0.1

entropy press

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

4 6

x 10 x 10

Plate 10 (Figure 17.1). A shock wave propagating through a

10

10

gas bounded on top and bottom by Lagrangian materials

with strength.

8 8

6 6

4

4

2

2

0

0

0 2 4 6 8 10 0 2 4 6 8 10

Plate 9 (Figure 15.12). This is the same calculation as in

Figure 15.11, except that 500 grid cells are used.

Plate 11 (Figure 18.1). A warm smoke plume injected from Plate 13. (Figure 19.1). A splash is generated as a sphere is

left to right rises under the influence of buoyancy. thrown into the water.

Plate 12. (Figure 18.2). Small-scale eddies are generated as Plate 14. (Figure 19.2). An interesting spray effect is gener-

smoke flows past a sphere. ated as a slightly submerged ellipse slips through the water.

Plate 15. (Figure 19.3). A thin water sheet is generated by

a sphere thrown into the water.

Plate 17. (Figure 19.5). Pouring water into a cylindrical

glass using the particle level set method.

Plate 16 (Figure 19.4). Pouring water into a cylindrical

glass using the particle level set method.

velocity field

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Plate 18. (Figure 20.1). An incompressible droplet traveling

to the right in a compressible gas flow. Note the lead shock

wave.

velocity field

1

Plate 20. (Figure 21.2). A water drop falls through the air

0.9

into the water. Surface tension forces cause the spherically

0.8 shaped region at the top of the water jet in the last frame.

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Plate 19. (Figure 20.2). A shock wave impinging on an

incompressible droplet producing a reflected wave and a

(very) weak transmitted wave.

Plate 23. (Figure 22.6). A flammable ball catches on fire as

it passes through a flame.

Plate 21. (Figure 22.4). Typical blue cores rendered using

the zero isocontour of the level set function.

Plate 22. (Figure 22.5). The density ratio of the unburnt to

burnt gas is increased from left to right, illustrating the Plate 24. (Figure 22.7). Campfire with realistic lighting of

effect of increased expansion. the surrounding rocks.

148

Copyright Design by: Sunlight webdesign