#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<png.h>

#define MAX_GRID_SIZE 50
#define MIN_GRID_SIZE 10
#define REGION_COUNT MAX_GRID_SIZE

png_image image;
png_bytep buffer[2];

int seed_point[REGION_COUNT][2],
    s, i, j, x, y, px, py, step, c_current, c_neighbor, *pixel,
    t, r, w, h,

/* Adjacency matrix between different color regions,
   or grid color assignments.                        */
    grid[MAX_GRID_SIZE][MAX_GRID_SIZE],

/* Color assignments for each region (1-based). */
    color_map[REGION_COUNT];

double transform[6], hx, hy;

double Rnd(double range)
{
   return (double)rand() / (RAND_MAX + 1u) * range;
}

/* Assign colors recursively.  Returns 1 on success. */
int RecursiveAssignColors()
{
   /* Set of colors that are adjacent to each region, one bit for each color. */
   int colors_used[REGION_COUNT],

   /* Lowest number of available colors for most constrained region. */
       lowest_degree_of_freedom = 5,

   /* Selected region index. */
       selected_region = -1,

       u, v, degree_of_freedom;

   /* Find region that don't have a color map entry yet, keeping
      one with the lowest degree of freedom.                     */
   for(memset(colors_used, u = 0, sizeof(colors_used)); u < REGION_COUNT; u++)
      if( !color_map[u] )
      {
         degree_of_freedom = 4;
         for(v = 0; v < REGION_COUNT; v++)
            if( u - v && grid[u][v] && color_map[v] &&
                !(colors_used[u] & (1 << color_map[v])) )
            {
               /* Neighbor used a color that is not used by current region,
                  mark that color as unavailable for use.                   */
               colors_used[u] |= 1 << color_map[v];
               degree_of_freedom--;
               if( degree_of_freedom == 0 )
                  return 0;  /* No colors available for current region. */
            }
         if( lowest_degree_of_freedom > degree_of_freedom )
         {
            lowest_degree_of_freedom = degree_of_freedom;
            selected_region = u;
         }
      }

   /* Return success if there are no uncolored regions left. */
   if( selected_region < 0 )
      return 1;

   /* Try different colors for the selected region. */
   for(u = 1; u <= 4; u++)
      if( !(colors_used[selected_region] & (1 << u)) )
      {
         color_map[selected_region] = u;
         if( RecursiveAssignColors() )
            return 1;

         /* Failed to make progress with this selected color, backtrack. */
         color_map[selected_region] = 0;
      }

   /* None of the color assignments worked. */
   return 0;
}

int DistanceSquared(int color)
{
   r = seed_point[--color][0] - px;
   t = seed_point[color][1] - py;
   return r * r + t * t;
}

int Mod(int u)
{
   return u > 0 ? u % MAX_GRID_SIZE
                : (MAX_GRID_SIZE - (-u % MAX_GRID_SIZE)) % MAX_GRID_SIZE;
}

