Implement convex hull for unsupported loudspeaker layouts

docs/conf.py

@@ -78,6 +78,9 @@
     html_theme = "sphinx_rtd_theme"
     html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
+# configure sphinxcontrib.bibtex
+bibtex_bibfiles = ['refs.bib']
+
 # make a pybtex style which uses the bibtex key for labels
 
 from pybtex.style.formatting.unsrt import Style as UnsrtStyle

docs/loudspeaker_layouts.rst

@@ -6,7 +6,19 @@ Loudspeaker Layouts
 An output loudspeaker layout is represented by :class:`Layout`, which contains
 a name, a list of :class:`Channel` objects, and the reference :type:`Screen`.
 
-Although :class:`Layout` objects could be constructed directly from scratch,
-only layouts listed in :cite:`bs2051` are supported. The function
-:func:`getLayout` is therefore provided to obtain a :class:`Layout` object
-given the layout name, e.g. ``4+5+0``.
+Loudspeaker layouts specified in :cite:`bs2051` are supported, including
+positions within the ranges specified. The function :func:`getLayout` is
+therefore provided to obtain a :class:`Layout` object given the layout name,
+e.g. ``4+5+0``.
+
+When using layouts with non-nominal positions, the
+:func:`Channel::polarPositionNominal` must match the position specified in Table
+1 in :cite:`bs2051`, and the :func:`Channel::polarPosition` must meet the
+specified constraints, including azimuth/elevation ranges and symmetry
+requirements.
+
+Non-standard loudspeaker layouts may be used, however their behaviour may
+change in future versions of the library. Loudspeaker layouts must be similar
+to those in :cite:`bs2051`, with 1, 2 or 3 layers and an optional T+000 or
+UH+180 loudspeaker. Using this functionality requires some understanding of the
+internals of the renderer, particularly section 6.1.3.1 of :cite:`bs2127`.

src/CMakeLists.txt

