/*  qdnlfilt -- quick and dirty nonlinear image smoothing

    Copyright (C) 1996  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
    -----------

    This is the third program in the quick & dirty image processing
    series.  This one does adaptive smoothing using a simpler algorithm
    than the one used by pnmnlfilt.

    qdnlfilt version 1 works ONLY on 8-bit rawbits pixmaps and
    graymaps with maxval = 255, and ONLY on stdin and stdout.

    Usage:  qdnlfilt factor < input > output

    Nonlinear sharpening is not supported because it's a bonehead
    thing to do.  Linear sharpening works better.

    Factor     Result
    ------     ------
     <0.0      Illegal
      0.0      No operation
      0.2      Light nonlinear smoothing
      0.8      Heavy nonlinear smoothing
      1.0      Maximum smoothing
     >1.0      Illegal

    To compile:  gcc -O3 -o qdnlfilt qdnlfilt.c

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

    Version 1 1996-11-10
*/

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

unsigned char *prev, *this, *next;
int linlen;

void
usage ()
{
  fprintf (stderr, "Usage:  qdnlfilt factor < input > output\n\
\n\
    Factor     Result\n\
    ------     ------\n\
     <0.0      Illegal\n\
      0.0      No operation\n\
      0.2      Light nonlinear smoothing\n\
      0.8      Heavy nonlinear smoothing\n\
      1.0      Maximum smoothing\n\
     >1.0      Illegal\n");
  exit (-1);
}

