Some work on grayscale blue-noise

.gitignore

@@ -2,3 +2,5 @@ src/*.o
 Dithering
 
 compile_commands.json
+
+.cache/

src/blue_noise.cpp

@@ -5,6 +5,7 @@
 #include <iostream>
 #include <fstream>
 #include <memory>
+#include <string>
 
 #include <CL/opencl.h>
 
@@ -101,6 +102,64 @@ image::Bl dither::blue_noise(int width, int height, int threads, bool use_opencl
     return {};
 }
 
+image::Bl dither::blue_noise_grayscale(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) / 4;
+    std::vector<float> precomputed(internal::precompute_gaussian(filter_size));
+
+    while(true) {
+        printf("Iteration %d\n", iterations);
+
+        internal::compute_filter_grayscale(image,
+                                           width, height, count,
+                                           filter_size, filter_out,
+                                           &precomputed, threads);
+
+        int min, max;
+        float tempPixel;
+        int prevmin = -1;
+        int prevmax = -1;
+        std::tie(min, max) = internal::filter_minmax(filter_out);
+        printf("min == %4d, max == %4d\n", min, max);
+        tempPixel = image[max];
+        image[max] = image[min];
+        image[min] = tempPixel;
+        if(prevmin >= 0 && prevmax >= 0
+                && (utility::dist(min, prevmin, width) < 1.5F
+                    || utility::dist(max, prevmax, width) < 1.5F)) {
+            break;
+        }
+        prevmin = min;
+        prevmax = max;
+
+//#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);
+        }
+//#endif
+        ++iterations;
+    }
+
+    // TODO
+
+    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;
@@ -117,7 +176,7 @@ std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int t
     fprintf(random_noise_image, "P1\n%d %d\n", width, height);
     for(int y = 0; y < height; ++y) {
         for(int x = 0; x < width; ++x) {
-            fprintf(random_noise_image, "%d ", pbp[utility::twoToOne(x, y, width)] ? 1 : 0);
+            fprintf(random_noise_image, "%d ", pbp[utility::twoToOne(x, y, width, height)] ? 1 : 0);
         }
         fputc('\n', random_noise_image);
     }
@@ -212,7 +271,7 @@ std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int t
             fprintf(blue_noise_image, "P1\n%d %d\n", width, height);
             for(int y = 0; y < height; ++y) {
                 for(int x = 0; x < width; ++x) {
-                    fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width)] ? 1 : 0);
+                    fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width, height)] ? 1 : 0);
                 }
                 fputc('\n', blue_noise_image);
             }
@@ -229,7 +288,7 @@ std::vector<bool> dither::internal::blue_noise_impl(int width, int height, int t
     fprintf(blue_noise_image, "P1\n%d %d\n", width, height);
     for(int y = 0; y < height; ++y) {
         for(int x = 0; x < width; ++x) {
-            fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width)] ? 1 : 0);
+            fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width, height)] ? 1 : 0);
         }
         fputc('\n', blue_noise_image);
     }
@@ -458,7 +517,7 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
         fprintf(random_noise_image, "P1\n%d %d\n", width, height);
         for(int y = 0; y < height; ++y) {
             for(int x = 0; x < width; ++x) {
-                fprintf(random_noise_image, "%d ", pbp[utility::twoToOne(x, y, width)] ? 1 : 0);
+                fprintf(random_noise_image, "%d ", pbp[utility::twoToOne(x, y, width, height)] ? 1 : 0);
             }
             fputc('\n', random_noise_image);
         }
@@ -540,7 +599,7 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
             fprintf(blue_noise_image, "P1\n%d %d\n", width, height);
             for(int y = 0; y < height; ++y) {
                 for(int x = 0; x < width; ++x) {
-                    fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width)] ? 1 : 0);
+                    fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width, height)] ? 1 : 0);
                 }
                 fputc('\n', blue_noise_image);
             }
@@ -556,7 +615,7 @@ std::vector<bool> dither::internal::blue_noise_cl_impl(
         fprintf(blue_noise_image, "P1\n%d %d\n", width, height);
         for(int y = 0; y < height; ++y) {
             for(int x = 0; x < width; ++x) {
-                fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width)] ? 1 : 0);
+                fprintf(blue_noise_image, "%d ", pbp[utility::twoToOne(x, y, width, height)] ? 1 : 0);
             }
             fputc('\n', blue_noise_image);
         }

src/blue_noise.hpp

@@ -22,6 +22,8 @@ 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<bool> blue_noise_cl_impl(
@@ -54,6 +56,27 @@ 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[i] = 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_squared = 1.5f * 1.5f;
 
     inline float gaussian(float x, float y) {
@@ -87,8 +110,9 @@ namespace internal {
             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)]) {
+                if(pbp[utility::twoToOne(p_prime, q_prime, width, height)]) {
-                    sum += gaussian((float)p - filter_size/2.0f, (float)q - filter_size/2.0f);
+                    sum += gaussian((float)p - filter_size/2.0f,
+                                    (float)q - filter_size/2.0f);
                 }
             }
         }
@@ -96,6 +120,24 @@ 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,
@@ -107,8 +149,8 @@ namespace internal {
             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)]) {
+                if(pbp[utility::twoToOne(p_prime, q_prime, width, height)]) {
-                    sum += precomputed[utility::twoToOne(p, q, filter_size)];
+                    sum += precomputed[utility::twoToOne(p, q, filter_size, filter_size)];
                 }
             }
         }
@@ -116,6 +158,25 @@ namespace internal {
         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)]
+                        * precomputed[utility::twoToOne(p, q, filter_size, filter_size)];
+            }
+        }
+
+        return sum;
+    }
+
     inline void compute_filter(
             const std::vector<bool> &pbp, int width, int height,
             int count, int filter_size, std::vector<float> &filter_out,
@@ -125,7 +186,7 @@ namespace internal {
             if(precomputed) {
                 for(int y = 0; y < height; ++y) {
                     for(int x = 0; x < width; ++x) {
-                        filter_out[utility::twoToOne(x, y, width)] =
+                        filter_out[utility::twoToOne(x, y, width, height)] =
                             internal::filter_with_precomputed(
                                 pbp, x, y, width, height, filter_size, *precomputed);
                     }
@@ -133,7 +194,7 @@ namespace internal {
             } else {
                 for(int y = 0; y < height; ++y) {
                     for(int x = 0; x < width; ++x) {
-                        filter_out[utility::twoToOne(x, y, width)] =
+                        filter_out[utility::twoToOne(x, y, width, height)] =
                             internal::filter(pbp, x, y, width, height, filter_size);
                     }
                 }
@@ -211,6 +272,36 @@ namespace internal {
 
     }
 
+    inline void compute_filter_grayscale(
+            const std::vector<float> &image, int width, int height,
+            int count, int filter_size, std::vector<float> &filter_out,
+            const std::vector<float> *precomputed = nullptr,
+            int threads = 1) {
+        if(precomputed) {
+            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 {
+            for(int y = 0; y < height; ++y) {
+                for(int x = 0; x < width; ++x) {
+                    filter_out[utility::twoToOne(x, y, width, height)] =
+                        internal::filter_grayscale(image,
+                                                   x, y,
+                                                   width, height,
+                                                   filter_size);
+                }
+            }
+        }
+    }
+
     inline std::pair<int, int> filter_minmax(const std::vector<float>& filter) {
         float min = std::numeric_limits<float>::infinity();
         float max = 0.0f;
@@ -287,7 +378,7 @@ namespace internal {
                     break;
                 }
             }
-            next = utility::twoToOne(xy.first, xy.second, width);
+            next = utility::twoToOne(xy.first, xy.second, width, height);
             if((get_one && pbp[next]) || (!get_one && !pbp[next])) {
                 return next;
             }

src/image.cpp

@@ -27,6 +27,16 @@ width(width),
 height(data.size() / width)
 {}
 
+image::Bl::Bl(const std::vector<float> &data, int width) :
+    data{},
+    width(width),
+    height(data.size() / width)
+{
+    for(float gspixel : data) {
+        this->data.push_back(static_cast<uint8_t>(255.0F * gspixel));
+    }
+}
+
 void image::Bl::randomize() {
     if(!isValid()) {
         return;
@@ -112,7 +122,7 @@ bool image::Bl::writeToFile(file_type type, bool canOverwrite, const char *filen
         return false;
     }
     for(unsigned int i = 0; i < data.size(); ++i) {
-        if(i % width == 0) {
+        if(type == file_type::PBM && i % width == 0) {
             fprintf(file, "\n");
         }
         switch(type) {
@@ -120,10 +130,14 @@ bool image::Bl::writeToFile(file_type type, bool canOverwrite, const char *filen
             fprintf(file, "%d ", data[i] == 0 ? 0 : 1);
             break;
         case file_type::PGM:
-            fprintf(file, "%c ", data[i]);
+            //fprintf(file, "%c ", data[i]);
+            fputc(data[i], file);
             break;
         case file_type::PPM:
-            fprintf(file, "%c %c %c ", data[i], data[i], data[i]);
+            //fprintf(file, "%c %c %c ", data[i], data[i], data[i]);
+            fputc(data[i], file);
+            fputc(data[i], file);
+            fputc(data[i], file);
             break;
         default:
             fclose(file);

src/image.hpp

@@ -53,6 +53,7 @@ namespace image {
         Bl(int width, int height);
         Bl(const std::vector<uint8_t> &data, int width);
         Bl(std::vector<uint8_t> &&data, int width);
+        Bl(const std::vector<float> &data, int width);
         virtual ~Bl() {}
 
         Bl(const Bl &other) = default;

src/main.cpp

@@ -4,10 +4,13 @@
 
 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(100, 100, 8, true);
-    bl.writeToFile(image::file_type::PBM, true, "blueNoiseOut.pbm");
+//    bl.writeToFile(image::file_type::PBM, true, "blueNoiseOut.pbm");
 //#endif
 
+    image::Bl bl = dither::blue_noise_grayscale(64, 64);
+    bl.writeToFile(image::file_type::PGM, true, "blueNoiseGrayscaleOut.pgm");
+
     return 0;
 }

src/utility.hpp

@@ -5,7 +5,9 @@
 #include <cmath>
 
 namespace utility {
-    inline int twoToOne(int x, int y, int width) {
+    inline int twoToOne(int x, int y, int width, int height) {
+        x = x % width;
+        y = y % height;
         return x + y * width;
     }