4.1.1 Closest Point Transformation

As mentioned in Section 3.1.3 the LS function is usually initialized as a signed distance function. In case of the sparse field LS method it is better to use the smallest (signed) Manhattan distance

$${\Phi}({\vec{p}}):=\pm\min_{{\vec{x}}\in{\mathcal{S}}}{\lVert{\vec{p}}-{\vec{x}}\rVert}_1,$$ (4.1)

rather than the smallest Euclidean distance, for the initialization. With the latter, the first time step of the sparse field LS method gives wrong results for the position of a (non-axis-parallel) plane moving with constant speed. However, if initialized with the Manhattan distance, trilinear interpolation of the LS function correctly describes the position of the plane after the first time step. The reason for this behavior is the update scheme (3.29) which also corresponds to a Manhattan distance calculation.

It is not necessary to calculate the initial LS values for the entire grid. Only grid points close to the boundary have to be considered. In particular, the sparse field method requires the determination of all active grid points together with their LS values only at the beginning. If the sign of the LS function is known for all other grid points, the additional layers can be determined using (3.28). The required LS values can then be computed using (3.29). The knowledge of the sign for the LS values at all points of the grid is also a prerequisite for the setup of the H-RLE data structure.

The signs of the LS values are unambiguously defined, if the set of all grid points with an opposite signed neighbor is known. This set of grid points clearly separates the positive and the negative grid points from each other. To satisfy all the requirements for the initialization of the sparse field method and the H-RLE data structure, the LS framework determines all grid points for which the LS value is in the range $\left[-1,1\right]$ .

This set of grid points contains all the required information, because due to (4.1) all grid points are included, which have an opposite signed neighbor. Before setting up the H-RLE data structure, the LS framework collects all these grid points together with their LS values and stores them in a list. In the following sections an efficient technique is described to obtain this list in linear time.

Since the Manhattan distance is only needed for grid points close to the boundary, the distance computation can be simplified. For these points the Manhattan distance can be approximated by the smallest distance to the boundary along any paraxial direction. In other words, the smallest unsigned distance of a grid point is approximated by

$$\lvert{\Phi}({\vec{p}})\rvert:=\min\lbrace \alpha\geq 0:\exists k... ...n{\mathcal{S}}\ \vee {\vec{p}}- \alpha\cdot{\vec{e}}_k\in{\mathcal{S}} \rbrace.$$ (4.2)

Here ${\vec{e}}_k$ denotes the unit vector pointing in the ${x}_k$ -direction ( $1\leq k\leq{D}$ ).

In practice the distance transform is carried out by an iteration over all segments of the discretized boundary representation. For each segment all intersecting grid lines are determined. The number of possible grid lines can be confined using the bounding box of the segment. In the two-dimensional case the test whether a grid line intersects a line segment or not is trivial. In the three-dimensional case the intersection test with a triangle is more difficult and is described in Appendix A. For each intersecting grid line the intersection point ${\vec{q}}=\left({q}_1,\ldots,{q}_{D}\right)$ is calculated and all grid points ${\vec{p}}$ on that grid line are considered for which

$$\left\vert{p}_k-{q}_k\right\vert\leq1+{\varepsilon}_1.$$ (4.3)

Here it is assumed that the grid line is parallel to the $k$ -th axis direction, hence parallel to ${\vec{e}}_k$ . ${\varepsilon}_1$ is a small positive constant ( $0<{\varepsilon}_1\ll1$ ), introduced for numerical reasons, which ensures the calculation of the initial LS values for all grid points directly neighbored to the boundary. Thus, ${\varepsilon}_1$ must be larger than the maximum numerical error for the calculation of the distance $\left\vert{p}_k-{q}_k\right\vert$ . In this work ${\varepsilon}_1=\num{e-4}$ is used.

The signed distance of a grid point ${\vec{p}}$ to the current segment is given by

$${\tau}({\vec{p}}, k) = \funcsgn ({n}_k) \cdot \left({p}_k-{q}_k\right).$$ (4.4)

The sign of the $k$ -th component of the normal vector ${\vec{n}}$ can be obtained by considering the orientation of the surface segment with respect to ${\vec{e}}_k$ . In two dimensions the distance is positive, if the vertex ${\vec{x}}_1$ is on the left-hand side. In case of three dimensions the sign is positive, if the vertices ${\vec{x}}_1$ , ${\vec{x}}_2$ , and ${\vec{x}}_3$ are seen counterclockwise.

Figure 4.2: (a) Two segments of the surface mesh (gray) meet on a grid line (black). Hence, for both grid points the distance is equal to both segments. As consequence the determination of the signed distance according to (4.4) is ambiguous. (b) The distance transform can produce inefficient sets of active grid points. The bottom left grid point does not have an opposite signed neighbor.

If for each grid point ${\vec{p}}$ the distance with the smallest absolute value $\left\vert{p}_k-{q}_k\right\vert$ is kept, at least all grid points with LS values in the range $\left[-1,1\right]$ will be properly initialized after iterating over all segments of the boundary mesh. However, this procedure can lead to problems as depicted in Figure 4.2a, where the wrong sign could be assigned to the LS value of both grid points, since both segments are equally distanced. To resolve this ambiguity without additional consideration of neighbor segments, another distance is measured using

$${\tau}'({\vec{p}}, k) = \funcsgn ({p}_k-{q}_k) \cdot ({p}_k-{q}_k-{\varepsilon}_2\cdot {t}_k)$$ (4.5)

in order to determine the closest segment for a certain grid point ${\vec{p}}$ . Here ${\vec{t}}=\left({t}_1,\ldots,{t}_{D}\right)$ denotes the unit vector pointing from ${\vec{q}}$ to the centroid $\vec{c}$ of the segment

$${\vec{t}}:=\frac{{\vec{c}}-{\vec{q}}}{\lVert{\vec{c}}-{\vec{q}}\rVert} \ {\vec{c}}:=\frac{1}{{D}}\sum_{i=1}^{{D}}{\vec{x}}_i.$$ (4.6)

${\varepsilon}_2$ is again a small positive constant ( $0<{\varepsilon}_2\ll1$ ). However, the distance which is finally assigned to ${\Phi}({\vec{p}})$ is still calculated using (4.4). ${\varepsilon}_2=\num{e-6}$ is used for all LS initializations in this work.

The initialization procedure is now as follows: For each segment and for all intersecting grid lines all grid points fulfilling (4.3) are determined. The indices ${\vec{p}}$ of these grid points are appended to a list along with their corresponding distances ${\tau}$ and ${\tau}'$ , defined in (4.4) and (4.5), respectively. Finally the list is lexicographically sorted with respect to the grid point indices. If there are more entries with the same ${\vec{p}}$ , (which is not very often the case,) ${\tau}$ of that entry with the smallest corresponding ${\tau}'$ is used for the initial LS value ${\Phi}({\vec{p}})$ .

The entire initialization algorithm has a complexity of ${\mathcal{O}}({{N}_\text{D}}\log {{N}_\text{D}}+ {{N}_\text{S}})$ , where ${{N}_\text{D}}$ is the final number of defined grid grid points in the H-RLE data structure. ${{N}_\text{S}}$ is the total number of segments of the boundary mesh. The logarithmic term is due to the sorting algorithm which cannot be avoided, since the setup of the H-RLE data structure requires a sorted list as well. This initialization algorithm can produce inefficient sets of active grid points as exemplified in Figure 4.2b, which can be avoided by appending a pruning procedure, as mentioned earlier (see Section 3.4.2).