diff --git a/docs/make.jl b/docs/make.jl index 3bc2dec41..a00209908 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -104,8 +104,8 @@ makedocs(sitename="TrixiParticles.jl", "Tutorial" => "tutorial.md", "Examples" => "examples.md", "Visualization" => "visualization.md", - "preprocessing" => [ - "Sample Geometries" => joinpath("preprocessing", "preprocessing.md"), + "Preprocessing" => [ + "Sampling of Geometries" => joinpath("preprocessing", "preprocessing.md"), ], "Components" => [ "Overview" => "overview.md", diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md index 8d83c9caf..5d367d015 100644 --- a/docs/src/preprocessing/preprocessing.md +++ b/docs/src/preprocessing/preprocessing.md @@ -1,6 +1,245 @@ +# Sampling of Geometries -TODO +Generating the initial configuration of a simulation requires filling volumes (3D) or surfaces (2D) of complex geometries with particles. +The algorithm to sample a complex geometry should be robust and fast, +since for large problems (large numbers of particles) or complex geometries (many geometry faces), +generating the initial configuration is not trivial and can be very expensive in terms of computational cost. +We therefore use a [winding number](https://en.wikipedia.org/wiki/Winding_number) approach for an inside-outside segmentation of an object. +The winding number ``w(\mathbf{p})`` is a signed integer-valued function of a point ``\mathbf{p}`` and is defined as +```math +w(\mathbf{p}) = \frac{1}{2 \pi} \sum^n_{i=1} \Theta_i. +``` + +Here, ``\Theta_i`` is the *signed* angle between ``\mathbf{c}_i - \mathbf{p}`` and ``\mathbf{c}_{i+1} - \mathbf{p}`` where ``\mathbf{c}_i`` and ``\mathbf{c}_{i+1}`` are two consecutive vertices on a curve. +In 3D, we refer to the solid angle of an *oriented* triangle with respect to ``\mathbf{p}``. + +We provide the following methods to calculate ``w(\mathbf{p})``: +- Horman et al. (2001) evaluate the winding number combined with an even-odd rule, but only for 2D polygons (see [WindingNumberHorman](@ref)). +- Naive winding: Jacobson et al. (2013) generalized the winding number so that the algorithm can be applied for both 2D and 3D geometries (see [WindingNumberJacobson](@ref)). +- Hierarchical winding: Jacobson et al. (2013) also introduced a fast hierarchical evaluation of the winding number. For further information see the description below. + +## [Hierarchical Winding](@id hierarchical_winding) +According to Jacobson et al. (2013) the winding number with respect to a polygon (2D) or triangle mesh (3D) is the sum of the winding numbers with respect to each edge (2D) or face (3D). +We can show this with the following example in which we determine the winding number for each edge of a triangle separately and sum them up: + +```julia +using TrixiParticles +using Plots + +triangle = [125.0 375.0 250.0 125.0; + 175.0 175.0 350.0 175.0] + +# Delete all edges but one +edge1 = deleteat!(TrixiParticles.Polygon(triangle), [2, 3]) +edge2 = deleteat!(TrixiParticles.Polygon(triangle), [1, 3]) +edge3 = deleteat!(TrixiParticles.Polygon(triangle), [1, 2]) + +algorithm = WindingNumberJacobson() + +grid = hcat(([x, y] for x in 1:500, y in 1:500)...) + +_, w1 = algorithm(edge1, grid; store_winding_number=true) +_, w2 = algorithm(edge2, grid; store_winding_number=true) +_, w3 = algorithm(edge3, grid; store_winding_number=true) + +w = w1 + w2 + w3 + +heatmap(1:500, 1:500, reshape(w1, 500, 500)', color=:coolwarm, showaxis=false, + tickfontsize=12, size=(570, 500), margin=6 * Plots.mm) +heatmap(1:500, 1:500, reshape(w2, 500, 500)', color=:coolwarm, showaxis=false, + tickfontsize=12, size=(570, 500), margin=6 * Plots.mm) +heatmap(1:500, 1:500, reshape(w3, 500, 500)', color=:coolwarm, showaxis=false, + tickfontsize=12, size=(570, 500), margin=6 * Plots.mm) +heatmap(1:500, 1:500, reshape(w, 500, 500)', color=:coolwarm, showaxis=false, + tickfontsize=12, size=(570, 500), margin=6 * Plots.mm, clims=(-1, 1)) + +``` + +```@raw html +
+ triangle +
+``` + +This summation property has some interesting consequences that we can utilize for an efficient computation of the winding number. +Let ``\mathcal{S}`` be an open surface and ``\bar{\mathcal{S}}`` an arbitrary closing surface, such that + +```math +\partial \bar{\mathcal{S}} = \partial \mathcal{S} +``` + +and ``\mathcal{B} = \bar{\mathcal{S}} \cup \mathcal{S}`` is some closed oriented surface. +For any query point ``\mathbf{p}`` outside of ``\mathcal{B}``, we know that + +```math +w_{\mathcal{S}}(\mathbf{p}) + w_{\bar{\mathcal{S}}}(\mathbf{p}) = w_{\mathcal{B}}(\mathbf{p}) = 0. +``` + +This means + +```math +w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p}), +``` + +regardless of how ``\bar{\mathcal{S}}`` is constructed (as long as ``\mathbf{p}`` is outside of ``\mathcal{B}``). + +We can use this property in the discrete case to efficiently compute the winding number of a query point +by partitioning the polygon or mesh in a "small" part (as in consisting of a small number of edges/faces) and a "large" part. +For the small part we just compute the winding number, and for the large part we construct a small closing and compute its winding number. +The partitioning is based on a hierarchical construction of bounding boxes. + +### Bounding volume hierarchy + +To efficiently find a "small part" and a "large part" as mentioned above, we construct a hierarchy of bounding boxes by starting with the whole domain and recursively splitting it in two equally sized boxes. +The resulting hierarchy is a binary tree. + +The algorithm by Jacobsen et al. (Algorithm 2, p. 5) traverses this binary tree recursively until we find the leaf in which the query point is located. +The recursion stops with the following criteria: + +- if the bounding box ``T`` is a leaf then ``T.\mathcal{S} = \mathcal{S} \cap T``, the part of ``\mathcal{S}`` + that lies inside ``T``, is the "small part" mentioned above, so evaluate the winding number naively as ``w(\mathbf{p}, T.\mathcal{S})``. +- else if ``\mathbf{p}`` is outside ``T`` then ``T.\mathcal{S}`` is the "large part", so evaluate the winding number naively + as ``-w(\mathbf{p}, T.\bar{\mathcal{S}})``, where ``T.\bar{\mathcal{S}}`` is the closing surface of ``T.\mathcal{S}``. + +#### Continuous example + +Now consider the following continuous (not discretized to a polygon) 2D example. +We compute the winding number of the point ``\mathbf{p}`` with respect to ``\mathcal{S}`` using the depicted hierarchy of bounding boxes. + +```@raw html +
+ continuous closing +
+``` + +(1): +- Recurse left: ``w_{\text{left}} = \text{\texttt{hierarchical\_winding}} (\mathbf{p}, T.\text{left})`` +- Recurse right: ``w_{\text{right}} = \text{\texttt{hierarchical\_winding}} (\mathbf{p},T.\text{right})`` + +(2): +- Query point ``\mathbf{p}`` is outside bounding box ``T``, so don't recurse deeper. +- Compute ``w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p})`` with the closure ``T.\bar{\mathcal{S}}``, which is generally much smaller (fewer edges in the discrete version) than ``T.\mathcal{S}``: + +```math +w_{\text{left}} = -\text{\texttt{naive\_winding}} (\mathbf{p}, T.\bar{\mathcal{S}}) +``` + +(3): +- Bounding box ``T`` is a leaf. Use open surface ``T.\mathcal{S}``: + +```math +w_{\text{right}} = \text{\texttt{naive\_winding}} (\mathbf{p}, T.\mathcal{S}) +``` + +The reconstructed surface will then look as in the following image. + +```@raw html +
+ reconstructed surface +
+``` + +We finally sum up the winding numbers + +```math +w = w_{\text{left}} + w_{\text{right} } = -w_{T_{\text{left}}.\bar{\mathcal{S}}} + w_{T_{\text{right}}.\mathcal{S}} +``` + +#### Discrete example + +We will now go through the discrete version of the example above. + +```@raw html +
+ discrete geometry +
+``` + +To construct the hierarchy for the discrete piecewise-linear example in (1), we have to do the following. + +(2): +Each edge is distributed to the child whose box contains the edge's barycenter (red dots in (2)). +Splitting stops when the number of a box's edges slips below a +threshold (usually ``\approx 100`` faces in 3D, here: 6 edges). + +(3): +For the closure, Jacobson et al. (2013) define *exterior vertices* (*exterior edges* in 3D) +as boundary vertices of such a segmentation (red dots in (3)). +To find them, we traverse around each edge (face in 3D) in order, and +increment or decrement for each vertex (edge) a specific counter. + +```julia +v1 = edge_vertices_ids[edge][1] +v2 = edge_vertices_ids[edge][2] + +vertex_count[v1] += 1 +vertex_count[v2] -= 1 +``` + +In 2D, a vertex is declared as exterior if `vertex_count(vertex) != 0`, so there is not the same amount of edges in this box going into versus out of the vertex. +To construct the closing surface, the exterior vertices are then connected to one arbitrary +exterior vertex using appropriately oriented line segments: + +```julia +edge = vertex_count[v] > 0 ? (closing_vertex, v) : (v, closing_vertex) +``` + +The resulting closed surface ``T.S \cup T.\bar{S}`` then has the same number of edges going into and out of each vertex. + +#### Incorrect evaluation + +If we follow the algorithm, we know that recursion stops if + +- the bounding box ``T`` is a leaf or +- the query point ``\mathbf{p}`` is outside the box. + +```@raw html +
+ incorrect evaluation +
+``` + +(1): The query point ``\mathbf{p}`` is outside the box, so we calculate the winding number with the (red) closure of the box. + +(2): The query point ``\mathbf{p}`` is inside the box, so we use the (blue) edges distributed to the box. + +(3): In this case, it leads to an incorrect evaluation of the winding number. +The query point is clearly inside the box, but not inside the reconstructed surface. +This is because the property ``w_{\mathcal{S}}(\mathbf{p}) = - w_{\bar{\mathcal{S}}}(\mathbf{p})`` +only holds when ``\mathbf{p}`` is outside of ``\mathcal{B}``, which is not the case here. + +#### Correct evaluation +Jacobson et al. (2013) don't mention this problem or provide a solution to it. +We contacted the authors and found that they know about this problem and solve it +by resizing the bounding box to fully include the closing surface +of the neighboring box, since it doesn't matter if the boxes overlap. + +```@raw html +
+ correct evaluation resizing +
+``` + +To avoid resizing, we take a different approach and calculate the closure of the bounding box differently: +- Exclude intersecting edges in the calculation of the exterior vertices. +- This way, all exterior vertices are inside the bounding box, and so will be the closing surface. +- The intersecting edges are later added with flipped orientation, + so that the closing is actually a closing of the exterior plus intersecting edges. + +```@raw html +
+ correct evaluation intersecting +
+``` + +The evaluation then looks as follows. + +```@raw html +
+ correct evaluation intersecting 2 +
+``` ```@autodocs Modules = [TrixiParticles] @@ -17,7 +256,6 @@ Modules = [TrixiParticles] Pages = [joinpath("preprocessing", "geometries", "io.jl")] ``` - ### [References](@id references_complex_shape) - Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung "Robust inside-outside segmentation using generalized winding numbers". In: ACM Transactions on Graphics, 32.4 (2013), pages 1--12. diff --git a/src/preprocessing/point_in_poly/winding_number_horman.jl b/src/preprocessing/point_in_poly/winding_number_horman.jl index b1b6f9cba..cae358049 100644 --- a/src/preprocessing/point_in_poly/winding_number_horman.jl +++ b/src/preprocessing/point_in_poly/winding_number_horman.jl @@ -4,8 +4,6 @@ Algorithm for inside-outside segmentation of a complex geometry proposed by [Horman et al. (2001)](@ref references_complex_shape). It is only supported for 2D geometries. [`WindingNumberHorman`](@ref) might handle edge cases a bit better, since the winding number is an integer value. -Also, it is faster than [`WindingNumberJacobson`](@ref) for 2D geometries with more than about 100 edges. - !!! warning "Experimental Implementation" This is an experimental feature and may change in any future releases. diff --git a/src/preprocessing/point_in_poly/winding_number_jacobson.jl b/src/preprocessing/point_in_poly/winding_number_jacobson.jl index 888df222c..fabf2f541 100644 --- a/src/preprocessing/point_in_poly/winding_number_jacobson.jl +++ b/src/preprocessing/point_in_poly/winding_number_jacobson.jl @@ -58,7 +58,7 @@ Algorithm for inside-outside segmentation of a complex geometry proposed by Jaco - `geometry`: Complex geometry returned by [`load_geometry`](@ref) and is only required when using `hierarchical_winding=true`. - `hierarchical_winding`: If set to `true`, an optimized hierarchical approach will be used, - which gives a significant speedup. + which gives a significant speedup. For further information see [Hierarchical Winding](@ref hierarchical_winding). - `winding_number_factor`: For leaky geometries, a factor of `0.4` will give a better inside-outside segmentation. !!! warning "Experimental Implementation"