Experimenting with space-filling curve access order for image convolution.

convolve.cpp

#include <array>
#include <cstdint>
#include <cstdlib>
#include <climits>
using namespace std;

uint32_t interleave(const uint16_t x, const uint16_t y) {
  uint32_t result;
  for (int i = 0; i < sizeof(uint16_t) * CHAR_BIT; ++i)
    result |= (x & 1u << i) << i | (y & 1u << i) << (i + 1);
  return result;
}

template<class T, uint16_t ROWS, uint16_t COLUMNS>
class mat {
private:

  union U {
    uint32_t whole;
    struct {
      uint16_t high;
      uint16_t low;
    } pair;
  };

public:

#ifdef INTERLEAVE

  inline T& operator()(const uint16_t row, const uint16_t column) {
    U index { .whole = interleave(row, column) };
    return data[index.pair.high * COLUMNS + index.pair.low];
  }

  inline T operator()(const uint16_t row, const uint16_t column) const {
    U index { .whole = interleave(row, column) };
    return data[index.pair.high * COLUMNS + index.pair.low];
  }

#else

  inline T& operator()(const uint16_t row, const uint16_t column) {
    return data[row * COLUMNS + column];
  }

  inline T operator()(const uint16_t row, const uint16_t column) const {
    return data[row * COLUMNS + column];
  }

#endif

  array<T, COLUMNS * ROWS> data;

  static constexpr uint16_t rows = ROWS;
  static constexpr uint16_t columns = COLUMNS;

  // auto begin() -> decltype(data.begin()) { return data.begin(); }
  // auto end() -> decltype(data.end()) { return data.end(); }

};

constexpr int size = 100;

int main() {
  const mat<int8_t, 3, 3> kernel {{{
    -1, -1, -1,
     0,  0,  0,
    +1, +1, +1
  }}};
  
  mat<uint8_t, size, size> image;
  srand(1000);

  for (int i = 0; i < 1e5; ++i) {

    for (uint16_t row = 0; row < image.rows; ++row)
      for (uint16_t column = 0; column < image.columns; ++column)
        image(row, column) = uint8_t(rand());

    decltype(image) output;
    for (uint16_t row = 0; row < size; ++row) {
      for (uint16_t column = 0; column < size; ++column) {
        output(row, column)
          = kernel(0, 0) * image(row-1, column-1)
          + kernel(0, 1) * image(row-1, column  )
          + kernel(0, 2) * image(row-1, column+1)
          + kernel(1, 0) * image(row  , column-1)
          + kernel(1, 1) * image(row  , column  )
          + kernel(1, 2) * image(row  , column+1)
          + kernel(2, 0) * image(row+1, column-1)
          + kernel(2, 1) * image(row+1, column  )
          + kernel(2, 2) * image(row+1, column+1);
      }
    }

  }

}