Shape Representation Using a Volume Coverage Model

(1)

Master of Science Thesis in Computer Graphics

Department of Electrical Engineering, Linköping University, 2019

Shape representation using

a volume coverage model

(2)

(3)

Abstract

Geometric shapes can be represented in a variety of different ways. A distance map is a map from points to distances. This can be used as a shape representation which can be created through a process known as a distance transform. This the-sis project tests a method for three-dimensional distance transforms using frac-tional volume coverage. This method produces distance maps with subvoxel dis-tance values. The result which is achieved is clearly better than what would be expected from a binary distance transform and similar to the one known from previous work. The resulting code has been published under a free and open source software license.

(6)

(7)

Abstract vii

Acknowledgements

First off I would like to thank my supervisor and my examiner Stefan and Inge-mar for all the help during this project. I would also like to thank my family and my friends for supporting me and giving advice when I was stuck.

(8)

(9)

1

Introduction

When working in the field of 3D graphics, some sort of shape representation is required. A common way to represent a geometric shape is as a polygonal mesh. Polygonal meshes are for example used in video games and movie production. Another way to represent a shape is as a collection of voxels. Voxels are perhaps most widely known as the building blocks in the video game Minecraft.

A distance map is a a map from points to distances represented as a grid of cells. It can be used as a shape representation with an arbitrary number of dimensions. The shape is represented by each cell containing the distance to the closest edge of the shape[5]. This is in contrast to a voxel representation where each cell only contains a value which says whether the cell is inside the shape or not. There are several possible applications of distance maps such as erosion, dilation, opening, closing, shape matching and skeletonization/thinning of shapes. Distance maps have for example been used in video games to render text on GPUs[2].

Voxel data can be used to create a distance map. The voxel data could be purely binary voxels where 1 means that the cell is filled and 0 means that the cell is empty. A better result can be achieved if the voxel data contains some information about partial coverage of elements, a form of antialiasing[3]. For example, a voxel with the value 0.5 means that the object’s edge passes through that voxel in such a way that exactly half of that voxel is covered by the shape. A voxel coverage method could potentially reduce the necessary resolution of the resulting distance map.

In this thesis the usage of distance maps generated from a volume coverage model will be investigated. The goal is to build on previous research on this sub-ject and make a new implementation of an already known algorithm. To create a coverage model from a polygonal mesh, either an analytical method or a numeric method can be used to calculate the coverage of partially covered elements. The

(10)

1.2 Research questions

The thesis tackles the following questions/tasks:

1. How well does the new implementation of a previously known distance transform algorithm perform?

2. Is a numeric approach to creating the coverage model good enough? 3. How can the generated distance map be visualized?

4. Can the software developed in this project be reused?

1.3 Limitations

(11)

2

Theory

This chapter will summarize some of the theory that will be used in the thesis. Coverage models are explained in Section 2.1 and distance transforms in Sec-tion 2.2.

2.1 Coverage model

In a pixel or voxel representation the value in each element can be used as a factor of how much that element is covered, see Figure 2.1. A value of 0 would mean that the element does not contain any part of the object. Similarly, a value of 1 means that the element is completely within the object. The elements with a value between 0 and 1 are those where the edge of the object is located some-where within. The accuracy of these elements are the ones of most significance to improve.

2.2 Distance transform

A distance map is a map from each image point to its smallest distance to a se-lected subset of image points. This subset can for example be the border of a shape. The process of creating a distance transform is called a “distance trans-form” (DT). Figure 2.2 shows a distance map generated from a part of the im-age in Figure 2.1. A distance map is similar to a distance function distance = f (x, y, z) but it is only defined in discrete points, though of course the values of these points could be interpolated to create a continous distance function. Note that Figure 2.2 shows a double sided distance map with distances on both the inside and outside of the shape. This is a special case of a distance map.

The simplest form of a distance transform would be applied to a pure binary 3

(12)

Figure 2.1: Area coverage of part of a circle

Figure 2.2:Distance map of part of a circle

image, one where each element is either completely filled or empty. This however produces a result with jagged edges in those cases where the edges are not aligned with the sampling grid. To get a better result the “grayscale” coverage data can be used.

Distance transforms can use different distance metrics. They can keep track of scalar distance or vectors to the nearest surface points. They can also use dif-ferent ways of propagating distances, one way being a sweeping scheme, where distances are propagated by scanning the grid from one side to the other. An alter-native is a wavefront scheme, where distances propagate in the order of distance from the surface. [5]

