Making a flat compound ray visualiser, plus some spherical projections

examples/CMakeLists.txt

@@ -537,3 +537,7 @@ endif()
 
 add_executable(show_boundingboxes show_boundingboxes.cpp)
 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)

examples/cray_eye.cpp

@@ -0,0 +1,80 @@
+/*
+ * A compound ray eye viewer
+ *
+ * Demonstrating use of mplot::compoundray::EyeVisual
+ */
+
+#include <iostream>
+#include <string>
+#include <memory>
+
+#include <sm/vec>
+#include <sm/vvec>
+
+#include <mplot/Visual.h>
+#include <mplot/ColourMap.h>
+#include <mplot/SphereVisual.h>
+#include <mplot/VectorVisual.h>
+#include <mplot/compoundray/EyeVisual.h>
+
+int main (int argc, char** argv)
+{
+    std::string eyefile = "";
+    if (argc < 2) {
+        std::cout << "Usage: " << argv[0] << " path/to/eyefile.eye [projection sphere radius]\n";
+        return -1;
+    } else {
+        eyefile = std::string (argv[1]);
+    }
+
+    float psrad = 0.6f;
+
+    [[maybe_unused]] float cx = 0.5f;
+    if (argc > 2) { cx = std::atof (argv[2]); }
+
+    auto v = mplot::Visual<>(1024, 768, "mplot::compoundray::EyeVisual");
+
+    // We read the information from the eye file into a vector of Ommatidium objects.  Ommatidium is
+    // defined in "cameras/CompoundEyeDataTypes.h" in compound ray, mplot::Ommatidium is a
+    // mplot/Seb's maths style equivalent. It contains 2 3D float vectors and two scalar floating point
+    // values.
+    auto ommatidia = std::make_unique<std::vector<mplot::compoundray::Ommatidium>>();
+    std::vector<std::array<float, 3>> ommatidiaColours;
+
+    // Read the eye file into ommatidia data structure. compound ray does this internally when we're
+    // using it, but for this example we instead make use of mplot::compoundray::readEye
+    if (mplot::compoundray::readEye (ommatidia.get(), eyefile) == nullptr) { std::cout << "Failed to read eye\n"; return -1; }
+
+    // Make some dummy data to demo the eye
+    sm::vvec<float> ommatidiaData;
+    ommatidiaData.linspace (0, 1, ommatidia->size());
+    // Colour it with a colour map
+    mplot::ColourMap cm (mplot::ColourMapType::Plasma);
+    ommatidiaColours.resize (ommatidia->size());
+    for (size_t i = 0; i < ommatidia->size(); ++i) {
+        ommatidiaColours[i] = cm.convert (ommatidiaData[i]);
+    }
+
+    auto eyevm = std::make_unique<mplot::compoundray::EyeVisual<>> (sm::vec<>{}, &ommatidiaColours, ommatidia.get());
+    v.bindmodel (eyevm);
+    eyevm->name = "Big Eye";
+    eyevm->show_cones = false;
+
+    auto ptype = mplot::compoundray::EyeVisual<>::projection_type::equirectangular; // mercator, equirectangular or cassini
+    sm::vec<> centre = { 0, 0, 0 };
+    sm::mat44<float> twod_tr;
+    twod_tr.translate (sm::vec<>{0,0,-2});
+    eyevm->add_spherical_projection (ptype, twod_tr, centre, psrad);
+
+    //centre = { cx, 0, 0 };
+    //eyevm->add_spherical_projection (ptype, twod_offset, centre, psrad, 780, 1560);
+
+    eyevm->show_sphere = true;
+    eyevm->show_rays = true;
+    eyevm->finalize();
+
+    [[maybe_unused]] auto ep = v.addVisualModel (eyevm);
+    ep->reinitColours();
+
+    v.keepOpen();
+}

examples/healpix.cpp

@@ -1,12 +1,16 @@
 /*
  * Make a healpix visual, showing the NEST index in a colour map
+ *
+ * This program also demonstrates 2D spherical projections that are available in sm::geometry and
+ * the Visual class mplot::SphereProjectionVisual
  */
 #include <cstdint>
 #include <cstdlib>
 #include <sstream>
 #include <sm/vec>
 #include <mplot/Visual.h>
 #include <mplot/HealpixVisual.h>
