diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile
index 284c1afe6..ec73932cc 100644
--- a/.devcontainer/Dockerfile
+++ b/.devcontainer/Dockerfile
@@ -29,7 +29,6 @@ RUN apt-get update && apt-get upgrade -y && \
 # Clone, build and clean up LASzip (3.x supports LAS 1.4 via laszip_api.h C API)
 RUN git clone https://github.com/LASzip/LASzip.git && \
     cd LASzip && \
-    git checkout tags/3.4.4 && \
     mkdir build && cd build && \
     cmake .. -DCMAKE_BUILD_TYPE=Release && \
     make -j$(nproc) && \
@@ -39,7 +38,6 @@ RUN git clone https://github.com/LASzip/LASzip.git && \
 # Clone, build and clean up Qhull
 RUN git clone https://github.com/qhull/qhull.git && \
     cd qhull && \
-    git checkout tags/v7.3.2 && \
     cd build && \
     cmake .. -DCMAKE_POSITION_INDEPENDENT_CODE:BOOL=true && \
     make -j$(nproc) && \
@@ -49,7 +47,6 @@ RUN git clone https://github.com/qhull/qhull.git && \
 # Clone, build and clean up libnabo
 RUN git clone https://github.com/ethz-asl/libnabo.git && \
     cd libnabo && \
-    git checkout tags/1.0.7 && \
     mkdir build && cd build && \
     cmake .. -DCMAKE_BUILD_TYPE=RelWithDebInfo && \
     make -j$(nproc) && \
@@ -67,4 +64,4 @@ RUN chmod +x /usr/local/bin/compile_install.sh
 WORKDIR /workspaces/raycloudtools
 
 # Set the default shell to bash
-SHELL ["/bin/bash", "-c"]
\ No newline at end of file
+SHELL ["/bin/bash", "-c"]
diff --git a/.devcontainer/commands/compile_install.sh b/.devcontainer/commands/compile_install.sh
old mode 100644
new mode 100755
index f2ae4d625..25a9dcc2e
--- a/.devcontainer/commands/compile_install.sh
+++ b/.devcontainer/commands/compile_install.sh
@@ -5,7 +5,7 @@ set -e
 echo "Building RayCloudTools..."
 cd /workspaces/raycloudtools
 mkdir -p build
-cd build 
+cd build
 
 # Construct the cmake command
 cmake .. \
@@ -41,4 +41,4 @@ make -j"$(nproc)"
 make install
 ldconfig /usr/local/lib
 
-echo "Build process completed successfully."
\ No newline at end of file
+echo "Build process completed successfully."
diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json
index 1768d359c..4c7553a0b 100644
--- a/.devcontainer/devcontainer.json
+++ b/.devcontainer/devcontainer.json
@@ -1,21 +1,30 @@
 {
-    "name": "RayCloudTools Development",
-    "build": {
-        "dockerfile": "Dockerfile",
-        "context": ".."
-    },
-    "customizations": {
-        "vscode": {
-            "extensions": [
-                "ms-vscode.cpptools",
-                "ms-vscode.cmake-tools",
-                "twxs.cmake",
-                "ms-vscode.cpptools-extension-pack"
-            ]
-        }
-    },
-    "remoteUser": "root",
-    "forwardPorts": [],
-    "postCreateCommand": "compile_install.sh",
-    "runArgs": ["--security-opt", "label=disable"]
-}
\ No newline at end of file
+  "name": "RayCloudTools Development",
+  "build": {
+    "dockerfile": "Dockerfile",
+    "context": "..",
+    "options": [
+      "--progress=plain"
+    ]
+  },
+  "customizations": {
+    "vscode": {
+      "extensions": [
+        "ms-vscode.cpptools",
+        "ms-vscode.cmake-tools",
+        "ms-vscode.cpptools-extension-pack"
+      ]
+    }
+  },
+  "mounts": [
+    "type=bind,source=${localWorkspaceFolder}/.devcontainer/commands/compile_install.sh,target=/usr/local/bin/compile_install.sh",
+  ],
+  "remoteUser": "root",
+  "forwardPorts": [],
+  "postCreateCommand": "compile_install.sh",
+  "runArgs": [
+    "--security-opt=label=disable"
+  ],
+  "workspaceFolder": "/workspaces/raycloudtools",
+  "workspaceMount": "type=bind,source=${localWorkspaceFolder},target=/workspaces/raycloudtools"
+}
diff --git a/.gitignore b/.gitignore
index 2ea21aab5..399bfaf18 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,4 @@
-.vscode/
 .idea/
 build*/
 cmake-build*/
+treetools/
diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 000000000..0637bc338
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1,21 @@
+{
+    "editor.formatOnSave": true,
+    "editor.rulers": [
+        120
+    ],
+    "files.trimFinalNewlines": true,
+    "files.insertFinalNewline": true,
+    "files.trimTrailingWhitespace": true,
+    "C_Cpp.default.cppStandard": "c++20",
+    "C_Cpp.default.cStandard": "c99",
+    "C_Cpp.default.compileCommands": "${workspaceFolder}/build/compile_commands.json",
+    "C_Cpp.files.exclude": {
+        "**/.vscode": true,
+        "build/**": true,
+    },
+    "C_Cpp.exclusionPolicy": "checkFolders",
+    "[cpp]": {
+        "editor.defaultFormatter": "ms-vscode.cpptools",
+        "editor.tabSize": 2
+    }
+}
diff --git a/docker/Dockerfile b/docker/Dockerfile
index 0df8349c2..e06cf5a41 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -34,7 +34,6 @@ RUN mkdir -p /deps && \
 # Clone, build and clean up LASzip (3.x supports LAS 1.4 via laszip_api.h C API)
 RUN git clone https://github.com/LASzip/LASzip.git && \
     cd LASzip && \
-    git checkout tags/3.4.4 && \
     mkdir build && cd build && \
     cmake .. -DCMAKE_BUILD_TYPE=Release && \
     make -j$(nproc) && \
@@ -97,7 +96,7 @@ RUN git clone https://github.com/csiro-robotics/treetools.git && \
     make install && \
     cd ../.. && \
     rm -rf treetools
-    
+
 # Update ldconfig
 RUN ldconfig /usr/local/lib
 
diff --git a/raycloudtools/raydiff/raydiff.cpp b/raycloudtools/raydiff/raydiff.cpp
index ee60d38c6..e4f29d6b9 100644
--- a/raycloudtools/raydiff/raydiff.cpp
+++ b/raycloudtools/raydiff/raydiff.cpp
@@ -17,7 +17,7 @@
 void usage(int exit_code = 1)
 {
   // clang-format off
-  std::cout << "Difference between two ray clouds output a coloured cloud, cloud1 differences in red, cloud2 differences in green, and similarity printed to screen." << std::endl;
+  std::cout << "Difference between two ray clouds output a coloured cloud, cloud1 differences in scarlet, cloud2 differences in green/cyan, and similarity printed to screen." << std::endl;
   std::cout << "usage (--view / -v to view results):" << std::endl;
   std::cout << "raydiff cloud1.ply cloud2.ply" << std::endl;
   std::cout << "                              --distance 0 - optional threshold in m for colouring differences. Default auto-detects distribution shoulder" << std::endl;
diff --git a/raycloudtools/rayimport/rayimport.cpp b/raycloudtools/rayimport/rayimport.cpp
index a93a6dae9..c04cc72d6 100644
--- a/raycloudtools/rayimport/rayimport.cpp
+++ b/raycloudtools/rayimport/rayimport.cpp
@@ -240,4 +240,4 @@ int main(int argc, char *argv[])
   if (save_file != "")
     ray::viewFile(save_file); 
   return res;
-}
\ No newline at end of file
+}
diff --git a/raycloudtools/rayinfo/rayinfo.cpp b/raycloudtools/rayinfo/rayinfo.cpp
index 1838e1425..806949cc2 100644
--- a/raycloudtools/rayinfo/rayinfo.cpp
+++ b/raycloudtools/rayinfo/rayinfo.cpp
@@ -74,6 +74,8 @@ int rayInfo(int argc, char *argv[])
   double min_ray_length = std::numeric_limits<double>::max();
   double max_ray_length = std::numeric_limits<double>::lowest();
   double max_bounded_ray_length = std::numeric_limits<double>::lowest();
