From 317e045e6c9ca844a49d35e4ca1304418136b4fa Mon Sep 17 00:00:00 2001
From: LasNikas <niklas.neher@web.de>
Date: Wed, 7 Aug 2024 15:17:04 +0200
Subject: [PATCH] add oriented bounding box

---
 src/TrixiParticles.jl                         |  1 +
 .../boundary/open_boundary/boundary_zones.jl  | 40 +++++++++++++++++++
 2 files changed, 41 insertions(+)

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 637a8bc08..3ef70b531 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -23,6 +23,7 @@ using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!
                  get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
+using Statistics: Statistics
 using StrideArrays: PtrArray, StaticInt
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
 using TrixiBase: trixi_include, @trixi_timeit, timer, timeit_debug_enabled,
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
index 586d565b1..fb095d7c2 100644
--- a/src/schemes/boundary/open_boundary/boundary_zones.jl
+++ b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -275,6 +275,22 @@ function calculate_spanning_vectors(plane, zone_width)
     return spanning_vectors(Tuple(plane), zone_width), SVector(plane[1]...)
 end
 
+function calculate_spanning_vectors(plane::TriangleMesh, zone_width)
+    plane_normal = normalize(sum(plane.face_normals) / nfaces(plane))
+
+    plane_points = oriented_bounding_box(stack(plane.vertices))
+
+    # Vectors spanning the plane
+    edge1 = plane_points[:, 2] - plane_points[:, 1]
+    edge2 = plane_points[:, 3] - plane_points[:, 1]
+
+    if !isapprox(abs.(normalize(cross(edge2, edge1))), abs.(plane_normal), atol=1e-2)
+        throw(ArgumentError("`plane` might be not planar"))
+    end
+
+    return hcat(plane_normal * zone_width, edge1, edge2), SVector(plane_points[:, 1]...)
+end
+
 function spanning_vectors(plane_points::NTuple{2}, zone_width)
     plane_size = plane_points[2] - plane_points[1]
 
@@ -338,3 +354,27 @@ function remove_outside_particles(initial_condition, spanning_set, zone_origin)
     return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density),
                             particle_spacing)
 end
+
+# According to:
+# https://logicatcore.github.io/scratchpad/lidar/sensor-fusion/jupyter/2021/04/20/3D-Oriented-Bounding-Box.html
+function oriented_bounding_box(point_cloud)
+    covariance_matrix = Statistics.cov(point_cloud; dims=2)
+    eigen_vectors = Statistics.eigvecs(covariance_matrix)
+    means = Statistics.mean(point_cloud, dims=2)
+
+    centered_data = point_cloud .- means
+
+    aligned_coords = eigen_vectors' * centered_data
+
+    min_corner = SVector(minimum(aligned_coords[1, :]),
+                         minimum(aligned_coords[2, :]),
+                         minimum(aligned_coords[3, :]))
+    max_corner = SVector(maximum(aligned_coords[1, :]),
+                         maximum(aligned_coords[2, :]),
+                         maximum(aligned_coords[3, :]))
+
+    plane_points = hcat(min_corner, max_corner,
+                        [min_corner[1], max_corner[2], min_corner[3]])
+
+    return eigen_vectors * plane_points .+ means
+end