Skip to content

Commit

Permalink
Improve oriented bounding box again
Browse files Browse the repository at this point in the history
  • Loading branch information
sebastiangrimberg committed Jun 26, 2024
1 parent 8e385c8 commit 36cb377
Showing 1 changed file with 36 additions and 19 deletions.
55 changes: 36 additions & 19 deletions palace/utils/geodata.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1004,39 +1004,56 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm,
// Compute the center as halfway along the main diagonal.
Vector3dMap(box.center.data()) = 0.5 * (v_000 + v_111);

// Compute the box axes. using the 4 extremal points, we find the first two axes as the
// Compute the box axes. Using the 4 extremal points, we find the first two axes as the
// edges which are closest to perpendicular. For a perfect rectangular prism point
// cloud, we could instead compute the axes and length in each direction using the
// found edges of the cuboid, but this does not work for non-rectangular cross-sections
// or pyramid shapes
// found edges of the cuboid, but this does not work for non-rectangular prism
// cross-sections or pyramid shapes.
{
std::array<const Eigen::Vector3d *, 4> verts = {&v_000, &v_001, &v_011, &v_111};
Eigen::Vector3d e_0 = Eigen::Vector3d::Zero(), e_1 = Eigen::Vector3d::Zero();
double dot_min = mfem::infinity();
for (int i = 0; i < 4; i++)
bool found_orthog = false;
for (int i_0 = 0; i_0 < 4; i_0++)
{
for (int j = 0; j < 4; j++)
for (int j_0 = i_0 + 1; j_0 < 4; j_0++)
{
if (j == i)
for (int i_1 = 0; i_1 < 4; i_1++)
{
continue;
}
for (int k = j + 1; k < 4; k++)
{
if (k == i)
for (int j_1 = i_1 + 1; j_1 < 4; j_1++)
{
continue;
if (i_1 == i_0 && j_1 == j_0)
{
continue;
}
const auto e_ij_0 = (*verts[j_0] - *verts[i_0]).normalized();
const auto e_ij_1 = (*verts[j_1] - *verts[i_1]).normalized();
const auto dot = std::abs(e_ij_0.dot(e_ij_1));
if (dot < dot_min)
{
dot_min = dot;
e_0 = e_ij_0;
e_1 = e_ij_1;
if (dot_min < rel_tol)
{
found_orthog = true;
break;
}
}
}
const auto e_ij = (*verts[j] - *verts[i]).normalized();
const auto e_ik = (*verts[k] - *verts[i]).normalized();
const auto dot = std::abs(e_ij.dot(e_ik));
if (dot < dot_min)
if (found_orthog)
{
dot_min = dot;
e_0 = e_ij;
e_1 = e_ik;
break;
}
}
if (found_orthog)
{
break;
}
}
if (found_orthog)
{
break;
}
}
Vector3dMap(box.axes[0].data()) = e_0;
Expand Down

0 comments on commit 36cb377

Please sign in to comment.