+  double mean_bounded_ray_length = 0.0;
+  int mean_count = 0;
   ray::RGBA min_col(255, 255, 255, 255), max_col(0, 0, 0, 0);
   int num_pixels_covered = 0;
   const double voxel_width = 0.5;
@@ -125,6 +127,8 @@ int rayInfo(int argc, char *argv[])
       if (colours[i].alpha > 0)
       {
         max_bounded_ray_length = std::max(max_bounded_ray_length, ray_length);
+        mean_bounded_ray_length += ray_length;
+        mean_count++;
 
         min_col.red = std::min(min_col.red, colours[i].red);
         min_col.green = std::min(min_col.green, colours[i].green);
@@ -145,6 +149,7 @@ int rayInfo(int argc, char *argv[])
   {
     usage();
   }
+  mean_bounded_ray_length /= (double)mean_count;
 
   // print the results to screen
   std::cout << std::endl;
@@ -196,7 +201,7 @@ int rayInfo(int argc, char *argv[])
       std::cout << num_minutes << " mins ";
     std::cout << seconds << " s \t(in seconds: " << path_period << ")" << std::endl;
   }
-  std::cout << "  ray length: \t\t" << min_ray_length << " to " << max_ray_length << " m, max end point ray length: " << max_bounded_ray_length << " m" << std::endl;
+  std::cout << "  ray length: \t\t" << min_ray_length << " to " << max_ray_length << " m, max end point ray length: " << max_bounded_ray_length << " m, mean: " << mean_bounded_ray_length << std::endl;
   std::cout << "  colour range (RGBA): \t" << (int)min_col.red << "," << (int)min_col.green << "," << (int)min_col.blue << "," << (int)min_col.alpha << " to " << 
                                             (int)max_col.red << "," << (int)max_col.green << "," << (int)max_col.blue << "," << (int)max_col.alpha << std::endl;
  
diff --git a/raylib/extraction/rayforest_watershed.cpp b/raylib/extraction/rayforest_watershed.cpp
index e31ba7360..fc346cb90 100644
--- a/raylib/extraction/rayforest_watershed.cpp
+++ b/raylib/extraction/rayforest_watershed.cpp
@@ -11,6 +11,7 @@
 #include <iostream>
 #include <queue>
 #include "../rayply.h"
