WIP work on generating blue-noise

CMakeLists.txt

@@ -18,12 +18,15 @@ endif()
 
 find_package(Threads REQUIRED)
 find_package(OpenCL REQUIRED)
+find_package(PNG REQUIRED)
 
 add_executable(Dithering ${Dithering_SOURCES})
 target_compile_features(Dithering PUBLIC cxx_std_17)
 target_include_directories(Dithering PUBLIC
     Threads::Threads
-    ${OpenCL_INCLUDE_DIRS})
+    ${OpenCL_INCLUDE_DIRS}
+    ${PNG_INCLUDE_DIRS})
 target_link_libraries(Dithering PUBLIC 
     Threads::Threads
-    ${OpenCL_LIBRARIES})
+    ${OpenCL_LIBRARIES}
+    ${PNG_LIBRARIES})

src/blue_noise.cl

@@ -1,6 +1,23 @@
+int twoToOne(x, y, width, height) {
+    while(x < 0) {
+        x += width;
+    }
+    while(y < 0) {
+        y += height;
+    }
+    x = x % width;
+    y = y % height;
+    return x + y * width;
+}
+
+//float gaussian(float x, float y) {
+//    return exp(-(x*x + y*y) / (1.5F * 1.5F * 2.0F));
+//}
+
 __kernel void do_filter(
-        __global float *filter_out, __global float *precomputed,
+        __global float *filter_out, __global const float *precomputed,
-        __global int *pbp, int width, int height, int filter_size) {
+        __global const int *pbp, const int width, const int height,
+        const int filter_size) {
     int i = get_global_id(0);
     if(i < 0 || i >= width * height) {
         return;
@@ -10,11 +27,12 @@ __kernel void do_filter(
 
     float sum = 0.0f;
     for(int q = 0; q < filter_size; ++q) {
-        int q_prime = (height + filter_size / 2 + y - q) % height;
+        int q_prime = height - filter_size / 2 + y + q;
         for(int p = 0; p < filter_size; ++p) {
-            int p_prime = (width + filter_size / 2 + x - p) % width;
+            int p_prime = width - filter_size / 2 + x + p;
-            if(pbp[p_prime + q_prime * width] != 0) {
+            if(pbp[twoToOne(p_prime, q_prime, width, height)] != 0) {
-                sum += precomputed[p + q * filter_size];
+                sum += precomputed[twoToOne(p, q, filter_size, filter_size)];
+                //sum += gaussian(p - filter_size / 2.0F + 0.5F, q - filter_size / 2.0F + 0.5F);
             }
         }
     }

src/blue_noise.cpp

@@ -6,6 +6,7 @@
 #include <fstream>
 #include <memory>
 #include <string>
+#include <unordered_set>
 
 #include <CL/opencl.h>
 
@@ -81,14 +82,14 @@ image::Bl dither::blue_noise(int width, int height, int threads, bool use_opencl
             }
 
             std::cout << "OpenCL: Initialized, trying cl_impl..." << std::endl;
-            std::vector<bool> result = internal::blue_noise_cl_impl(
+            std::vector<unsigned int> result = internal::blue_noise_cl_impl(
                 width, height, filter_size, context, device, program);
 
             clReleaseProgram(program);
             clReleaseContext(context);
 
             if(!result.empty()) {
-                return internal::toBl(result, width);
+                return internal::rangeToBl(result, width);
             }
         } while (false);
     }
@@ -96,160 +97,13 @@ image::Bl dither::blue_noise(int width, int height, int threads, bool use_opencl
     if(!using_opencl) {
         std::cout << "OpenCL: Failed to setup/use or is not enabled, using regular impl..."
             << std::endl;
-        return internal::toBl(internal::blue_noise_impl(width, height, threads), width);
+        return internal::rangeToBl(internal::blue_noise_impl(width, height, threads), width);
     }
 
     return {};
 }
 
-image::Bl dither::blue_noise_grayscale(int width, int height, int threads) {
+std::vector<unsigned int> dither::internal::blue_noise_impl(int width, int height, int threads) {
-    int count = width * height;
-    std::vector<float> filter_out;
-    filter_out.resize(count);
-
-    std::vector<float> image = internal::random_noise_grayscale(count);
-
-    int iterations = 0;
-    int filter_size = (width + height) / 2;
-    std::vector<float> precomputed(internal::precompute_gaussian(filter_size));
-
-    // TODO DEBUG
-    //float pmax = 0.01F;
-    //for(float value : precomputed) {
-    //    if(value > pmax) {
-    //        pmax = value;
-    //    }
-    //}
-    //for(float &value: precomputed) {
-    //    value /= pmax;
-    //}
-    //return image::Bl(precomputed, filter_size);
-
-    int min, max, min2, max2;
-    int prevmax = -1;
-    int prevmax2 = -1;
-    float tempPixel;
-    while(true) {
-        printf("Iteration %d\n", iterations);
-
-        internal::compute_filter_grayscale(image,
-                                           width, height, count,
-                                           filter_size, filter_out,
-                                           &precomputed, 0);
-
-        std::tie(min, max) = internal::filter_minmax(filter_out);
-        //std::tie(std::ignore, max) = internal::filter_minmax_in_range(max,
-        //                                                              width,
-        //                                                              height,
-        //                                                              7,
-        //                                                              filter_out);
-        printf("min == %4d, max == %4d", min, max);
-        tempPixel = image[max];
-        image[max] = 0.0F;
-
-        internal::compute_filter_grayscale(image,
-                                           width, height, count,
-                                           filter_size, filter_out,
-                                           &precomputed, 0);
-
-        std::tie(min2, max2) = internal::filter_minmax(filter_out);
-        //std::tie(min2, std::ignore) = internal::filter_minmax_in_range(min2,
-        //                                                               width,
-        //                                                               height,
-        //                                                               7,
-        //                                                               filter_out);
-        printf(", min2 == %4d, max2 == %4d\n", min2, max2);
-
-        if(min2 != min) {
-            image[max] = tempPixel;
-            break;
-        } else {
-            image[max] = image[min];
-            image[min] = tempPixel;
-        }
-        //if(prevmax == max && prevmax2 == max2) {
-        //    image[max] = tempPixel;
-        //    break;
-        //} else {
-        //    image[max] = image[min2];
-        //    image[min2] = tempPixel;
-        //}
-
-        prevmax = max;
-        prevmax2 = max2;
-
-//#ifndef NDEBUG
-        if(iterations % 20 == 0) {
-            std::string name;
-            name.append("tempGrayscale");
-            if(iterations < 10) {
-                name.append("00");
-            } else if(iterations < 100) {
-                name.append("0");
-            }
-            name.append(std::to_string(iterations));
-            name.append(".pgm");
-            image::Bl(image, width).writeToFile(image::file_type::PGM, true, name);
-
-            name.clear();
-            name.append("tempFilter");
-            if(iterations < 10) {
-                name.append("00");
-            } else if(iterations < 100) {
-                name.append("0");
-            }
-            name.append(std::to_string(iterations));
-            name.append(".pgm");
-
-            internal::compute_filter_grayscale(image,
-                                               width, height, count,
-                                               filter_size, filter_out,
-                                               &precomputed, 0);
-            std::vector<float> normalizedFilter(filter_out);
-            float fmax = -std::numeric_limits<float>::infinity();
-            float fmin = std::numeric_limits<float>::infinity();
-            for(float value : normalizedFilter) {
-                if(value > fmax) {
-                    fmax = value;
-                }
-                if(value < fmin) {
-                    fmin = value;
-                }
-            }
-            fmax -= fmin;
-            for(float &value : normalizedFilter) {
-                value = (value - fmin) / fmax;
-            }
-            image::Bl(normalizedFilter, width).writeToFile(image::file_type::PGM, true, name);
-        }
-//#endif
-        ++iterations;
-    }
-
-    internal::compute_filter_grayscale(image,
-                                       width, height, count,
-                                       filter_size, filter_out,
-                                       &precomputed, 0);
-    std::vector<float> normalizedFilter(filter_out);
-    float fmax = -std::numeric_limits<float>::infinity();
-    float fmin = std::numeric_limits<float>::infinity();
-    for(float value : normalizedFilter) {
-        if(value > fmax) {
-            fmax = value;
-        }
-        if(value < fmin) {
-            fmin = value;
-        }
-    }
-    fmax -= fmin;
-    for(float &value : normalizedFilter) {
-        value = (value - fmin) / fmax;
-    }
-    image::Bl(normalizedFilter, width).writeToFile(image::file_type::PGM, true, "filterOut.pgm");
-    return image::Bl(image, width);
-}
-
-std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int threads) {
     int count = width * height;
     std::vector<float> filter_out;
     filter_out.resize(count);
@@ -301,36 +155,11 @@ std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int t
 //        }
 #endif
 
-        int min, max, min_zero, max_one;
+        int min, max;
-        std::tie(min, max) = internal::filter_minmax(filter_out);
+        std::tie(min, max) = internal::filter_minmax(filter_out, pbp);
-        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;
+        pbp[max] = false;
 
         // get filter values again
         internal::compute_filter(pbp, width, height, count, filter_size,
@@ -338,20 +167,13 @@ std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int t
 
         // get second buffer's min
         int second_min;
-        std::tie(second_min, std::ignore) = internal::filter_minmax(filter_out);
+        std::tie(second_min, std::ignore) = internal::filter_minmax(filter_out, pbp);
-        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(utility::dist(max_one, second_min, width) < 1.5f) {
+        if(second_min == max) {
-            pbp[max_one] = true;
+            pbp[max] = true;
             break;
         } else {
-            pbp[min_zero] = true;
+            pbp[second_min] = true;
         }
 
         if(iterations % 100 == 0) {
@@ -384,10 +206,45 @@ std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int t
     fclose(blue_noise_image);
 //#endif
 
-    return pbp;
+    std::cout << "Generating dither_array...\n";
+    std::vector<unsigned int> dither_array(count);
+    int min, max;
+    {
+        std::vector<bool> pbp_copy(pbp);
+        std::cout << "Ranking minority pixels...\n";
+        for (unsigned int i = pixel_count; i-- > 0;) {
+            std::cout << i << ' ';
+            internal::compute_filter(pbp, width, height, count, filter_size,
+                    filter_out, precomputed.get(), threads);
+            std::tie(std::ignore, max) = internal::filter_minmax(filter_out, pbp);
+            pbp[max] = false;
+            dither_array[max] = i;
+        }
+        pbp = pbp_copy;
+    }
+    std::cout << "\nRanking remainder of first half of pixels...\n";
+    for (unsigned int i = pixel_count; i < (unsigned int)((count + 1) / 2); ++i) {
+        std::cout << i << ' ';
+        internal::compute_filter(pbp, width, height, count, filter_size,
+                filter_out, precomputed.get(), threads);
+        std::tie(min, std::ignore) = internal::filter_minmax(filter_out, pbp);
+        pbp[min] = true;
+        dither_array[min] = i;
+    }
+    std::cout << "\nRanking last half of pixels...\n";
+    for (unsigned int i = (count + 1) / 2; i < (unsigned int)count; ++i) {
+        std::cout << i << ' ';
+        internal::compute_filter(pbp, width, height, count, filter_size,
+                filter_out, precomputed.get(), threads);
+        std::tie(std::ignore, max) = internal::filter_minmax(filter_out, pbp);
+        pbp[max] = true;
+        dither_array[max] = i;
+    }
+
+    return dither_array;
 }
 
-std::vector<bool> dither::internal::blue_noise_cl_impl(
+std::vector<unsigned int> dither::internal::blue_noise_cl_impl(
         int width, int height, int filter_size, cl_context context, cl_device_id device, cl_program program) {
     cl_int err;
     cl_kernel kernel;
@@ -418,18 +275,6 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
         return {};
     }
 
-    /*
-    err = clEnqueueWriteBuffer(queue, d_pbp, CL_TRUE, 0, count * sizeof(int), &pbp_i[0], 0, nullptr, nullptr);
-    if(err != CL_SUCCESS) {
-        std::cerr << "OpenCL: Failed to write to d_pbp buffer\n";
-        clReleaseMemObject(d_pbp);
-        clReleaseMemObject(d_precomputed);
-        clReleaseMemObject(d_filter_out);
-        clReleaseCommandQueue(queue);
-        return {};
-    }
-    */
-
     kernel = clCreateKernel(program, "do_filter", &err);
     if(err != CL_SUCCESS) {
         std::cerr << "OpenCL: Failed to create kernel: ";
@@ -635,29 +480,10 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
             break;
         }
 
-        int min, max, min_zero, max_one;
+        int min, max;
-        std::tie(min, max) = internal::filter_minmax(filter);
+        std::tie(min, max) = internal::filter_minmax(filter, pbp);
-        if(!pbp[max]) {
-            max_one = internal::get_one_or_zero(pbp, true, max, width, height);
-        } else {
-            max_one = max;
-        }
-        if(!pbp[max_one]) {
-            std::cerr << "ERROR: Failed to find pbp[max] one" << std::endl;
-            break;
-        }
 
-        if(pbp[min]) {
+        pbp[max] = false;
-            min_zero = internal::get_one_or_zero(pbp, false, min, width, height);
-        } else {
-            min_zero = min;
-        }
-        if(pbp[min_zero]) {
-            std::cerr << "ERROR: Failed to find pbp[min] zero" << std::endl;
-            break;
-        }
-
-        pbp[max_one] = false;
 
         if(!get_filter()) {
             std::cerr << "OpenCL: Failed to execute do_filter\n";
@@ -666,20 +492,13 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
 
         // get second buffer's min
         int second_min;
-        std::tie(second_min, std::ignore) = internal::filter_minmax(filter);
+        std::tie(second_min, std::ignore) = internal::filter_minmax(filter, pbp);
-        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(utility::dist(max_one, second_min, width) < 1.5f) {
+        if(second_min == max) {
-            pbp[max_one] = true;
+            pbp[max] = true;
             break;
         } else {
-            pbp[min_zero] = true;
+            pbp[second_min] = true;
         }
 
         if(iterations % 100 == 0) {
@@ -711,10 +530,83 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
         fclose(blue_noise_image);
     }
 
+#ifndef NDEBUG
+    {
+        image::Bl pbp_image = toBl(pbp, width);
+        pbp_image.writeToFile(image::file_type::PNG, true, "debug_pbp_before.png");
+    }
+#endif
+
+    std::cout << "Generating dither_array...\n";
+    std::unordered_set<unsigned int> set;
+    std::vector<unsigned int> dither_array(count);
+    int min, max;
+    {
+        std::vector<bool> pbp_copy(pbp);
+        std::cout << "Ranking minority pixels...\n";
+        for (unsigned int i = pixel_count; i-- > 0;) {
+            std::cout << i << ' ';
+            get_filter();
+            std::tie(std::ignore, max) = internal::filter_minmax(filter, pbp);
+            pbp.at(max) = false;
+            dither_array.at(max) = i;
+            if (set.find(max) != set.end()) {
+                std::cout << "\nWARNING: Reusing index " << max << '\n';
+            } else {
+                set.insert(max);
+            }
+        }
+        pbp = pbp_copy;
+    }
+    std::cout << "\nRanking remainder of first half of pixels...\n";
+    for (unsigned int i = pixel_count; i < (unsigned int)((count + 1) / 2); ++i) {
+        std::cout << i << ' ';
+        get_filter();
+        std::tie(min, std::ignore) = internal::filter_minmax(filter, pbp);
+        pbp.at(min) = true;
+        dither_array.at(min) = i;
+        if (set.find(min) != set.end()) {
+            std::cout << "\nWARNING: Reusing index " << min << '\n';
+        } else {
+            set.insert(min);
+        }
+    }
+#ifndef NDEBUG
+    {
+        get_filter();
+        internal::write_filter(filter, width, "filter_mid.pgm");
+        image::Bl pbp_image = toBl(pbp, width);
+        pbp_image.writeToFile(image::file_type::PNG, true, "debug_pbp_mid.png");
+    }
+#endif
+    std::cout << "\nRanking last half of pixels...\n";
+    for (unsigned int i = (count + 1) / 2; i < (unsigned int)count; ++i) {
+        std::cout << i << ' ';
+        get_filter();
+        std::tie(std::ignore, max) = internal::filter_minmax(filter, pbp);
+        pbp.at(max) = true;
+        dither_array.at(max) = i;
+        if (set.find(max) != set.end()) {
+            std::cout << "\nWARNING: Reusing index " << max << '\n';
+        } else {
+            set.insert(max);
+        }
+    }
+    std::cout << std::endl;
+
+#ifndef NDEBUG
+    {
+        get_filter();
+        internal::write_filter(filter, width, "filter_after.pgm");
+        image::Bl pbp_image = toBl(pbp, width);
+        pbp_image.writeToFile(image::file_type::PNG, true, "debug_pbp_after.png");
+    }
+#endif
+
     clReleaseKernel(kernel);
     clReleaseMemObject(d_pbp);
     clReleaseMemObject(d_precomputed);
     clReleaseMemObject(d_filter_out);
     clReleaseCommandQueue(queue);
-    return pbp;
+    return dither_array;
 }

src/blue_noise.hpp

@@ -1,6 +1,7 @@
 #ifndef BLUE_NOISE_HPP
 #define BLUE_NOISE_HPP
 
+#include <limits>
 #include <vector>
 #include <functional>
 #include <unordered_set>
@@ -13,6 +14,7 @@
 #include <random>
 #include <cassert>
 #include <stdexcept>
+#include <iostream>
 
 #include <sys/sysinfo.h>
 
@@ -25,11 +27,9 @@ namespace dither {
 
 image::Bl blue_noise(int width, int height, int threads = 1, bool use_opencl = true);
 
-image::Bl blue_noise_grayscale(int width, int height, int threads = 1);
-
 namespace internal {
-    std::vector<bool> blue_noise_impl(int width, int height, int threads = 1);
+    std::vector<unsigned int> blue_noise_impl(int width, int height, int threads = 1);
-    std::vector<bool> blue_noise_cl_impl(
+    std::vector<unsigned int> blue_noise_cl_impl(
         int width, int height, int filter_size,
         cl_context context, cl_device_id device, cl_program program);
 
@@ -59,27 +59,6 @@ namespace internal {
         return pbp;
     }
 
-    inline std::vector<float> random_noise_grayscale(unsigned int size) {
-        std::vector<float> graynoise;
-        graynoise.reserve(size);
-        std::default_random_engine re(std::random_device{}());
-        std::uniform_real_distribution<float> dist(0.0F, 1.0F);
-
-        for(unsigned int i = 0; i < size; ++i) {
-            graynoise.push_back(static_cast<float>(i) / static_cast<float>(size - 1));
-            //graynoise.push_back(dist(re));
-        }
-        for(unsigned int i = 0; i < size - 1; ++i) {
-            std::uniform_int_distribution<unsigned int> range(i + 1, size - 1);
-            unsigned int ridx = range(re);
-            float temp = graynoise[i];
-            graynoise[i] = graynoise[ridx];
-            graynoise[ridx] = temp;
-        }
-
-        return graynoise;
-    }
-
     constexpr float mu = 1.5F;
     constexpr float mu_squared = mu * mu;
     constexpr float double_mu_squared = 2.0F * mu * mu;
@@ -95,8 +74,8 @@ namespace internal {
         for(int i = 0; i < size * size; ++i) {
             auto xy = utility::oneToTwo(i, size);
             precomputed.push_back(gaussian(
-                (float)xy.first - (float)size / 2.0f,
+                (float)xy.first - (float)size / 2.0F + 0.5F,
-                (float)xy.second - (float)size / 2.0f));
+                (float)xy.second - (float)size / 2.0F + 0.5F));
         }
 
         return precomputed;
@@ -113,12 +92,12 @@ namespace internal {
         // p' = (M + x - (p - M/2)) % M = (3M/2 + x - p) % M
         // q' = (N + y - (q - M/2)) % N = (N + M/2 + y - q) % N
         for(int q = 0; q < filter_size; ++q) {
-            int q_prime = (height + filter_size / 2 + y - q) % height;
+            int q_prime = (height - filter_size / 2 + y + q) % height;
             for(int p = 0; p < filter_size; ++p) {
-                int p_prime = (width + filter_size / 2 + x - p) % width;
+                int p_prime = (width - filter_size / 2 + x + p) % width;
                 if(pbp[utility::twoToOne(p_prime, q_prime, width, height)]) {
-                    sum += gaussian((float)p - filter_size/2.0f,
+                    sum += gaussian((float)p - filter_size/2.0F + 0.5F,
-                                    (float)q - filter_size/2.0f);
+                                    (float)q - filter_size/2.0F + 0.5F);
                 }
             }
         }
@@ -126,24 +105,6 @@ namespace internal {
         return sum;
     }
 
-    inline float filter_grayscale(
-            const std::vector<float> &image,
-            int x, int y,
-            int width, int height, int filter_size) {
-        float sum = 0.0F;
-        for(int q = 0; q < filter_size; ++q) {
-            int q_prime = (height + filter_size / 2 + y - q) % height;
-            for(int p = 0; p < filter_size; ++p) {
-                int p_prime = (width + filter_size / 2 + x - p) % width;
-                sum += image[utility::twoToOne(p_prime, q_prime, width, height)]
-                        * gaussian((float)p - filter_size/2.0F,
-                                   (float)q - filter_size/2.0F);
-            }
-        }
-
-        return sum;
-    }
-
     inline float filter_with_precomputed(
             const std::vector<bool>& pbp,
             int x, int y,
@@ -151,35 +112,13 @@ namespace internal {
             const std::vector<float> &precomputed) {
         float sum = 0.0f;
 
-        for(int q = 0; q < filter_size; ++q) {
-            int q_prime = (height + filter_size / 2 + y - q) % height;
-            for(int p = 0; p < filter_size; ++p) {
-                int p_prime = (width + filter_size / 2 + x - p) % width;
-                if(pbp[utility::twoToOne(p_prime, q_prime, width, height)]) {
-                    sum += precomputed[utility::twoToOne(p, q, filter_size, filter_size)];
-                }
-            }
-        }
-
-        return sum;
-    }
-
-    inline float filter_with_precomputed_grayscale(
-            const std::vector<float>& image,
-            int x, int y,
-            int width, int height, int filter_size,
-            const std::vector<float> &precomputed) {
-        float sum = 0.0F;
-
         for(int q = 0; q < filter_size; ++q) {
             int q_prime = (height - filter_size / 2 + y + q) % height;
             for(int p = 0; p < filter_size; ++p) {
                 int p_prime = (width - filter_size / 2 + x + p) % width;
-                sum += image[utility::twoToOne(p_prime, q_prime, width, height)]
+                if(pbp[utility::twoToOne(p_prime, q_prime, width, height)]) {
-                        * precomputed[utility::twoToOne(p,
+                    sum += precomputed[utility::twoToOne(p, q, filter_size, filter_size)];
-                                                        q,
+                }
-                                                        filter_size,
-                                                        filter_size)];
             }
         }
 
@@ -281,95 +220,61 @@ namespace internal {
 
     }
 
-    inline void compute_filter_grayscale(
+    inline std::pair<int, int> filter_minmax(const std::vector<float> &filter,
-            const std::vector<float> &image, int width, int height,
+                                             std::vector<bool> pbp) {
-            int count, int filter_size, std::vector<float> &filter_out,
+        // ensure minority pixel is "true"
-            const std::vector<float> *precomputed = nullptr,
+        unsigned int count = 0;
-            int threads = 1) {
+        for (bool value : pbp) {
-        if(threads == 1) {
+            if(value) {
-            if(precomputed) {
+                ++count;
-                for(int y = 0; y < height; ++y) {
-                    for(int x = 0; x < width; ++x) {
-                        filter_out[utility::twoToOne(x, y, width, height)] =
-                            internal::filter_with_precomputed_grayscale(
-                                image,
-                                x, y,
-                                width, height,
-                                filter_size,
-                                *precomputed);
             }
         }
-            } else {
+        if (count * 2 >= pbp.size()) {
-                for(int y = 0; y < height; ++y) {
+            for (unsigned int i = 0; i < pbp.size(); ++i) {
-                    for(int x = 0; x < width; ++x) {
+                pbp[i] = !pbp[i];
-                        filter_out[utility::twoToOne(x, y, width, height)] =
-                            internal::filter_grayscale(image,
-                                                       x, y,
-                                                       width, height,
-                                                       filter_size);
-                    }
-                }
-            }
-        } else {
-            if(threads == 0) {
-                threads = get_nprocs();
-                if(threads == 0) {
-                    throw std::runtime_error("0 threads detected, "
-                            "should be impossible");
             }
         }
 
-            if(precomputed) {
-                const auto tfn = [] (unsigned int ymin, unsigned int ymax,
-                                     unsigned int width, unsigned int height,
-                                     unsigned int filter_size,
-                                     const std::vector<float> *const image,
-                                     std::vector<float> *const filter_out,
-                                     const std::vector<float> *const precomputed) {
-                    for(unsigned int y = ymin; y < ymax; ++y) {
-                        for(unsigned int x = 0; x < width; ++x) {
-                            (*filter_out)[utility::twoToOne(x, y, width, height)] =
-                                internal::filter_with_precomputed_grayscale(
-                                    *image,
-                                    x, y,
-                                    width, height,
-                                    filter_size,
-                                    *precomputed);
-                        }
-                    }
-                };
-                unsigned int step = height / threads;
-                std::vector<std::thread> threadHandles;
-                for(int i = 0; i < threads; ++i) {
-                    unsigned int starty = i * step;
-                    unsigned int endy = (i + 1) * step;
-                    if(i + 1 == threads) {
-                        endy = height;
-                    }
-                    threadHandles.emplace_back(tfn, starty, endy,
-                                               width, height,
-                                               filter_size,
-                                               &image,
-                                               &filter_out,
-                                               precomputed);
-                }
-                for(int i = 0; i < threads; ++i) {
-                    threadHandles[i].join();
-                }
-            } else {
-                // TODO unimplemented
-                throw std::runtime_error("Unimplemented");
-            }
-        }
-    }
-
-    inline std::pair<int, int> filter_minmax(const std::vector<float>& filter) {
         float min = std::numeric_limits<float>::infinity();
         float max = -std::numeric_limits<float>::infinity();
-        int min_index = 0;
+        int min_index = -1;
-        int max_index = 0;
+        int max_index = -1;
 
         for(std::vector<float>::size_type i = 0; i < filter.size(); ++i) {
+            if(!pbp[i] && filter[i] < min) {
+                min_index = i;
+                min = filter[i];
+            }
+            if(pbp[i] && filter[i] > max) {
+                max_index = i;
+                max = filter[i];
+            }
+        }
+
+        return {min_index, max_index};
+    }
+
+    inline std::pair<int, int> filter_abs_minmax(
+            const std::vector<float> &filter) {
+        float min = std::numeric_limits<float>::infinity();
+        float max = -std::numeric_limits<float>::infinity();
+        int min_index = -1;
+        int max_index = -1;
+
+        std::default_random_engine re(std::random_device{}());
+        std::size_t startIdx = std::uniform_int_distribution<std::size_t>(0, filter.size() - 1)(re);
+
+        for(std::vector<float>::size_type i = startIdx; i < filter.size(); ++i) {
+            if(filter[i] < min) {
+                min_index = i;
+                min = filter[i];
+            }
+            if(filter[i] > max) {
+                max_index = i;
+                max = filter[i];
+            }
+        }
+        for(std::vector<float>::size_type i = 0; i < startIdx; ++i) {
             if(filter[i] < min) {
                 min_index = i;
                 min = filter[i];
@@ -449,7 +354,7 @@ namespace internal {
 
     inline void write_filter(const std::vector<float> &filter, int width, const char *filename) {
         int min, max;
-        std::tie(min, max) = filter_minmax(filter);
+        std::tie(min, max) = filter_abs_minmax(filter);
 
         printf("Writing to %s, min is %.3f, max is %.3f\n", filename, filter[min], filter[max]);
         FILE *filter_image = fopen(filename, "w");
@@ -472,12 +377,40 @@ namespace internal {
                 && "New image::Bl size too small (pbp's size is not a multiple of width)");
 
         for(unsigned int i = 0; i < pbp.size(); ++i) {
-            bwImage.getData()[i] = pbp[i] ? 1 : 0;
+            bwImage.getData()[i] = pbp[i] ? 255 : 0;
         }
 
         return bwImage;
     }
 
+    inline image::Bl rangeToBl(const std::vector<unsigned int> &values, int width) {
+        int min = std::numeric_limits<int>::max();
+        int max = std::numeric_limits<int>::min();
+
+        for (int value : values) {
+            if (value < min) {
+                min = value;
+            }
+            if (value > max) {
+                max = value;
+            }
+        }
+
+        std::cout << "rangeToBl: Got min == " << min << " and max == " << max << std::endl;
+
+        max -= min;
+
+        image::Bl grImage(width, values.size() / width);
+        assert((unsigned long)grImage.getSize() >= values.size()
+                && "New image::Bl size too small (values' size is not a multiple of width)");
+
+        for(unsigned int i = 0; i < values.size(); ++i) {
+            grImage.getData()[i] = ((float)((int)(values[i]) - min) / (float)max) * 255.0F;
+        }
+
+        return grImage;
+    }
+
     inline std::pair<int, int> filter_minmax_in_range(int start, int width,
                                                    int height,
                                                    int range,

src/image.cpp

@@ -2,6 +2,9 @@
 
 #include <cstdio>
 #include <random>
+#include <iostream>
+
+#include <png.h>
 
 image::Bl::Bl() :
 data(),
@@ -83,6 +86,7 @@ bool image::Bl::canWriteFile(file_type type) {
     case file_type::PBM:
     case file_type::PGM:
     case file_type::PPM:
+    case file_type::PNG:
         return true;
     default:
         return false;
@@ -104,6 +108,66 @@ bool image::Bl::writeToFile(file_type type, bool canOverwrite, const char *filen
         fclose(file);
     }
 
+    if(type == file_type::PNG) {
+        FILE *outfile = fopen(filename, "wb");
+        if (outfile == nullptr) {
+            return false;
+        }
+        const static auto pngErrorLFn = [] (png_structp /* unused */,
+                                            png_const_charp message) {
+            std::cerr << "WARNING [libpng]: " << message << std::endl;
+        };
+        const static auto pngWarnLFn = [] (png_structp /* unused */,
+                                           png_const_charp message) {
+            std::cerr << "ERROR [libpng]: " << message << std::endl;
+        };
+
+        png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
+                                                      nullptr,
+                                                      pngErrorLFn,
+                                                      pngWarnLFn);
+
+        if (png_ptr == nullptr) {
+            fclose(outfile);
+            return false;
+        }
+
+        png_infop info_ptr = png_create_info_struct(png_ptr);
+        if (info_ptr == nullptr) {
+            png_destroy_write_struct(&png_ptr, nullptr);
+            fclose(outfile);
+            return false;
+        }
+
+        if (setjmp(png_jmpbuf(png_ptr))) {
+            png_destroy_write_struct(&png_ptr, &info_ptr);
+            fclose(outfile);
+            return false;
+        }
+
+        png_init_io(png_ptr, outfile);
+
+        png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_GRAY,
+                     PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT,
+                     PNG_FILTER_TYPE_DEFAULT);
+
+        png_write_info(png_ptr, info_ptr);
+
+        //png_set_filler(png_ptr, 0, PNG_FILLER_AFTER);
+
+        for (unsigned int j = 0; j < this->data.size() / this->width; ++j) {
+            unsigned char *dataPtr = &this->data.at(j * this->width);
+            png_write_rows(png_ptr, &dataPtr, 1);
+        }
+
+        png_write_end(png_ptr, nullptr);
+
+        png_destroy_write_struct(&png_ptr, &info_ptr);
+
+        fclose(outfile);
+        return true;
+    }
+
     switch(type) {
     case file_type::PBM:
         file = fopen(filename, "w");

src/image.hpp

@@ -18,7 +18,7 @@ namespace image {
         PBM,
         PGM,
         PPM,
-        // TODO PNG support
+        PNG,
     };
 
     class Base {

src/main.cpp

@@ -4,13 +4,10 @@
 
 int main(int argc, char **argv) {
 //#ifndef NDEBUG
-//    std::cout << "Trying blue_noise..." << std::endl;
+    std::cout << "Trying blue_noise..." << std::endl;
-//    image::Bl bl = dither::blue_noise(100, 100, 8, true);
+    image::Bl bl = dither::blue_noise(32, 32, 15, true);
-//    bl.writeToFile(image::file_type::PBM, true, "blueNoiseOut.pbm");
+    bl.writeToFile(image::file_type::PNG, true, "blueNoiseOut.png");
 //#endif
 
-    image::Bl bl = dither::blue_noise_grayscale(64, 64);
-    bl.writeToFile(image::file_type::PGM, true, "blueNoiseGrayscaleOut.pgm");
-
     return 0;
 }