Polar plotting

examples/CMakeLists.txt

@@ -472,3 +472,14 @@ target_link_libraries(coordarrows OpenGL::GL glfw Freetype::Freetype)
 
 add_executable(sph_to_cart sph_to_cart.cpp)
 target_link_libraries(sph_to_cart OpenGL::GL glfw Freetype::Freetype)
+
+add_executable(polar polar.cpp)
+target_link_libraries(polar OpenGL::GL glfw Freetype::Freetype)
+
+if(HAVE_STD_FORMAT)
+  add_executable(zernike_radial zernike_radial.cpp)
+  target_link_libraries(zernike_radial OpenGL::GL glfw Freetype::Freetype)
+
+  add_executable(zernike zernike.cpp)
+  target_link_libraries(zernike OpenGL::GL glfw Freetype::Freetype)
+endif()

examples/polar.cpp

@@ -0,0 +1,55 @@
+/*
+ * A simple polar plot example.
+ */
+
+#include <sm/mathconst>
+#include <sm/vvec>
+#include <mplot/Visual.h>
+#include <mplot/PolarVisual.h>
+
+int main()
+{
+    using mc = sm::mathconst<double>;
+
+    mplot::Visual<> v (1024, 768, "Polar plot");
+
+    constexpr std::size_t N = 40;
+    sm::vvec<double> rho;
+    rho.linspace (0.0, 1.0, N);
+    sm::vvec<double> theta;
+    theta.linspace (0.0, mc::two_pi, N);
+
+    auto pv = std::make_unique<mplot::PolarVisual<double>> (sm::vec<float>{0.0f});
+    v.bindmodel (pv);
+    pv->cm.setType (mplot::ColourMapType::Cork);
+    pv->setFrameColour (mplot::colour::black);
+    pv->setTextColour (mplot::colour::black);
+    pv->radius = 1.0f;
+    pv->tf.fontsize = 0.08f;
+    pv->ticklabelgap = 0.05f;
+    pv->numrings = N;
+    pv->numsegs = N;
+    pv->twodimensional = false;
+
+    // Make a function that increases radially near theta == 0 and theta == pi/2 and is 0 everywhere else
+    sm::vvec<double> data;
+    for (auto rh : rho) {
+        for (auto th : theta) {
+            if (std::abs (th) < 0.1 || std::abs (th - mc::pi_over_2) < 0.1) {
+                data.push_back (rh);
+            } else {
+                data.push_back (0.0);
+            }
+        }
+    }
+
+    pv->setScalarData (&data);
+
+    pv->zScale.output_range = {-1.0f, 1.0f};
+    pv->zScale.compute_scaling (-2, 2);
+
+    pv->finalize();
+    v.addVisualModel (pv);
+
+    v.keepOpen();
+}

examples/zernike.cpp

@@ -0,0 +1,77 @@
+// Demo Zernike polynomials
+
+#include <format>
+#include <complex>
+#include <algorithm>
+#include <sm/mathconst>
+#include <sm/algo>
+#include <sm/vvec>
+#include <mplot/Visual.h>
+#include <mplot/PolarVisual.h>
+
+int main()
+{
+    // You can require that n - |m| is even to "ensure the rotational invariant property is met"
+    // (see https://doi.org/10.7554/eLife.54026, p 18)
+    constexpr bool require_n_minus_abs_m_even = false;
+
+    // Choose flat or 3D plots
+    constexpr bool flat_plots = true;
+
+    using mc = sm::mathconst<double>;
+
+    mplot::Visual<> v (1024, 768, "Zernike Polynomials");
+
+    constexpr std::size_t N = 100;
+    sm::vvec<double> rho;
+    rho.linspace (0.0, 1.0, N);
+
+    sm::vvec<double> theta;
+    theta.linspace (0.0, mc::two_pi, N);
+
+    for (unsigned int n = 0; n <= 16; ++n) {
+
+        for (int m = 0; m <= static_cast<int>(n); ++m) {
+
+            if constexpr (require_n_minus_abs_m_even) { if ((n - std::abs(m)) % 2 != 0) { continue; } }
+
+            auto pv = std::make_unique<mplot::PolarVisual<double>> (sm::vec<float>{1.25f * n, 1.25f * m, 0});
+            v.bindmodel (pv);
+            pv->cm.setType (mplot::ColourMapType::Cork);
+            pv->setFrameColour (mplot::colour::goldenrod1);
+            pv->setTextColour (mplot::colour::black);
+            pv->radius = 0.5f;
+            pv->tf.fontsize = 0.05f;
+            pv->numrings = N;
+            pv->numsegs = N;
+            pv->addLabel (std::format ("n{}, m{}", n, m), sm::vec<float>{0.2,-0.58,0}, mplot::TextFeatures(0.08f));
+            pv->twodimensional = false;
+
+            sm::vvec<double> Vnm_real;
+            for (auto rh : rho) {
+                double r_nm = std::clamp (sm::algo::zern_radial_poly (n, m, rh), -10.0, 10.0);
+                for (auto th : theta) {
+                    Vnm_real.push_back (std::clamp (std::real(sm::algo::zern_polynomial (m, r_nm, th)), -10.0, 10.0));
+                }
+            }
+
+            // Replace NaN with 0
+            for (auto& v : Vnm_real) { if (std::isnan (v)) { v = 0.0; } }
+
+            pv->setScalarData (&Vnm_real);
+
+            // Flat plot if necessary (above clamps and isnan testing should make this unnecessary)
+            if constexpr (flat_plots) {
+                pv->zScale.null_scaling();
+            } else {
+                pv->zScale.output_range = {-1.0f, 1.0f};
+                pv->zScale.compute_scaling (-10, 10);
+            }
+
+            pv->finalize();
+            v.addVisualModel (pv);
+        }
+    }
+
+    v.keepOpen();
+}