A simple way to use the coverage value to improve the distance value is to assume that the object border is aligned with the grid axes. For a more accurate result, a more complex method that estimates the border direction can be used.

2.2.1 Sweeping distance propagation

When sweeping an image, a so called “mask” is used. A mask consists of the neighboring grid cells that should be compared with the current cell to find the shortest total distance. There are different masks that can be chosen with differ-ent characteristics, for example how suitable they are for parallellization. For this thesis a simple mask is shown in Figure 2.3. This is the mask that would be used when scanning along the Y axis in the positive direction, assuming that the Y axis is pointing up. The masks used for different axes don’t have to be symmetrical with each other but in this case they are. This particular mask is simply rotated so that the gray cube points along the primary scanning axis. The gray cube marks the current cell whose value will be overwritten with the lowest possible distance value calculated from any of the cells in the mask.

An example of a scanning distance transform are shown in Figure 2.4. This scan uses a mask like the one in Figure 2.3. It can either be seen as a two-dimensional case or as a slice of the three-two-dimensional case. The black cell is

(13)

2.2 Distance transform 5

Figure 2.3:Distance propagation mask1

Figure 2.4: Scanning dis-tance transform1

Figure 2.5: Parallell scan-ning distance transform1

the one which is currently being processed. Dark grey cells are those that are being compared with the current cell to possibly find a better value. Light grey areas are those that have already been processed in this scan. Finally the white area is simply cells that remain to be processed. A parallell variant of the scan can be seen in Figure 2.5. With the chosen mask, cell processing only depends on cells in the previous layer and can therefore be parallellized within a layer in any arbitrary way[7].

See Figure 2.6 for an example of one step of the distance propagation through the grid. The black cell containing the value 2.6 is the one which is currently being processed. New potential values for the current cell will be calulated by adding the distance to each cell marked by the mask to the value in that cell. The smallest of these potential values will then be chosen as the new value for the current cell. This gives us four alternatives:

• 2.6 + 0 = 2.6

(14)

Figure 2.6:Distance propagation • 2.2 + √ 12_{+ 1}2≈_3.6 • 1.7 + √ 12_{+ 0}2≈_2.7 • 1.1 + √ 12_{+ 1}2≈_2.5

Comparing these values, 2.5 is the smallest and so it will be the new distance value written to the cell before the algorithm moves to the next cell.

2.3 Implicit meshes

One way that distance maps can be visualized is as 2-dimensional slices, for ex-ample in a grayscale format. To be able to see a projected image of the repre-sented volume it can for example either be directly rendered using ray marching or preprocessed into a polygonal mesh. Both of these methods would render the surface at the location where the distance function f (x, y, z) = 0. Defining a sur-face where a function equals some constant is known as an isosursur-face. “Marching cubes” is a well known algorithm for extracting a polygonal mesh from an isosur-face.

2.4 Related work

Under this section a number of research papers are collected and discussed in relation to this thesis.

(15)

2.4 Related work 7

2.4.1 Euclidean distance mapping

Danielsson [1] presents an early algorithm for computing an approximate eu-clidean distance transform based on vector propagation. This algorithm uses a pure binary image as its input. Many of the later algorithms are based on this one.

2.4.2 Anti-aliased Euclidean distance transform

Gustavson and Strand [3] present a distance measure that has a better accuracy for lower resolution images than earlier methods. The method uses area coverage to calculate sub-pixel distances and improve the precision.

2.4.3 Precise Euclidean transforms in 3D from voxel coverage

representation

Ilić et al. [4] present two variants of an algorithm for estimating Euclidean Dis-tance Transform (EDT) with sub-voxel precision. This algorithm will be referred to as the “Voxel Coverage Distance Algorithm” in the rest of this document and it is further described in Section 2.5. The variants differ in how they estimate the di-rection of the object border. One uses only the integer vector didi-rection calculated from a binary EDT and the other uses the binary EDT vector direction for long distances and a gradient vector of the coverage data close to the object border.

2.4.4 Exact Linear Time Euclidean Distance Transforms of Grid

Line Sampled Shapes

Lindblad and Sladoje [6] propose a method for computing the Euclidean distance transform. This method does not use a coverage representation as input however but rather a “Grid line sampling”. This is supposed to have a similar result to the method found by Ilić et al. [4] while being faster and easier to control/evalu-ate. The algorithm doesn’t take coverage data as its input but the paper gives an explanation of how grid line sampling can be approximated using the coverage data.

