/* Copyright (C) 1995 Aladdin Enterprises.  All rights reserved.
  
  This file is part of Aladdin Ghostscript.
  
  Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  or distributor accepts any responsibility for the consequences of using it,
  or for whether it serves any particular purpose or works at all, unless he
  or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  License (the "License") for full details.
  
  Every copy of Aladdin Ghostscript must include a copy of the License,
  normally in a plain ASCII text file named PUBLIC.  The License grants you
  the right to copy, modify and redistribute Aladdin Ghostscript, but only
  under certain conditions described in the License.  Among other things, the
  License requires that the copyright notice and this notice be preserved on
  all copies.
*/

/* siscale.c */
/* Image scaling filters */
#include "math_.h"
#include "memory_.h"
#include "stdio_.h"
#include "strimpl.h"
#include "siscale.h"

/*
 *	Image scaling code is based on public domain code from
 *	Graphics Gems III (pp. 414-424), Academic Press, 1992.
 */

/* ---------------- ImageScaleEncode/Decode ---------------- */

public_st_IScale_state();	/* public so clients can allocate */

#define ss ((stream_IScale_state *)st)

/* ------ Digital filter definition ------ */

/* Mitchell filter definition */
#define Mitchell_support 2.0
#define B (1.0 / 3.0)
#define C (1.0 / 3.0)
private double
Mitchell_filter(double t)
{	double t2 = t * t;

	if ( t < 0 )
	  t = -t;

	if ( t < 1 )
	  return
	    ((12 - 9 * B - 6 * C) * (t * t2) +
	     (-18 + 12 * B + 6 * C) * t2 +
	     (6 - 2 * B)) / 6;
	else if ( t < 2 )
	  return
	    ((-1 * B - 6 * C) * (t * t2) +
	     (6 * B + 30 * C) * t2 +
	     (-12 * B - 48 * C) * t +
	     (8 * B + 24 * C)) / 6;
	else
	  return 0;
}

#define filter_support Mitchell_support
#define filter_proc Mitchell_filter
#define fproc(t) filter_proc(t)
#define fWidthIn filter_support

/*
 * The environment provides the following definitions:
 *	typedef PixelTmp
 *	double fproc(double t)
 *	double fWidthIn
 *	PixelTmp {min,max,unit}PixelTmp
 */
#define CLAMP(v, mn, mx)\
  (v < mn ? mn : v > mx ? mx : v)

/* ------ Auxiliary procedures ------ */

/* Pre-calculate filter contributions for a row or a column. */
private int
contrib_pixels(double scale)
{	return (int)(fWidthIn / min(scale, 1.0) * 2 + 1);
}
private void
calculate_contrib(CLIST *contrib, CONTRIB *items, double scale,
  int size, int limit, int modulus, int stride,
  double rescale_factor)
{
	double WidthIn, fscale;
	bool squeeze;
	int npixels;
	int i, j;

	if ( scale < 1.0 )
	  {	WidthIn = fWidthIn / scale;
		fscale = 1.0 / scale;
		squeeze = true;
	  }
	else
	  {	WidthIn = fWidthIn;
		fscale = 1.0;
		squeeze = false;
	  }
	npixels = (int)(WidthIn * 2 + 1);

	for ( i = 0; i < size; ++i )
	  {	double center = i / scale;
		int left = (int)ceil(center - WidthIn);
		int right = (int)floor(center + WidthIn);
		int max_n = -1;
		contrib[i].n = 0;
		contrib[i].p = items + i * npixels;
		if ( squeeze )
		  {	for ( j = left; j <= right; ++j )
			  {	double weight =
				  fproc((center - j) / fscale) / fscale;
				int n =
				  (j < 0 ? -j :
				   j >= limit ? (limit - j) + limit - 1 :
				   j);
				int k = contrib[i].n++;
				if ( n > max_n )
				  max_n = n;
				contrib[i].p[k].pixel = (n % modulus) * stride;
				contrib[i].p[k].weight =
				  weight * rescale_factor;
			  }
		  }
		else
		  {	for ( j = left; j <= right; ++j )
			  {	double weight = fproc(center - j);
				int n =
				  (j < 0 ? -j :
				   j >= limit ? (limit - j) + limit - 1 :
				   j);
				int k = contrib[i].n++;
				if ( n > max_n )
				  max_n = n;
				contrib[i].p[k].pixel = (n % modulus) * stride;
				contrib[i].p[k].weight =
				  weight * rescale_factor;
			  }
		  }
		contrib[i].max_index = max_n;
	  }
}