examples/zernike_radial.cpp

@@ -0,0 +1,68 @@
+// Demo Zernike polynomials
+
+#include <format>
+#include <sm/algo>
+#include <sm/vvec>
+#include <mplot/Visual.h>
+#include <mplot/GraphVisual.h>
+
+int main()
+{
+    mplot::Visual<> v (1024, 768, "Zernike Radial Polynomials (Rnm)");
+
+    sm::vvec<double> rho;
+    rho.linspace (0.001, 1, 200);
+
+    mplot::DatasetStyle ds (mplot::stylepolicy::lines);
+
+    // You can require that n - |m| is even to "ensure the rotational invariant property is met"
+    // (see https://doi.org/10.7554/eLife.54026, p 18)
+    constexpr bool require_n_minus_abs_m_even = false;
+
+    for (unsigned int n = 0; n < 10; ++n) {
+        sm::vec<float> offset = {1.4f * n, 2.4f, 0};
+        sm::vec<float> offset2 = {1.4f * n, 0, 0};
+        auto gv = std::make_unique<mplot::GraphVisual<double>> (offset);
+        v.bindmodel (gv);
+        gv->xlabel = "rho";
+        gv->ylabel = "Rnm";
+        gv->setlimits (0.0, 1.0, -10.0, 10.0);
+
+        auto gv2 = std::make_unique<mplot::GraphVisual<double>> (offset2);
+        v.bindmodel (gv2);
+        gv2->xlabel = "rho";
+        gv2->ylabel = "Rnm";
+        gv2->setlimits (0.0, 1.0, -10.0, 10.0);
+
+        bool have_data = false;
+
+        for (int m = 0; m <= static_cast<int>(n); ++m) {
+
+            if constexpr (require_n_minus_abs_m_even) {
+                if ((n - std::abs(m)) % 2 != 0) { continue; }
+            }
+
+            sm::vvec<double> Rnm;
+            sm::vvec<double> Rnm2;
+            for (auto rh : rho) {
+                Rnm.push_back (sm::algo::zern_radial_poly (n, m, rh));
+                Rnm2.push_back (sm::algo::zern_radial_poly (n, -m, rh));
+            }
+            ds.datalabel = std::format ("n{}, m{}", n, m);
+            ds.linecolour = mplot::DatasetStyle::datacolour (m + n);
+            gv->setdata (rho, Rnm, ds);
+            ds.datalabel = std::format ("n{}, m{}", n, -m);
+            gv2->setdata (rho, Rnm2, ds);
+            have_data = true;
+        }
+
+        if (have_data) {
+            gv->finalize();
+            gv2->finalize();
+            v.addVisualModel (gv);
+            v.addVisualModel (gv2);
+        }
+    }
+
+    v.keepOpen();
+}

mplot/CMakeLists.txt