@@ -13,6 +13,7 @@ add_custom_command(
 add_library(ear
   bs2051.cpp
   bs2051_layouts.cpp
+  common/convex_hull.cpp
   common/facets.cpp
   common/geom.cpp
   common/point_source_panner.cpp

src/common/convex_hull.cpp

@@ -0,0 +1,102 @@
+#include "convex_hull.hpp"
+#include <Eigen/Geometry>
+#include "ear/helpers/assert.hpp"
+
+namespace ear {
+
+  std::vector<Facet> convex_hull(const std::vector<Eigen::Vector3d> &positions,
+                                 double tolerance) {
+    using Vec = Eigen::Vector3d;
+    using Tri = std::array<size_t, 3>;
+
+    // the mean point, guaranteed to be inside the convex hull
+    Vec inside_point = Vec::Zero();
+    for (auto &pos : positions) inside_point += pos;
+    inside_point /= positions.size();
+
+    // find all triangles on the convex hull:
+    // - iterate through all possible triangles
+    // - find the triangle normal, checking for collinearity
+    // - it's on the convex hull if:
+    //    - the inside point is not on the plane of the triangle
+    //    - all points are on the same side of the plane of the triangle as the
+    //    inside point
+    std::vector<Tri> hull_tris;
+    std::vector<Vec> hull_tri_normals;
+
+    for (size_t i = 0; i < positions.size(); i++)
+      for (size_t j = i + 1; j < positions.size(); j++)
+        for (size_t k = j + 1; k < positions.size(); k++) {
+          Vec normal =
+              (positions[j] - positions[i]).cross(positions[k] - positions[i]);
+          ear_assert(normal.squaredNorm() > tolerance,
+                     "collinear points in convex hull");
+          normal.normalize();  // XXX: not really required, but helps make
+                               // tolerances consistent
+
+          double dot_inside = normal.dot(inside_point - positions[i]);
+
+          if (std::abs(dot_inside) < tolerance)
+            continue;  // tri coplanar with inside point
+
+          bool points_on_same_side_as_inside = true;
+          for (auto &pos : positions) {
+            double dot_point = normal.dot(pos - positions[i]);
+
+            // check if signs are equal, with tolerance
+            if (!(dot_inside > 0 ? dot_point > -tolerance
+                                 : dot_point < tolerance)) {
+              points_on_same_side_as_inside = false;
+              break;
+            }
+          }
+
+          if (points_on_same_side_as_inside) {
+            hull_tris.push_back({i, j, k});
+            hull_tri_normals.push_back(normal);
+          }
+        }
+
+    // merge coplanar triangles into facets
+    std::vector<Facet> hull_facets;
+    std::vector<Vec> hull_facet_normals;
+
+    for (size_t tri_i = 0; tri_i < hull_tris.size(); tri_i++) {
+      const Vec &tri_point = positions[hull_tris[tri_i][0]];
+      const Vec &tri_norm = hull_tri_normals[tri_i];
+
+      // check each facet to see if this triangle is on the same plane.
+      // if it is, add the points to the facet. if no facet is found, make a
+      // new facet.
+      bool found_facet = false;
+      for (size_t facet_i = 0; facet_i < hull_facets.size(); facet_i++) {
+        const Vec &facet_norm = hull_facet_normals[facet_i];
+        const Vec &facet_point = positions[*hull_facets[facet_i].begin()];
+
+        // they are on the same plane if the line between the test points is
+        // not along the normal, and if the normals point in the same or
+        // opposite direction. alternatively, we could just test all tri
+        // vertices.
+        if (std::abs((tri_point - facet_point).dot(facet_norm)) < tolerance &&
+            facet_norm.cross(tri_norm).squaredNorm() < tolerance) {
+          for (size_t corner_idx : hull_tris[tri_i]) {
+            hull_facets[facet_i].insert((Facet::value_type)corner_idx);
+          }
+
+          found_facet = true;
+          break;
+        }
+      }
+
+      if (!found_facet) {
+        Facet f;
+        for (size_t corner_idx : hull_tris[tri_i])
+          f.insert((Facet::value_type)corner_idx);
+        hull_facets.push_back(std::move(f));
+        hull_facet_normals.push_back(tri_norm);
+      }
+    }
+
+    return hull_facets;
+  }
+}  // namespace ear

src/common/convex_hull.hpp

@@ -0,0 +1,17 @@
+#pragma once
+#include <Eigen/Core>
+#include "facets.hpp"
+
+namespace ear {
+  /// very simple convex hull implementation
+  ///
+  /// limitations:
+  /// - doesn't deal with degenerate cases (less than 4 points, or all points
+  ///   being coplanar)
+  /// - doesn't deal with collinear points at all (we only use points a unit
+  ///   distance from the origin, so this never occurs)
+  /// - doesn't check for valid geometry, so the tolerance has to be high
+  ///   enough
+  std::vector<Facet> convex_hull(const std::vector<Eigen::Vector3d> &positions,
+                                 double tolerance = 1e-5);
+}  // namespace ear

src/common/point_source_panner.cpp

@@ -3,6 +3,7 @@
 #include <Eigen/Dense>
 #include <boost/algorithm/clamp.hpp>
 #include <boost/make_unique.hpp>
+#include "convex_hull.hpp"
 #include "ear/bs2051.hpp"
 #include "ear/helpers/assert.hpp"
 #include "facets.hpp"
@@ -419,16 +420,18 @@ namespace ear {
     return std::make_shared<PolarPointSourcePanner>(std::move(regions));
   }
 
-  std::shared_ptr<PointSourcePanner> configureFullPolarPanner(
-      const Layout& layout) {
+  std::tuple<std::vector<Eigen::Vector3d>, std::vector<Eigen::Vector3d>,
+             std::set<int>, Eigen::MatrixXd>
+  getAugmentedLayout(const Layout& layout) {
     // add some extra height speakers that are treated as real speakers until
     // the downmix in PointSourcePannerDownmix
+    std::vector<Channel> allChannels(layout.channels());
+
     std::vector<Channel> extraChannels;
     Eigen::MatrixXd downmix;
-    Layout layoutExtra(layout);
     std::tie(extraChannels, downmix) = extraPosVerticalNominal(layout);
     for (const auto& extraChannel : extraChannels) {
-      layoutExtra.channels().push_back(extraChannel);
+      allChannels.push_back(extraChannel);
     }
 
     // add some virtual speakers above and below that will be used as the
@@ -447,25 +450,39 @@ namespace ear {
       virtualPositions.push_back(Eigen::Vector3d{0.0, 0.0, 1.0});
     }
 
-    std::vector<PolarPosition> positionsRealVec = layoutExtra.positions();
     std::vector<Eigen::Vector3d> positionsReal;
-    std::transform(positionsRealVec.begin(), positionsRealVec.end(),
-                   std::back_inserter(positionsReal),
-                   [](PolarPosition pos) -> Eigen::Vector3d {
-                     return toNormalisedVector3d(pos);
-                   });
-    positionsReal.insert(positionsReal.end(), virtualPositions.begin(),
-                         virtualPositions.end());
-
+    std::vector<Eigen::Vector3d> positionsNominal;
     std::set<int> virtualVerts;
-    virtualVerts.insert(static_cast<int>(layoutExtra.channels().size()));
-    if (virtualPositions.size() == 2) {
-      virtualVerts.insert(static_cast<int>(layoutExtra.channels().size() + 1));
+
+    for (const Channel& channel : allChannels) {
+      positionsReal.push_back(toNormalisedVector3d(channel.polarPosition()));
+      positionsNominal.push_back(
+          toNormalisedVector3d(channel.polarPositionNominal()));
+    }
+    for (const Eigen::Vector3d& pos : virtualPositions) {
+      virtualVerts.insert(static_cast<int>(positionsReal.size()));
+      positionsReal.push_back(pos);
+      positionsNominal.push_back(pos);
     }
 
+    return {positionsReal, positionsNominal, virtualVerts, downmix};
+  }
+
+  std::shared_ptr<PointSourcePanner> configureFullPolarPanner(
+      const Layout& layout) {
+    std::vector<Eigen::Vector3d> positionsReal;
+    std::vector<Eigen::Vector3d> positionsNominal;
+    std::set<int> virtualVerts;
+    Eigen::MatrixXd downmix;
+    std::tie(positionsReal, positionsNominal, virtualVerts, downmix) =
+        getAugmentedLayout(layout);
+
     // Facets of the convex hull; each set represents a facet and contains the
     // indices of its corners in positions.
-    std::vector<Facet> facets = FACETS.at(layout.name());
+    auto facets_it = FACETS.find(layout.name());
+    std::vector<Facet> facets = facets_it != FACETS.end()
+                                    ? facets_it->second
+                                    : convex_hull(positionsNominal);
 
     // Turn the facets into regions for the point source panner.
     std::vector<std::unique_ptr<RegionHandler>> regions;

src/common/point_source_panner.hpp

@@ -2,6 +2,7 @@
 #include <Eigen/Core>
 #include <memory>
 #include <numeric>
+#include <set>
 #include <string>
 #include "ear/layout.hpp"
 #include "ear/metadata.hpp"
@@ -174,6 +175,17 @@ namespace ear {
   std::pair<std::vector<Channel>, Eigen::MatrixXd> extraPosVerticalNominal(
       Layout layout);
 
+  // given a layout, determine the full set of loudspeaker positions used for
+  // panning. returns:
+  //
+  // - the real position of real and virtual loudspeakers
+  // - the nominal position of the real and virtual loudspeakers
+  // - the indices of the virtual loudspeakers in the two position lists
+  // - a downmix matrix to be applied to the output of the real loudspeakers
+  std::tuple<std::vector<Eigen::Vector3d>, std::vector<Eigen::Vector3d>,
+             std::set<int>, Eigen::MatrixXd>
+  getAugmentedLayout(const Layout& layout);
+
   std::shared_ptr<PointSourcePanner> configureFullPolarPanner(
       const Layout& layout);
 

tests/point_source_panner_tests.cpp

@@ -5,6 +5,8 @@
 #include <catch2/catch.hpp>
 #include <random>
 #include <string>
+#include "common/convex_hull.hpp"
+#include "common/facets.hpp"
 #include "common/geom.hpp"
 #include "common/point_source_panner.hpp"
 #include "ear/bs2051.hpp"
@@ -532,3 +534,40 @@ TEST_CASE("screen_loudspeaker_positions") {
     }
   }
 }
+
+TEST_CASE("hull") {
+  for (auto& layoutFull : loadLayouts()) {
+    if (layoutFull.name() == "0+2+0") continue;
+    auto layout = layoutFull.withoutLfe();
+
+    SECTION(layout.name()) {
+      std::vector<Eigen::Vector3d> positionsReal;
+      std::vector<Eigen::Vector3d> positionsNominal;
+      std::set<int> virtualVerts;
+      Eigen::MatrixXd downmix;
+      std::tie(positionsReal, positionsNominal, virtualVerts, downmix) =
+          getAugmentedLayout(layout);
+
+      std::vector<Facet> facets_precomputed = FACETS.at(layout.name());
+
+      // check that we're not close to a tolerance that is too big or small
+      for (double tolerance : {1e-6, 1e-5, 1e-4}) {
+        SECTION("tol = " + std::to_string(tolerance)) {
+          std::vector<Facet> facets_calculated =
+              convex_hull(positionsNominal, tolerance);
+
+          std::sort(facets_precomputed.begin(), facets_precomputed.end());
+          std::sort(facets_calculated.begin(), facets_calculated.end());
+
+          REQUIRE(facets_precomputed == facets_calculated);
+        }
+      }
+
+#ifdef CATCH_CONFIG_ENABLE_BENCHMARKING
+      if (layout.name() == "9+10+3") {
+        BENCHMARK("hull") { return convex_hull(positionsNominal); };
+      }
+#endif
+    }
+  }
+}