+#include "../rayrandom.h"
 
 namespace ray
 {
@@ -39,9 +40,9 @@ void Forest::renderWatershed(const std::string &cloud_name_stub, std::vector<Tre
         continue;
       }
       srand(1 + ind);
-      colour.red = (uint8_t)(rand() % 256);
-      colour.green = (uint8_t)(rand() % 256);
-      colour.blue = (uint8_t)(rand() % 256);
+      colour.red = (uint8_t)(::ray::rand() % 256);
+      colour.green = (uint8_t)(::ray::rand() % 256);
+      colour.blue = (uint8_t)(::ray::rand() % 256);
       Eigen::Vector3d pos =
         min_bounds_ + voxel_width_ * Eigen::Vector3d(0.5 + static_cast<double>(x), 0.5 + static_cast<double>(y), 0);
       pos[2] = original_heightfield_(x, y);
@@ -54,10 +55,10 @@ void Forest::renderWatershed(const std::string &cloud_name_stub, std::vector<Tre
   for (auto &ind : indices)
   {
     Eigen::Vector3d tip;
-    srand(1 + ind);
-    colour.red = (uint8_t)(rand() % 256);
-    colour.green = (uint8_t)(rand() % 256);
-    colour.blue = (uint8_t)(rand() % 256);
+    ::ray::srand(1 + ind);
+    colour.red = (uint8_t)(::ray::rand() % 256);
+    colour.green = (uint8_t)(::ray::rand() % 256);
+    colour.blue = (uint8_t)(::ray::rand() % 256);
     if (trees[ind].area < min_area_)
     {
       continue;
diff --git a/raylib/extraction/raygrid2d.cpp b/raylib/extraction/raygrid2d.cpp
index 36efceb9d..fb5536064 100644
--- a/raylib/extraction/raygrid2d.cpp
+++ b/raylib/extraction/raygrid2d.cpp
@@ -203,7 +203,7 @@ void RayIndexGrid2D::init(const Eigen::Vector3d &min_bound, const Eigen::Vector3
             << std::endl;
 
   pixels_.resize(dims_[0] * dims_[1]);
-  memset(&pixels_[0], 0, sizeof(Pixel) * pixels_.size());
+  std::fill(pixels_.begin(), pixels_.end(), Pixel{});
 }
 
 void RayIndexGrid2D::fillRays(const Cloud &cloud)
diff --git a/raylib/extraction/rayleaves.cpp b/raylib/extraction/rayleaves.cpp
index 896bfdd59..733075d1c 100644
--- a/raylib/extraction/rayleaves.cpp
+++ b/raylib/extraction/rayleaves.cpp
@@ -10,6 +10,7 @@
 #include "../rayply.h"
 #include "../raymesh.h"
 #include "../rayforeststructure.h"
+#include "../rayrandom.h"
 #define STB_IMAGE_IMPLEMENTATION
 #include "raylib/imageread.h"
 
@@ -144,10 +145,10 @@ bool generateLeaves(const std::string &cloud_stub, const std::string &trees_file
   };
   std::vector<Leaf> leaves;
   std::vector<double> leaf_counter(grid.voxels().size());
-  std::srand(1);
+  ::ray::srand(1);
   for (size_t i = 0; i<grid.voxels().size(); i++)
   {
-    leaf_counter[i] = (double)(std::rand()%10000) / 10000.0; // a random start stops regions of low density have 0 leaves
+    leaf_counter[i] = (double)(::ray::rand()%10000) / 10000.0; // a random start stops regions of low density have 0 leaves
   }
 
 
@@ -170,7 +171,7 @@ bool generateLeaves(const std::string &cloud_stub, const std::string &trees_file
       double num_leaves_d = desired_leaf_area / leaf_area;
       double num_points = (double)voxel.numHits();
       double &count = leaf_counter[index];
-      count += num_leaves_d / num_points;
+      count += num_leaves_d / num_points; 
       bool add_leaf = false;
       if (count >= 1.0)
       {
diff --git a/raylib/extraction/rayterrain.cpp b/raylib/extraction/rayterrain.cpp
index 718bf8004..e4425a26d 100644
--- a/raylib/extraction/rayterrain.cpp
+++ b/raylib/extraction/rayterrain.cpp
@@ -9,6 +9,7 @@
 #include "../rayply.h"
 #include "../rayprogress.h"
 #include "../rayprogressthread.h"
+#include "../rayrandom.h"
 
 #if RAYLIB_WITH_TBB
 #include <tbb/enumerable_thread_specific.h>
@@ -49,8 +50,8 @@ struct Node
   {
     num_visits++;
 // This checks in a cone rather than just the corner of a cube shape that you would get
-// from a raw Pareto front calculation    
-#define CONE_CHECK  
+// from a raw Pareto front calculation
+#define CONE_CHECK
     Vector4d dif = corner - pos;
     if (dif == Vector4d(0, 0, 0, 0))
     {
@@ -123,7 +124,7 @@ void constructOctalSpacePartition(std::vector<Node> &nodes, std::vector<Vector4d
   int i = 0;
   while (points.size() > 0)
   {
-    const int ind = rand() % static_cast<int>(points.size());
+    const int ind = ::ray::rand() % static_cast<int>(points.size());
     nodes[i++].pos = points[ind];
     points[ind] = points.back();
     points.pop_back();
@@ -177,7 +178,7 @@ void Terrain::getParetoFront(const std::vector<Vector4d> &points, std::vector<Ve
 #if RAYLIB_WITH_TBB
   tbb::parallel_for<size_t>(0, nodes.size(), process_rays);
 #else
-  #pragma omp parallel for
+#pragma omp parallel for
   for (size_t n = 0; n < nodes.size(); n++)
   {
     process_rays(n);
@@ -200,7 +201,7 @@ void Terrain::getParetoFront(const std::vector<Vector4d> &points, std::vector<Ve
 void Terrain::growUpwards(const std::vector<Eigen::Vector3d> &positions, double gradient)
 {
 #if RAYLIB_WITH_QHULL
-  // The idea behind ground extraction is to tilt the upwards vector to the (1,1,1) direction then 
+  // The idea behind ground extraction is to tilt the upwards vector to the (1,1,1) direction then
   // find the Pareto front in the three principle axes. https://en.wikipedia.org/wiki/Pareto_front
   //
   // Efficient Pareto front calculation is based on: Algorithms and Analyses for Maximal Vector Computation. Godfrey
@@ -255,6 +256,10 @@ void Terrain::growUpwards(const std::vector<Eigen::Vector3d> &positions, double
 
   mesh_.indexList() = hull.mesh().indexList();
   mesh_.vertices() = vecs;
+#else
+  RAYLIB_UNUSED(positions);
+  RAYLIB_UNUSED(gradient);
+  std::cerr << "growUpwards: extracting terrain requires QHull, see README instructions for installation" << std::endl;
 #endif
 }
 
@@ -349,11 +354,20 @@ void Terrain::growUpwardsFast(const std::vector<Eigen::Vector3d> &ends, double p
   std::cout << "size before: " << ends.size() << ", size after: " << points.size() << std::endl;
 
   growUpwards(points, gradient);
+#else
+  RAYLIB_UNUSED(ends);
+  RAYLIB_UNUSED(pixel_width);
+  RAYLIB_UNUSED(min_bound);
+  RAYLIB_UNUSED(max_bound);
+  RAYLIB_UNUSED(gradient);
+  std::cerr << "growUpwardsFast: extracting terrain requires QHull, see README instructions for installation"
+            << std::endl;
 #endif
 }
 
-// Convert the @c cloud input to the mesh_ member variable. 
-void Terrain::extract(const Cloud &cloud, const Eigen::Vector3d &offset, const std::string &file_prefix, double gradient, bool verbose)
+// Convert the @c cloud input to the mesh_ member variable.
+void Terrain::extract(const Cloud &cloud, const Eigen::Vector3d &offset, const std::string &file_prefix,
+                      double gradient, bool verbose)
 {
 #if RAYLIB_WITH_QHULL
   // preprocessing to make the cloud smaller.
@@ -386,6 +400,11 @@ void Terrain::extract(const Cloud &cloud, const Eigen::Vector3d &offset, const s
     local_cloud.save(file_prefix + "_terrain.ply");
   }
 #else
+  RAYLIB_UNUSED(cloud);
+  RAYLIB_UNUSED(offset);
+  RAYLIB_UNUSED(file_prefix);
+  RAYLIB_UNUSED(gradient);
+  RAYLIB_UNUSED(verbose);
   std::cerr << "Error: extracting terrain requires QHull, see README instructions for installation" << std::endl;
 #endif
 }
diff --git a/raylib/extraction/raytrunks.cpp b/raylib/extraction/raytrunks.cpp
index 0b83acf94..1040e87be 100644
--- a/raylib/extraction/raytrunks.cpp
+++ b/raylib/extraction/raytrunks.cpp
@@ -9,6 +9,7 @@
 #include "../raycuboid.h"
 #include "../raygrid.h"
 #include "../rayply.h"
+#include "../rayrandom.h"
 #include "raygrid2d.h"
 
 namespace ray
@@ -652,9 +653,9 @@ void Trunks::saveDebugTrunks(const std::string &filename, bool verbose, const st
     for (auto &id : lowest_trunk_ids)
     {
       const Trunk &trunk = trunks[id];
-      colour.red = uint8_t(rand() % 255);
-      colour.green = uint8_t(rand() % 255);
-      colour.blue = uint8_t(rand() % 255);
+      colour.red = uint8_t(::ray::rand() % 255);
+      colour.green = uint8_t(::ray::rand() % 255);
+      colour.blue = uint8_t(::ray::rand() % 255);
 
       const Eigen::Vector3d side1 = trunk.dir.cross(Eigen::Vector3d(1, 2, 3)).normalized();
       const Eigen::Vector3d side2 = side1.cross(trunk.dir);
@@ -725,6 +726,7 @@ std::vector<std::pair<Eigen::Vector3d, double>> Trunks::load(const std::string &
         if (i < 3)
         {
           base[i] = std::stod(token.c_str());
+          radius = 0.0;
         }
         else
         {
diff --git a/raylib/imageread.h b/raylib/imageread.h
index bafd7117a..8a5e7fdb8 100644
--- a/raylib/imageread.h
+++ b/raylib/imageread.h
@@ -3038,9 +3038,9 @@ static int stbi__parse_entropy_coded_data(stbi__jpeg *z)
 
 static void stbi__jpeg_dequantize(short *data, stbi__uint16 *dequant)
 {
-   int i;
+   std::size_t i;
    for (i=0; i < 64; ++i)
-      data[i] *= dequant[i];
+      data[i] *= static_cast<short>(dequant[i]);
 }
 
 static void stbi__jpeg_finish(stbi__jpeg *z)
@@ -4996,10 +4996,10 @@ static void stbi__de_iphone(stbi__png *z)
             stbi_uc a = p[3];
             stbi_uc t = p[0];
             if (a) {
-               stbi_uc half = a / 2;
-               p[0] = (p[2] * 255 + half) / a;
-               p[1] = (p[1] * 255 + half) / a;
-               p[2] = ( t   * 255 + half) / a;
+               stbi_uc half = a >> 1; // divided by 2.
+               p[0] = static_cast<stbi_uc>((p[2] * 255 + half) / a);
+               p[1] = static_cast<stbi_uc>((p[1] * 255 + half) / a);
+               p[2] = static_cast<stbi_uc>(( t   * 255 + half) / a);
             } else {
                p[0] = p[2];
                p[2] = t;
@@ -7055,9 +7055,10 @@ static void stbi__hdr_convert(float *output, stbi_uc *input, int req_comp)
    if ( input[3] != 0 ) {
       float f1;
       // Exponent
-      f1 = (float) ldexp(1.0f, input[3] - (int)(128 + 8));
+      f1 = ldexp(1.0f, input[3] - (int)(128 + 8));
       if (req_comp <= 2)
-         output[0] = (input[0] + input[1] + input[2]) * f1 / 3;
+         // Promote stbi_uc int to prevent overflow. Result always fits in range of float.
+         output[0] = static_cast<float>(static_cast<int>(input[0]) + input[1] + input[2]) * f1 / 3.0f;
       else {
          output[0] = input[0] * f1;
          output[1] = input[1] * f1;
diff --git a/raylib/rayalignment.cpp b/raylib/rayalignment.cpp
index c94421011..fae1a506e 100644
--- a/raylib/rayalignment.cpp
+++ b/raylib/rayalignment.cpp
@@ -56,7 +56,7 @@ void Array3D::init(const Eigen::Vector3d &box_min, double voxel_width, const Eig
   voxel_width_ = voxel_width;
   dims_ = dimensions;
   cells_.resize(dims_[0] * dims_[1] * dims_[2]);
-  memset(&cells_[0], 0, cells_.size() * sizeof(Complex));
+  std::fill(cells_.begin(), cells_.end(), Complex{});
   null_cell_ = 0;
 }
 
@@ -150,7 +150,7 @@ void Array3D::fillWithRays(const Cloud &cloud)
 void Array1D::init(int length)
 {
   cells_.resize(length);
-  memset(&cells_[0], 0, cells_.size() * sizeof(Complex));
+  std::fill(cells_.begin(), cells_.end(), Complex{});
 }
 
 void Array1D::operator*=(const Array1D &other)
diff --git a/raylib/raybuildinggen.cpp b/raylib/raybuildinggen.cpp
index 86b0e9d0a..3c7c4ffd9 100644
--- a/raylib/raybuildinggen.cpp
+++ b/raylib/raybuildinggen.cpp
@@ -67,7 +67,7 @@ void BuildingGen::splitRoom(const Cuboid &cuboid, std::vector<Cuboid> &cuboids,
       double table_width = 1.0;
       double table_length = 1.5;
       const double table_height = 1.1;
-      if (ray::rand() % 2)
+      if (::ray::rand() % 2)
         std::swap(table_width, table_length);
       Eigen::Vector3d table_rad = 0.5 * Eigen::Vector3d(table_width, table_length, table_height);
       Eigen::Vector3d table_pos(
@@ -87,8 +87,8 @@ void BuildingGen::splitRoom(const Cuboid &cuboid, std::vector<Cuboid> &cuboids,
       double cupboard_width = 0.6;
       double cupboard_length = 2.0;
       double cupboard_height = random(1.0, 3.0);
-      int wall_id = ray::rand() % 2;
-      int side = ray::rand() % 2;
+      int wall_id = ::ray::rand() % 2;
+      int side = ::ray::rand() % 2;
       if (wall_id == 1)
         std::swap(cupboard_width, cupboard_length);
       Eigen::Vector3d cupboard_rad = 0.5 * Eigen::Vector3d(cupboard_width, cupboard_length, cupboard_height);
@@ -201,7 +201,7 @@ void BuildingGen::generate()
   left.max_bound_[0] = building.min_bound_[0] - params_.outer_wall_width;
   left.min_bound_[0] = -big;
   cuboids.push_back(left);
-  int num_windows = ray::rand() % 11;
+  int num_windows = ::ray::rand() % 11;
   for (int i = 0; i < num_windows; i++)
   {
     Eigen::Vector3d pos =
diff --git a/raylib/raydiffer.cpp b/raylib/raydiffer.cpp
index c2b17ee14..5d6e149ef 100644
--- a/raylib/raydiffer.cpp
+++ b/raylib/raydiffer.cpp
@@ -13,7 +13,7 @@
 //  if the points are in a plane  then k=2
 //  if the points are in a volume then k=3
 // this macro finds the k that best fit the uniform distribution. Otherwise we default to k=2
-#define FIND_CORRELATION_DIMENSION 
+#define FIND_CORRELATION_DIMENSION
 
 namespace ray
 {
@@ -30,7 +30,7 @@ void calcNearestNeighbourDistances(const std::vector<Eigen::Vector3f> &cloud1, c
       num_bounded2++;
   Eigen::MatrixXf points_p(3, num_bounded);
   int j = 0;
-  for (unsigned int i = 0; i < cloud1.size(); i++) 
+  for (unsigned int i = 0; i < cloud1.size(); i++)
   {
     if (cloud1[i][0] != std::numeric_limits<float>::max())
       points_p.col(j++) = cloud1[i];
@@ -39,7 +39,7 @@ void calcNearestNeighbourDistances(const std::vector<Eigen::Vector3f> &cloud1, c
 
   Eigen::MatrixXf points_q(3, num_bounded2);
   j = 0;
-  for (unsigned int i = 0; i < cloud2.size(); i++) 
+  for (unsigned int i = 0; i < cloud2.size(); i++)
   {
     if (cloud2[i][0] != std::numeric_limits<float>::max())
       points_q.col(j++) = cloud2[i];
@@ -74,7 +74,7 @@ double getShoulder(double k, std::vector<float> sorted_dists, double &min_error_
   // 1. transform the result to make the distances uniform if their 3D points are uniformly distributed
   for (int i = 0; i<num; i++)
   {
-    sorted_dists[i] = (float)std::pow((double)sorted_dists[i], k); 
+    sorted_dists[i] = (float)std::pow((double)sorted_dists[i], k);
   }
 
   // 2. accumulate outside term backwards
@@ -88,7 +88,7 @@ double getShoulder(double k, std::vector<float> sorted_dists, double &min_error_
     double I = (float)i - yN;
     outside_const[i] = (float)((double)outside_const[i+1] + I*I);
     outside_linear[i] = (float)((double)outside_linear[i+1] + 2.0*I*d);
-    outside_square[i] =(float)((double)outside_square[i+1] + d*d);      
+    outside_square[i] =(float)((double)outside_square[i+1] + d*d);
   }
 
   // 3. accumulate linear term forwards, but store only best results
@@ -96,7 +96,7 @@ double getShoulder(double k, std::vector<float> sorted_dists, double &min_error_
   double inside_linear = 0.0;
   double inside_square = 0.0;
   double min_error_sqr = 0.0;
-  double min_error_i = 0.0; 
+  double min_error_i = 0.0;
   min_error_dist = 0.0;
   for (int i = 0; i<num; i++)
   {
@@ -111,7 +111,7 @@ double getShoulder(double k, std::vector<float> sorted_dists, double &min_error_
     double linear = outside_linear[i]/d;
     double a = square;
     double b = -linear - 2.0*yN*square;
-    double c = outside_const[i] + linear*yN + square*yN*yN; 
+    double c = outside_const[i] + linear*yN + square*yN*yN;
 
     // add the inside part
     a += inside_square/sqr(sorted_dists[i]); // the division is to match the gradient to y
@@ -227,10 +227,10 @@ double printDistanceStatistics(const std::vector<float> &dists_to_cloud1, const
   return similarity;
 }
 
-bool writeDifferencesToRayClouds(const std::string &cloud1_namestub, const std::string &cloud2_namestub, std::vector<float> &dists_to_cloud1, 
+bool writeDifferencesToRayClouds(const std::string &cloud1_namestub, const std::string &cloud2_namestub, std::vector<float> &dists_to_cloud1,
    std::vector<float> &dists_to_cloud2, double dist_threshold, bool individual_files, bool visualise)
 {
-  std::cout << "saving out differences, coloured red for differences to " << cloud1_namestub << ".ply and green for differences in " << cloud2_namestub << ".ply" << std::endl;
+  std::cout << "saving out differences, coloured scarlet for differences to " << cloud1_namestub << ".ply and green/cyan for differences in " << cloud2_namestub << ".ply" << std::endl;
 
   // now render visuals
   CloudWriter writer;
@@ -243,11 +243,11 @@ bool writeDifferencesToRayClouds(const std::string &cloud1_namestub, const std::
   std::vector<Eigen::Vector3i> samples;
 
   int j = 0;
-  Eigen::Vector3d diff_col(255,0,0);
+  Eigen::Vector3d diff_col(255,0,127);
   std::vector<float> *dists = &dists_to_cloud1;
   const float eps = 1e-8f; // in case there is inaccuracy in the KNN distance estimation for co-located point pairs
   auto colour = [&](std::vector<Eigen::Vector3d> &starts, std::vector<Eigen::Vector3d> &ends,
-                    std::vector<double> &times, std::vector<RGBA> &colours) 
+                    std::vector<double> &times, std::vector<RGBA> &colours)
   {
     // firstly we store a count per cell
     chunk.clear();
@@ -277,13 +277,13 @@ bool writeDifferencesToRayClouds(const std::string &cloud1_namestub, const std::
   j = 0;
   dists_to_cloud1.clear();
   dists_to_cloud1.shrink_to_fit();
-  Eigen::Vector3d diff2_col(0,255,0);
+  Eigen::Vector3d diff2_col(0,255,127);
 
   dists = &dists_to_cloud2;
   if (!individual_files)
   {
     auto colour2 = [&](std::vector<Eigen::Vector3d> &starts, std::vector<Eigen::Vector3d> &ends,
-                      std::vector<double> &times, std::vector<RGBA> &colours) 
+                      std::vector<double> &times, std::vector<RGBA> &colours)
     {
       // firstly we store a count per cell
       chunk.clear();
@@ -313,7 +313,9 @@ bool writeDifferencesToRayClouds(const std::string &cloud1_namestub, const std::
     writer.end();
 
     if (visualise)
+    {
       viewFile(cloud1_namestub + "_diff.ply");
+    }
   }
   else
   {
@@ -326,7 +328,9 @@ bool writeDifferencesToRayClouds(const std::string &cloud1_namestub, const std::
       return false;
     writer.end();
     if (visualise)
+    {
       viewFile(cloud1_namestub + "_diff.ply", cloud2_namestub + "_diff.ply");
+    }
   }
   return true;
 }
diff --git a/raylib/raylaz.cpp b/raylib/raylaz.cpp
index 439c9ce75..e429c0017 100644
--- a/raylib/raylaz.cpp
+++ b/raylib/raylaz.cpp
@@ -153,6 +153,7 @@ bool readLas(const std::string &file_name,
   std::cout << "loaded " << file_name << " with " << number_of_points << " points" << std::endl;
   return true;
 #else   // RAYLIB_WITH_LAS
+  RAYLIB_UNUSED(offset_to_remove);
   RAYLIB_UNUSED(max_intensity);
   RAYLIB_UNUSED(file_name);
   RAYLIB_UNUSED(apply);
@@ -169,8 +170,7 @@ bool readLas(std::string file_name, std::vector<Eigen::Vector3d> &positions, std
 {
   std::vector<Eigen::Vector3d> starts;  // dummy as lax just reads in point clouds, not ray clouds
   auto apply = [&](std::vector<Eigen::Vector3d> &start_points, std::vector<Eigen::Vector3d> &end_points,
-                   std::vector<double> &time_points, std::vector<RGBA> &colour_values)
-  {
+                   std::vector<double> &time_points, std::vector<RGBA> &colour_values) {
     starts.insert(starts.end(), start_points.begin(), start_points.end());
     positions.insert(positions.end(), end_points.begin(), end_points.end());
     times.insert(times.end(), time_points.begin(), time_points.end());
diff --git a/raylib/rayrandom.h b/raylib/rayrandom.h
index 7c676c6f8..5e1fbf7a7 100644
--- a/raylib/rayrandom.h
+++ b/raylib/rayrandom.h
@@ -18,12 +18,13 @@ namespace ray
 /// This is about ~5x faster then using the standard std::rand() method with no additional drawbacks.
 class RAYLIB_EXPORT PCGRandomGenerator
 {
-public:
+private:
   PCGRandomGenerator()
   {
     seed(12748, 3147, 792751, 14992);
   }  // we can't initialise using std::rand, as it isn't platform independent
 
+public:
   /// The minimum possible result value (inclusive). For compatibility with @c std::shuffle()
   static constexpr unsigned min() { return std::numeric_limits<unsigned int>::min(); }
   /// The maximum possible result value (inclusive). For compatibility with @c std::shuffle()
@@ -62,14 +63,14 @@ class RAYLIB_EXPORT PCGRandomGenerator
 };
 
 /// Return a random number using some magic engine behind the hood.
-/// Use this over std::rand() as there's a good chance this will be faster/better.
+/// Use this over std::rand() as it will be consistent across platforms.
 unsigned int RAYLIB_EXPORT rand();
 void RAYLIB_EXPORT srand(unsigned int seed);
 
 /// Return a uniformed number between [0,1).
 inline double randUniformDouble()
 {
-  return static_cast<double>(ray::rand() % std::numeric_limits<int>::max()) /
+  return static_cast<double>(::ray::rand() % std::numeric_limits<int>::max()) /
          static_cast<double>(std::numeric_limits<int>::max());
 }
 }  // namespace ray
diff --git a/raylib/rayutils.h b/raylib/rayutils.h
index 470c5e27a..c2e016442 100644
--- a/raylib/rayutils.h
+++ b/raylib/rayutils.h
@@ -9,7 +9,10 @@
 #include "raylib/raylibconfig.h"
 #include "rayrandom.h"
 
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
 #include <Eigen/Dense>
+#pragma GCC diagnostic pop
 #include <algorithm>
 #include <cassert>
 #include <chrono>
@@ -45,13 +48,13 @@ inline int runWithMemoryCheck(std::function<bool(int argc, char *argv[])> main_f
     std::cerr << "Error: Not enough memory to process the input file," << std::endl;
     std::cerr << "consider using raydecimate or raysplit grid to operate on a smaller file." << std::endl;
     return 1;
-  }  
+  }
   catch (std::length_error const &)  // catch any memory allocation problems in generating large images
   {
     std::cerr << "Error: Not enough memory to process the input file," << std::endl;
     std::cerr << "consider using raydecimate or raysplit grid to operate on a smaller file." << std::endl;
     return 1;
-  }    
+  }
 }
 
 inline std::vector<std::string> split(const std::string &s, char delim)
@@ -165,7 +168,7 @@ inline T median(std::vector<T> list)
   {
     typename std::vector<T>::iterator middle2 = middle - 1;
     nth_element(first, middle2, last);
-    return (*middle + *middle2) / static_cast<T>(2);
+    return static_cast<T>(*middle + *middle2) / static_cast<T>(2);
   }
 }
 
@@ -198,16 +201,21 @@ inline std::vector<T> componentList(const std::vector<U> &list, const T &compone
 
 struct RGBA
 {
-  RGBA(){}
-  RGBA(uint8_t r, uint8_t g, uint8_t b, uint8_t a) : red(r), green(g), blue(b), alpha(a) {}
+  RGBA() {}
+  RGBA(uint8_t r, uint8_t g, uint8_t b, uint8_t a)
+    : red(r)
+    , green(g)
+    , blue(b)
+    , alpha(a)
+  {}
   uint8_t red;
   uint8_t green;
   uint8_t blue;
   uint8_t alpha;
-  static RGBA white(){ return RGBA(255, 255, 255, 255); }
-  static RGBA terrain(){ return RGBA(149,105,72, 255); }
-  static RGBA treetrunk(){ return RGBA(192,166,141, 255); }
-  static RGBA leaves(){ return RGBA(60,102,44, 255); }
+  static RGBA white() { return RGBA(255, 255, 255, 255); }
+  static RGBA terrain() { return RGBA(149, 105, 72, 255); }
+  static RGBA treetrunk() { return RGBA(192, 166, 141, 255); }
+  static RGBA leaves() { return RGBA(60, 102, 44, 255); }
 };
 
 /// Converts a value from 0 to 1 into a RGBA structure
@@ -370,33 +378,34 @@ inline int sign(double x)
 }
 
 // for similar appraoch see: https://github.com/StrandedKitty/tiles-intersect/blob/master/src/index.js
-template<class T> 
+template <class T>
 void walkGrid(const Eigen::Vector3d &start, const Eigen::Vector3d &end, T &object)
 {
   Eigen::Vector3d direction = end - start;
   double max_length = direction.norm();
   Eigen::Vector3i p = Eigen::Vector3d(std::floor(start[0]), std::floor(start[1]), std::floor(start[2])).cast<int>();
-  const Eigen::Vector3i target = Eigen::Vector3d(std::floor(end[0]), std::floor(end[1]), std::floor(end[2])).cast<int>();
-  
+  const Eigen::Vector3i target =
+    Eigen::Vector3d(std::floor(end[0]), std::floor(end[1]), std::floor(end[2])).cast<int>();
+
   const Eigen::Vector3i step(sign(direction[0]), sign(direction[1]), sign(direction[2]));
   direction /= max_length;
-  float eps = 1e-10; // remove tiny about so grid walking doesn't exceed its boundary
+  float eps = 1e-10f;  // remove tiny about so grid walking doesn't exceed its boundary
   max_length -= eps;
   Eigen::Vector3d lengths, length_delta;
-  for (int j = 0; j<3; j++)
+  for (int j = 0; j < 3; j++)
   {
-    const double to = step[j] > 0 ? (double)p[j] + 1.0 - start[j] : start[j] - (double)p[j];       
+    const double to = step[j] > 0 ? (double)p[j] + 1.0 - start[j] : start[j] - (double)p[j];
     const double dir = std::max(std::numeric_limits<double>::epsilon(), std::abs(direction[j]));
     lengths[j] = to / dir;
     length_delta[j] = 1.0 / dir;
   }
   int ax = lengths[0] < lengths[1] && lengths[0] < lengths[2] ? 0 : (lengths[1] < lengths[2] ? 1 : 2);
   if (object(p, target, 0.0, lengths[ax], max_length))
-  {      
-    return; // only adding to one cell
+  {
+    return;  // only adding to one cell
   }
-  
-  while (lengths[ax] < max_length) 
+
+  while (lengths[ax] < max_length)
   {
     p[ax] += step[ax];
     const double in_length = lengths[ax];
@@ -404,19 +413,23 @@ void walkGrid(const Eigen::Vector3d &start, const Eigen::Vector3d &end, T &objec
     ax = lengths[0] < lengths[1] && lengths[0] < lengths[2] ? 0 : (lengths[1] < lengths[2] ? 1 : 2);
     if (object(p, target, in_length, lengths[ax], max_length))
     {
-      break; // only adding to one cell
-    }          
-  }     
+      break;  // only adding to one cell
+    }
+  }
 }
 
 inline int viewFile(const std::string &file_name, const std::string &file2_name = "")
-{ 
+{
   // Force Qt to use X11 platform (not Wayland). Can remove the first string if wish to use default render platform
   std::string command = "QT_QPA_PLATFORM=xcb " + std::string(R_VISTOOL) + " " + file_name;
   if (file2_name != "")
     command += " " + file2_name;
-  return system(command.c_str());
-}  
+  const auto status = system(command.c_str());
+  if (status != 0) {
+    std::cerr << "failed to view file with command: " << command << std::endl;
+  }
+  return status;
+}
 
 }  // namespace ray
 
diff --git a/tests/raytest/raytests.cpp b/tests/raytest/raytests.cpp
index ef3b99960..e7a3a3b79 100644
--- a/tests/raytest/raytests.cpp
+++ b/tests/raytest/raytests.cpp
@@ -38,25 +38,35 @@ namespace raytest
 
   /// Compare the statistical (1st and 2nd order) moments of the two ray clouds. This almost surely
   /// detects differing clouds, and always equal clouds, given a tolerance @c eps.
-  void compareMoments(const Eigen::ArrayXd &m1, const std::vector<double> &m2, double eps = 0.1)
+  std::size_t compareMoments(const Eigen::ArrayXd &m1, const std::vector<double> &m2, double eps = 0.1)
   {
-    for (size_t i = 0; i<m2.size(); i++)
+    size_t i;
+    for (i = 0; i<m2.size(); i++)
     {
-      EXPECT_GT(m1[i], m2[i]-eps);
-      EXPECT_LT(m1[i], m2[i]+eps);
+      EXPECT_NEAR(m1[i], m2[i], eps);
+      if (std::abs(m1[i] - m2[i]) > eps)
+      {
+        return i;
+      }
     }
+    return i;
   }
 
   /// Compare the statistical (1st and 2nd order) moments of the two ray clouds. This almost surely
   /// detects differing clouds, and always equal clouds, given a tolerance @c eps.
-  void compareMomentsPercentageError(const Eigen::ArrayXd &m1, const std::vector<double> &m2, double percentage = 5.0)
+  std::size_t compareMomentsPercentageError(const Eigen::ArrayXd &m1, const std::vector<double> &m2, double percentage = 5.0)
   {
-    double eps = 0.01 * percentage;
-    for (size_t i = 0; i<m2.size(); i++)
+    size_t i;
+    for (i = 0; i<m2.size(); i++)
     {
-      EXPECT_GE(m1[i], m2[i] * (1.0 - eps));
-      EXPECT_LE(m1[i], m2[i] * (1.0 + eps));
+      double eps = std::abs(m2[i]) * percentage * 0.01;
+      EXPECT_NEAR(m1[i], m2[i], eps);
+      if (std::abs(m1[i] - m2[i]) > eps)
+      {
+        return i;
+      }
     }
+    return i;
   }
 
   /// Creates two copies of the same room with a rotational difference, then aligns the first onto the second 
@@ -68,8 +78,8 @@ namespace raytest
     EXPECT_EQ(command("rayalign room.ply room2.ply"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_aligned.ply"));
-    compareMoments(cloud.getMoments(), {-0.0618268, -0.077552, 0.0531072, 7.58334e-08, 7.97642e-08, 1.93877e-08, -0.180532, -0.219257, 0.0654452, 2.47241, 2.08183, 1.28226, 17.539, 10.1994, 0.304682, 0.761892, 0.429502, 0.987362, 0.318932, 0.225742, 0.389901, 0.111705});  }
-
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.0618268, -0.077552, 0.0531072, 7.58334e-08, 7.97642e-08, 1.93877e-08, -0.180532, -0.219257, 0.0654452, 2.47241, 2.08183, 1.28226, 17.539, 10.1994, 0.304682, 0.761892, 0.429502, 0.987362, 0.318932, 0.225742, 0.389901, 0.111705}), cloud.getMoments().size());
+  }
   /// Colours a room according to the normal direction of the surfaces, comparing to the expected results
   TEST(Basic, RayColour)
   {
@@ -77,7 +87,7 @@ namespace raytest
     EXPECT_EQ(command("raycolour room.ply normal"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_coloured.ply"));
-    compareMoments(cloud.getMoments(), {-0.108066, -0.0410134, 0.052168, 7.05134e-08, 8.45038e-08, 1.93877e-08, -0.276144, -0.0760758, 0.065631, 2.42455, 2.13738, 1.28226, 17.539, 10.1994, 0.497919, 0.496369, 0.490293, 0.987362, 0.248361, 0.203648, 0.385192, 0.111705});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.108066, -0.0410134, 0.052168, 7.05134e-08, 8.45038e-08, 1.93877e-08, -0.276144, -0.0760758, 0.065631, 2.42455, 2.13738, 1.28226, 17.539, 10.1994, 0.497919, 0.496369, 0.490293, 0.987362, 0.248361, 0.203648, 0.385192, 0.111705}), cloud.getMoments().size());
   }
   
   /// Creates two rooms, with different transformations, then combines them, and compares to the expected result.
@@ -90,7 +100,7 @@ namespace raytest
     EXPECT_EQ(command("./raycombine min room.ply room2.ply 1 rays"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_combined.ply"));
-    compareMoments(cloud.getMoments(), {-0.0867714, -0.0679941, 0.546619, 0.0215326, 0.0272819, 0.499969, -0.305657, -0.186353, 0.582642, 2.95777, 2.47531, 1.63323, 17.4967, 10.1789, 0.305355, 0.763356, 0.427376, 0.979005, 0.318409, 0.225661, 0.389366, 0.143369});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.0867714, -0.0679941, 0.546619, 0.0215326, 0.0272819, 0.499969, -0.305657, -0.186353, 0.582642, 2.95777, 2.47531, 1.63323, 17.4967, 10.1789, 0.305355, 0.763356, 0.427376, 0.979005, 0.318409, 0.225661, 0.389366, 0.143369}), cloud.getMoments().size());
   }
   
   /// Creates a building with random seed 1, and compares to the expected results
@@ -99,7 +109,7 @@ namespace raytest
     EXPECT_EQ(command("raycreate building 1"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("building.ply"));
-    compareMoments(cloud.getMoments(), {-3.168, 16.5472, 7.04175, 3.28539, 18.0871, 2.96654, -3.19551, 16.564, 7.30826, 4.14359, 18.2463, 3.34431, 935.715, 540.236, 0.499997, 0.500408, 0.429416, 0.998894, 0.372134, 0.372389, 0.390499, 0.0332398});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-3.168, 16.5472, 7.04175, 3.28539, 18.0871, 2.96654, -3.19551, 16.564, 7.30826, 4.14359, 18.2463, 3.34431, 935.715, 540.236, 0.499997, 0.500408, 0.429416, 0.998894, 0.372134, 0.372389, 0.390499, 0.0332398}), cloud.getMoments().size());
   }
   
   /// Creates a forest and decimates it, comparing to the expected result 
@@ -111,7 +121,7 @@ namespace raytest
     EXPECT_TRUE(cloud.load("forest_decimated.ply"));
     // Below does not compare the time values (or the colour values, which are based on time here)
     // because spatial decimation does not constraint which time it picks points from.
-    compareMoments(cloud.getMoments(), {-0.222571, 1.08156, 1.67264, 6.00755, 5.78731, 0.508713, -0.202668, 1.09517, 2.6238, 6.0285, 5.85715, 3.22093, 69.0574, 35.2775, 0.48969, 0.498403, 0.443549, 1, 0.379062, 0.366963, 0.389535, 0});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.222571, 1.08156, 1.67264, 6.00755, 5.78731, 0.508713, -0.202668, 1.09517, 2.6238, 6.0285, 5.85715, 3.22093, 69.0574, 35.2775, 0.48969, 0.498403, 0.443549, 1, 0.379062, 0.366963, 0.389535, 0}), cloud.getMoments().size());
   }
 
   /// Creates a room, and calls denoise using a fixed distance threshols, and compares to expected result
@@ -121,7 +131,7 @@ namespace raytest
     EXPECT_EQ(command("raydenoise room.ply 3 cm"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_denoised.ply"));
-    compareMoments(cloud.getMoments(), {-0.108066, -0.0410134, 0.052168, 8.67026e-08, 8.81787e-08, 2.24394e-08, -0.464107, -0.113806, 0.161496, 2.82122, 2.34281, 1.35279, 17.81, 10.2005, 0.297047, 0.758802, 0.440232, 0.975166, 0.317215, 0.226682, 0.390971, 0.155618});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.108066, -0.0410134, 0.052168, 8.67026e-08, 8.81787e-08, 2.24394e-08, -0.464107, -0.113806, 0.161496, 2.82122, 2.34281, 1.35279, 17.81, 10.2005, 0.297047, 0.758802, 0.440232, 0.975166, 0.317215, 0.226682, 0.390971, 0.155618}), cloud.getMoments().size());
   }
 
   /// Creates two rooms, the second is decimated and transformed, then rayrestore is called to apply this transformation to
@@ -136,7 +146,7 @@ namespace raytest
     EXPECT_EQ(command("rayrestore room2_decimated.ply 10 cm room.ply"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_restored.ply"));
-    compareMoments(cloud.getMoments(), {2.07399, 0.575952, 3.05217, 7.85442e-08, 7.70963e-08, 1.93877e-08, 1.9391, 0.682169, 3.06563, 2.10068, 2.45642, 1.28226, 17.539, 10.1994, 0.304682, 0.761892, 0.429502, 0.987362, 0.318932, 0.225742, 0.389901, 0.111705});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {2.07399, 0.575952, 3.05217, 7.85442e-08, 7.70963e-08, 1.93877e-08, 1.9391, 0.682169, 3.06563, 2.10068, 2.45642, 1.28226, 17.539, 10.1994, 0.304682, 0.761892, 0.429502, 0.987362, 0.318932, 0.225742, 0.389901, 0.111705}), cloud.getMoments().size());
   }  
 
   /// Creates a forest and rotates it in all three axes, comparing to the expected result
@@ -146,7 +156,7 @@ namespace raytest
     EXPECT_EQ(command("rayrotate forest.ply 10,20,30"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("forest.ply"));
-    compareMoments(cloud.getMoments(), {-0.254879, 0.846076, 2.02322, 5.62648, 5.68622, 2.56306, 0.266873, 0.772016, 3.2888, 5.54595, 5.69021, 4.28394, 62.683, 36.1903, 0.514327, 0.504407, 0.413534, 1, 0.372377, 0.365965, 0.391709, 0});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.254879, 0.846076, 2.02322, 5.62648, 5.68622, 2.56306, 0.266873, 0.772016, 3.2888, 5.54595, 5.69021, 4.28394, 62.683, 36.1903, 0.514327, 0.504407, 0.413534, 1, 0.372377, 0.365965, 0.391709, 0}), cloud.getMoments().size());
   }  
 
   /// Creates a room and smooths this ray cloud, comparing to the expected result
@@ -156,7 +166,7 @@ namespace raytest
     EXPECT_EQ(command("raysmooth room.ply"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_smooth.ply"));
-    compareMoments(cloud.getMoments(), {-0.108066, -0.0410134, 0.052168, 7.05134e-08, 8.45038e-08, 1.93877e-08, -0.27615, -0.0761079, 0.0656267, 2.42413, 2.13691, 1.28163, 17.539, 10.1994, 0.304682, 0.761892, 0.429502, 0.987362, 0.318932, 0.225742, 0.389901, 0.111705});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.108066, -0.0410134, 0.052168, 7.05134e-08, 8.45038e-08, 1.93877e-08, -0.27615, -0.0761079, 0.0656267, 2.42413, 2.13691, 1.28163, 17.539, 10.1994, 0.304682, 0.761892, 0.429502, 0.987362, 0.318932, 0.225742, 0.389901, 0.111705}), cloud.getMoments().size());
   }  
 
   /// Creates a room, then splits it around a plane, comparing agaisnt the expected result
@@ -166,7 +176,7 @@ namespace raytest
     EXPECT_EQ(command("raysplit room.ply plane 0,0.1,1.5"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_outside.ply"));
-    compareMoments(cloud.getMoments(), {-0.467731, 1.05075, 1.43662, 2.20441, 1.60162, 0.106775, -0.77974, 1.03139, 1.57353, 3.67521, 2.64766, 0.485084, 17.3995, 10.279, 0.311066, 0.759795, 0.425206, 0.951355, 0.321609, 0.226785, 0.39073, 0.215125});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-0.467731, 1.05075, 1.43662, 2.20441, 1.60162, 0.106775, -0.77974, 1.03139, 1.57353, 3.67521, 2.64766, 0.485084, 17.3995, 10.279, 0.311066, 0.759795, 0.425206, 0.951355, 0.321609, 0.226785, 0.39073, 0.215125}), cloud.getMoments().size());
   }  
 
   /// Creates a room and runs raytransients, comparing the identified transients ray cloud to the expected results
@@ -176,7 +186,7 @@ namespace raytest
     EXPECT_EQ(command("raytransients min room.ply 1 rays"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("room_transient.ply"));
-    compareMoments(cloud.getMoments(), {-1.05406, -0.240721, -0.0629182, 5.05649e-08, 3.32941e-08, 2.54759e-08, 0.268724, -0.136746, -0.596782, 1.04798, 0.921776, 0.527205, 32.1452, 6.7491, 0.205871, 0.395641, 0.884296, 1, 0.225501, 0.296487, 0.153923, 0});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {-1.05406, -0.240721, -0.0629182, 5.05649e-08, 3.32941e-08, 2.54759e-08, 0.268724, -0.136746, -0.596782, 1.04798, 0.921776, 0.527205, 32.1452, 6.7491, 0.205871, 0.395641, 0.884296, 1, 0.225501, 0.296487, 0.153923, 0}), cloud.getMoments().size());
   }  
 
   /// Creates a forest and translates it in all three axes, comparing to the expected result
@@ -186,7 +196,7 @@ namespace raytest
     EXPECT_EQ(command("raytranslate forest.ply 10,20,30"), 0);
     ray::Cloud cloud;
     EXPECT_TRUE(cloud.load("forest.ply"));
-    compareMoments(cloud.getMoments(), {9.66298, 21.3454, 31.7177, 6.0926, 5.75511, 0.56438, 9.69155, 21.3605, 33.0883, 6.10555, 5.82564, 3.20507, 62.683, 36.1903, 0.514327, 0.504407, 0.413534, 1, 0.372377, 0.365965, 0.391709, 0});
+    EXPECT_EQ(compareMoments(cloud.getMoments(), {9.66298, 21.3454, 31.7177, 6.0926, 5.75511, 0.56438, 9.69155, 21.3605, 33.0883, 6.10555, 5.82564, 3.20507, 62.683, 36.1903, 0.514327, 0.504407, 0.413534, 1, 0.372377, 0.365965, 0.391709, 0}), cloud.getMoments().size());
   }
 
 #if RAYLIB_WITH_QHULL
@@ -197,7 +207,7 @@ namespace raytest
     EXPECT_EQ(command("raywrap terrain.ply upwards 1.0"), 0);
     ray::Mesh mesh;
     EXPECT_TRUE(ray::readPlyMesh("terrain_mesh.ply", mesh));
-    compareMoments(mesh.getMoments(), {0.0386662, -1.52168, -0.139079, 3.30621, 3.35391, 0.705937});
+    EXPECT_EQ(compareMoments(mesh.getMoments(), {0.0386662, -1.52168, -0.139079, 3.30621, 3.35391, 0.705937}), mesh.getMoments().size());
   }  
 
   /// Tests extraction of terrain and extraction of trees
@@ -207,24 +217,24 @@ namespace raytest
     EXPECT_EQ(command("rayextract terrain forest.ply"), 0);
     ray::Mesh mesh;
     EXPECT_TRUE(ray::readPlyMesh("forest_mesh.ply", mesh));
-    compareMoments(mesh.getMoments(), {-0.00147491, -0.00191917, -0.0617946, 5.77995, 5.80266, 0.0426993});
+    EXPECT_EQ(compareMoments(mesh.getMoments(), {-0.00147491, -0.00191917, -0.0617946, 5.77995, 5.80266, 0.0426993}), mesh.getMoments().size());
     EXPECT_EQ(command("rayextract trees forest.ply forest_mesh.ply"), 0);
 
     ray::ForestStructure forest;
     EXPECT_TRUE(forest.load("forest_trees.txt"));
-    compareMomentsPercentageError(forest.getMoments(), {20, 22.2172, 1058.42, 1.39416, 0.112794, 1.6403, 21918, 0, 75.75});
+    EXPECT_EQ(compareMomentsPercentageError(forest.getMoments(), {20, 22.2172, 1058.42, 1.39416, 0.112794, 2.2863, 21918, 0, 105.4}), forest.getMoments().size());
 
     EXPECT_EQ(command("rayextract forest forest.ply --ground forest_mesh.ply"), 0);
 
     ray::ForestStructure forest2;
     EXPECT_TRUE(forest2.load("forest_forest.txt"));
-    compareMoments(forest2.getMoments(), {11, 8.43829, 586.427, 1.40054, 0.200644, 0, 16697, 8.49419, 3.09917});
+    EXPECT_EQ(compareMomentsPercentageError(forest2.getMoments(), {11, 8.43829, 586.427, 1.40054, 0.200644, 0, 16697, 8.49419, 3.09917}), forest2.getMoments().size());
 
     EXPECT_EQ(command("rayextract trunks forest.ply"), 0);
 
     ray::ForestStructure forest3;
     EXPECT_TRUE(forest3.load("forest_trunks.txt"));
-    compareMoments(forest3.getMoments(), {21, 20.0797, 1124.61, 1.60427, 0.135159, 0, 0, 0, 0});
+    EXPECT_EQ(compareMomentsPercentageError(forest3.getMoments(), {21, 20.0797, 1124.61, 1.60427, 0.135159, 0, 0, 0, 0}), forest3.getMoments().size());
   }  
 #endif  // RAYLIB_WITH_QHULL
 } // raytest