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

#define PI 3.14159265358979323846264338327950288419716939937510
#define SQUARE_CORNER 0.25

png_image image;
png_bytep buffer[2];
double fx, fy, rx, ry,

/* Transformations. */
       g0, g1, g2, g5, h0, h1, h2,

/* Projected grid coordinates. */
       grid_x, grid_y,

/* Computed dot sizes. */
       min_dot, max_dot;

int s, i, j, k, x, y, d, w, h,

/* Distortion parameters. */
    phase[2][2], period[2][2];

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

/* Apply transformation from screen coordinates to grid coordinates,
   storing results in (grid_x, grid_y).                              */
void ApplyGridTransform(double u, double v)
{
   grid_x = g0 * u - g1 * v + g2;
   grid_y = g1 * u + g0 * v + g5;
}

/* Check if a grid point is within radius of a single nearby halftone dot,
   returns 1 if so.                                                        */
int WithinDotRadius(double px, double py, double target_x, double target_y)
{
   px -= target_x;
   py -= target_y;

   fx = h0 * target_x + h1 * target_y + h2;
   return px * px + py * py < fx * fx;
}

/* Determine if a grid point is within the radius of nearby halftone dots,
   returns 1 if so.                                                        */
int WithinRadiusOfNearbyDots()
{
   return WithinDotRadius(grid_x, grid_y, floor(grid_x), floor(grid_y)) ||
          WithinDotRadius(grid_x, grid_y, floor(grid_x), ceil(grid_y)) ||
          WithinDotRadius(grid_x, grid_y, ceil(grid_x), floor(grid_y)) ||
          WithinDotRadius(grid_x, grid_y, ceil(grid_x), ceil(grid_y));
}

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]) )
   {
      w = image.width;
      h = image.height;
      for(image.format = PNG_FORMAT_RGBA; i < 2;)
         buffer[i++] = (png_bytep)malloc(s = w * h * 4);
      d = *buffer && buffer[1] &&
          png_image_finish_read(&image, NULL, *buffer, 0, NULL);
   }
   if( !d )
      return printf("Error reading %s\n", argv[1]);

   memcpy(buffer[1], *buffer, s);

   srand(time(NULL));

   /* Generate transformation from screen coordinates to grid coordinates.

      Halftone dots will be centered around integer intervals in grid
      coordinates, maximum distance between the center of any two dots
      will be sqrt(2).

      grid_x = g[0] * x + g[1] * y + g[2]
      grid_y = g[3] * x + g[4] * y + g[5]                                  */
   fx = Rnd(90) + 30;
   fx /= w > h ? w : h;

   g0 = cos(fy = Rnd(PI)) * fx;
   g1 = sin(fy) * fx;
   g2 = Rnd(2);
   g5 = Rnd(2);

   /* Generate a second transformation from grid coordinates to
      gradient coordinates.  The X value here will be used to
      compute halftone dot sizes.

      dot_size = h[0] * grid_x + h[1] * grid_y + h[2]

      Only the rotation matter for this transformation, we will
      compute scale and offset such that the dot radii are zero
      in one corner and 0.5*sqrt(2) in the other corner.        */
   h0 = cos(fy = Rnd(2.0 * PI));
   h1 = sin(fy);
   h2 = 0;

   /* Apply rotation to the 4 corners of the rectangle enclosed by
      (0,0) .. (width,height)

      Note that it goes up to (width,height) as opposed to
      (width-1,height-1), thus the rectangle is always at least one
      pixel wide.                                                   */
   for(; j < 4; j++)
   {
      ApplyGridTransform(j % 2 * w, (j & 2) / 2 * h);

      fx = h0 * grid_x + h1 * grid_y + h2;
      min_dot = !j || min_dot > fx ? fx : min_dot;
      max_dot = !j || max_dot < fx ? fx : max_dot;
   }

   /* Update gradient transform so that corners of the screen cover
      the full size range of [0 .. 0.5*sqrt(2)].                    */
   h0 *= fx = 0.5 * sqrt(2) / (max_dot - min_dot);
   h1 *= fx;

   /* Generate distortion parameters. */
   for(h2 = -min_dot * fx; x < 2; x++)
      for(y = 0; y < 2; period[x][y++] = Rnd(6) + 2)
         phase[x][y] = Rnd(8);

   /* Render pixels. */
   for(y = i = 0; y < h; y++)
      for(x = 0; x < w; x++, i++)
      {
         ApplyGridTransform(x, y);
         if( sizeof('c') > 1 )
         {
            /* C mode: halftone dots. */
            s = WithinRadiusOfNearbyDots();
         }
         else
         {
            /* C++ mode: distortion squares. */

            /* Check if pixel is near corner of the grid cells.
               If not, set parity based on square coordinates.  */
            fx = floor(grid_x);
            fy = floor(grid_y);
            rx = grid_x - fx;
            ry = grid_y - fy;
            if( rx + ry < SQUARE_CORNER ||
                rx + ry > 2.0 - SQUARE_CORNER ||
                ry - rx > 1.0 - SQUARE_CORNER ||
                rx - ry > 1.0 - SQUARE_CORNER )
            {
               /* Near corner, determine which diagonal line the
                  pixels are on.                                 */
               rx = round(grid_x);
               ry = round(grid_y);
               d = (int)(rx + ry) & 1;

               /* Apply wave function specific to that diagonal line,
                  such that parity alternates every few corners along
                  that line.                                          */
               for(j = s = 0; j < 2; s ^= k > period[j++][d] / 2)
                  if( (k = fmod(phase[j][d] + rx - j * ry, period[j][d])) < 0 )
                     k += period[j][d];

               /* Render corner as a 2x2 checkerboard. */
               s ^= grid_x - rx > fabs(grid_y - ry) ||
                    grid_x - rx < -fabs(grid_y - ry);
            }
            else
               /* Away from corners. */
               s = (int)(fx + fy) & 1;
         }
         memset(buffer[s] + i * 4, 0, 4);
      }

   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;
}