Initial commit

1 files changed, 246 insertions, 0 deletions

diff --git a/imagealgos.cpp b/imagealgos.cpp
new file mode 100644
index 0000000..23902d1
--- /dev/null
+++ b/imagealgos.cpp
@@ -0,0 +1,246 @@
+#include "imagealgos.h"
+
+#include <cassert>
+#include <cstdint>
+#include <cstring>
+
+#include <algorithm>
+#include <chrono>
+#include <iostream>
+#include <limits>
+#include <mutex>
+#include <typeinfo>
+#include <utility>
+
+// #include <arm_neon.h>
+
+#include "genetic_algos.h"
+#include "macro.h"
+
+uint8_t pgm_image[64 + img_width * img_height * sizeof(uint16_t)];
+size_t pgm_image_size = 0;
+std::mutex pgm_image_mtx;
+
+size_t pgm_save(Image *img, FILE *outfile, bool really_save) {
+    std::lock_guard<std::mutex> lg(pgm_image_mtx);
+
+    size_t n = 0;
+
+    // n += fprintf(outfile, "P5\n%d %d\n%d\n",
+    // img->width, img->height, 0xFF);
+    n += sprintf((char*)pgm_image, "P5\n%d %d\n%d\n",
+                 img->width, img->height, 0xFF);
+
+    for (size_t i = 0; i < img->width * img->height; ++i)
+    {
+        uint16_t *pixels = (uint16_t*)img->data;
+        const auto p = pixels[i];
+        uint8_t value = (pixels[i] & 0xFF00) >> 8;
+
+        // n += fwrite(&value, 1, 1, outfile);
+        memcpy((void*)(pgm_image + n), &value, sizeof(value));
+        n += sizeof(value);
+    }
+
+    pgm_image_size = n;
+
+    // std::cout << "size is " << n << std::endl;
+
+    if (really_save)
+    {
+        fwrite(pgm_image, 1, pgm_image_size, outfile);
+        fflush(outfile);
+    }
+
+    return n;
+}
+
+void unpack_10bit(uint8_t const *src, Image const &image, uint16_t *dest)
+{
+    unsigned int w_align = image.width & ~3;
+    for (unsigned int y = 0; y < image.height; y++, src += image.stride)
+    {
+        uint8_t const *ptr = src;
+        unsigned int x;
+        for (x = 0; x < w_align; x += 4, ptr += 5)
+        {
+            *dest++ = (ptr[0] << 2) | ((ptr[4] >> 0) & 3);
+            *dest++ = (ptr[1] << 2) | ((ptr[4] >> 2) & 3);
+            *dest++ = (ptr[2] << 2) | ((ptr[4] >> 4) & 3);
+            *dest++ = (ptr[3] << 2) | ((ptr[4] >> 6) & 3);
+        }
+        for (; x < image.width; x++)
+            *dest++ = (ptr[x & 3] << 2) | ((ptr[4] >> ((x & 3) << 1)) & 3);
+    }
+}
+
+void unpack_16bit(uint8_t const *src, Image const &image, uint16_t *dest)
+{
+    start_timer(unpack_16bit);
+    /* Assume the pixels in memory are already in native byte order */
+    unsigned int w = image.width;
+
+    for (unsigned int y = 0; y < image.height; y++)
+    {
+        memcpy(dest, src, 2 * w);
+        dest += w;
+        src += image.stride;
+    }
+    stop_timer(unpack_16bit);
+}
+
+void rotate(Image &image)
+{
+    start_timer(rotate);
+
+    using namespace std;
+
+    for (size_t i = 0; i < image.height; ++i)
+    {
+        for (size_t j = 0; j < image.width; ++j)
+        {
+            image.rotated_cw[j][i] = image.data[image.height - i][j];
+        }
+    }
+
+    stop_timer(rotate);
+}
+
+template<class T, size_t N>
+constexpr size_t mysize(T (&)[N]) { return N; }
+
+float process_column(uint16_t (&column)[])
+{
+    float result = std::numeric_limits<float>::quiet_NaN();
+
+    constexpr uint32_t patternSize = 16; // good
+    constexpr uint32_t signalThreshold = 450; // = SKO * sqrt(patternSize)
+    static constexpr uint32_t patternOffset = patternSize - ((patternSize % 2 == 1) ? 1 : 0);
+    const uint32_t correlationSize = img_height - patternSize +
+                                     ((patternSize % 2 == 1) ? 1 : 0);
+    uint32_t correlation[img_height];
+    uint32_t integralSum[img_height];
+    uint32_t maxSum = signalThreshold * 50;
+    uint32_t x1 = 0;
+    int32_t y1 = 0;
+    int32_t y2 = 0;
+
+    memset(correlation, 0, img_height * sizeof(uint32_t));
+    integralSum[0] = 0;
+
+    for(uint32_t i = 1; i < img_height; ++i) {
+        if (column[i] < 100) {
+            column[i] = 1;
+        }
+        integralSum[i] = column[i] / 256 + integralSum[i - 1];
+    }
+
+    // maxSum = 0 ;
+    // size_t maxIdx { 0 };
+
+    // for (size_t i = 0; i < img_height - patternSize; ++i) {
+    //     const auto sum = integralSum[i + patternSize] - integralSum[i];
+    //     if (sum > maxSum) {
+    //         maxSum = sum;
+    //         maxIdx = i;
+    //     }
+    // }
+
+    // Algo genetic(column + maxIdx);
+    // // std::cout << "maxIdx " << maxIdx << std::endl;
+
+    // // return maxIdx + genetic.run().a;
+    // return 500;
+    // return img_height - maxIdx - genetic.run().a;
+
+
+    for(uint32_t i = 0; i < correlationSize; ++i)
+        correlation[i + patternSize / 2] =
+            column[i + patternSize / 2] / 256 *
+            (integralSum[i + patternOffset] - integralSum[i]);
+
+    for(uint32_t i = 3; i < img_height - 2; ++i)
+    {
+        const auto sum = correlation[i - 1] +
+                         correlation[i] +
+                         correlation[i + 1];
+
+        if(sum > maxSum)
+        {
+            const int32_t rioux0 = int32_t(correlation[i - 2 - 1] + correlation[i - 1 - 1]) -
+                                   int32_t(correlation[i + 1 - 1] + correlation[i + 2 - 1]);
+
+            if(rioux0 < 0)
+            {
+                const int32_t rioux1 = int32_t(correlation[i - 2] + correlation[i - 1]) -
+                                       int32_t(correlation[i + 1] + correlation[i + 2]);
+
+                if(rioux1 >= 0)
+                {
+                    x1 = i - 1;
+                    y1 = rioux0;
+                    y2 = rioux1;
+                    maxSum = sum;
+                }
+            }
+        }
+    }
+
+    result = (y2 != y1) ?
+                 (img_height - (float(x1) - (float(y1) / (y2 - y1)))) :
+                 std::numeric_limits<float>::quiet_NaN();
+
+
+    static bool result_done = false;
+    if (!result_done) {
+        std::cout << "result " << result << std::endl;
+        result_done = true;
+    }
+    // std::cout << "result is '" << result << "'\n";
+
+    return result;
+
+// center of mass
+#if 0
+    auto max_el = std::max_element(std::begin(accumulated_sum),
+                                   std::end(accumulated_sum) - window_size);
+
+    size_t max_sum_idx = max_el - std::begin(accumulated_sum) + window_size;
+
+    double sum_w = 0;
+    double prod_wx = 0;
+    double wmc = 0;
+
+    for(int i = max_sum_idx - window_size; i < max_sum_idx; ++i)
+    {
+        prod_wx += column[i] * i;
+        sum_w += column[i];
+    }
+
+    wmc = float(prod_wx) / float(sum_w);
+
+    result = img_height - wmc;
+
+    return result;
+#endif
+}
+
+void process_columns(Image &image)
+{
+    std::cout << "here\n";
+    start_timer(process_columns);
+
+    for (size_t i = 0; i < image.width; i++)
+    {
+        image.pixels[i] = process_column(image.rotated_cw[i]);
+        // Algo genetic(image.rotated_cw[i]);
+
+        // image.pixels[i] = genetic.run().a;
+
+        // if (i == 0) {
+            // std::cout << "pixel: " << image.pixels[i] << std::endl;
+          // }
+    }
+
+    stop_timer(process_columns);
+}