/*
uncomposite-raw.c 0.0.1 (2004-Feb-12-Thu)
http://www.nicemice.net/amc/utils/#ppmuncomposite
Adam M. Costello
http://www.nicemice.net/amc/

uncomposite-raw <onblack> <onwhite> <uncomp> <alpha>

The first two files are inputs, the second two are outputs.  The first
two contain RGB triplets already composited onto black and white,
respectively.  The second two will contain the reverse-engineered
uncomposited (non-premultiplied) RGB triplets and alpha samples,
respectively.  There is no guarantee that the two input files can
be reproduced by compositing the output onto black and white;
approximations may be necessary in general.  But if the two input
files were created by compositing an RGBA image onto black and white
backgrounds, then this program will reconstruct that RGBA image, to
within roundoff error.

All samples are 16-bit unsigned, most significant byte first
(big-endian).  The RGB samples are intensity samples, not gamma-encoded.
This program does not depend on the image dimensions, nor on the order
of the red, green, and blue samples within each triplet.

This program should perhaps take more general inputs and have more
graceful error handling (rather than just calling assert() everywhere),
but it was written quickly and was only designed to be called by a
particular script (ppmuncomposite).

*/

#include <assert.h>
#include <stdio.h>

#define readbyte(x) ((unsigned int) getc(x))

int main(int argc, char **argv)
{
  FILE *onblack, *onwhite, *uncomp, *alpha;
  unsigned int b[3], w[3], u, a;
    /* b[i] and w[i] are samples read from onblack and onwhite. */
    /* u and a are samples to be written to uncomp and alpha.   */
  int i, c;
  double backweight, foreweight, x;

  assert(argc == 5);
  onblack = fopen(argv[1], "rb");
  assert(onblack != 0);
  onwhite = fopen(argv[2], "rb");
  assert(onwhite != 0);
  uncomp = fopen(argv[3], "wb");
  assert(uncomp != 0);
  alpha = fopen(argv[4], "wb");
  assert(alpha != 0);

  for (;;) {
    /* Check for end-of-file: */
    c = getc(onblack);
    if (c == EOF) break;
    ungetc(c,onblack);

    backweight = 0.0;  /* see below */

    for (i = 0;  i < 3;  ++i) {
      /* Read a sample from onblack and onwhite: */
      b[i] = (readbyte(onblack) << 8) + readbyte(onblack);
      w[i] = (readbyte(onwhite) << 8) + readbyte(onwhite);

      /* The formula for compositing a foreground sample and a        */
      /* background sample is:                                        */
      /*                                                              */
      /* foreweight = alpha                                           */
      /* backweight = 1 - alpha                                       */
      /* compsamp = foreweight * foresamp + backweight * backsamp     */
      /*                                                              */
      /* When backsamp is 0 (black) or 1 (white), the formula         */
      /* simplifies to:                                               */
      /*                                                              */
      /* onblacksamp = foreweight * foresamp                          */
      /* onwhitesamp = foreweight * foresamp + backweight             */
      /*                                                              */
      /* Subtracting yields:                                          */
      /*                                                              */
      /* backweight = onwhitesamp - onblacksamp                       */
      /*                                                              */
      /* We can compute that difference for the red, green, and blue  */
      /* samples, and we'll get three separate values for backweight, */
      /* but the output must have a single value for alpha, which is  */
      /* 1 - backweight.  So we'll take the mean of the three values, */
      /* by summing them here and dividing by 3 after the loop.       */

      backweight += (double) w[i] - b[i];

      /* Since w[i] and b[i] are scaled with 65535 representing 1.0,  */
      /* backweight is similarly scaled.                              */
    }

    backweight /= 3.0;

    /* Verify that we didn't reach end-of-file */
    /* in the middle of reading the triplets:  */
    assert(!feof(onblack));
    assert(!feof(onwhite));

    /* The foreground and background weights sum */
    /* to one, which is represented by 65535.    */
    foreweight = 65535.0 - backweight;

    for (i = 0;  i < 3;  ++i) {
      /* If the foreground weight is zero, the uncomposited sample is */
      /* irrelevant, and we need to avoid a division by zero, so just */
      /* set it to zero:                                              */

      if (foreweight == 0.0) u = 0;
      else {
        /* Taking into account the scale factor of 65535, */
        /* the compositing formulas are (ideally):        */
        /*                                                */
        /* b[i] = u[i]*foreweight/65535                   */
        /* w[i] = u[i]*foreweight/65535 + backweight      */
        /*                                                */
        /* We'd like to find a u[i] that satisfies both,  */
        /* but that might not be possible, so we solve    */
        /* both and take the mean, and round:             */

        x = (-backweight + w[i] + b[i]) * 32767.5 / foreweight;
        u = x <= 0.0 ? 0 : x >= 65535.0 ? 65535 : x + 0.5;
      }

      putc(u >> 8, uncomp);
      putc(u & 255, uncomp);
    }

    a = foreweight <=     0.0 ?     0 :
        foreweight >= 65535.0 ? 65535 :
        foreweight + 0.5 ;
    putc(a >> 8, alpha);
    putc(a & 255, alpha);
  }

  /* Verify that the other input is also exhausted: */
  c = getc(onwhite);
  assert(c == EOF);
  return 0;
}