Allows use of polygonal domain boundaries for Voronoi diagrams

examples/CMakeLists.txt

@@ -411,6 +411,9 @@ target_link_libraries(triangle_intersect OpenGL::GL glfw Freetype::Freetype)
 add_executable(voronoi_random voronoi_random.cpp)
 target_link_libraries(voronoi_random OpenGL::GL glfw Freetype::Freetype)
 
+add_executable(voronoi_boundary voronoi_boundary.cpp)
+target_link_libraries(voronoi_boundary OpenGL::GL glfw Freetype::Freetype)
+
 add_executable(voronoi_rectangular voronoi_rectangular.cpp)
 target_link_libraries(voronoi_rectangular OpenGL::GL glfw Freetype::Freetype)
 
@@ -558,3 +561,12 @@ target_link_libraries(show_boundingboxes OpenGL::GL glfw Freetype::Freetype)
 # if have compound-ray header
 add_executable(cray_eye cray_eye.cpp)
 target_link_libraries(cray_eye OpenGL::GL glfw Freetype::Freetype)
+
+if(ARMADILLO_FOUND)
+  # Make use of principle component analysis from arma in this example
+  add_executable(ellipse_pca ellipse_pca.cpp)
+  target_link_libraries(ellipse_pca ${ARMADILLO_LIBRARY} ${ARMADILLO_LIBRARIES} OpenGL::GL glfw Freetype::Freetype)
+endif()
+
+add_executable(trace_boundary trace_boundary.cpp)
+target_link_libraries(trace_boundary OpenGL::GL glfw Freetype::Freetype)

examples/ellipse_pca.cpp

@@ -0,0 +1,81 @@
+/*
+ * Applying principle component analysis (using arma) to some random data
+ */
+#include <iostream>
+
+#include <sm/mathconst>
+#include <sm/random>
+#include <sm/mat22>
+
+#include <armadillo>
+
+#include <mplot/Visual.h>
+#include <mplot/GraphVisual.h>
+
+int main()
+{
+    constexpr unsigned int n_samp = 100;
+
+    // Create data: Draw some samples from 2D Gaussian and rotate
+    sm::rand_normal<float> rn1 (0.0f, 2.0f);
+    sm::rand_normal<float> rn2 (0.0f, 0.5f);
+    sm::vvec<sm::vec<float, 2>> _x (n_samp, {0});
+    sm::mat22<float> rotn;
+    rotn.rotate (sm::mathconst<float>::pi_over_8);
+    for (unsigned int i = 0; i < n_samp; ++i) {
+        _x[i] = rotn * sm::vec<float, 2>{ rn1.get(), rn2.get() };
+    }
+
+    // Graph the data
+    mplot::Visual v(1024, 768, "Principle component analysis with armadillo");
+    // Create a GraphVisual object (obtaining a unique_ptr to the object) with a spatial offset within the scene of 0,0,0
+    auto gv = std::make_unique<mplot::GraphVisual<float>> (sm::vec<float>{-0.5f,-0.5f,0.0f});
+    mplot::DatasetStyle ds (mplot::stylepolicy::markers);
+    ds.datalabel = std::string("data");
+    v.bindmodel (gv);
+    gv->setlimits (-8, 8, -8, 8);
+    gv->setdata (_x, ds);
+
+    std::cout << "\narma gives:\n";
+    // Place data in arma::Mat
+    arma::Mat<float> x(_x.size(), 2);
+    for (unsigned int i = 0; i < _x.size(); ++i) {
+        x(i, 0) = _x[i][0];
+        x(i, 1) = _x[i][1];
+    }
+
+    // Apply arma::princomp
+    arma::Mat<float> co, sc;
+    arma::Col<float> lat, tsq;
+    arma::princomp (co, sc, lat, tsq, x);
+
+    std::cout << "coeff: " << co << std::endl;
+    //std::cout << "scores: " << sc << std::endl;
+    std::cout << "latent: " << lat << std::endl;
+
+    // Mat access is (r, c)
+    sm::vec<float, 2> pc1vec = { co(0, 0), co(1, 0) };
+    float angle1 = pc1vec.angle();
+    std::cout << "Angle of first component " << pc1vec << " is " << angle1 * sm::mathconst<float>::rad2deg
+              << " and length is " << lat(0) << std::endl;
+    sm::vec<float, 2> pc2vec = { co(0, 1), co(1, 1) };
+    float angle2 = pc2vec.angle();
+    std::cout << "Angle of 2nd component " << pc2vec << " is " << angle2 * sm::mathconst<float>::rad2deg
+              << " and length is " << lat(1) << std::endl;
+
+    mplot::DatasetStyle ds2 (mplot::stylepolicy::lines);
+    // Show 3 sigma axes
+    ds2.setcolour (mplot::colour::crimson);
+    ds2.datalabel = std::string("PC1 3") + mplot::unicode::toUtf8 (mplot::unicode::sigma);
+    sm::vvec<sm::vec<float, 2>> pc1vv = { {0,0}, pc1vec * 3.0f * std::sqrt(lat(0))};
+    gv->setdata (pc1vv, ds2);
+
+    ds2.setcolour (mplot::colour::orange);
+    ds2.datalabel = std::string("PC2 3") + mplot::unicode::toUtf8 (mplot::unicode::sigma);
+    sm::vvec<sm::vec<float, 2>> pc2vv = { {0,0}, pc2vec * 3.0f * std::sqrt(lat(1)) };
+    gv->setdata (pc2vv, ds2);
+
+    gv->finalize();
+    v.addVisualModel (gv);
+    v.keepOpen();
+}