@@ -79,6 +79,7 @@ install(
   LengthscaleVisual.h
   PointRowsMeshVisual.h
   PointRowsVisual.h
+  PolarVisual.h
   PolygonVisual.h
   QuadsMeshVisual.h
   QuadsVisual.h

mplot/CartGridVisual.h

@@ -68,6 +68,9 @@ namespace mplot {
         //! Do the computations to initialize the vertices that will represent the HexGrid.
         void initializeVertices()
         {
+            this->determine_datasize();
+            if (this->datasize == 0) { return; }
+
             // Optionally compute an offset to ensure that the cartgrid is centred about the mv_offset.
             if (this->centralize == true) {
                 float left_lim = -this->cg->width()/2.0f;
@@ -119,38 +122,12 @@ namespace mplot {
             this->idx = 0;
             unsigned int nrect = this->cg->num();
 
-            if (this->scalarData != nullptr) {
-                this->dcopy.resize (this->scalarData->size());
-                this->zScale.transform (*(this->scalarData), dcopy);
-                this->dcolour.resize (this->scalarData->size());
-                this->colourScale.transform (*(this->scalarData), dcolour);
-            } else if (this->vectorData != nullptr) {
-                this->dcopy.resize (this->vectorData->size());
-                this->dcolour.resize (this->vectorData->size());
-                this->dcolour2.resize (this->vectorData->size());
-                this->dcolour3.resize (this->vectorData->size());
-                std::vector<float> veclens(dcopy);
-                for (unsigned int i = 0; i < this->vectorData->size(); ++i) {
-                    veclens[i] = (*this->vectorData)[i].length();
-                    this->dcolour[i] = (*this->vectorData)[i][0];
-                    this->dcolour2[i] = (*this->vectorData)[i][1];
-                    // Could also extract a third colour for Trichrome vs Duochrome
-                    this->dcolour3[i] = (*this->vectorData)[i][2];
-                }
-                this->zScale.transform (veclens, this->dcopy);
-                if (this->cm.getType() != mplot::ColourMapType::RGB
-                    && this->cm.getType() != mplot::ColourMapType::RGBMono
-                    && this->cm.getType() != mplot::ColourMapType::RGBGrey) {
-                    this->colourScale.transform (this->dcolour, this->dcolour);
-                    this->colourScale2.transform (this->dcolour2, this->dcolour2);
-                    this->colourScale3.transform (this->dcolour3, this->dcolour3);
-                }
-            }
+            this->setupScaling();
 
             for (unsigned int ri = 0; ri < nrect; ++ri) {
                 std::array<float, 3> clr = this->setColour (ri);
                 this->vertex_push (this->cg->d_x[ri]+centering_offset[0],
-                                   this->cg->d_y[ri]+centering_offset[1], dcopy[ri], this->vertexPositions);
+                                   this->cg->d_y[ri]+centering_offset[1], this->dcopy[ri], this->vertexPositions);
                 this->vertex_push (clr, this->vertexColors);
                 this->vertex_push (0.0f, 0.0f, 1.0f, this->vertexNormals);
             }
@@ -189,40 +166,8 @@ namespace mplot {
             unsigned int nrect = this->cg->num();
             this->idx = 0;
 
-            if (this->scalarData != nullptr) {
-                this->dcopy.resize (this->scalarData->size());
-                this->zScale.transform (*(this->scalarData), dcopy);
-                this->dcolour.resize (this->scalarData->size());
-                this->colourScale.transform (*(this->scalarData), dcolour);
-            } else if (this->vectorData != nullptr) {
-                this->dcopy.resize (this->vectorData->size());
-                this->dcolour.resize (this->vectorData->size());
-                this->dcolour2.resize (this->vectorData->size());
-                this->dcolour3.resize (this->vectorData->size());
-                std::vector<float> veclens(dcopy);
-                for (unsigned int i = 0; i < this->vectorData->size(); ++i) {
-                    veclens[i] = (*this->vectorData)[i].length();
-                    this->dcolour[i] = (*this->vectorData)[i][0];
-                    this->dcolour2[i] = (*this->vectorData)[i][1];
-                    // Could also extract a third colour for Trichrome vs Duochrome (or for raw RGB signal)
-                    this->dcolour3[i] = (*this->vectorData)[i][2];
-                }
-                this->zScale.transform (veclens, this->dcopy);
-
-                // Handle case where this->cm.getType() == mplot::ColourMapType::RGB and there is
-                // exactly one colour. ColourMapType::RGB assumes R/G/B data all in range 0->1
-                // ALREADY and therefore they don't need to be re-scaled with this->colourScale.
-                if (this->cm.getType() != mplot::ColourMapType::RGB
-                    && this->cm.getType() != mplot::ColourMapType::RGBMono
-                    && this->cm.getType() != mplot::ColourMapType::RGBGrey) {
-                    this->colourScale.transform (this->dcolour, this->dcolour);
-                    // Dual axis colour maps like Duochrome and HSV will need to use colourScale2 to
-                    // transform their second colour/axis,
-                    this->colourScale2.transform (this->dcolour2, this->dcolour2);
-                    // Similarly for Triple axis maps
-                    this->colourScale3.transform (this->dcolour3, this->dcolour3);
-                } // else assume dcolour/dcolour2/dcolour3 are all in range 0->1 (or 0-255) already
-            }
+            this->setupScaling();
+
             float datumC = 0.0f;   // datum at the centre
             float datumNE = 0.0f;  // datum at the hex to the east.
             float datumNNE = 0.0f;
@@ -239,17 +184,17 @@ namespace mplot {
             for (unsigned int ri = 0; ri < nrect; ++ri) {
 
                 // Use the linear scaled copy of the data, dcopy.
-                datumC  = dcopy[ri];
-                datumNE =  R_HAS_NE(ri)  ? dcopy[R_NE(ri)] : datumC;
+                datumC  = this->dcopy[ri];
+                datumNE =  R_HAS_NE(ri)  ? this->dcopy[R_NE(ri)] : datumC;
                 //std::cout << "NE? " << (R_HAS_NE(ri) ? "yes\n" : "no\n");
-                datumNN =  R_HAS_NN(ri)  ? dcopy[R_NN(ri)] : datumC;
-                datumNW =  R_HAS_NW(ri)  ? dcopy[R_NW(ri)] : datumC;
+                datumNN =  R_HAS_NN(ri)  ? this->dcopy[R_NN(ri)] : datumC;
+                datumNW =  R_HAS_NW(ri)  ? this->dcopy[R_NW(ri)] : datumC;
                 //std::cout << "NW? " << (R_HAS_NW(ri) ? "yes\n" : "no\n");
-                datumNS =  R_HAS_NS(ri)  ? dcopy[R_NS(ri)] : datumC;
-                datumNNE = R_HAS_NNE(ri) ? dcopy[R_NNE(ri)] : datumC;
-                datumNNW = R_HAS_NNW(ri) ? dcopy[R_NNW(ri)] : datumC;
-                datumNSW = R_HAS_NSW(ri) ? dcopy[R_NSW(ri)] : datumC;
-                datumNSE = R_HAS_NSE(ri) ? dcopy[R_NSE(ri)] : datumC;
+                datumNS =  R_HAS_NS(ri)  ? this->dcopy[R_NS(ri)] : datumC;
+                datumNNE = R_HAS_NNE(ri) ? this->dcopy[R_NNE(ri)] : datumC;
+                datumNNW = R_HAS_NNW(ri) ? this->dcopy[R_NNW(ri)] : datumC;
+                datumNSW = R_HAS_NSW(ri) ? this->dcopy[R_NSW(ri)] : datumC;
+                datumNSE = R_HAS_NSE(ri) ? this->dcopy[R_NSE(ri)] : datumC;
 
                 // Use a single colour for each rect, even though rectangle's z
                 // positions are interpolated. Do the _colour_ scaling:
@@ -471,37 +416,9 @@ namespace mplot {
         float border_thickness_fixed = 0.0f;
 
     protected:
-        //! An overridable function to set the colour of rect ri
-        std::array<float, 3> setColour (unsigned int ri)
-        {
-            std::array<float, 3> clr = { 0.0f, 0.0f, 0.0f };
-            if (this->cm.numDatums() == 3) {
-                //if constexpr (std::is_same<std::decay_t<T>, unsigned char>::value == true) {
-                if constexpr (std::is_integral<std::decay_t<T>>::value) {
-                    // Differs from above as we divide by 255 to get value in range 0-1
-                    clr = this->cm.convert (this->dcolour[ri]/255.0f, this->dcolour2[ri]/255.0f, this->dcolour3[ri]/255.0f);
-                } else {
-                    clr = this->cm.convert (this->dcolour[ri], this->dcolour2[ri], this->dcolour3[ri]);
-                }
-            } else if (this->cm.numDatums() == 2) {
-                // Use vectorData
-                clr = this->cm.convert (this->dcolour[ri], this->dcolour2[ri]);
-            } else {
-                clr = this->cm.convert (this->dcolour[ri]);
-            }
-            return clr;
-        }
-
         //! The cartgrid to visualize
         const sm::cartgrid* cg;
 
-        //! A copy of the scalarData which can be transformed suitably to be the z value of the surface
-        std::vector<float> dcopy;
-        //! A copy of the scalarData (or first field of vectorData), scaled to be a colour value
-        std::vector<float> dcolour;
-        std::vector<float> dcolour2;
-        std::vector<float> dcolour3;
-
         // A centering offset to make sure that the Cartgrid is centred on
         // this->mv_offset. This is computed so that you *add* centering_offset to each
         // computed x/y/z position for a rectangle, and this means that the rectangle