2.4.5 3D distance fields: a survey of techniques and

applications

Jones et al. [5] present a survey comparing different algorithms that compute three dimensional euclidean distance transforms. They explain, among other things, the difference between sweeping and a wavefront scheme and the differ-ence between a scalar and vector distance transforms.

(16)

Figure 2.7:Volume covered by object (v) and distance to center (d)

2.5 Voxel Coverage Distance Algorithm

The Voxel Coverage Distance Algorithm[4] is the algorithm used in this project to calculate the fractional component of the edge distance from voxel coverage data. The “edge distance” is the distance d from the surface of the object to the center of the cell. The voxel coverage is the fraction v of the cell’s volume which is covered by the object. See Figure 2.7 which can be seen as a cross section of a three-dimensional voxel. The object’s surface is approximated as a flat plane. Section 2.5.1 will show how the volume of the object’s intersection with a grid cell can be calculated from a known distance between the object’s surface and the grid cell’s corner. Section 2.5.2 will then show how the inverse of these equations are used in the algorithm.

2.5.1 Distance to volume

This calculation has been split into multiple cases. Depending on the number of corners that are assumed to be covered by the volume. See Figure 2.8 for an example. The first one has at most a single corner, the next one has two corners and so on up until the last which covers all but one corner. These volumes are

(17)

2.5 Voxel Coverage Distance Algorithm 9

Figure 2.8:Different cases of volume coverage

x

e

_x

e

_y y

z

e

_z

Figure 2.9:The case of a single corner covered by the object2 going to look different for every value of the surface’s normal n.

Figure 2.9 shows the case where a single corner is covered by the object. The distances from this corner to the points where the object intersect with the cell’s edges are denoted as ex, eyand ez. The values (ex, ey, ez) are related to the surface

normal n in the following way: ex= dc 1 nx , ey= dc 1 ny , ez = dc 1 nz

where dcis the shortest distance from the covered corner to the object surface.

The volume V of the intersecting part can be calculated as: V = 1

6(exeyez).

Figure 2.10 shows the case where two corners are covered by the object. The volume can be calculated with this formula:

V = 1 6(exeyez−(ex−1)ey ex−1 ex !2 ez).

Figure 2.11 shows the case where three corners are covered by the object. The volume can be calculated with this formula:

V = 1 6(exeyez−(ex−1)ey ex−1 ex !2 ez−(ey−1)ex ey−1 ey !2 ez).

(18)

Figure 2.10:The case of two corners covered by the object2

e

_x

e

_x

_z

(19)

2.5 Voxel Coverage Distance Algorithm 11

Figure 2.12:A case of four corners covered by the object2

When four corners are covered by the object, there are actually two different possibilities with different volume calculations. Figure 2.12 shows one case. The measurements have not been written out since they are similar to the three corner case in Figure 2.11. The volume can be calculated with this formula:

V = 1 6(exeyez−(ex−1)ey ex−1 ex !2 ez−(ey−1)ex ey−1 ey !2 ez−(ez−1)ex ez −₁ ez !2 ey).

Figure 2.13 shows the other possibility when four corners are covered by the object. The volume can be calculated with this formula:

V = 1 6(exeyez−(ex−1)ey ex−1 ex !2 ez−(ey−1)ex ey−1 ey !2 ez+ (ex ey−1 ey −_1)(e_yex−1 ex −_1)e_zex−1 ex (ey ex−1 ex −_1)/(e_yex−1 ex )).

2.5.2 Volume to distance

The previous section showed the calculation from distance to volume. This sec-tion will show the inverse, distance calculated from volume. This is the Voxel Coverage Distance Algorithm.

To ease the splitting of the calculation into different cases, a number of values can be defined. Each vnis defined as the largest possible volume when n corners

_z

e

x

_e

–1

x

/ ( e

_y

)

e

_x

–1

e

_x

( e

_y

– 1)

(21)

2.5 Voxel Coverage Distance Algorithm 13 v3=              1 2 ny+ nz nx if ny+ nz ≤nx 1 6 n2_z−_3n_x_n_z_{+ 3n}2 x nxny − 1 6 (nx−ny)3 nxnynz otherwise v4= 1 − v3, v5 = 1 − v2, v6= 1 − v1