+#include <mplot/SphericalProjectionVisual.h>
 
 int main (int argc, char** argv)
 {
@@ -25,7 +29,17 @@ int main (int argc, char** argv)
     // scheme. If we fill it with sequential values, then the colour map will show the
     // hierarchical nature of the HEALPix.
     for (int64_t p = 0; p < hpv->n_pixels(); ++p) {
-        hpv->pixeldata[p] = static_cast<float>(p) / hpv->n_pixels();
+        hpv->pixeldata[p] = 0.5f * static_cast<float>(p) / hpv->n_pixels();
+        // Mark circles for the x, y and z axes to check they line up
+        hp::t_vec pv = hp::nest2vec (hpv->get_nside(), p);
+        auto v = sm::vec<double>({pv.x, pv.y, pv.z});
+        if (v.angle (sm::vec<double>::ux()) < 0.05) {
+            hpv->pixeldata[p] = 0.6;
+        } else if (v.angle (sm::vec<double>::uy()) < 0.05) {
+            hpv->pixeldata[p] = 0.8;
+        } else if (v.angle (sm::vec<double>::uz()) < 0.05) {
+            hpv->pixeldata[p] = 1.0;
+        }
     }
 
     std::stringstream ss;
@@ -37,7 +51,112 @@ int main (int argc, char** argv)
 
     // Finalize and add
     hpv->finalize();
-    v.addVisualModel (hpv);
+    auto hpvp = v.addVisualModel (hpv);
+
+    // Show some 2D projections, as well
+    sm::vvec<std::array<float, 3>> hpvcolours (hpvp->pixeldata.size(), mplot::colour::crimson);
+    sm::vvec<sm::vec<float, 2>> latlong (hpvp->pixeldata.size());
+    for (uint32_t i = 0; i < hpvp->pixeldata.size(); ++i) {
+
+        // ang.theta is co-latitude (and needs re-casting to be longitude), ang.phi is a longitude
+        hp::t_ang ang = hp::nest2ang (hpvp->get_nside(), i);
+
+        // ang.theta is pi for the South pole and 0 for the North pole
+        // latitude should be in range -pi/2 (S) < lat <= pi/2 (N)
+        float _lat = (sm::mathconst<float>::pi - ang.theta) - sm::mathconst<float>::pi_over_2;
+
+        latlong[i] = { _lat, static_cast<float>(ang.phi) };
+        hpvcolours[i] = hpvp->cm.convert (hpvp->pixeldata[i]);
+    }
+
+    // Add two-dimensional projections
+    auto spv = std::make_unique<mplot::SphericalProjectionVisual<float>> (sm::vec<float>{5,0,0});
+    v.bindmodel (spv);
+    spv->twodimensional (true);
+    spv->proj_type = sm::geometry::spherical_projection::type::mercator;
+    spv->latlong = latlong;
+    spv->lambda0 = 0.0f;
+    spv->colour = hpvcolours;
+    spv->finalize();
+    auto spvp = v.addVisualModel (spv);
+    auto ext = spvp->extents(); // Use VisualModel::extents() to help place the label
+    spvp->addLabel ("Mercator projection", sm::vec<>{ ext[0].min, ext[1].min - 0.16f, 0.0f }, mplot::TextFeatures(0.08f));
+
+    spv = std::make_unique<mplot::SphericalProjectionVisual<float>> (sm::vec<float>{13,0,0});
+    v.bindmodel (spv);
+    spv->twodimensional (true);
+    spv->proj_type = sm::geometry::spherical_projection::type::mercator;
+    spv->latlong = latlong;
+    spv->lambda0 = sm::mathconst<float>::pi_over_4;
+    spv->colour = hpvcolours;
+    spv->finalize();
+    spvp = v.addVisualModel (spv);
+    ext = spvp->extents(); // Use VisualModel::extents() to help place the label
+    std::stringstream ss1;
+    ss1 << "Mercator projection with "
+        << mplot::unicode::toUtf8(mplot::unicode::lambda)
+        << mplot::unicode::toUtf8(mplot::unicode::subs0)
+        << " = " << mplot::unicode::toUtf8(mplot::unicode::pi) << "/4";
+    spvp->addLabel (ss1.str(),
+                    sm::vec<>{ ext[0].min, ext[1].min - 0.16f, 0.0f }, mplot::TextFeatures(0.08f));
+
+    spv = std::make_unique<mplot::SphericalProjectionVisual<float>> (sm::vec<float>{-5,-4,0});
+    v.bindmodel (spv);
+    spv->twodimensional (true);
+    spv->proj_type = sm::geometry::spherical_projection::type::equirectangular;
+    spv->latlong = latlong;
+    spv->colour = hpvcolours;
+    spv->finalize();
+    spvp = v.addVisualModel (spv);
+    ext = spvp->extents();
+    spvp->addLabel ("Equirectangular projection", sm::vec<>{ ext[0].min, ext[1].min - 0.16f, 0.0f }, mplot::TextFeatures(0.08f));
+
+    spv = std::make_unique<mplot::SphericalProjectionVisual<float>> (sm::vec<float>{-5,-8,0});
+    v.bindmodel (spv);
+    spv->twodimensional (true);
+    spv->lambda0 = sm::mathconst<float>::pi_over_4;
+    //spv->phi0 = sm::mathconst<float>::pi_over_4; // Latitude offset
+    //spv->phi1 = sm::mathconst<float>::pi_over_4; // Longitude scaling
+    spv->proj_type = sm::geometry::spherical_projection::type::equirectangular;
+    spv->latlong = latlong;
+    spv->colour = hpvcolours;
+    spv->finalize();
+    spvp = v.addVisualModel (spv);
+    ext = spvp->extents();
+    std::stringstream ss2;
+    ss2 << "Equirectangular projection with "
+        << mplot::unicode::toUtf8(mplot::unicode::lambda)
+        << mplot::unicode::toUtf8(mplot::unicode::subs0)
+        << " = " << mplot::unicode::toUtf8(mplot::unicode::pi) << "/4";
+    spvp->addLabel (ss2.str(), sm::vec<>{ ext[0].min, ext[1].min - 0.16f, 0.0f }, mplot::TextFeatures(0.08f));
+
+    spv = std::make_unique<mplot::SphericalProjectionVisual<float>> (sm::vec<float>{-9,3,0});
+    v.bindmodel (spv);
+    spv->twodimensional (true);
+    spv->proj_type = sm::geometry::spherical_projection::type::cassini;
+    spv->latlong = latlong;
+    spv->colour = hpvcolours;
+    spv->finalize();
+    spvp = v.addVisualModel (spv);
+    ext = spvp->extents();
+    spvp->addLabel ("Cassini projection", sm::vec<>{ ext[0].min, ext[1].min - 0.16f, 0.0f }, mplot::TextFeatures(0.08f));
+
+    spv = std::make_unique<mplot::SphericalProjectionVisual<float>> (sm::vec<float>{-4,3,0});
+    v.bindmodel (spv);
+    spv->twodimensional (true);
+    spv->lambda0 = sm::mathconst<float>::pi_over_4;
+    spv->proj_type = sm::geometry::spherical_projection::type::cassini;
+    spv->latlong = latlong;
+    spv->colour = hpvcolours;
+    spv->finalize();
+    spvp = v.addVisualModel (spv);
+    ext = spvp->extents();
+    std::stringstream ss3;
+    ss3 << "Cassini projection with "
+        << mplot::unicode::toUtf8(mplot::unicode::lambda)
+        << mplot::unicode::toUtf8(mplot::unicode::subs0)
+        << " = " << mplot::unicode::toUtf8(mplot::unicode::pi) << "/4";
+    spvp->addLabel (ss3.str(), sm::vec<>{ ext[0].min, ext[1].min - 0.16f, 0.0f }, mplot::TextFeatures(0.08f));
 
     v.keepOpen();
     return 0;