/*****************************************************************
 * flblue.c: FBM Release 1.0 25-Feb-90 Michael Mauldin
 *
 * Copyright (C) 1989,1990 by Michael Mauldin.  Permission is granted
 * to use this file in whole or in part for any purpose, educational,
 * recreational or commercial, provided that this copyright notice
 * is retained unchanged.  This software is available to all free of
 * charge by anonymous FTP and in the UUNET archives.
 *
 * flblue.c: Blue noise dithering
 *
 * CONTENTS
 *	bluenoise_fbm (input, output, noiselevel)
 *
 * EDITLOG
 *	LastEditDate = Mon Jun 25 00:04:47 1990 - Michael Mauldin
 *	LastFileName = /usr2/mlm/src/misc/fbm/flblue.c
 *
 * HISTORY
 * 25-Jun-90  Michael Mauldin (mlm@cs.cmu.edu) Carnegie Mellon
 *	Package for Release 1.0
 *
 * 20-May-89  Michael Mauldin (mlm) at Carnegie Mellon University
 *	Bug fix from Dave Cohrs <dave@cs.wisc.edu>
 *
 * 07-Mar-89  Michael Mauldin (mlm) at Carnegie Mellon University
 *	Beta release (version 0.9) mlm@cs.cmu.edu
 *
 * 12-Nov-88  Michael Mauldin (mlm) at Carnegie-Mellon University
 *	Created.
 *****************************************************************/

# include <stdio.h>
# include <math.h>
# include <ctype.h>
# include "fbm.h"

/*****************************************************************
 * bluenoise_fbm:  Do Floyd-Steinberg halftoning with serpentine
 *		   raster and 'noiselevel' random weights.
 *		   (noise level runs from 0 to 100 percent)
 *
 * REFERENCES
 *	Digital Halftoning, by Robert Ulichney (1986 MIT Press)
 *****************************************************************/

# define RAND(RN) (((seed = 1103515245 * seed + 12345) >> 12) % (RN))
# define INITERR(X,Y) \
	(((int) X) - (((int) Y)?WHITE:BLACK) + ((WHITE/2)-((int) X))/2)

#ifndef lint
static char *fbmid =
"$FBM flblue.c <1.0> 25-Jun-90  (C) 1989,1990 by Michael Mauldin, source \
code available free from MLM@CS.CMU.EDU and from UUNET archives$";
#endif