/* Apply filter to zoom horizontally from src to tmp. */
private void
zoom_x(PixelTmp *tmp, const void /*PixelIn*/ *src, int sizeofPixelIn,
  int tmp_width, int WidthIn, int Colors, const CLIST *contrib)
{	int c, i, j;
	for ( c = 0; c < Colors; ++c )
	  {
#define zoom_x_loop(PixelIn)\
		const PixelIn *raster = &((PixelIn *)src)[c];\
		for ( i = 0; i < tmp_width; ++i )\
		  {	double weight = 0.0;\
			for ( j = 0; j < contrib[i].n; ++j )\
				weight +=\
				  raster[contrib[i].p[j].pixel]\
				    * contrib[i].p[j].weight;\
			tmp[i * Colors + c] =\
			  (PixelTmp)CLAMP(weight, minPixelTmp, maxPixelTmp);\
		  }

		if ( sizeofPixelIn == 1 )
		  {	zoom_x_loop(byte)
		  }
		else			/* sizeofPixelIn == 2 */
		  {	zoom_x_loop(bits16)
		  }
	  }
}

/* Apply filter to zoom vertically from tmp to dst. */
/* This is simpler because we can treat all columns identically */
/* without regard to the number of samples per pixel. */
private void
zoom_y(void /*PixelOut*/ *dst, int sizeofPixelOut, uint MaxValueOut,
  const PixelTmp *tmp, int WidthOut, int tmp_width,
  int Colors, const CLIST *contrib)
{	int kn = WidthOut * Colors;
	int cn = contrib->n;
	const CONTRIB *cp = contrib->p;
	int kc;
	double max_weight = (double)MaxValueOut;

#define zoom_y_loop(PixelOut)\
	for ( kc = 0; kc < kn; ++kc )\
	  {	const PixelTmp *raster = &tmp[kc];\
		double weight = 0.0;\
		int j;\
		for ( j = 0; j < cn; ++j )\
		  weight += raster[cp[j].pixel] * cp[j].weight;\
		((PixelOut *)dst)[kc] =\
		  (PixelOut)CLAMP(weight, 0, max_weight);\
	  }

	if ( sizeofPixelOut == 1 )
	  {	zoom_y_loop(byte)
	  }
	else			/* sizeofPixelOut == 2 */
	  {	zoom_y_loop(bits16)
	  }
}

/* ------ Stream implementation ------ */

#define tmp_width WidthOut
#define tmp_height HeightIn

/* Forward references */
private void s_IScale_release(P1(stream_state *st));

/* Initialize the filter. */
private int
s_IScale_init(stream_state *st)
{	gs_memory_t *mem = ss->memory;

	ss->sizeofPixelIn = ss->BitsPerComponentIn / 8;
	ss->sizeofPixelOut = ss->BitsPerComponentOut / 8;
	ss->xscale = (double)ss->WidthOut / (double)ss->WidthIn;
	ss->yscale = (double)ss->HeightOut / (double)ss->HeightIn;

	ss->src_y = 0;
	ss->src_size = ss->WidthIn * ss->sizeofPixelIn * ss->Colors;
	ss->src_offset = 0;
	ss->dst_y = -1;
	ss->dst_size = ss->WidthOut * ss->sizeofPixelOut * ss->Colors;
	ss->dst_offset = ss->dst_size;

	/* create intermediate image to hold horizontal zoom */
	ss->tmp = (PixelTmp *)gs_alloc_byte_array(mem,
				ss->tmp_height,
				ss->tmp_width * ss->Colors * sizeof(PixelTmp),
				"image_scale tmp");
	ss->contrib = (CLIST *)gs_alloc_byte_array(mem,
				max(ss->WidthOut, ss->HeightOut),
				sizeof(CLIST), "image_scale contrib");
	ss->items = (CONTRIB *)gs_alloc_byte_array(mem,
				max(contrib_pixels(ss->xscale) * ss->WidthOut,
				    contrib_pixels(ss->yscale) * ss->HeightOut),
				sizeof(CONTRIB), "image_scale contrib[*]");
	/* Allocate buffers for 1 row of source and destination. */
	ss->dst = gs_alloc_byte_array(mem, ss->WidthOut * ss->Colors,
				ss->sizeofPixelOut, "image_scale dst");
	ss->src = gs_alloc_byte_array(mem, ss->WidthIn * ss->Colors,
				ss->sizeofPixelIn, "image_scale src");
	if ( ss->tmp == 0 || ss->contrib == 0 || ss->items == 0 ||
	     ss->dst == 0 || ss->src == 0
	   )
	  {	s_IScale_release(st);
		return ERRC;			/****** WRONG ******/
	  }

	/* Pre-calculate filter contributions for a row. */
	calculate_contrib(ss->contrib, ss->items, ss->xscale,
			  ss->WidthOut, ss->WidthIn, ss->WidthIn, ss->Colors,
			  (double)unitPixelTmp / ss->MaxValueIn);

	return 0;

}

