C++ sample code for The Barnsley Fern

contents/barnsley/barnsley.md

@@ -125,6 +125,8 @@ The biggest differences between the two code implementations is that the Barnsle
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/barnsley.jl)
+{% sample lang="cpp" %}
+[import, lang:"cpp"](code/c++/barnsley.cpp)
 {% endmethod %}
 
 ### Bibliography

contents/barnsley/code/c++/barnsley.cpp

@@ -0,0 +1,118 @@
+#include <array>
+#include <cassert>
+#include <fstream>
+#include <random>
+
+using Vec2 = std::array<double, 2>;
+using Vec3 = std::array<double, 3>;
+using Row = std::array<double, 3>;
+using Op = std::array<Row, 3>;
+
+constexpr auto OpN = 4;
+
+template <int N>
+auto operator+(std::array<double, N> x, std::array<double, N> y) {
+  for (auto i = 0U; i < N; ++i)
+    x[i] += y[i];
+  return x;
+}
+
+template <int N>
+auto operator*(double k, std::array<double, N> v) {
+  for (auto i = 0U; i < N; ++i)
+    v[i] *= k;
+  return v;
+}
+
+template <int N>
+auto operator*(std::array<double, N> v, double k) {
+  return k * v;
+}
+
+auto operator*(const Op& x, const Vec3& y) {
+  auto ret = Vec3{};
+  for (auto i = 0; i < 3; ++i) {
+    ret[i] = 0;
+    for (auto j = 0; j < 3; ++j)
+      ret[i] += y[j] * x[i][j];
+  }
+  return ret;
+}
+
+// Returns a pseudo-random number generator
+std::default_random_engine& rng() {
+  // Initialize static pseudo-random engine with non-deterministic random seed
+  static std::default_random_engine randEngine(std::random_device{}());
+  return randEngine;
+}
+
+// Returns a random double in [0, 1)
+double drand() {
+  return std::uniform_real_distribution<double>(0.0, 1.0)(rng());
+}
+
+// This is a function that reads in the Hutchinson operator and
+// corresponding
+//     probabilities and outputs a randomly selected transform
+// This works by choosing a random number and then iterating through all
+//     probabilities until it finds an appropriate bin
+auto select_array(
+    const std::array<Op, OpN>& hutchinson_op,
+    const std::array<double, OpN>& probabilities) {
+
+  // random number to be binned
+  auto rnd = drand();
+
+  // This checks to see if a random number is in a bin, if not, that
+  //     probability is subtracted from the random number and we check the
+  //     next bin in the list
+  for (auto i = 0U; i < probabilities.size(); ++i) {
+    if (rnd < probabilities[i])
+      return hutchinson_op[i];
+    rnd -= probabilities[i];
+  }
+  assert(!"check if probabilities adding up to 1");
+}
+
+// This is a general function to simulate a chaos game
+// n is the number of iterations
+// initial_location is the the starting point of the chaos game
+// hutchinson_op is the set of functions to iterate through
+// probabilities is the set of probabilities corresponding to the likelihood
+//     of choosing their corresponding function in hutchinson_op
+auto chaos_game(
+    size_t n,
+    Vec2 initial_location,
+    const std::array<Op, OpN>& hutchinson_op,
+    const std::array<double, OpN>& probabilities) {
+
+  // Initializing the output array and the initial point
+  auto output_points = std::vector<Vec2>{};
+
+  // extending point to 3D for affine transform
+  auto point = Vec3{initial_location[0], initial_location[1], 1};
+
+  for (auto i = 0U; i < n; ++i) {
+    output_points.push_back(Vec2{point[0], point[1]});
+    point = select_array(hutchinson_op, probabilities) * point;
+  }
+
+  return output_points;
+}
+
+int main() {
+
+  const std::array barnsley_hutchinson = {
+      Op{Row{0.0, 0.0, 0.0}, Row{0.0, 0.16, 0.0}, Row{0.0, 0.0, 1.0}},
+      Op{Row{0.85, 0.04, 0.0}, Row{-0.04, 0.85, 1.60}, Row{0.0, 0.0, 1.0}},
+      Op{Row{0.20, -0.26, 0.0}, Row{0.23, 0.22, 1.60}, Row{0.0, 0.0, 1.0}},
+      Op{Row{-0.15, 0.28, 0.0}, Row{0.26, 0.24, 0.44}, Row{0.0, 0.0, 1.0}}};
+
+  const std::array barnsley_probabilities = {0.01, 0.85, 0.07, 0.07};
+  auto output_points = chaos_game(
+      10'000, Vec2{0, 0}, barnsley_hutchinson, barnsley_probabilities);
+
+  std::ofstream ofs("out.dat");
+  for (auto pt : output_points)
+    ofs << pt[0] << '\t' << pt[1] << '\n';
+}