The inputs to the Voxel Coverage Distance Algorithm are the voxel coverage value v where 0 ≤ v ≤ 0.5 and the direction of the object boundary n = (nx, ny, nz),

where nx, nyand nz are expected to be sorted such that nx ≥ny≥ nz. If v > 0.5,

the algorithm can still be applied by calculating −dv(n, 1 − v) instead because dv

is antisymmetric in v = 0.5 and dv(0.5) = 0. dv(n, v) =                  0.5 − v if ny= 0 (implies nz = 0) 1 2(nx+ ny+ nz) −p6vn3 xnynz if ny, 0 and 0 ≤ v ≤ v1 1 2(nx+ ny) − q 2vnxny−121n2z if ny, 0 and v1≤v ≤ v2 1 2(nx+ ny+ nz) − xnx otherwise (2.1)

The first case of Equation 2.1 is for when the normal of the object’s surface aligns with one of the grid’s axes (ny= nz = 0). This is the simplest one to

calcu-late the distance for since only one dimension needs to be taken into account. The last case of Equation 2.1 contains a value called x. This is calculated by solving one of two third degree equations. If v ≤ v3 and nz , 0 then x is calculated by solving Equation 2.2.

0 = x3n 3 z n3y −_3x2(ny+ nz)n 3 z nxn3y + 3xn 3 z(n2y+ n2z) n2xn3y + x      6v n4_z n2xn2y −n 3_(n3 y+ n3z) n3xn3y       (2.2) select x where            ny nx ≤_{x ≤} ny+ nz nx if ny+ nz≤nx ny nx ≤_{x ≤ 1} _otherwise Otherwise x is calculated by solving Equation 2.3.

0 = 2x3n 3 z n3y −_3x2(nx+ ny+ nz)n 3 z nxn3y +3xn 3 z(n2x+ n2y+ n2z) n2xn3y +x      6v n4z n2xn2y − n 3_(n3 x+ n3y+ n3z) n3xn3y       (2.3) select x where            ny nx ≤_{x ≤} ny+ nz nx if ny+ nz≤nx ny nx ≤_{x ≤ 1} _otherwise

Because all expected x values are real, the equations do not need to be solved for any complex solutions.

(22)

(23)

3

Method

The work was split in several parts: Preparation by reading previous studies and making a test implementation for a two-dimensional distance transform, an im-plementation of a three-dimensional distance transform presented by Ilić et al. [4] and the evaluation of that three-dimensional distance transform.

3.1 Preparation

Several papers were read to study distance transforms in general, binary eu-clidean distance transforms, multi-dimensional distance transforms and distance transforms that utilize coverage data.

A simple implementation of a two-dimensional distance transform based on the algorithm presented by Gustavson and Strand [3] was developed and tested on grayscale images.

3.2 Implementation

The two-dimensional program was modified to use three-dimensional data and to use the Voxel Coverage Distance Algorithm [4] from Section 2.5.

3.2.1 Volume sampling

During testing, volume data was calculated as the overlap between different ob-jects and a grid. Both perfect spheres and polygonal 3d meshes were used.

A sphere’s coverage of a cell was calculated recursively splitting a voxel and comparing the sphere’s radius against the corners of each “subvoxel”. The recur-sive process can be described like this: if a voxel is outside the sphere that voxel’s

(24)

# c o u n t i n g a s e n t i r e l y on one s i d e , o t h e r w i s e t h i s would # r e c u r s e i n f i n i t e l y .

SUBDIVISION_EPSILON = 1 . 0 / 128.0

def s u b v o x e l _ c o v e r a g e ( pos1 , s u b v o x e l _ s i d e ) : pos2 = pos1 + vec3 ( s u b v o x e l _ s i d e )

c e n t e r _ d i s t 1 = pos1 − s p h e r e _ c e n t e r c e n t e r _ d i s t 2 = pos2 − s p h e r e _ c e n t e r

min_dist = m i n _ a b s o l u t e _ c o o r d i n a t e s ( pos1 , pos2 ) max_dist = m a x _ a b s o l u t e _ c o o r d i n a t e s ( pos1 , pos2 )

i f magnitude ( min_dist ) − r a d i u s > −SUBDIVISION_EPSILON : # T h i s s u b v o x e l i s o u t s i d e i f t h e c l o s e s t c o r n e r i s return 0 . 0