/* Process a buffer.  Note that this handles both Encode and Decode. */
private int
s_IScale_process(stream_state *st, stream_cursor_read *pr,
  stream_cursor_write *pw, bool last)
{
	/* Read input data and scale horizontally into tmp. */

	while ( ss->src_y < ss->HeightIn )
	  { uint rleft = pr->limit - pr->ptr;
	    uint rcount = ss->src_size - ss->src_offset;
	    if ( rleft == 0 )
	      return 0;			/* need more input */
	    if ( rleft >= rcount )
	      { /* We're going to fill up a row. */
		const byte *row;
		if ( ss->src_offset == 0 )
		  { /* We have a complete row.  Read the data */
		    /* directly from the input. */
		    row = pr->ptr + 1;
		  }
		else
		  { /* We're buffering a row in src. */
		    row = ss->src;
		    memcpy((byte *)ss->src + ss->src_offset, pr->ptr + 1,
			   rcount);
		    ss->src_offset = 0;
		  }
		/* Apply filter to zoom horizontally from src to tmp. */
		zoom_x(ss->tmp + ss->src_y * ss->tmp_width * ss->Colors, row,
		       ss->sizeofPixelIn, ss->tmp_width, ss->WidthIn,
		       ss->Colors, ss->contrib);
		pr->ptr += rcount;
		if ( ++(ss->src_y) == ss->HeightIn )
		  { /* We've read all the input.  Prepare for output. */
		    /* Pre-calculate filter contributions for a column. */
		    calculate_contrib(ss->contrib, ss->items, ss->yscale,
				      ss->HeightOut, ss->HeightIn,
				      ss->tmp_height, ss->WidthOut * ss->Colors,
				      (double)ss->MaxValueOut / unitPixelTmp);
		  }
	      }
	    else
	      { /* We don't have a complete row.  Copy data to src buffer. */
		memcpy((byte *)ss->src + ss->src_offset, pr->ptr + 1, rleft);
		ss->src_offset += rleft;
		pr->ptr += rleft;
		break;
	      }
	  }

	/* Scale tmp vertically and write into output. */

	while ( ss->dst_y < ss->HeightOut )
	  { uint wleft = pw->limit - pw->ptr;
	    uint wcount = ss->dst_size - ss->dst_offset;
	    if ( wcount != 0 )
	      { /* We're buffering an output row. */
		uint ncopy = min(wleft, wcount);
		if ( wleft == 0 )
		  return 1;			/* need more output space */
		memcpy(pw->ptr + 1, (byte *)ss->dst + ss->dst_offset, ncopy);
		pw->ptr += ncopy;
		ss->dst_offset += ncopy;
	      }
	    else
	      { /* We're going to fill up an output row. */
		byte *row;
		if ( ++(ss->dst_y) == ss->HeightOut )
		  break;
		if ( wleft >= ss->dst_size )
		  { /* We can scale the row directly into the output. */
		    row = pw->ptr + 1;
		    pw->ptr += ss->dst_size;
		  }
		else if ( wleft == 0 )
		  { --(ss->dst_y);
		    return 1;			/* need more output space */
		  }
		else
		  { /* We'll have to buffer the row. */
		    row = ss->dst;
		    ss->dst_offset = 0;
		  }
		/* Apply filter to zoom vertically from tmp to dst. */
		zoom_y(row, ss->sizeofPixelOut, ss->MaxValueOut, ss->tmp,
		       ss->WidthOut, ss->tmp_width, ss->Colors,
		       ss->contrib + ss->dst_y);
	      }
	  }

	return EOFC;
}

/* Release the filter's storage. */
private void
s_IScale_release(stream_state *st)
{	gs_memory_t *mem = ss->memory;
	gs_free_object(mem, (void *)ss->src, "image_scale src"); /* no longer const */
	  ss->src = 0;
	gs_free_object(mem, ss->dst, "image_scale dst");
	  ss->dst = 0;
	gs_free_object(mem, ss->items, "image_scale contrib[*]");
	  ss->items = 0;
	gs_free_object(mem, ss->contrib, "image_scale contrib");
	  ss->contrib = 0;
	gs_free_object(mem, ss->tmp, "image_scale tmp");
	  ss->tmp = 0;
}

#undef ss

/* Stream template */
const stream_template s_IScale_template =
{	&st_IScale_state, s_IScale_init, s_IScale_process, 1, 1,
	s_IScale_release
};