examples/trace_boundary.cpp

@@ -0,0 +1,47 @@
+/*
+ * Trace around a boundary of points in 2D. This code uses Graham's scan to compute the convex hull
+ * of a set of randomly generated points.
+ */
+#include <iostream>
+
+#include <sm/mathconst>
+#include <sm/random>
+#include <sm/vvec>
+#include <sm/vec>
+#include <sm/algo>
+
+#include <mplot/Visual.h>
+#include <mplot/GraphVisual.h>
+
+int main()
+{
+    // Get ready to graph the data
+    mplot::Visual v(1024, 768, "Graham's scan computes the convex hull");
+    mplot::DatasetStyle ds (mplot::stylepolicy::markers);
+    auto gv = std::make_unique<mplot::GraphVisual<float>> (sm::vec<float>{-0.5f,-0.5f,0.0f});
+    v.bindmodel (gv);
+
+    constexpr unsigned int n_points = 200;
+
+    // Create data
+    sm::rand_uniform<float> rngxy(-0.8f, 0.8f);
+    sm::vvec<sm::vec<float, 2>> points(n_points);
+    for (unsigned int i = 0; i < n_points; ++i) { points[i] = { rngxy.get(), rngxy.get() }; }
+
+    // Trace algorithm
+    sm::vvec<sm::vec<float, 2>> boundary = sm::geometry::graham_scan (points);
+    // Close the boundary
+    boundary.push_back (boundary.front());
+
+    gv->setlimits (-1, 1, -1, 1);
+
+    gv->setdata (points, ds);
+
+    mplot::DatasetStyle ds2 (mplot::stylepolicy::lines);
+    ds2.linecolour = mplot::colour::crimson;
+    gv->setdata (boundary, ds2);
+
+    gv->finalize();
+    v.addVisualModel (gv);
+    v.keepOpen();
+}

examples/voronoi_boundary.cpp

