Add OrientedBoundingBox for spatial region analysis and particle selection

docs/src/general/initial_condition.md

@@ -11,3 +11,8 @@ Pages = [joinpath("general", "initial_condition.jl")]
 Modules = [TrixiParticles]
 Pages = map(file -> joinpath("setups", file), readdir(joinpath("..", "src", "setups")))
 ```
+
+```@autodocs
+Modules = [TrixiParticles]
+Filter = t -> t === TrixiParticles.OrientedBoundingBox
+```

src/TrixiParticles.jl

@@ -96,7 +96,7 @@ export RectangularTank, RectangularShape, SphereShape, ComplexShape
 export ParticlePackingSystem, SignedDistanceField
 export WindingNumberHormann, WindingNumberJacobson
 export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry,
-       sample_boundary, planar_geometry_to_face
+       sample_boundary, planar_geometry_to_face, OrientedBoundingBox, is_in_oriented_box
 export SourceTermDamping
 export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
        GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection

src/preprocessing/geometries/geometries.jl

@@ -92,7 +92,7 @@ face, face_normal = planar_geometry_to_face(planar_geometry)
 function planar_geometry_to_face(planar_geometry::TriangleMesh)
     face_normal = normalize(sum(planar_geometry.face_normals) / nfaces(planar_geometry))
 
-    face_vertices = oriented_bounding_box(stack(planar_geometry.vertices))
+    face_vertices, _, _ = oriented_bounding_box(stack(planar_geometry.vertices))
 
     # Vectors spanning the face
     edge1 = face_vertices[:, 2] - face_vertices[:, 1]
@@ -122,11 +122,247 @@ function oriented_bounding_box(point_cloud)
 
     aligned_coords = eigen_vectors' * centered_data
 
-    min_corner = minimum(aligned_coords, dims=2)
-    max_corner = maximum(aligned_coords, dims=2)
+    min_corner = minimum(aligned_coords, dims=2) .- eps()
+    max_corner = maximum(aligned_coords, dims=2) .+ eps()
 
-    face_vertices = hcat(min_corner, max_corner,
-                         [min_corner[1], max_corner[2], min_corner[3]])
+    if length(min_corner) == 2
+        rect_coords = hcat(min_corner, max_corner, [min_corner[1], max_corner[2]])
+    else
+        rect_coords = hcat(min_corner, max_corner,
+                           [min_corner[1], max_corner[2], min_corner[3]])
+    end
+
+    rotated_rect_coords = eigen_vectors * rect_coords .+ means
+
+    return rotated_rect_coords, eigen_vectors, (min_corner, max_corner)
+end
+
+"""
+    OrientedBoundingBox(; box_origin, orientation_matrix, edge_lengths, local_axis_scale::Tuple))
+    OrientedBoundingBox(point_cloud::AbstractMatrix; local_axis_scale::Tuple)
+    OrientedBoundingBox(geometry::Union{Polygon, TriangleMesh}; local_axis_scale::Tuple)
+
+Creates an oriented bounding box (rectangle in 2D or cuboid in 3D) that can be
+rotated and positioned arbitrarily in space.
+
+The box is defined either by explicit parameters
+or by automatically fitting it around an existing geometry or a point cloud with optional scaling.
+
+# Arguments
+- `point_cloud`: An array where the ``i``-th column holds the coordinates of a point ``i``.
+- `geometry`: Geometry returned by [`load_geometry`](@ref).
+
+# Keywords
+- `box_origin`: The corner point from which the box is constructed.
+- `orientation_matrix`: A matrix defining the orientation of the box in space.
+    Each column of the matrix represents one of the local axes of the box (e.g., x, y, z in 3D).
+    The matrix must be orthogonal, ensuring that the local axes are perpendicular to each other
+    and normalized.
+- `edge_lengths`: The lengths of the edges of the box:
+    - In 2D: `(width, height)`
+    - In 3D: `(width, height, depth)`
+- `local_axis_scale`: Allows for anisotropic scaling along the oriented axes of the `OrientedBoundingBox`
+    (the eigenvectors of the geometry's covariance matrix). Default is no scaling.
+    The tuple components correspond to:
+    - first element: scaling along the first eigenvector (local x-axis),
+    - second element: scaling along the second eigenvector (local y-axis),
+    - third element (only in 3D): scaling along the third eigenvector (local z-axis).
+
+    **Note:** Scaling is always applied in the local `OrientedBoundingBox`
+    coordinate system, i.e. along its oriented axes.
+    Scaling along arbitrary world directions is not supported,
+    as this would break the orthogonality of the spanning vectors.
+
+# Examples
+```jldoctest; output=false, setup = :(using LinearAlgebra: I)
+# 2D axis aligned
+OrientedBoundingBox(box_origin=[0.0, 0.0], orientation_matrix=I(2), edge_lengths=[2.0, 1.0])
+
+# 2D rotated
+orientation_matrix = [cos(π/4)  -sin(π/4);
+                      sin(π/4)   cos(π/4)]
+OrientedBoundingBox(; box_origin=[0.0, 0.0], orientation_matrix, edge_lengths=[2.0, 1.0])
+
+# 2D point cloud
+# Create a point cloud from a spherical shape (quarter circle sector)
+shape = SphereShape(0.1, 1.5, (0.2, 0.4), 1.0, n_layers=4,
+                    sphere_type=RoundSphere(; start_angle=0, end_angle=π/4))
+OrientedBoundingBox(shape.coordinates)
+
+# 3D rotated
+orientation_matrix = [cos(π/4)  -sin(π/4)  0.0;  # x-axis after rotation
+                      sin(π/4)   cos(π/4)  0.0;  # y-axis after rotation
+                      0.0        0.0       1.0]  # z-axis remains unchanged
+OrientedBoundingBox(; box_origin=[0.5, -0.2, 0.0], orientation_matrix,
+                    edge_lengths=[1.0, 2.0, 3.0])
+
+# output
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ OrientedBoundingBox (3D)                                                                         │
+│ ════════════════════════                                                                         │
+│ box origin: ………………………………………………… [0.5, -0.2, 0.0]                                                 │
+│ edge lengths: …………………………………………… (1.0, 2.0, 3.0)                                                  │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+```
+"""
+struct OrientedBoundingBox{NDIMS, ELTYPE <: Real, SV}
+    box_origin       :: SVector{NDIMS, ELTYPE}
+    spanning_vectors :: SV
+end
+
+function OrientedBoundingBox(; box_origin, orientation_matrix, edge_lengths,
+                             local_axis_scale::Tuple=ntuple(_ -> 0, length(box_origin)))
+    box_origin_ = SVector(box_origin...)
+    NDIMS = length(box_origin_)
+
+    @assert size(orientation_matrix) == (NDIMS, NDIMS)
+    @assert length(edge_lengths) == NDIMS
+    @assert length(local_axis_scale) == NDIMS
+
+    spanning_vectors = ntuple(i -> SVector{NDIMS}(orientation_matrix[:, i] *
+                                                  edge_lengths[i]), NDIMS)
+
+    prod(local_axis_scale) > 0 || return OrientedBoundingBox(box_origin_, spanning_vectors)
+
+    return scaled_bbox(SVector(local_axis_scale), box_origin_, spanning_vectors)
+end
+
+function OrientedBoundingBox(point_cloud::AbstractMatrix;
+                             local_axis_scale::Tuple=ntuple(_ -> 0, size(point_cloud, 1)))
+    NDIMS = size(point_cloud, 1)
+    vertices, eigen_vectors, (min_corner, max_corner) = oriented_bounding_box(point_cloud)
+
+    # Use the first vertex as box origin (bottom-left-front corner)
+    box_origin = SVector{NDIMS}(vertices[:, 1])
+
+    # Calculate edge lengths from min/max corners
+    edge_lengths = max_corner - min_corner
+
+    return OrientedBoundingBox(; box_origin, orientation_matrix=eigen_vectors, edge_lengths,
+                               local_axis_scale)
+end
+
+function OrientedBoundingBox(geometry::Union{Polygon, TriangleMesh};
+                             local_axis_scale::Tuple=ntuple(_ -> 0, ndims(geometry)))
+    point_cloud = stack(geometry.vertices)
+
+    return OrientedBoundingBox(point_cloud; local_axis_scale)
+end
+
+@inline Base.ndims(::OrientedBoundingBox{NDIMS}) where {NDIMS} = NDIMS
+
+function Base.show(io::IO, ::MIME"text/plain", box::OrientedBoundingBox)
+    @nospecialize box # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, box)
+    else
+        summary_header(io, "OrientedBoundingBox ($(ndims(box))D)")
+        summary_line(io, "box origin", box.box_origin)
+        summary_line(io, "edge lengths", norm.(box.spanning_vectors))
+        summary_footer(io)
+    end
+end
+
+function scaled_bbox(local_axis_scale::SVector{2}, box_origin, spanning_vectors)
+    # Uniform scaling about the center, center remains unchanged
+    v1, v2 = spanning_vectors
+    center = box_origin + (v1 + v2) / 2
+
+    # Scaling factor per oriented axis
+    s1, s2 = local_axis_scale
+
+    # New spanning vectors
+    v1p, v2p = s1 * v1, s2 * v2
+
+    new_origin = center - (v1p + v2p) / 2
+
+    return OrientedBoundingBox(new_origin, (v1p, v2p))
+end
+
+function scaled_bbox(local_axis_scale::SVector{3}, box_origin, spanning_vectors)
+    # Uniform scaling about the center, center remains unchanged
+    v1, v2, v3 = spanning_vectors
+    center = box_origin + (v1 + v2 + v3) / 2
+
+    # Scaling factor per oriented axis
+    s1, s2, s3 = local_axis_scale
+
+    # New spanning vectors
+    v1p, v2p, v3p = s1 * v1, s2 * v2, s3 * v3
+
+    new_origin = center - (v1p + v2p + v3p) / 2
+
+    return OrientedBoundingBox(new_origin, (v1p, v2p, v3p))
+end
+
+function is_in_oriented_box(coordinates::AbstractArray, box)
+    is_in_box = fill(false, size(coordinates, 2))
+    @threaded default_backend(coordinates) for particle in axes(coordinates, 2)
+        particle_coords = current_coords(coordinates, box, particle)
+        is_in_box[particle] = is_in_oriented_box(particle_coords, box)
+    end
+
+    return findall(is_in_box)
+end
+
+@inline function is_in_oriented_box(particle_coords::SVector{NDIMS}, box) where {NDIMS}
+    (; spanning_vectors, box_origin) = box
+    relative_position = particle_coords - box_origin
+
+    for dim in 1:NDIMS
+        span_dim = spanning_vectors[dim]
+        # Checks whether the projection of the particle position
+        # falls within the range of the zone
+        if !(0 <= dot(relative_position, span_dim) <= dot(span_dim, span_dim))
+
+            # Particle is not in box
+            return false
+        end
+    end
+
+    # Particle is in box
+    return true
+end
+
+@inline function Base.intersect(initial_condition::InitialCondition,
+                                boxes::OrientedBoundingBox...)
+    (; coordinates, density, mass, velocity, pressure, particle_spacing) = initial_condition
+    box = first(boxes)
+
+    if ndims(box) != ndims(initial_condition)
+        throw(ArgumentError("all passed `OrientedBoundingBox`s must have the same dimensionality as the initial condition"))
+    end
+
+    keep_indices = is_in_oriented_box(coordinates, box)
+
+    result = InitialCondition(; coordinates=coordinates[:, keep_indices],
+                              density=density[keep_indices], mass=mass[keep_indices],
+                              velocity=velocity[:, keep_indices],
+                              pressure=pressure[keep_indices],
+                              particle_spacing)
+
+    return intersect(result, Base.tail(boxes)...)
+end
+
+@inline function Base.setdiff(initial_condition::InitialCondition,
+                              boxes::OrientedBoundingBox...)
+    (; coordinates, density, mass, velocity, pressure, particle_spacing) = initial_condition
+    box = first(boxes)
+
+    if ndims(box) != ndims(initial_condition)
+        throw(ArgumentError("all passed `OrientedBoundingBox`s must have the same dimensionality as the initial condition"))
+    end
+
+    delete_indices = is_in_oriented_box(coordinates, box)
+    keep_indices = setdiff(eachparticle(initial_condition), delete_indices)
+
+    result = InitialCondition(; coordinates=coordinates[:, keep_indices],
+                              density=density[keep_indices],
+                              mass=mass[keep_indices],
+                              velocity=velocity[:, keep_indices],
+                              pressure=pressure[keep_indices],
+                              particle_spacing)
 
-    return eigen_vectors * face_vertices .+ means
+    return setdiff(result, Base.tail(boxes)...)
 end

src/visualization/recipes_plots.jl

@@ -170,3 +170,41 @@ RecipesBase.@recipe function f(geometry::Polygon)
         return (x, y)
     end
 end
+
+RecipesBase.@recipe function f(box::OrientedBoundingBox{2})
+    # Get the box origin and spanning vectors
+    origin = box.box_origin
+    v1, v2 = box.spanning_vectors
+
+    # Calculate all 4 corners of the rectangle
+    # Corner 1: origin
+    # Corner 2: origin + v1
+    # Corner 3: origin + v1 + v2 (diagonal corner)
+    # Corner 4: origin + v2
+    corners = hcat(origin,
+                   origin + v1,
+                   origin + v1 + v2,
+                   origin + v2,
+                   origin)  # Close the rectangle by returning to start
+
+    x = corners[1, :]
+    y = corners[2, :]
+
+    aspect_ratio --> :equal
+    grid --> false
+
+    # First plot the vertices as a scatter plot
+    @series begin
+        seriestype --> :scatter
+        return (x, y)
+    end
+
+    # Now plot the edges as a line plot
+    @series begin
+        # Ignore series in legend and color cycling. Note that `:=` forces the attribute,
+        # whereas `-->` would only set it if it is not already set.
+        primary := false
+
+        return (x, y)
+    end
+end