/* chinatsu.c - Don Yang (uguu.org)

   01/02/12
                                                                           */
#include<stdio.h>
#include<stdlib.h>
#include<string.h>

typedef struct
{
   int min_p, max_p;  /* Primary pixel ranges */
   int min_q, max_q;  /* Secondary pixel ranges */
} Ranges;

/* Parse PPM header, returning offset to pixel data */
static int ParseHeader(FILE *infile, int *width, int *height)
{
   char buffer[256];
   int i, c;

   fseek(infile, 0, SEEK_SET);
   if( (int)fread(buffer, 1, sizeof(buffer), infile) <= 0 )
      return 0;
   if( sscanf(buffer, "P6 %d %d %d", width, height, &i) != 3 || i != 255 )
      return 0;

   for(i = c = 0; i < sizeof(buffer); i++)
   {
      if( buffer[i] == '\n' )
      {
         if( ++c == 3 )
            return i + 1;
      }
   }
   return 0;
}

/* Scale an input value from [0,256) to new range */
static int ScaleValue(int input, int range_min, int range_max)
{
   return input * (range_max - range_min) / 256 + range_min;
}
static float ScaleValueF(int input, int range_min, int range_max)
{
   return (float)input * (range_max - range_min) / 256.f + range_min;
}

/* Compute area of the selected range */
static int ComputeArea(const Ranges *range)
{
   return (range->max_p - range->min_p) * (range->max_q - range->min_q);
}

/* Check if all pixel combinations will be valid after applying
   the specified scaling parameters, returns 0 if not.          */
static int Fit(const Ranges *range, char pixels[256][256])
{
   int x, y, x0, y0, k;

   for(x = 0; x < 256; x++)
   {
      x0 = ScaleValue(x, range->min_p, range->max_p);
      for(y = 0; y < 256; y++)
      {
         y0 = ScaleValue(y, range->min_q, range->max_q);
         if( pixels[x][y] == 0 )
            continue;
         k = (3 * x0 - y0) / 2;
         if( k < 0 || k > 255 )
            return 0;
      }
   }
   return 1;
}

/* Compute scaling parameters */
static void GetScaledRange(FILE *primary, FILE *secondary,
                           int width, int height,
                           unsigned char *buffer1, unsigned char *buffer2,
                           Ranges *ranges)
{
   char pixels[256][256];  /* [primary][secondary] */
   int x, y, step, ip0, ip1, iq0, iq1;
   Ranges best, z;
   int best_area, z_area;

   /* Mark all pixel combinations as unused */
   memset(pixels, 0, sizeof(pixels));

   /* Record all pixel combinations */
   for(y = 0; y < height; y++)
   {
      fread(buffer1, width * 3, 1, primary);
      fread(buffer2, width * 3, 1, secondary);
      for(x = 0; x < width * 3; x++)
         pixels[buffer1[x]][buffer2[x]] = 1;
   }

   /* Set initial range and try to expand outward.  This initial range
      guarantees that all pixel combinations will be valid, but it may
      waste some pixel resolution if the two images are similar.  We
      will try to expand this in smaller increments to find the
      biggest range where the pixel combinations are still valid.      */
   ranges->min_p = ranges->min_q = 64;
   ranges->max_p = ranges->max_q = 256 - 64;
   memcpy(&best, ranges, sizeof(best));
   best_area = ComputeArea(&best);
   for(step = 64; step > 0; step /= 2)
   {
      for(ip0 = -1; ip0 <= 0; ip0++)
      {
         z.min_p = ranges->min_p + ip0 * step;
         if( z.min_p < 0 )
            continue;
         for(iq0 = -1; iq0 <= 0; iq0++)
         {
            z.min_q = ranges->min_q + iq0 * step;
            if( z.min_q < 0 )
               continue;
            for(ip1 = 1; ip1 >= 0; ip1--)
            {
               z.max_p = ranges->max_p + ip1 * step;
               if( z.max_p > 256 )
                  continue;
               for(iq1 = 1; iq1 >= 0; iq1--)
               {
                  z.max_q = ranges->max_q + iq1 * step;
                  if( z.max_q > 256 )
                     continue;

                  /* Ignore combinations that yields areas that are
                     no better than current best.                   */
                  z_area = ComputeArea(&z);
                  if( z_area <= best_area )
                     continue;

                  /* Keep new ranges if all pixels will fit */
                  if( Fit(&z, pixels) )
                  {
                     memcpy(&best, &z, sizeof(z));
                     best_area = z_area;
                  }
               }
            }
         }
      }
      memcpy(ranges, &best, sizeof(best));
   }
}

/* Convert floating point pixel value to integer */
static unsigned char ConvertToInt(float p)
{
   int p0 = (int)p;
   float r = p - (float)p0;

   /* Use random number generator for dithering.  This random number
      generator is not seeded, so the output may be deterministic.
      We don't really care either way for our purpose.              */
   if( ((float)rand() / (float)RAND_MAX) < r )
      p0++;
   if( p0 < 0 )
      return 0;
   if( p0 > 255 )
      return 255;
   return p0;
}