@@ -0,0 +1,75 @@
+/*
+ * This example generates a number (n_points) of random (but bounded) coordinates and
+ * uses the VoronoiVisual to display the coordinates as a map, with the order of
+ * random-choice being used to colourize the Voronoi cells.
+ *
+ * Test harness for drawing an arbitrary boundary around the Voronoi region
+ *
+ * Author Seb James
+ * Date 2026
+ */
+#include <iostream>
+#include <sm/vec>
+#include <sm/random>
+#include <mplot/Visual.h>
+#include <mplot/VoronoiVisual.h>
+
+int main (int argc, char** argv)
+{
+    int rtn = -1;
+
+    float border_width = 0.2f;
+    unsigned int n_points = 80;
+    int domshape = 0; // 0 for rectangular, 1 for circ, 2 for ellipse, 3 for traced
+
+    if (argc > 1) { domshape = std::stoi (argv[1]); }
+    if (argc > 2) { border_width = std::stof (argv[2]); }
+    if (argc > 3) { n_points = std::stoi (argv[3]); }
+
+    mplot::Visual v(1024, 768, "VoronoiVisual");
+
+    sm::rand_uniform<float> rngxy(-0.8f, 0.8f, n_points);
+    sm::rand_uniform<float> rngz(0.8f, 1.0f, n_points);
+
+    // make n_points random coordinates
+    std::vector<sm::vec<float>> points(n_points);
+    std::vector<float> data(n_points);
+
+    for (unsigned int i = 0; i < n_points; ++i) {
+        points[i] = { rngxy.get(), rngxy.get(), rngz.get() };
+        data[i] = static_cast<float>(i) / n_points;
+    }
+
+    mplot::ColourMapType cmap_t = mplot::ColourMapType::Plasma;
+
+    sm::vec<float, 3> offset = { 0.0f };
+    auto vorv = std::make_unique<mplot::VoronoiVisual<float>> (offset);
+    v.bindmodel (vorv);
+    vorv->show_voronoi2d = true; // true to show the 2D voronoi edges
+    vorv->debug_dataCoords = true; // true to show coordinate spheres
+
+    // traced, circular, ellipsoid, rectangular:
+    if (domshape == 1) {
+        vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::circular;
+    } else if (domshape == 2) {
+        vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::ellipsoid;
+    } else if (domshape == 3) {
+        vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::traced;
+    } else {
+        vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::rectangular;
+    }
+    // A border to add to our domain boundary
+    vorv->border_width = border_width;
+    // Some domain shapes need to know a number of points:
+    vorv->num_boundary_points = 12;
+
+    vorv->cm.setType (cmap_t);
+    vorv->setDataCoords (&points);
+    vorv->setScalarData (&data);
+    vorv->finalize();
+    v.addVisualModel (vorv);
+
+    v.keepOpen();
+
+    return rtn;
+}

examples/voronoi_function.cpp

@@ -44,7 +44,10 @@ int main()
     vorv->show_voronoi2d = false; // true to show the 2D voronoi edges
     vorv->debug_dataCoords = false; // true to show coordinate spheres
     float length_scale = 4.0f / std::sqrt (n_points);
-    vorv->border_width  = length_scale;
+    vorv->border_width = length_scale;
+    //vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::traced;
+    //vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::circular;
+    //vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::rectangular; // default
     vorv->cm.setType (cmap_t);
     vorv->setDataCoords (&points);
     vorv->setScalarData (&data);

examples/voronoi_function_flat.cpp

@@ -43,7 +43,10 @@ int main()
     vorv->show_voronoi2d = true; // true to show the 2D voronoi edges
     vorv->debug_dataCoords = false; // true to show coordinate spheres
     float length_scale = 4.0f / std::sqrt (n_points);
-    vorv->border_width  = length_scale;
+    vorv->border_width = length_scale;
+    //vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::traced;
+    //vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::circular;
+    //vorv->dom_shape = mplot::VoronoiVisual<float>::domain_shape::rectangular; // default
     vorv->cm.setType (cmap_t);
     vorv->setDataCoords (&points);
     vorv->setScalarData (&data);
@@ -61,7 +64,7 @@ int main()
         if (fcount++% 600 == 0) {
             vorvp->cm.setType (++cmap_t);
         }
-        vorvp->reinitColours(); // Not quite working when I change the colourmap
+        vorvp->reinitColours();
         v.waitevents(0.018);
         v.render();
         k += 0.01f;

mplot/VisualDataModel.h

@@ -8,6 +8,7 @@
 #include <sm/vec>
 #include <sm/vvec>
 #include <sm/scale>
+#include <sm/centroid>
 #include <mplot/VisualModel.h>
 #include <mplot/ColourMap.h>
 
@@ -141,6 +142,8 @@ namespace mplot
             this->reinit();
         }
 
+        sm::vec<float> coordsCentroid() const { return sm::algo::centroid (*this->dataCoords); }
+
         //! An overridable function to set the colour of rect ri
         std::array<float, 3> setColour (uint64_t ri)
         {