/*  qdvignet -- correct vignetting / falloff in digitized pictures
    (quick & dirty)

    Copyright (C) 1998  David Flater.

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    If you cannot find a copy of the GNU General Public License,
    write to the Free Software Foundation, Inc., 675 Mass Ave,
    Cambridge, MA 02139, USA.


    Description
    -----------

    Cheap lenses tend to produce pictures with dark corners.  The
    darkening gets worse the farther you get from the center of the
    picture.  qdvignet is a special-purpose filter to compensate for
    this problem in digital images.

    According to the Photographic Lenses FAQ v1.18, falloff naturally
    occurs at the rate of cos^4 of the angle off axis, but a good
    optical design will compensate for this somewhat.  So, this
    program just uses a linear interpolation that yields subjectively
    good results for whatever optical design my cheap P&S cameras
    have.  The only objectionable visual artifact produced by this
    method is the discontinuity in the first derivative of the
    correction function around the cutoff radius.  In plain English,
    that means that you might perceive a circle around the image at
    the point where the correction starts to take effect.

    The qd* series of filters are "quick and dirty" standalone
    programs that use the PPM and PGM image formats.  qdvignet version
    1 works ONLY on 8-bit rawbits pixmaps and graymaps with maxval =
    255, and ONLY on stdin and stdout.

    Usage:  qdvignet cor [radius] [-l | -r] < input > output

    cor is the maximum brightness correction that will be applied in
    the farthest corners.  cor above 1 will brighten the corners; cor
    between 0 and 1 will darken them.

    The radius selects the cutoff for where the brightening or
    darkening should begin.  Specifying -l or -r restricts the
    filtering to only the left or right side of the image
    respectively.

    To compile: gcc -O3 -o qdvignet qdvignet.c -lm


    Version 2 2000-02-21
      Hastily added code to avoid choking on comments in PNM files.

    Version 1 1998-05-23
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#ifdef MSDOS
#include <io.h>
#include <fcntl.h>
#endif

void
usage ()
{
  fprintf (stderr, "\
    Usage:  qdvignet cor [radius] [-l | -r] < input > output\n\
\n\
    cor is the maximum brightness correction that will be applied in\n\
    the farthest corners.  cor above 1 will brighten the corners; cor\n\
    between 0 and 1 will darken them.\n\
\n\
    The radius selects the cutoff for where the brightening or\n\
    darkening should begin.  Specifying -l or -r restricts the\n\
    filtering to only the left or right side of the image\n\
    respectively.\n");
  exit (-1);
}

void
bogus ()
{
  fprintf (stderr,
  "qdvignet works ONLY on 8-bit rawbits pixmaps and\n");
  fprintf (stderr,
  "graymaps with maxval = 255, and ONLY on stdin and stdout.\n");
  exit (-1);
}

void
getnoncommentline (char buf[1000]) {
  do
    assert (fgets (buf, 1000, stdin));
  while (buf[0] == '#');
}

int
main (int argc, char **argv)
{
  int looper, tval, width, height, maxval, bpp, x, y, whichside = 0, gotminr;
  double cor, minr, maxr, midx, midy, r, dx, dy;
  char magicnum[3], junk, buf[1000], *lr;

#ifdef MSDOS
  setmode (fileno (stdin), O_BINARY);
  setmode (fileno (stdout), O_BINARY);
#endif

  /* Parse command line */
  if (argc < 2 || argc > 4)
    usage ();
  if (sscanf (argv[1], "%lf", &cor) != 1)
    usage ();
  if (cor <= 0.0)
    usage ();
  lr = NULL;
  gotminr = 0;
  if (argc > 2) {
    gotminr = sscanf (argv[2], "%lf", &minr);
    if (argc == 4) {
      if (gotminr)
        lr = argv[3];
      else
        usage();
    } else if (!gotminr)
      lr = argv[2];
  }
  if (lr) {
    if (!strcmp (lr, "-l"))
      whichside = 1;
    else if (!strcmp (lr, "-r"))
      whichside = 2;
    else
      usage ();
  }

  getnoncommentline (buf);
  if (sscanf (buf, "%2s", magicnum) != 1) {
    fprintf (stderr, "Error getting header data\n");
    bogus();
  }
  getnoncommentline (buf);
  if (sscanf (buf, "%d %d", &width, &height) != 2) {
    fprintf (stderr, "Error getting header data\n");
    bogus();
  }
  getnoncommentline (buf);
  if (sscanf (buf, "%d", &maxval) != 1) {
    fprintf (stderr, "Error getting header data\n");
    bogus();
  }
  if (strcmp (magicnum, "P5") != 0 && strcmp (magicnum, "P6") != 0) {
    fprintf (stderr, "Bad magic number\n");
    bogus ();
  }
  if (maxval != 255) {
    fprintf (stderr, "Maxval != 255\n");
    bogus ();
  }

  if (magicnum[1] == '5')
    bpp = 1;
  else
    bpp = 3;
  assert (width > 0 && height > 0);
  printf ("%s\n%d %d\n%d\n", magicnum, width, height, maxval);

  /* More figuring */
  if (!gotminr)
    minr = (double)(width < height ? width : height) / 2.0;
  midx = (double)width / 2.0;
  midy = (double)height / 2.0;
  maxr = sqrt (midx * midx + midy * midy);

  /* Main loop */
  for (y=0; y<height; y++) {
    dy = (double)y - midy;
    for (x=0; x<width; x++) {
      dx = (double)x - midx;
      r = sqrt (dx * dx + dy * dy);
      for (looper=0; looper < bpp; looper++) {
        tval = getchar();
        if (!((whichside == 1 && dx > 0.0) ||
              (whichside == 2 && dx < 0.0) ||
              (r < minr))) {
          /* Linear interp */
          tval = (int)((double)tval + ((cor - 1.0) * (double)tval) *
           ((r - minr) / (maxr - minr)) + 0.5);
          if (tval > 255)
            tval = 255;
        }
        assert (putchar (tval) != EOF);
      }
    }
  }
  fflush (stdout);
  exit (0);
}