bluenoise_fbm (input, output, noiselevel)
FBM *input, *output;
double noiselevel;
{ register unsigned char *bmp, *obm;
  register unsigned seed = 0;
  register int i, j, rowlen, gray, error, w, h, den, outrow;
  int w1, w3, w5, w7, smrange, lgrange, smnoise, lgnoise;
  int *eerr, *oerr;

  if (input->hdr.planes != 1)
  { fprintf (stderr, "bluenoise_fbm: can't halftone color images\n");
    return (0);
  }

  fprintf (stderr, "Blue noise, %1.2lf%% weights\n", noiselevel);

  /* Allocate output */
  free_fbm (output);
  output->hdr = input->hdr;
  output->hdr.bits = 1;
  output->hdr.physbits = 8;
  outrow = 16 * ((input->hdr.cols + 15) / 16); /* Pad to even byte boundary */
  output->hdr.rowlen = outrow;
  output->hdr.plnlen = outrow*output->hdr.rows;
  alloc_fbm (output);

  w = input->hdr.cols;
  h = input->hdr.rows;
  rowlen = input->hdr.rowlen;
  
  /* Allocate space for error arrays */
  eerr = (int *) malloc ((unsigned) w * sizeof (*eerr));
  oerr = (int *) malloc ((unsigned) w * sizeof (*oerr));
  for (i=0; i<w; i++) eerr[i] = oerr[i] = 0;

    /* The left border */
  error = 0;
  for (j=0; j<h; j++)
  { register int thresh = (WHITE/2 + RAND (129) - 64);

    gray = input->bm[j*rowlen] + error;
    den = gray > thresh ? WHITE : BLACK;
    error = gray - den;
    output->bm[j*outrow] = den;
  }

  /* The right border */
  error = 0;
  for (j=0; j<h; j++)
  { register int thresh = (WHITE/2 + RAND (129) - 64);

    gray = input->bm[j*rowlen + (w-1)] + error;
    den = gray > thresh ? WHITE : BLACK;
    error = gray - den;
    output->bm[j*outrow + (w-1)] = den;
  }

  /* The top border */
  error = 0;
  for (i=0; i<w; i++)
  { register int inp = input->bm[i], thresh = (WHITE/2 + RAND (129) - 64);

    gray = inp + error;
    den = gray > thresh ? WHITE : BLACK;
    error = gray - den;
    output->bm[i] = den;
    eerr[i] = INITERR (inp, den);
  }

  /*
   * Now process the interior bits
   *
   *  Weights:			1+n1  5+n5  3-n1
   *				7-n5    *
   *
   * n1 and n5 are random noise from -0.5 to 0.5 and -2.5 to 2.5
   */
   
  smrange = 2000 * noiselevel/100.0; smrange += 1;
  lgrange = 10000 * noiselevel/100.0; lgrange += 1;
  
# ifdef DEBUG
  fprintf (stderr, "Blue noise level %6.2lf (small %d..%d, large %d..%d)\n",
	  noiselevel, -(smrange/2), smrange/2, -(lgrange/2), lgrange/2);
# endif

  for (j=1; j<h; j++)
  { bmp = &input->bm[j*rowlen];
    obm = &output->bm[j*outrow];

    if (j&1)				/* Odd rows */
    { oerr[0] = INITERR (bmp[0], obm[0]);

      for (i=1; i<w-1; i++)
      { /* Set random weights */
	w1 = 1000; w3 = 3000; w5 = 5000; w7 = 7000;
  
	smnoise = RAND (smrange) - smrange/2;
	w1 += smnoise; w3 -= smnoise;
  
	lgnoise = RAND (lgrange) - lgrange/2;
	w5 += lgnoise; w7 -= lgnoise;
  
	error =  (w1 * eerr[i-1] +
		  w5 * eerr[i] +
		  w3 * eerr[i+1] +
		  w7 * oerr[i-1]) / 16000;
	gray = bmp[i] + error;
  
# ifdef DEBUG
	if (j>10 && j<14 && i>40 && i<44)
	{ fprintf (stderr,
		   "\n<%3d,%3d> input %d, error %d\n",
		   i, j, bmp[i], error);
	  fprintf (stderr, "Noise {%d,%d} Weights {%d,%d,%d,%d}\n",
		   smnoise, lgnoise, w1, w3, w5, w7);
	  fprintf (stderr,
		   "Errs:\t%5d  %5d  %5d\n\t%5d      *\n",
		   eerr[i-1], eerr[i], eerr[i+1], oerr[i-1]);
	}
# endif

	if (gray > (WHITE/2))
	{ obm[i] = 1;	oerr[i] = gray - WHITE; }
	else
	{ obm[i] = 0;	oerr[i] = gray; }
      }
      
      /* Set errors for ends of this row */
      oerr[0]   = INITERR (bmp[0], obm[0]);
      oerr[w-1] = INITERR (bmp[w-1], obm[w-1]);
    }
    else
    { eerr[w-1] = INITERR (bmp[w-1], obm[w-1]);

      for (i=w-2; i>0; i--)
      { /* Set random weights */
	w1 = 1000; w3 = 3000; w5 = 5000; w7 = 7000;
  
	smnoise = RAND (smrange) - smrange/2;
	w1 += smnoise; w3 -= smnoise;
  
	lgnoise = RAND (lgrange) - lgrange/2;
	w5 += lgnoise; w7 -= lgnoise;
  
	error =  (w1 * oerr[i+1] +
		  w5 * oerr[i] +
		  w3 * oerr[i-1] +
		  w7 * eerr[i+1]) / 16000;
	gray = bmp[i] + error;
  
	if (gray > (WHITE/2))
	{ obm[i] = 1;	eerr[i] = gray - WHITE; }
	else
	{ obm[i] = 0;	eerr[i] = gray; }
      }

      /* Set errors for ends of this row */
      eerr[0]   = INITERR (bmp[0], obm[0]);
      eerr[w-1] = INITERR (bmp[w-1], obm[w-1]);
    }
  }
  return (1);
}