diff --git a/Makefile b/Makefile
index 9e06ccc..449bec9 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,5 @@
 
-COMMON_FLAGS=-Wall -Wextra -Wpedantic -std=c++17
+COMMON_FLAGS=-Wall -Wextra -Wpedantic -std=c++17 -lpthread
 ifdef DEBUG
 	CPPFLAGS=${COMMON_FLAGS} -g -O0
 else
diff --git a/src/blue_noise.cpp b/src/blue_noise.cpp
index 63a2488..8c8dbcb 100644
--- a/src/blue_noise.cpp
+++ b/src/blue_noise.cpp
@@ -1,11 +1,130 @@
 #include "blue_noise.hpp"
 
 #include <random>
+#include <cassert>
+#include <iostream>
+#include <condition_variable>
+#include <mutex>
+#include <thread>
+#include <chrono>
 
-std::vector<bool> dither::blue_noise(std::size_t width, std::size_t height) {
+#ifndef NDEBUG
+# include <cstdio>
+#endif
+
+void dither::internal::recursive_apply_radius(
+        std::size_t idx, std::size_t width, std::size_t height,
+        std::size_t radius, const std::function<bool(std::size_t)>& fn) {
+    std::unordered_set<std::size_t> visited;
+#ifndef NDEBUG
+    if(recursive_apply_radius_impl(idx, width, height, radius, fn, visited)) {
+        puts("recursive_apply_radius_impl found result");
+    } else {
+        puts("recursive_apply_radius_impl did NOT find result");
+    }
+#else
+    recursive_apply_radius_impl(idx, width, height, radius, fn, visited);
+#endif
+}
+
+bool dither::internal::recursive_apply_radius_impl(
+        std::size_t idx, std::size_t width, std::size_t height,
+        std::size_t radius, const std::function<bool(std::size_t)>& fn,
+        std::unordered_set<std::size_t>& visited) {
+    if(fn(idx)) {
+        return true;
+    }
+    std::size_t x, y, temp;
+    std::tie(x, y) = oneToTwo(idx, width);
+
+    if(x + 1 < width) {
+        temp = idx + 1;
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    temp, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    } else {
+        temp = twoToOne(0, y, width);
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    twoToOne(0, y, width),
+                    width, height, radius - 1,
+                    fn, visited)) {
+                return true;
+            }
+        }
+    }
+
+    if(x > 0) {
+        temp = idx - 1;
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    idx - 1, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    } else {
+        temp = twoToOne(width - 1, y, width);
+        if(visited.find(temp) == visited.end()) {
+            if(recursive_apply_radius_impl(
+                    temp, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    }
+
+    if(y + 1 < height) {
+        temp = idx + width;
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    temp, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    } else {
+        temp = twoToOne(x, 0, width);
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    temp, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    }
+
+    if(y > 0) {
+        temp = idx - width;
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    temp, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    } else {
+        temp = twoToOne(x, height - 1, width);
+        if(visited.find(temp) == visited.end()) {
+            visited.insert(temp);
+            if(recursive_apply_radius_impl(
+                    temp, width, height, radius - 1, fn, visited)) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+
+std::vector<bool> dither::blue_noise(std::size_t width, std::size_t height, std::size_t threads) {
     std::size_t count = width * height;
-    std::vector<std::size_t> dither_array;
-    dither_array.resize(count);
+    std::vector<double> filter_out;
+    filter_out.resize(count);
 
     std::vector<bool> pbp; // Prototype Binary Pattern
     pbp.resize(count);
@@ -13,9 +132,11 @@ std::vector<bool> dither::blue_noise(std::size_t width, std::size_t height) {
     std::default_random_engine re(std::random_device{}());
     std::uniform_int_distribution<std::size_t> dist(0, count - 1);
 
+    const std::size_t pixel_count = count * 4 / 10;
+
     // initialize pbp
     for(std::size_t i = 0; i < count; ++i) {
-        if(i < count / 2) {
+        if(i < pixel_count) {
             pbp[i] = true;
         } else {
             pbp[i] = false;
@@ -30,8 +151,203 @@ std::vector<bool> dither::blue_noise(std::size_t width, std::size_t height) {
         pbp[i] = pbp[ridx];
         pbp[ridx] = temp;
     }
+//#ifndef NDEBUG
+    printf("Inserting %ld pixels into image of max count %ld\n", pixel_count, count);
+    // generate image from randomized pbp
+    FILE *random_noise_image = fopen("random_noise.pbm", "w");
+    fprintf(random_noise_image, "P1\n%ld %ld\n", width, height);
+    for(std::size_t y = 0; y < height; ++y) {
+        for(std::size_t x = 0; x < width; ++x) {
+            fprintf(random_noise_image, "%d ", pbp[internal::twoToOne(x, y, width)] ? 1 : 0);
+        }
+        fputc('\n', random_noise_image);
+    }
+    fclose(random_noise_image);
+//#endif
 
+//#ifndef NDEBUG
+    std::size_t iterations = 0;
+//#endif
 
+    while(true) {
+//#ifndef NDEBUG
+//        if(++iterations % 10 == 0) {
+            printf("Iteration %ld\n", ++iterations);
+//        }
+//#endif
+        // get filter values
+        if(threads == 1) {
+            for(std::size_t y = 0; y < height; ++y) {
+                for(std::size_t x = 0; x < width; ++x) {
+                    filter_out[internal::twoToOne(x, y, width)] =
+                        internal::filter(pbp, x, y, width, height);
+                }
+            }
+        } else {
+            if(threads == 0) {
+                threads = 10;
+            }
+            std::size_t active_count = 0;
+            std::mutex cv_mutex;
+            std::condition_variable cv;
+            for(std::size_t i = 0; i < count; ++i) {
+                {
+                    std::unique_lock lock(cv_mutex);
+                    active_count += 1;
+                }
+                std::thread t([] (std::size_t *ac, std::mutex *cvm,
+                            std::condition_variable *cv, std::size_t i,
+                            const std::vector<bool> *pbp, std::size_t width,
+                            std::size_t height, std::vector<double> *fout) {
+                        std::size_t x, y;
+                        std::tie(x, y) = internal::oneToTwo(i, width);
+                        (*fout)[i] = internal::filter(*pbp, x, y, width, height);
+                        std::unique_lock lock(*cvm);
+                        *ac -= 1;
+                        cv->notify_all();
+                    },
+                    &active_count, &cv_mutex, &cv, i, &pbp, width, height, &filter_out);
+                t.detach();
 
-    return {};
+                std::unique_lock lock(cv_mutex);
+                while(active_count >= threads) {
+#ifndef NDEBUG
+//                    std::cout << "0, active_count = " << active_count
+//                        << ", pre wait_for" << std::endl;
+#endif
+                    cv.wait_for(lock, std::chrono::seconds(1));
+#ifndef NDEBUG
+//                    std::cout << "0, active_count = " << active_count
+//                        << ", post wait_for" << std::endl;
+#endif
+                }
+            }
+            std::unique_lock lock(cv_mutex);
+            while(active_count > 0) {
+                cv.wait_for(lock, std::chrono::seconds(1));
+            }
+        }
+
+#ifndef NDEBUG
+//        for(std::size_t i = 0; i < count; ++i) {
+//            std::size_t x, y;
+//            std::tie(x, y) = internal::oneToTwo(i, width);
+//            printf("%ld (%ld, %ld): %f\n", i, x, y, filter_out[i]);
+//        }
+#endif
+
+        std::size_t min, max, min_zero, max_one;
+        std::tie(min, max) = internal::filter_minmax(filter_out);
+        if(!pbp[max]) {
+            max_one = internal::get_one_or_zero(pbp, true, max, width, height);
+#ifndef NDEBUG
+            std::cout << "Post get_one(...)" << std::endl;
+#endif
+        } else {
+            max_one = max;
+        }
+        if(!pbp[max_one]) {
+            std::cerr << "ERROR: Failed to find pbp[max] one" << std::endl;
+            break;
+        }
+
+        if(pbp[min]) {
+            min_zero = internal::get_one_or_zero(pbp, false, min, width, height);
+#ifndef NDEBUG
+            std::cout << "Post get_zero(...)" << std::endl;
+#endif
+        } else {
+            min_zero = min;
+        }
+        if(pbp[min_zero]) {
+            std::cerr << "ERROR: Failed to find pbp[min] zero" << std::endl;
+            break;
+        }
+
+        // remove 1
+        pbp[max_one] = false;
+
+        // get filter values again
+        if(threads == 1) {
+            for(std::size_t y = 0; y < height; ++y) {
+                for(std::size_t x = 0; x < width; ++x) {
+                    filter_out[internal::twoToOne(x, y, width)] =
+                        internal::filter(pbp, x, y, width, height);
+                }
+            }
+        } else {
+            std::size_t active_count = 0;
+            std::mutex cv_mutex;
+            std::condition_variable cv;
+            for(std::size_t i = 0; i < count; ++i) {
+                {
+                    std::unique_lock lock(cv_mutex);
+                    active_count += 1;
+                }
+                std::thread t([] (std::size_t *ac, std::mutex *cvm,
+                            std::condition_variable *cv, std::size_t i,
+                            const std::vector<bool> *pbp, std::size_t width,
+                            std::size_t height, std::vector<double> *fout) {
+                        std::size_t x, y;
+                        std::tie(x, y) = internal::oneToTwo(i, width);
+                        (*fout)[i] = internal::filter(*pbp, x, y, width, height);
+                        std::unique_lock lock(*cvm);
+                        *ac -= 1;
+                        cv->notify_all();
+                    },
+                    &active_count, &cv_mutex, &cv, i, &pbp, width, height, &filter_out);
+                t.detach();
+
+                std::unique_lock lock(cv_mutex);
+                while(active_count >= threads) {
+#ifndef NDEBUG
+//                    std::cout << "1, active_count = " << active_count
+//                        << ", pre wait_for" << std::endl;
+#endif
+                    cv.wait_for(lock, std::chrono::seconds(1));
+#ifndef NDEBUG
+//                    std::cout << "1, active_count = " << active_count
+//                        << ", post wait_for" << std::endl;
+#endif
+                }
+            }
+            std::unique_lock lock(cv_mutex);
+            while(active_count > 0) {
+                cv.wait_for(lock, std::chrono::seconds(1));
+            }
+        }
+
+        // get second buffer's min
+        std::size_t second_min;
+        std::tie(second_min, std::ignore) = internal::filter_minmax(filter_out);
+        if(pbp[second_min]) {
+            second_min = internal::get_one_or_zero(pbp, false, second_min, width, height);
+            if(pbp[second_min]) {
+                std::cerr << "ERROR: Failed to find pbp[second_min] zero" << std::endl;
+                break;
+            }
+        }
+
+        if(min_zero != second_min) {
+            pbp[max_one] = true;
+            break;
+        } else {
+            pbp[min_zero] = true;
+        }
+    }
+
+//#ifndef NDEBUG
+    // generate blue_noise image from pbp
+    FILE *blue_noise_image = fopen("blue_noise.pbm", "w");
+    fprintf(blue_noise_image, "P1\n%ld %ld\n", width, height);
+    for(std::size_t y = 0; y < height; ++y) {
+        for(std::size_t x = 0; x < width; ++x) {
+            fprintf(blue_noise_image, "%d ", pbp[internal::twoToOne(x, y, width)] ? 1 : 0);
+        }
+        fputc('\n', blue_noise_image);
+    }
+    fclose(blue_noise_image);
+//#endif
+
+    return pbp;
 }
diff --git a/src/blue_noise.hpp b/src/blue_noise.hpp
index 7f587d7..acaba0b 100644
--- a/src/blue_noise.hpp
+++ b/src/blue_noise.hpp
@@ -4,10 +4,12 @@
 #include <vector>
 #include <tuple>
 #include <cmath>
+#include <functional>
+#include <unordered_set>
 
 namespace dither {
 
-std::vector<bool> blue_noise(std::size_t width, std::size_t height);
+std::vector<bool> blue_noise(std::size_t width, std::size_t height, std::size_t threads = 1);
 
 namespace internal {
     inline std::size_t twoToOne(std::size_t x, std::size_t y, std::size_t width) {
@@ -21,22 +23,22 @@ namespace internal {
     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));
+        return std::exp(-(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) {
+            std::size_t width, std::size_t height) {
         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) {
+        // q' = (N + y - (q - M/2)) % N = (N + M/2 + y - q) % N
+        for(std::size_t q = 0; q < width; ++q) {
+            std::size_t q_prime = (height + width / 2 + y - q) % height;
+            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) {
@@ -47,6 +49,61 @@ namespace internal {
 
         return sum;
     }
+
+    inline std::tuple<std::size_t, std::size_t> filter_minmax(const std::vector<double>& filter) {
+        double min = std::numeric_limits<double>::infinity();
+        double max = 0.0;
+        std::size_t min_index = 0;
+        std::size_t max_index = 0;
+
+        for(std::size_t i = 0; i < filter.size(); ++i) {
+            if(filter[i] < min) {
+                min_index = i;
+                min = filter[i];
+            }
+            if(filter[i] > max) {
+                max_index = i;
+                max = filter[i];
+            }
+        }
+
+        return {min_index, max_index};
+    }
+
+    void recursive_apply_radius(
+        std::size_t idx, std::size_t width,
+        std::size_t height, std::size_t radius,
+        const std::function<bool(std::size_t)>& fn);
+
+    bool recursive_apply_radius_impl(
+        std::size_t idx, std::size_t width,
+        std::size_t height, std::size_t radius,
+        const std::function<bool(std::size_t)>& fn,
+        std::unordered_set<std::size_t>& visited);
+
+    inline std::size_t get_one_or_zero(
+            const std::vector<bool>& pbp, bool get_one,
+            std::size_t idx, std::size_t width, std::size_t height) {
+        std::size_t found_idx;
+        bool found = false;
+        for(std::size_t radius = 1; radius <= 12; ++radius) {
+            recursive_apply_radius(
+                idx, width, height, radius,
+                [&found_idx, &found, &pbp, &get_one] (std::size_t idx) {
+                    if((get_one && pbp[idx]) || (!get_one && !pbp[idx])) {
+                        found_idx = idx;
+                        found = true;
+                        return true;
+                    } else {
+                        return false;
+                    }
+                });
+            if(found) {
+                return found_idx;
+            }
+        }
+        return idx;
+    }
 } // namespace dither::internal
 
 } // namespace dither
diff --git a/src/main.cpp b/src/main.cpp
index 11b7fad..d86e9c7 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,3 +1,12 @@
+#include "blue_noise.hpp"
+
+#include <iostream>
+
 int main(int argc, char **argv) {
+//#ifndef NDEBUG
+    std::cout << "Trying blue_noise..." << std::endl;
+    dither::blue_noise(70, 70, 8);
+//#endif
+
     return 0;
 }