Fix camera/rendering + apparent physical bias (#16)

* Introduce verlet.h (Verlet integration)

To investigate velocity bias issue.

Integrator does not seem to be cause of problem.

Why: Both Yoshida and Verlet have same problem.

* WIP: Add routine to test circle-rectangle collision

I'm doing this because Raylib's built-in routine for collision detection seems to be erroneous. (Additionally, it's never performed well. In my 1,000 particle test runs, circle testing has been one of the more consistent bottlenecks despite the O(n^2) vs. O(n) time complexities [n = number of particles]. According to my inexpert reading of the profiler's reports, because of the time the CPU spends in converting between integers and floating-points as part of the implementation detail of the routine, which may be causing certain unknown effects that propagate outward.)

Or, more precisely, disk-rectangle collision (two-dimensional figures), which means that, if one figure is contained in the other, they are said to be colliding (intersecting).

(As opposed to the one-dimensional case where it would be considered not intersecting because no point of intersection is created on the respective curves.)

Remaining:
1. At least one affirmative case.
2. Coverage of every branch.
3. Rotated version of each test.

Later, I will compute the branch coverage for all tests.

* Remove hypot from disk-rectangle intersection detection

The hypot calls were estimated to be O(n) where n is the number of particles in the previous implementation.

I will add in the tests soon when I find the time to do it.

* Revise disk and (area) rectangle testing routine

* Fix "flickering near boundaries" issue

The erroneous "CheckCollisionCircleRec" (Raylib) routine was replaced with a custom routine.

* circle.h: Remove unintended math function calls

The signbit calls were observed to be calling a system function that would extract the sign bit from the floating-point number.

* demo: Remove collision bias

Refresh the Halton disk at every frame to remove statistical bias when particles are partially colliding

* Fix last remaining camera issue

I ignored the screen offset, which is now taken into account.

* Table.h: Revert to Yoshida integration

Verlet integration was only used for exposition of a problem.

demo/Table.h

@@ -6,13 +6,14 @@
 #include <circle.h>
 #include <kahan.h>
 #include <newton.h>
+#include <verlet.h>
 #include <yoshida.h>
 
 namespace phy {
 
 struct Particle {
-  /// @brief Kinematic properties (position: y0, velocity: y1).
-  dyn::Yoshida<float> kin;
+  /// @brief Kinematic properties.
+  std::complex<float> xy, v;
 
   /// @brief Mass and radius.
   float mass = 1.0f, radius = 1.0f;
@@ -28,11 +29,11 @@ struct Particle {
   /// @param r Radius, positive.
   constexpr Particle(std::complex<float> xy, std::complex<float> v, float m = 1,
                      float r = 1)
-      : kin{xy, v}, mass{m}, radius{r} {}
+      : xy{xy}, v{v}, mass{m}, radius{r} {}
 
   /// @brief Construct a circle that represents this particle.
   [[nodiscard]] constexpr dyn::Circle<float> circle() const {
-    return {kin.y0, radius};
+    return {xy, radius};
   }
 };
 
@@ -68,7 +69,10 @@ class Table : public std::vector<Particle> {
       };
       // step(): Integrate.
       // (Calls accel() above a few times with slightly different xy.)
-      copy[i].kin.step(dt, accel);
+      auto step = dyn::Yoshida<>{p.xy, p.v};
+      step.step(dt, accel);
+      auto &q = copy[i];
+      q.xy = step.y0, q.v = step.y1;
     }
     *this = copy;
   }
@@ -79,10 +83,25 @@ class Table : public std::vector<Particle> {
     dyn::Kahan<C<>> cm;
     dyn::Kahan<double> m;
     for (auto &&p : *this)
-      cm += double(p.mass) * C<>(p.kin.y0), m += double(p.mass);
+      cm += double(p.mass) * C<>(p.xy), m += double(p.mass);
     auto c = C<float>(cm() / m());
     for (auto &&p : *this)
-      p.kin.y0 -= c;
+      p.xy -= c;
+  }
+
+  /// @brief Refresh the "disk" used for parts of the calculation.
+  void refresh_disk() noexcept { gr.refresh_disk(); }
+
+  /// @brief Test whether the simulation is in "good state."
+  bool good() noexcept {
+    auto constexpr goodf = [](float f) { return std::isfinite(f); };
+    auto constexpr goodc = [](std::complex<float> f) {
+      return std::isfinite(f.real()) && std::isfinite(f.imag());
+    };
+    for (auto &&p : *this)
+      if (!goodc(p.xy) || !goodc(p.v))
+        return false;
+    return true;
   }
 };
 