/* Convert input pixel combinations to output pixel values */
static void MergePixels(const Ranges *ranges,
                        unsigned char p_in,
                        unsigned char q_in,
                        unsigned char *p_out,
                        unsigned char *q_out,
                        unsigned char *c1_out,
                        unsigned char *c2_out)
{
   float p = ScaleValueF(p_in, ranges->min_p, ranges->max_p);
   float q = ScaleValueF(q_in, ranges->min_q, ranges->max_q);
   float c = (3 * p - q) / 2.0f;
   *p_out = ConvertToInt(p);
   *q_out = ConvertToInt(q);
   *c1_out = *c2_out = ConvertToInt(c);
}

/* Merge two images to output */
static void MergeImages(const Ranges *ranges, int width, int height, int even,
                        unsigned char *buffer1, unsigned char *buffer2,
                        FILE *primary, FILE *secondary, FILE *output)
{
   /* Pixel schedules */
   static int index[2][12] =
   {
      /* Rc Gp Bc Rc Gc Bp
         Rp Gc Bc Rq Gq Bq */
      {6, 9, 0, 3,
       1, 10, 4, 7,
       5, 11, 2, 8},
      /* Rq Gq Bq Rp Gc Bc
         Rc Gc Bp Rc Gp Bc */
      {3, 0, 6, 9,
       10, 1, 4, 7,
       8, 2, 5, 11}
   };

   unsigned char pixel[12];
   int x, y, scanline;

   for(y = 0; y < height; y++)
   {
      fread(buffer1, width * 3, 1, primary);
      fread(buffer2, width * 3, 1, secondary);

      for(scanline = 0; scanline <= 6; scanline += 6)
      {
         for(x = 0; x < width; x++)
         {
            /* R */
            MergePixels(ranges, buffer1[x * 3], buffer2[x * 3],
                        &pixel[index[even][0]],
                        &pixel[index[even][1]],
                        &pixel[index[even][2]],
                        &pixel[index[even][3]]);
            /* G */
            MergePixels(ranges, buffer1[x * 3 + 1], buffer2[x * 3 + 1],
                        &pixel[index[even][4]],
                        &pixel[index[even][5]],
                        &pixel[index[even][6]],
                        &pixel[index[even][7]]);
            /* B */
            MergePixels(ranges, buffer1[x * 3 + 2], buffer2[x * 3 + 2],
                        &pixel[index[even][8]],
                        &pixel[index[even][9]],
                        &pixel[index[even][10]],
                        &pixel[index[even][11]]);

            fwrite(pixel + scanline, 6, 1, output);
         }
      }
   }
}

int main(int argc, char **argv)
{
   FILE *primary, *secondary, *output;
   int status, width, height, offset1, width2, height2, offset2;
   Ranges ranges;
   unsigned char *buffer1, *buffer2;

   /* Check command line arguments */
   if( argc < 4 )
   {
      printf("%s <primary.ppm> <secondary.ppm> <output.ppm> [odd]\n", *argv);
      return 1;
   }

   /* Check primary image */
   if( (primary = fopen(argv[1], "rb")) == NULL )
   {
      printf("Error opening %s for reading\n", argv[1]);
      return 1;
   }
   do
   {
      status = 1;
      offset1 = ParseHeader(primary, &width, &height);
      if( offset1 == 0 )
      {
         printf("Error reading header from %s\n", argv[1]);
         break;
      }

      /* Check secondary image */
      if( (secondary = fopen(argv[2], "rb")) == NULL )
      {
         printf("Error opening %s for reading\n", argv[2]);
         break;
      }
      do
      {
         offset2 = ParseHeader(secondary, &width2, &height2);
         if( offset2 == 0 )
         {
            printf("Error reading header from %s\n", argv[2]);
            break;
         }
         if( width != width2 || height != height2 )
         {
            puts("Image dimensions mismatched");
            break;
         }

         /* Allocate scanline buffers.  These usually succeeds unless
            image dimensions were bad.                                */
         buffer1 = (unsigned char*)malloc(width * 3);
         buffer2 = (unsigned char*)malloc(width * 3);
         if( buffer1 == NULL || buffer2 == NULL )
         {
            printf("Not enough memory for %d pixels\n", width);
            break;
         }

         /* Compute scaling parameters */
         fseek(primary, offset1, SEEK_SET);
         fseek(secondary, offset2, SEEK_SET);
         GetScaledRange(primary, secondary, width, height, buffer1, buffer2,
                        &ranges);

         do
         {
            /* Open output */
            if( (output = fopen(argv[3], "wb+")) == NULL )
            {
               printf("Can not open %s for writing\n", argv[3]);
               break;
            }

            /* Write output header */
            fprintf(output, "P6\n%d %d\n255\n", width * 2, height * 2);

            /* Merge images */
            fseek(primary, offset1, SEEK_SET);
            fseek(secondary, offset2, SEEK_SET);
            MergeImages(&ranges, width, height, (argc > 4) ? 1 : 0,
                        buffer1, buffer2,
                        primary, secondary, output);

            /* Cleanup */
            status = fclose(output);
         } while(0);
         free(buffer1);
         free(buffer2);
      } while(0);
      fclose(secondary);
   } while(0);
   fclose(primary);
   return status;
}