WIP blue_noise impl

Makefile

@@ -1,12 +1,14 @@
 
-COMMON_FLAGS=-Wall -Wextra -Wpedantic
+COMMON_FLAGS=-Wall -Wextra -Wpedantic -std=c++17
 ifdef DEBUG
 	CPPFLAGS=${COMMON_FLAGS} -g -O0
 else
 	CPPFLAGS=${COMMON_FLAGS} -O3 -DNDEBUG
 endif
 
-SOURCES=src/main.cpp
+SOURCES= \
+	src/main.cpp \
+	src/blue_noise.cpp
 OBJECTS=${subst .cpp,.o,${SOURCES}}
 
 all: Dithering

src/blue_noise.cpp

@@ -0,0 +1,37 @@
+#include "blue_noise.hpp"
+
+#include <random>
+
+std::vector<bool> dither::blue_noise(std::size_t width, std::size_t height) {
+    std::size_t count = width * height;
+    std::vector<std::size_t> dither_array;
+    dither_array.resize(count);
+
+    std::vector<bool> pbp; // Prototype Binary Pattern
+    pbp.resize(count);
+
+    std::default_random_engine re(std::random_device{}());
+    std::uniform_int_distribution<std::size_t> dist(0, count - 1);
+
+    // initialize pbp
+    for(std::size_t i = 0; i < count; ++i) {
+        if(i < count / 2) {
+            pbp[i] = true;
+        } else {
+            pbp[i] = false;
+        }
+    }
+    // randomize pbp
+    for(std::size_t i = 0; i < count-1; ++i) {
+        decltype(dist)::param_type range{i+1, count-1};
+        std::size_t ridx = dist(re, range);
+        // probably can't use std::swap since using std::vector<bool>
+        bool temp = pbp[i];
+        pbp[i] = pbp[ridx];
+        pbp[ridx] = temp;
+    }
+
+
+
+    return {};
+}

src/blue_noise.hpp

@@ -0,0 +1,54 @@
+#ifndef BLUE_NOISE_HPP
+#define BLUE_NOISE_HPP
+
+#include <vector>
+#include <tuple>
+#include <cmath>
+
+namespace dither {
+
+std::vector<bool> blue_noise(std::size_t width, std::size_t height);
+
+namespace internal {
+    inline std::size_t twoToOne(std::size_t x, std::size_t y, std::size_t width) {
+        return x + y * width;
+    }
+
+    inline std::tuple<std::size_t, std::size_t> oneToTwo(std::size_t i, std::size_t width) {
+        return {i % width, i / width};
+    }
+
+    constexpr double mu_squared = 1.5 * 1.5;
+
+    inline double gaussian(double x, double y) {
+        return std::exp(-std::sqrt(x*x + y*y)/(2*mu_squared));
+    }
+
+    inline double filter(
+            const std::vector<bool>& pbp,
+            std::size_t x, std::size_t y,
+            std::size_t width) {
+        double sum = 0.0;
+
+        // Should be range -M/2 to M/2, but size_t cannot be negative, so range
+        // is 0 to M.
+        // p' = (M + x - (p - M/2)) % M = (3M/2 + x - p) % M
+        // q' = (M + y - (q - M/2)) % M = (3M/2 + y - q) % M
+        for(std::size_t q = 0; q <= width; ++q) {
+            std::size_t q_prime = (3 * width / 2 + y - q) % width;
+            for(std::size_t p = 0; p <= width; ++p) {
+                std::size_t p_prime = (3 * width / 2 + x - p) % width;
+                bool pbp_value = pbp[twoToOne(p_prime, q_prime, width)];
+                if(pbp_value) {
+                    sum += gaussian((double)p - width/2.0, (double)q - width/2.0);
+                }
+            }
+        }
+
+        return sum;
+    }
+} // namespace dither::internal
+
+} // namespace dither
+
+#endif