void
bogus ()
{
  fprintf (stderr,
  "qdnlfilt 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
getline () {
  unsigned char *temp;
  temp = prev;
  prev = this;
  this = next;
  next = temp;
  assert (fread (next, 1, linlen, stdin) == linlen);
}

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

int
main (int argc, char **argv)
{
  int width, height, maxval, grayflag = 0, disp, x, y;
  float factor, noise, variance;
  char magicnum[3], junk, buf[1000];
  if (argc != 2)
    usage ();
  if (sscanf (argv[1], "%f", &factor) != 1)
    usage ();
  if (factor < 0.0 || factor > 1.0)
    usage ();
  noise = factor * 255.0;
  noise = noise * noise / 9.0;
#ifdef MSDOS
  setmode (fileno (stdin), O_BINARY);
  setmode (fileno (stdout), O_BINARY);
#endif

  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')
    grayflag = 1;
  printf ("%s\n%d %d\n%d\n", magicnum, width, height, maxval);

  if (width < 3 || height < 3) {
    /* Too small to do anything with; echo stdin to stdout */
    int tval;
    while ((tval = getchar ()) != EOF)
      assert (putchar (tval) != EOF);
    fflush (stdout);
    exit (0);
  }

  disp = (grayflag? 1 : 3);
  linlen = width * disp;
  assert (prev = malloc (linlen));
  assert (this = malloc (linlen));
  assert (next = malloc (linlen));

  getline();
  getline();

  assert (fwrite (this, 1, linlen, stdout) == linlen);
  for (y=1; y<height-1; y++) {
    getline();
    for (x=0; x<disp; x++)
      assert (putchar ((int)(this[x])) != EOF);
    for (x=1; x<width-1; x++) {
      if (grayflag) {
        float mean, temp;
        mean = prev[x-disp] + this[x-disp] + next[x-disp]
             + prev[x]      + this[x]      + next[x]
             + prev[x+disp] + this[x+disp] + next[x+disp];
        mean /= 9.0;
        variance = 0.0;
        temp = (float)(prev[x-disp]) - mean;
        variance += temp * temp;
        temp = (float)(this[x-disp]) - mean;
        variance += temp * temp;
        temp = (float)(next[x-disp]) - mean;
        variance += temp * temp;
        temp = (float)(prev[x]) - mean;
        variance += temp * temp;
        temp = (float)(this[x]) - mean;
        variance += temp * temp;
        temp = (float)(next[x]) - mean;
        variance += temp * temp;
        temp = (float)(prev[x+disp]) - mean;
        variance += temp * temp;
        temp = (float)(this[x+disp]) - mean;
        variance += temp * temp;
        temp = (float)(next[x+disp]) - mean;
        variance += temp * temp;
        variance /= 9.0;
        if (variance > 0.0) {
          temp = mean + (variance * ((float)(this[x]) - mean)) /
                 (variance + noise) + 0.5;
          assert (putchar ((int)temp) != EOF);
        } else
          assert (putchar ((int)(this[x])) != EOF);
      } else {
        float rmean, gmean, bmean, temp;
        int xx = x * 3;
        variance = 0.0;
        rmean = prev[xx-disp] + this[xx-disp] + next[xx-disp]
              + prev[xx]      + this[xx]      + next[xx]
              + prev[xx+disp] + this[xx+disp] + next[xx+disp];
        rmean /= 9.0;
        temp = (float)(prev[xx-disp]) - rmean;
        variance += temp * temp;
        temp = (float)(this[xx-disp]) - rmean;
        variance += temp * temp;
        temp = (float)(next[xx-disp]) - rmean;
        variance += temp * temp;
        temp = (float)(prev[xx]) - rmean;
        variance += temp * temp;
        temp = (float)(this[xx]) - rmean;
        variance += temp * temp;
        temp = (float)(next[xx]) - rmean;
        variance += temp * temp;
        temp = (float)(prev[xx+disp]) - rmean;
        variance += temp * temp;
        temp = (float)(this[xx+disp]) - rmean;
        variance += temp * temp;
        temp = (float)(next[xx+disp]) - rmean;
        variance += temp * temp;
        xx++;
        gmean = prev[xx-disp] + this[xx-disp] + next[xx-disp]
              + prev[xx]      + this[xx]      + next[xx]
              + prev[xx+disp] + this[xx+disp] + next[xx+disp];
        gmean /= 9.0;
        temp = (float)(prev[xx-disp]) - gmean;
        variance += temp * temp;
        temp = (float)(this[xx-disp]) - gmean;
        variance += temp * temp;
        temp = (float)(next[xx-disp]) - gmean;
        variance += temp * temp;
        temp = (float)(prev[xx]) - gmean;
        variance += temp * temp;
        temp = (float)(this[xx]) - gmean;
        variance += temp * temp;
        temp = (float)(next[xx]) - gmean;
        variance += temp * temp;
        temp = (float)(prev[xx+disp]) - gmean;
        variance += temp * temp;
        temp = (float)(this[xx+disp]) - gmean;
        variance += temp * temp;
        temp = (float)(next[xx+disp]) - gmean;
        variance += temp * temp;
        xx++;
        bmean = prev[xx-disp] + this[xx-disp] + next[xx-disp]
              + prev[xx]      + this[xx]      + next[xx]
              + prev[xx+disp] + this[xx+disp] + next[xx+disp];
        bmean /= 9.0;
        temp = (float)(prev[xx-disp]) - bmean;
        variance += temp * temp;
        temp = (float)(this[xx-disp]) - bmean;
        variance += temp * temp;
        temp = (float)(next[xx-disp]) - bmean;
        variance += temp * temp;
        temp = (float)(prev[xx]) - bmean;
        variance += temp * temp;
        temp = (float)(this[xx]) - bmean;
        variance += temp * temp;
        temp = (float)(next[xx]) - bmean;
        variance += temp * temp;
        temp = (float)(prev[xx+disp]) - bmean;
        variance += temp * temp;
        temp = (float)(this[xx+disp]) - bmean;
        variance += temp * temp;
        temp = (float)(next[xx+disp]) - bmean;
        variance += temp * temp;
        variance /= 27.0;
        if (variance > 0.0) {
          temp = rmean + (variance * ((float)(this[3*x]) - rmean)) /
                 (variance + noise) + 0.5;
          assert (putchar ((int)temp) != EOF);
          temp = gmean + (variance * ((float)(this[3*x+1]) - gmean)) /
                 (variance + noise) + 0.5;
          assert (putchar ((int)temp) != EOF);
          temp = bmean + (variance * ((float)(this[3*x+2]) - bmean)) /
                 (variance + noise) + 0.5;
          assert (putchar ((int)temp) != EOF);
        } else {
          assert (putchar ((int)(this[3*x])) != EOF);
          assert (putchar ((int)(this[3*x+1])) != EOF);
          assert (putchar ((int)(this[3*x+2])) != EOF);
        }
      }
    }
    for (x=linlen-disp; x<linlen; x++)
      assert (putchar ((int)(this[x])) != EOF);
  }
  assert (fwrite (next, 1, linlen, stdout) == linlen);

  fflush (stdout);
  exit (0);
}