int main(int argc, char **argv)
{
   if( argc - 4 )
      return printf("%s <input.png> <output1.png> <output2.png>\n", *argv);

   image.version = PNG_IMAGE_VERSION;
   if( png_image_begin_read_from_file(&image, argv[1]) )
   {
      image.format = PNG_FORMAT_RGBA;
      w = image.width;
      h = image.height;
      for(s = w * h; i < 2; i++)
         buffer[i] = (png_bytep)calloc(s, 4);
      t = *buffer && buffer[1] &&
          png_image_finish_read(&image, NULL, *buffer, 0, NULL);
   }
   if( !t )
      return printf("Error reading %s\n", argv[1]);

   srand(time(NULL));

   pixel = (int*)buffer[1];

   if( sizeof('c') > 1 )
   {
      /* C mode: hexagons. */

      /* Initialize transformation from screen coordinates to grid coordinates.

         px = t[0] * x + t[1] * y + t[2]
         py = t[3] * x + t[4] * y + t[5]                                     */
      hx = Rnd(MAX_GRID_SIZE - MIN_GRID_SIZE) + MIN_GRID_SIZE;
      hx /= w > h ? w : h;
      transform[0] = transform[4] = cos(hy = Rnd(3.14159265358979323846264338327950288)) * hx;
      transform[3] = -(transform[1] = sin(hy) * hx);
      transform[2] = Rnd(4.0);
      transform[5] = Rnd(4.0);

      /* Stretch in X axis to compensate for hexagonal coordinates. */
      transform[0] /= hx = sqrt(3.0);
      transform[1] /= hx;

      /* Initialize random colors for each grid point. */
      for(y = 0; y < MAX_GRID_SIZE; y++)
         for(x = 0; x < MAX_GRID_SIZE;)
            grid[y][x++] = Rnd(4);

      /* Render hexagons. */
      for(y = 0; y < h; y++)
         for(x = 0; x < w;)
         {
            /* Convert to grid coordinates. */
            hx = transform[0] * x + transform[1] * y + transform[2];
            hy = transform[3] * x + transform[4] * y + transform[5];

            /* Convert to hexagonal coordinates.
            https://www.redblobgames.com/grids/hexagons/more-pixel-to-hex.html
            */
            t = floor(hx + hy);
            r = floor((floor(hy - hx) + t) / 3);

            pixel[y * w + x++] =
               grid[Mod(floor((floor(2 * hx + 1) + t) / 3) - r)][Mod(r)];
         }
   }
   else
   {
      /* C++ mode: voronoi diagrams. */

      /* Initialize seed points.

         Given a large enough image, the seed points will likely be
         distinct, but we don't enforce this since it saves a bit
         of code space.                                             */
      for(i = 0; i < REGION_COUNT; pixel[y * w + x] = ++i)
      {
         seed_point[i][0] = x = Rnd(w);
         seed_point[i][1] = y = Rnd(h);
      }

      /* Generate voronoi coloring using jumping flood algorithm.
         https://en.wikipedia.org/wiki/Jump_flooding_algorithm    */
      for(step = 1; step * 2 < w && step * 2 < h; step *= 2);
      for(; step > 0; step /= 2)
         for(y = 0; y < h; y++)
            for(x = 0; x < w; x++)
               /* Only process colored pixels. */
               if( (c_current = pixel[y * w + x]) )
                  for(i = -step; i <= step; i += step)
                  {
                     py = y + i;
                     if( py >= 0 && py < h )
                        for(j = -step; j <= step; j += step)
                        {
                           px = x + j;
                           if( px >= 0 && px < w )
                           {
                              c_neighbor = pixel[py * image.width + px];

                              /* Overwrite current pixel color if it's
                                 not yet assigned, or if distance to
                                 current seed pixel is further away.   */
                              pixel[py * image.width + px] =
                                 !c_neighbor ||
                                 DistanceSquared(c_neighbor) >
                                 DistanceSquared(c_current)
                                    ? c_current
                                    : c_neighbor;
                           }
                        }
                  }

      /* Build adjacency matrix for all regions. */
      for(i = 0; i < s; i++)
      {
         c_current = pixel[i] - 1;
         if( i % w )
         {
            /* Add edge for pixel to the left. */
            c_neighbor = pixel[i - 1] - 1;
            grid[c_current][c_neighbor] =
            grid[c_neighbor][c_current] = 1;
         }
         if( i >= w )
         {
            /* Add edge for pixel above. */
            c_neighbor = pixel[i - w] - 1;
            grid[c_current][c_neighbor] =
            grid[c_neighbor][c_current] = 1;
         }
      }

      /* Assign colors. */
      RecursiveAssignColors();

      /* Remap region colors to reduced color set. */
      for(i = 0; i < s; i++)
         pixel[i] = color_map[pixel[i] - 1] - 1;
   }

   /* Move pixels. */
   r = Rnd(8) + 2;
   t = Rnd(8) + 2;
   px = 2 * t;
   py = 2 * r;
   for(y = i = 0; y < h; y++)
      for(x = 0; x < w; x++, i++)
         if( (s = pixel[i])
                ? s > 2
                   ? 1
                   : s < 2
                      ? Rnd(2) < ((x % r * t > y % t * r) ^
                                  ((x / r + y / t) % 2))
                      : Rnd(2) > ((x % px * py > (py - 1 - y % py) * px) ^
                                  ((x / px + y / py) % 2))
               : 0 )
         {
            memcpy(buffer[1] + i * 4, *buffer + i * 4, 4);
            memset(*buffer + i * 4, 0, 4);
         }
         else
            pixel[i] = 0;

   for(i = 0; i < 2; i++)
      if( !png_image_write_to_file(&image, argv[i + 2], 0, buffer[i], 0, NULL) )
         return printf("Error writing %s\n", argv[i + 2]);

   return 0;
}