e l i f magnitude ( max_dist ) − r a d i u s < SUBDIVISION_EPSILON : # T h i s s u b v o x e l i s i n s i d e i f t h e f u r t h e s t c o r n e r i s return s u b v o x e l _ s i d e ^3 e l s e: # Sum volume r e c u r s i v e l y f o r s u b d i v i d e d s u b v o x e l return sum(map( λ λ λ p : s u b v o x e l _ c o v e r a g e ( p , s u b v o x e l _ s i d e / 2 . 0 ) , s u b d i v i d e ( pos , s u b v o x e l _ s i d e ) ) )

For 3d models ray cast sampling was used instead. This method could likely have been used for the spheres as well but it is not symmetrical in all dimensions as implemented in this project. Rays were cast from the yz-plane along the x axis. This means that the x axis resolution will only be limited by floating point precision while the resolution of the y and z axis are limited to the resolution of the plane that the rays are cast from. To make the method symmetrical, rays would have to be cast in other directions.

(25)

3.3 Program interface 17 −_1.0 _0.0 _1.0 −_2.0 _0.0 _2.0 −_1.0 _0.0 _1.0 −_2.0 _0.0 _2.0 −_4.0 _0.0 _4.0 −_2.0 _0.0 _2.0 −_1.0 _0.0 _1.0 −_2.0 _0.0 _2.0 −_1.0 _0.0 _1.0

Figure 3.1:Kernel to calculate the x component of the gradient

3.2.2 Calculating coverage gradient

The variant of the Voxel Coverage Distance Algorithm used in this thesis requires gradient vectors as edge direction data. This was calculated by convolving the coverage data with a kernel for each axis, specifically a Sobel operator. See Fig-ure 3.1 for the kernel used to calculate the gradient’s x component. For y and z similar kernels were used, but with swapped axes.

3.2.3 Distance transform

The coverage data that was calculated was used as input for the implemented algorithm. The distance transform tested in this thesis is a sweeping DT that uses scalar values for sub-voxel distances but vectors for distances between voxels. The mask used for the scan is the one seen in Figure 2.3.

3.2.4 Polygonal mesh recreation

A variant of the marching cubes algorithm was used to recreate a polygonal mesh from the distance data. This was done to see how high the resolution of the sampling grid would have to be to get a satisfactory result. The implementation of this algorithm was not part of this thesis and a prepackaged software library was used3.

3.3 Program interface

The developed program4 has a command line interface that can process data in hdf5 files. There is also a library interface to be linked with other programs. The

