Convert an image into one made out of coloured circles with mean pixel values

circles.cpp

#include <algorithm>
#include <cmath>
#include <iostream>
#include <random>
#include <string>
#include <vector>

// RGB values
// Add arithmetic operations to take means squared deviations
struct rgb {
  long long r, g, b;
  rgb operator+(const rgb& p) const { return {r + p.r, g + p.g, b + p.b}; }
  rgb operator-(const rgb& p) const { return {r - p.r, g - p.g, b - p.b}; }
  rgb operator*(const rgb& p) const { return {r * p.r, g * p.g, b * p.b}; }
  rgb operator/(int n) const { return {r / n, g / n, b / n}; }
  rgb operator*(int n) const { return {r * n, g * n, b * n}; }
  long long norm() { return r * r + g * g + b * b; }
};

// PPM image
struct image {
  std::string header;
  int W, H, C;
  std::vector<std::vector<rgb>> pixel;
};

// Read PPM image
std::istream& operator>>(std::istream& in, image& img) {
  in >> img.header;
  in >> img.W >> img.H;
  in >> img.C;
  img.pixel = std::vector<std::vector<rgb>>(img.H, std::vector<rgb>(img.W));
  for (int i = 0; i < img.H; ++i) {
    for (int j = 0; j < img.W; ++j) {
      in >> img.pixel[i][j].r >> img.pixel[i][j].g >> img.pixel[i][j].b;
    }
  }
  return in;
}

// Write PPM image
std::ostream& operator<<(std::ostream& out, const image& img) {
  out << img.header << "\n";
  out << img.W << " " << img.H << "\n";
  out << img.C << "\n";
  for (int i = 0; i < img.H; ++i) {
    for (int j = 0; j < img.W; ++j) {
      out << img.pixel[i][j].r << " " << img.pixel[i][j].g << " "
          << img.pixel[i][j].b << "\n";
    }
  }
  return out;
}

int main() {
  std::ios::sync_with_stdio(false);

  // Read the image
  image img;
  std::cin >> img;

  // Find global mean to colour unfilled pixels with in the end
  rgb globalmean = {0, 0, 0};
  for (auto row : img.pixel) {
    globalmean = std::accumulate(row.begin(), row.end(), globalmean);
  }
  globalmean = globalmean / (img.H * img.W);

  // Make a copy of the image and set all pixels to invalid values
  image conv = img;
  for (auto& i : conv.pixel) {
    for (auto& j : i) {
      j = {-1, -1, -1};
    }
  }

  // Bin (dx,dy) pairs into bins so that increasing radii can be scanned easily
  int lim = std::min(img.H / 2, img.W / 2);
  std::vector<std::vector<std::pair<int, int>>> points(lim + 1);
  for (int i = -lim; i <= lim; ++i) {
    for (int j = -lim; j <= lim; ++j) {
      int r = sqrt(i * i + j * j);
      if (r <= lim) points[r].push_back({i, j});
    }
  }

  // Generate and shuffle the pixels in the image to go through them randomly
  std::vector<std::pair<int, int>> pvec;
  for (int i = 0; i < conv.H; ++i) {
    for (int j = 0; j < conv.W; ++j) {
      pvec.push_back({i, j});
    }
  }
  std::random_device rd;
  std::mt19937 g(rd());
  std::shuffle(pvec.begin(), pvec.end(), g);

  // Circles will be generated in three rounds
  // In each round, from each pixel, move as far as possible till we hit
  // the end of the board / a coloured pixel or the border pixels are too
  // different from the mean
  // Colour the circle with the mean if the radius is larger than threshold
  // In the rounds, the colour difference threshold increases and the radius
  // threshold decreases
  int norms[] = {3000, 10000, 25000};
  int maxrads[] = {20, 10, 5};
  for (int round = 0; round < 3; ++round) {
    for (auto p : pvec) {
      int maxrad, area = 0;
      rgb tot = {0, 0, 0}, mean = tot;
      for (maxrad = 0; maxrad <= lim; ++maxrad) {
        rgb curtot = {0, 0, 0};
        bool flag = true;
        for (auto dp : points[maxrad]) {
          int x = p.first + dp.first, y = p.second + dp.second;
          if (x < 0 or x >= conv.H or y < 0 or y >= conv.W or
              conv.pixel[x][y].r >= 0 or
              (maxrad > 3 and (img.pixel[x][y] - mean).norm() > norms[round])) {
            flag = false;
            break;
          }
          curtot = curtot + img.pixel[x][y];
        }
        if (!flag) break;
        tot = tot + curtot;
        area += points[maxrad].size();
        mean = tot / area;
      }

      if (maxrad <= maxrads[round]) continue;
      for (int i = 0; i < maxrad; ++i) {
        for (auto dp : points[i]) {
          int x = p.first + dp.first, y = p.second + dp.second;
          conv.pixel[x][y] = mean;
        }
      }
    }
  }

  // If a pixel is still uncoloured in the end, colour it with the mean
  // Probably better to colour it with some kind of a running average?
  for (auto& i : conv.pixel) {
    for (auto& j : i) {
      if (j.r == -1) j = globalmean;
    }
  }
  std::cout << conv << std::endl;
  return 0;
}