demo/main.cpp

@@ -1,12 +1,12 @@
 #include <raylib.h>
 
 #include <algorithm>
-#include <cfenv>
 #include <cmath>
 #include <complex>
 #include <sstream>
 #include <vector>
 
+#include <circle.h>
 #include <halton.h>
 
 #include "Table.h"
@@ -19,7 +19,7 @@ static int do_main() {
   auto constexpr RADIUS = 0.05f;
   auto constexpr MASS = 1.0f;
   auto constexpr too_far = [](Particle &p) {
-    return std::abs(p.kin.y0) > 5'000.0f;
+    return std::abs(p.xy) > 5'000.0f;
   };
 
   // Quasi-random number generator.
@@ -183,7 +183,7 @@ static int do_main() {
       auto o = std::complex<float>(n.x, n.y);
       auto p = random_particle();
       // Set location
-      p.kin.y0 = o;
+      p.xy = o;
       table.push_back(p);
 
       interactive.spawned_last_frame = true;
@@ -195,19 +195,18 @@ static int do_main() {
     }
 
   simulate:
-    // Clear floating point exceptions.
-    std::feclearexcept(FE_ALL_EXCEPT);
-
     // Do the simulation!
     if (!interactive.freeze)
       table.step(dt);
 
+    table.refresh_disk();
+
     // Don't be moving the particles while it's being manipulated.
-    if (!interactive.spawned_last_frame && !interactive.other_user_manip)
-      table.center();
+    // if (!interactive.spawned_last_frame && !interactive.other_user_manip)
+    // table.center();
 
     // Inspect floating point exceptions.
-    if (std::fetestexcept(FE_DIVBYZERO | FE_INVALID))
+    if (!table.good())
       // NaN or infinity somewhere. Reset the simulation.
       goto reset_sim;
 
@@ -216,16 +215,21 @@ static int do_main() {
 
     // Draw all particles in the frame.
     BeginMode2D(cam);
-    for (auto &&p : table) {
-      auto cp = p.circle();
-      Vector2 xy{cp.real(), cp.imag()};
-      // Should I issue the draw call or not? Tell it here:
-      // co stores screen dimensions calculated earlier.
-      auto co = cam.offset;
-      // x, y, width, then height.
-      Rectangle screen{-co.x / 2, -co.y / 2, co.x, co.y};
-      if (CheckCollisionCircleRec(xy, cp.radius, screen))
-        DrawCircleV(xy, cp.radius, WHITE);
+    {
+      // Compute the screen's less-less (ll) and greater-greater (gg) corners in
+      // the world coordinate system so that draw calls would only be issued for
+      // certainly visible particles.
+      std::complex<float> scwh{float(GetScreenWidth()),
+                               float(GetScreenHeight())};
+      std::complex<float> scof{cam.offset.x, cam.offset.y};
+      auto scsc = (scwh - scof) / cam.zoom;
+      std::complex<float> sctg{cam.target.x, cam.target.y};
+      auto ll = sctg - scof / cam.zoom, gg = sctg + scsc;
+      for (auto &&p : table) {
+        auto cp = p.circle();
+        if (dyn::disk_arrect_isct(cp, ll, gg))
+          DrawCircleV({cp.real(), cp.imag()}, cp.radius, WHITE);
+      }
     }
     EndMode2D();
 

dyn/CMakeLists.txt

@@ -5,5 +5,6 @@ add_library(dyn INTERFACE
         kahan.h
         halton.h
         newton.h
-        circle.h)
+        circle.h
+        verlet.h)
 target_include_directories(dyn INTERFACE .)

dyn/circle.h

@@ -5,6 +5,7 @@
 /// @brief Represent a circle by the center and radius.
 
 #include <complex>
+#include <numbers>
 
 namespace dyn {
 
@@ -22,6 +23,55 @@ template <typename F = float> struct Circle : public std::complex<F> {
       : std::complex<F>(center), radius(radius) {}
 };
 
+/// @brief Decide whether at least one intersection (point) exists between the
+/// area of the disk centered at the origin with the given radius and the area
+/// of the given rectangle (degenerate cases are unspecified).
+/// @param ll The less-less corner of the rectangle.
+/// @param gg The greater-greater corner of the rectangle.
+template <typename F = float>
+bool origin_disk_arrect_isct(F radius, std::complex<F> ll,
+                             std::complex<F> gg) noexcept {
+  // sin(45 deg) [also cos(45 deg)].
+  constexpr auto sc45 = std::numbers::sqrt2_v<F> / F(2);
+
+  // Decide whether the given rectangle and a square centered at (0,0) having a
+  // "radius" (half side length) of `radius` (thus a full side length of twice
+  // the `radius`).
+  auto arrec_arsq_isct = [gg, ll](auto radius) {
+    return ll.real() < radius && -radius < gg.real() && ll.imag() < radius &&
+           -radius < gg.imag();
+  };
+
+  // The rectangle does not touch the bounding square of the circle.
+  if (!arrec_arsq_isct(radius))
+    return false;
+
+  // The rectangle touches the interior square of the circle.
+  if (arrec_arsq_isct(radius * sc45))
+    return true;
+
+  // At least one of the four corners of the given rectangle is in the disk.
+  F cc[4] = {std::abs(ll), std::abs(gg),
+             std::abs(std::complex<F>{ll.real(), gg.imag()}),
+             std::abs(std::complex<F>{gg.real(), ll.imag()})};
+  for (auto &&c : cc)
+    if (c < radius)
+      return true;
+
+  // They aren't touching.
+  return false;
+}
+
+/// @brief Decide whether at least one intersection (point) exists between a
+/// circular disk and the area of a rectangle (degenerate cases are
+/// unspecified).
+template <typename F = float>
+bool disk_arrect_isct(Circle<F> circ, std::complex<F> ll,
+                      std::complex<F> gg) noexcept {
+  ll -= circ, gg -= circ;
+  return origin_disk_arrect_isct(circ.radius, ll, gg);
+}
+
 } // namespace dyn
 
 #endif // GRASS_CIRCLE_H

dyn/newton.h

@@ -23,6 +23,14 @@ template <typename F = float, unsigned short N_MONTE = 30> class Gravity {
   /// for Monte Carlo integration in the case of overlapping circles.
   std::array<std::complex<F>, N_MONTE> disk{};
 
+  // Each Halton sequence (a kind of low-discrepancy sequence) creates an
+  // evenly spaced set of points on the unit interval (0, 1); unlike the
+  // uniform distribution, however, the points look "uniformly distributed"
+  // (number of points being mostly proportional to length of any subset)
+  // even for a finite sample of points.
+  Halton<F, 2> h2;
+  Halton<F, 3> h3;
+
 public:
   /// @brief Create an instance with a quasi-random internal state.
   constexpr Gravity() { refresh_disk(); }
@@ -54,14 +62,6 @@ template <typename F = float, unsigned short N_MONTE = 30> class Gravity {
   /// case of intersecting circles) with new evenly distributed points on the
   /// unit disk centered about the origin. Call time to time.
   constexpr void refresh_disk() {
-    // Each Halton sequence (a kind of low-discrepancy sequence) creates an
-    // evenly spaced set of points on the unit interval (0, 1); unlike the
-    // uniform distribution, however, the points look "uniformly distributed"
-    // (number of points being mostly proportional to length of any subset)
-    // even for a finite sample of points.
-    Halton<F, 2> h2;
-    Halton<F, 3> h3;
-
     // Fill `disk` with random points on unit disk centered about origin
     // by rejection sampling.
     for (auto &&p : disk) {

dyn/verlet.h

@@ -0,0 +1,21 @@
+#ifndef GRASS_VERLET_H
+#define GRASS_VERLET_H
+
+#include <complex>
+
+namespace dyn {
+
+template <typename F = float> struct Verlet {
+  std::complex<F> y0, y1;
+
+  template <typename A> void step(F h, A y2) {
+    auto a = y2(y0);
+    y0 += h * y1 + h * h * F(0.5) * a;
+    auto b = y2(y0);
+    y1 += h * F(0.5) * (a + b);
+  }
+};
+
+} // namespace dyn
+
+#endif // GRASS_VERLET_H

tests/CMakeLists.txt

@@ -11,7 +11,8 @@ include(GoogleTest)
 
 # Now simply link against gtest or gtest_main as needed. Eg
 add_executable(units yoshida_test.cpp
-        newton_test.cpp)
+        newton_test.cpp
+        circle_test.cpp)
 target_precompile_headers(units INTERFACE "gtest/gtest.h")
 target_link_libraries(units gtest_main dyn)
 gtest_discover_tests(units)

tests/circle_test.cpp

@@ -0,0 +1,58 @@
+#include "gtest/gtest.h"
+
+#include <complex>
+
+#include <circle.h>
+
+/// @brief See `CircleTest0.png` for reference.
+class CircleTest0 : public testing::Test {
+protected:
+  // xy, r; rectangles (2 corners)
+  dyn::Circle<float> const circle{0, 2.4f};
+  std::complex<float> const less_less{-4.0f, -4.0f};
+  std::complex<float> const greater_greater{-2.0f, -2.0f};
+};
+
+TEST_F(CircleTest0, Disjoint0) {
+  ASSERT_FALSE(dyn::disk_arrect_isct(circle, less_less, greater_greater));
+}
+
+class CircleTest1 : public testing::Test {
+protected:
+  dyn::Circle<float> const circle{0, 2.0f};
+  std::complex<float> const less_less{-5.0f, 1.0f},
+      greater_greater{-1.0f, 5.0f};
+};
+
+TEST_F(CircleTest1, In0) {
+  ASSERT_TRUE(dyn::disk_arrect_isct(circle, less_less, greater_greater));
+}
+
+class CircleTest2 : public testing::Test {
+protected:
+  dyn::Circle<float> const circle{0, 3.0f};
+  std::complex<float> const less_less{0.6f, 2.7f}, greater_greater{1.8f, 3.6f};
+};
+
+TEST_F(CircleTest2, In0) {
+  ASSERT_TRUE(dyn::disk_arrect_isct(circle, less_less, greater_greater));
+}
+
+class CircleTest3 : public testing::Test {
+protected:
+  dyn::Circle<float> const circle{0, 4.0f};
+  std::complex<float> const c{-2.0f, -2.0f}, e{0, -1.0f}, l{-2.0f, -5.0f},
+      n{5.0f, 3.0f}, s{7.0f, 2.0f}, u{9.0f, 3.0f};
+};
+
+TEST_F(CircleTest3, In0) {
+  ASSERT_TRUE(dyn::disk_arrect_isct(circle, c, e));
+}
+
+TEST_F(CircleTest3, In1) {
+  ASSERT_TRUE(dyn::disk_arrect_isct(circle, l, n));
+}
+
+TEST_F(CircleTest3, Out0) {
+  ASSERT_FALSE(dyn::disk_arrect_isct(circle, s, u));
+}