3_{https://github.com/swiftcoder/isosurface}

(26)

command line program has subcommands for several tasks listed in Table 3.1. Every subcommand has parameters to control its execution. The full documenta-tion for the command line interface is located in Appendix A.

3.3.1 C library interface

The library interface contains a C function signature listed in Listing 3.2. The size_x, size_y and size_z arguments contain the size of the grid that the data was sampled in. The argument coverage_data contains the input cover-age data to be processed. It should have a length of size_x · size_y · size_z elements. The data layout in memory is expected to be designed such that it can be indexed like in (3.1).

Listing 3.2:C language interface void distance_transform_3d_volume_coverage ( s i z e _ t s i z e _ x , s i z e _ t s i z e _ y , s i z e _ t s i z e _ z , f l o a t (* coverage_to_distance ) ( s t r u c t V e c t o r 3 i i n t e g e r _ d i s t a n c e , s t r u c t V e c t o r 3 f gradient_normal , f l o a t coverage , f l o a t o l d _ d i s t a n c e ) , const f l o a t * coverage_data , f l o a t * distance_buffer )

index= size_x · size_y · z + size_x · y + x (3.1) The argument coverage_to_distance is a callback that the user provides. The callback receives, for each cell, an integer vector which is the current estima-tion to the closest edge, the normalized gradient in that cell of the coverage data,

(27)

3.4 Evaluation 19

the coverage fraction for the cell and the old distance value. The callback should return an estimate of the entire distance, including the length of the integer vec-tor. There is an included callback function called default_coverage_to_di stancewhich implements the Voxel Coverage Distance Algorithm described in Section 2.5.

3.4 Evaluation

To test the implementation both a perfect sphere and 3d model data was used as input data. There was a numerical evaluation that calculated the accuracy of dis-tance maps from volume sampled spheres. There was also a subjective evaluation of reconstructed 3d models.

3.4.1 Numeric evaluation

The finished implementation was tested in a way similar to how Ilić et al. [4] tested their implementation. 50 positions were randomly generated; these were used as the center point for spheres with radii between 1 and 30 voxels. Each sphere’s volume coverage values was sampled on a grid and this data was then transformed to a distance map. The calculated distance was then compared to the true distance to the randomized sphere’s surface and averaged between all the different center points.

3.4.2 Subjective evaluation

A survey was conducted to evaluate the quality of the method’s visual result when a model was reconstructed. People were asked how well they thought re-constructed models, made from distance maps of varying resolutions, matched the original.

A Utah teapot model, see Figure 3.2, was scanned using raycasts. Then the volume data was transformed into a distance map and from the distance map a new 3d model was generated. The rendered images can be seen in Figure 3.3 A number of people were asked the question: “Which resolution do you think leads to a sufficiently accurate recreation of the original model?”.

The grid resolutions that were available as choses were: • 23 × 12 × 15, labeled “10 voxels”

• 33 × 17 × 22, labeled “15 voxels” • 43 × 22 × 28, labeled “20 voxels” • 64 × 32 × 41, labeled “30 voxels”

Each picture was labeled based on the number entered as the resolution on the command line. This number is used as the length of the shortest side of the model with the rest automatically sized based on the proportions of the model.

(28)

Figure 3.2:The original Utah teapot 3d model.

During the calculation each grid dimension also had a margin of one voxel on each side. This is the reason all of the grids were two voxels larger than the model. The results of the survey are discussed in Section 4.3.1.

(29)

3.4 Evaluation 21

Figure 3.3:The images used in the survey. The label under each image is the width of the teapot in voxels.

(30)

(31)

4

Results

This chapter presents the results of the project. Section 4.1 presents the results of the preparation experiment, Section 4.2 presents the results of the implemen-tation and Section 4.3 contains the results of the evaluation.

4.1 Preparation

The image used to test the experimental application is shown in Figure 4.1. This was turned into the distance map in Figure 4.2. The calculated distance vectors can be seen in Figure 4.3.

4.2 Implementation

The final implementation can be used as either a command line application, a Rust library or a C library. The code written for this project is released under the GPLv3 license. It is written in the programming language Rust but bindings for C also exist for easier interoperability with other software.

The implementation did not end up using any parallellization to speed up the process. Making the code fully optimized was not a top priority of the project.

There was an experiment with a ray marching renderer of the distance map but this was later discarded due to the increased workload.

4.3 Evaluation

Figure 4.7 shows the error range when testing the implementation on the ran-domly positioned spheres. The error range was calculated as the difference be-tween the maximum and the minimum error. Each datapoint is the average value

(32)

Figure 4.1:Image used as input for the experimental 2d distance transform.

Figure 4.2: Distances of the resulting 2d distance map. Distances outside the object are shown in the left image and distances inside are shown in the right one.

(33)

4.3 Evaluation 25

Figure 4.3:A downsampled version of the vector field calculated when gen-erating the 2d distance map.

(34)

Figure 4.5: Recreated polygonal mesh. The input model was sampled on a 23 × 12 × 15 grid.

Figure 4.6: Recreated polygonal mesh. The input model was sampled on a 84 × 42 × 53 grid.

(35)

4.3 Evaluation 27

Figure 4.7:Error range for spheres

for all positions. Figure 4.8 shows the average absolute value error of all posi-tions.

4.3.1 Polygonal mesh recreation

Different people were asked which of the different images from Section 3.4.2 had a sufficient resolution. The responses from the 11 people who were asked are shown in Figure 4.9. Some people gave more detailed responses which can be seen in Table 4.1.

(36)

Figure 4.8:Mean absolute error (MAE) for spheres 0 10 voxels 0 15 voxels 6 20 voxels 5 30 voxels 0 1 2 3 4 5 6

(37)

4.3 Evaluation 29

You can see subtle differences between 20 and 30 but I think it looks good at 20. What does “sufficiently accurate” mean? It’s pretty obviously a reasonable teapot at 15 but at 20 it’s a good quality teapot.

I’d say 30 for identical, 20 is ok.

There’s a much bigger difference between 10 and 20 than there is between 20 and 30. So “the rate of improvement” has decreased inbetween.

I would say that 20 is minimum, of course it’s more exact the higher resolution you have.

(38)

(39)

5

Discussion

In this chapter the method from Chapter 3 is discussed in Section 5.1 and the results from Chapter 4 are discussed in Section 5.2.

5.1 Method

There was a condition in the described algorithm that lead to some problems. This is Equation 9 from Ilić et al. [4]:

dV CEDT 3=        da+ d grad v , for |dvgrad| ≤0.5 ∧ db≤0.5

dBEDT + dvort, otherwise.

da= di·n, db= |di×n|

where di is the direction of the edge as an integer vector and n is the normal

of the edge. dvgrad is the result of the Voxel Coverage Distance Algorithm using

the gradient as the direction of the normal and dvortis the result when using the

integer edge direction vector as the direction. The condition for the first alterna-tive is |dvgrad| ≤0.5. This will not be true in corners of cells where the distance to

the edge is higher than 0.5, see Figure 5.1, which means that dvortwill have to be

used. However, the integer distance vector is (0, 0, 0) inside the cell and therefore does not have a direction. 0.75 was used instead of 0.5 which gave a result closer to the one in the original report.

Perhaps the code could be optimized by first running a purely binary distance transform and when it finishes, switching to the variant that uses coverage data to improve the result.

(40)

Figure 5.1:Areas of a cell where the distance to the center is greater than 0.5

5.2 Result

It was discovered that edges mostly parallell to the axes of the sampling space actually ended up with larger errors than those with a steeper angle. This is probably because the subvoxel quantization is more coarse in that case. Because the change in distance is lower when the angle is small, all subvoxels quantize to the same level.

5.2.1 Subjective evaluation

The answers on the survey were pretty evenly split between saying that the 20 voxel resolution and that the 30 voxel resolution were sufficient. It was some-what clear when picking the images that the choice would be between those two. The question about “sufficient resolution” was a bit vague, something that was noted by the people who responded. The answer depends on what you define as sufficient.

(41)

6

Conclusion

The result of this thesis project is a command line application to transform vol-umetric coverage data into a distance map. A library was also created to let the code be used in other applications.

The implementation of the Volume Coverage Algorithm recreates the result of the original report well. The result when testing with sampled spheres is not quite as accurate but this might be due to differences in the sampling method or something else.

Using a marching cubes algorithm turned out to be an acceptable way to visu-alize the resulting distance map.

6.1 Future work

The method presented in [6] could potentially be better in the case of the ray casting method that was used for the 3d model because the originally sampled data is distances on a plane.

The numerice approach used when sampling the spheres seemed to be suf-ficient to be close to a recreation of the original result. However, an analytical method for calculating the volume coverage of the spheres would have given more accurate test data.

The resulting data could possibly have been rendered using a ray marching method instead to get a more accurate result for the end visualization. An at-tempt at this was made but it was not completed due to problems with overesti-mated distances. The visual result was not the central point of the project so it was not the most important thing to focus on.

The model that was raycasted was treated as if it consisted of square columns along the ray axis. The result could probably be improved by interpolating sur-face locations, possibly using the sursur-face normal, or by ray casting in multiple

(42)

(43)

A

Command line documentation

This appendix contains a detailed documentation of the command line interface of the developed application. All data that is written to files is stored in the Hdf5 format1. Datasets that already exist in a chosen file will be overwritten without warning. By default each subcommand uses its own unique dataset name to pre-vent conflicts. Any output file that does not already exist will be created. All sub-commands that output data write a single dataset containing a three-dimensional array of 32-bit floating point numbers.

A.1 Sample perfect spheres

The subcommand for calculating the volume coverage of a perfect sphere is sample-sphere <radius> <output-file>.

The parameter radius is the radius of the sphere to sample measured in voxel lengths. The radius must currently be an integer. The sampling grid will sized automatically to fit the sphere. The output-file should be an Hdf5 file. The dataset name defaults to “coverage”.

There are optional parameters to specify the name of the dataset that will be written to in the output file and to offset the sphere from the center of the grid by a certain distance. The parameters are described in Table A.1.

The offset is specified in the format x,y,z, for example: --offset 0.5,-0.3,0.

1_{https://en.wikipedia.org/wiki/Hierarchical_Data_Format#HDF5}

(44)

can be used. The obj-file should be an .obj file which is used as input. The output-fileworks the exact same way as the identically named parameter for sample-sphere, it should be an Hd5 file and the default dataset name is “cov-erage”. Note that this default dataset name is the same as for sample-sphere.

There are two optional parameters to specify the dataset in the output file and the resolution of the sampling grid, they are described in Table A.2. The specified resolution number will be used as the length of the shortest side and the others will be sized automatically to keep the size proportional to the original model. The default resolution is 16.

Option Description

--output-dataset dataset name for output data --resolution resolution of the sampling grid Table A.2:Optional flags for the subcommand sample-obj

A.3 Apply distance transform

To actually apply the distance transform, the subcommand transform <input-file>

is used. The input-file should be an Hdf5 file containing volume coverage samples. By default the same file is used for output and the default dataset names for input and output data are “coverage” and “distance” respectively. Flags to choose input dataset name, input file and output file are described in Table A.3.

A.4 Calculate error of sphere distance map

The subcommand

(45)

A.5 Generate and render polygonal model from distance map 37

Option Description

--input-dataset dataset name for input data --output-dataset dataset name for output data

--output-files file to write output data to Table A.3:Optional flags for the subcommand transform

is used to evaluate the accuracy of a distance map calculated from a sphere using sample-sphere. Naturally, the input file is chosen as input-file and sphere radius using radius. Remember to use the same radius as when generating the sphere’s volume coverage. There is an optional parameter called --offset that sets the sphere’s center offset. The center should obviously also be the same as when sampling the sphere.

A.5 Generate and render polygonal model from

distance map

The subcommand

render-polygonized <file>

is used to visualize a distance map in the provided file by generating a polygonal model and rendering it in an OpenGL window. A variant of Marching Cubes is used to generate the mesh. The resolution of the generated mesh and the dataset which is used from the input file can be customized using flags described in Ta-ble A.4. The default mesh resolution is 6 and “distance” is the default dataset name.

Option Description

--marching-cubes-resolution resolution of the generated mesh --dataset dataset name for input data Table A.4:Optional flags for the subcommand render-polygonized

(46)

(47)

Bibliography

[1] Per-Erik Danielsson. Euclidean distance mapping. Graphical Models /graph-ical Models and Image Processing /computer Vision, Graphics, and Image Processing - CVGIP, 14:227–248, 11 1980. doi: 10.1016/0146-664X(80) 90054-4.

[2] Chris Green. Improved alpha-tested magnification for vector textures and special effects. In ACM SIGGRAPH 2007 Courses, SIGGRAPH ’07, page 9–18, New York, NY, USA, 2007. Association for Computing Machinery. ISBN 9781450318235. doi: 10.1145/1281500.1281665. URL https://doi.org/ 10.1145/1281500.1281665.

[3] Stefan Gustavson and Robin Strand. Anti-aliased euclidean distance trans-form. Pattern Recognition Letters, 32(2):252 – 257, 2011. ISSN 0167-8655. doi: https://doi.org/10.1016/j.patrec.2010.08.010. URL http://www. sciencedirect.com/science/article/pii/S0167865510002953. [4] Vladimir Ilić, Joakim Lindblad, and Nataša Sladoje. Precise euclidean

dis-tance transforms in 3d from voxel coverage representation. Pattern Recog-nition Letters, 65:184 – 191, 2015. ISSN 0167-8655. doi: https://doi.org/ 10.1016/j.patrec.2015.07.035. URL http://www.sciencedirect.com/ science/article/pii/S0167865515002445.

[5] Mark W. Jones, J. Andreas Bærentzen, and Milos Sramek. 3d distance fields: A survey of techniques and applications. IEEE TRANSACTIONS ON VISU-ALIZATION AND COMPUTER GRAPHICS, 12:581–599, 2006.

[6] Joakim Lindblad and Nataša Sladoje. Exact linear time euclidean distance transforms of grid line sampled shapes. In Jón Atli Benediktsson, Jocelyn Chanussot, Laurent Najman, and Hugues Talbot, editors, Mathematical Mor-phology and Its Applications to Signal and Image Processing, pages 645–656, Cham, 2015. Springer International Publishing. ISBN 978-3-319-18720-4. [7] I Ragnemalm. The euclidean distance transform and its implementation on

simd architectures. In Proceedings of the 6th Scandinavian Conference on Image Analysis, pages 379–384, 1989.

Shape Representation Using a Volume Coverage Model

Master of Science Thesis in Computer Graphics

Department of Electrical Engineering, Linköping University, 2019