[imager.git] / scale.im

#include "imager.h"
#include "imageri.h"

/*
 * i_scale_mixing() is based on code contained in pnmscale.c, part of
 * the netpbm distribution.  No code was copied from pnmscale but
 * the algorthm was and for this I thank the netpbm crew.
 *
 * Tony
 */

/* pnmscale.c - read a portable anymap and scale it
**
** Copyright (C) 1989, 1991 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
**
*/


static void
zero_row(i_fcolor *row, int width, int channels);

#code
static void
IM_SUFFIX(accum_output_row)(i_fcolor *accum, double fraction, IM_COLOR const *in,
		 int width, int channels);
static void
IM_SUFFIX(horizontal_scale)(IM_COLOR *out, int out_width, 
                            i_fcolor const *in, int in_width,
                            int channels);
#/code

/*
=item i_scale_mixing

Returns a new image scaled to the given size.

Unlike i_scale_axis() this does a simple coverage of pixels from
source to target and doesn't resample.

Adapted from pnmscale.

=cut
*/
i_img *
i_scale_mixing(i_img *src, int x_out, int y_out) {
  i_img *result;
  i_fcolor *accum_row = NULL;
  int x, y, ch;
  int accum_row_bytes;
  double rowsleft, fracrowtofill;
  int rowsread;
  double y_scale;

  mm_log((1, "i_scale_mixing(src %p, x_out %d, y_out %d)\n", 
	  src, x_out, y_out));

  i_clear_error();

  if (x_out <= 0) {
    i_push_errorf(0, "output width %d invalid", x_out);
    return NULL;
  }
  if (y_out <= 0) {
    i_push_errorf(0, "output height %d invalid", y_out);
    return NULL;
  }

  if (x_out == src->xsize && y_out == src->ysize) {
    return i_copy(src);
  }

  y_scale = y_out / (double)src->ysize;

  result = i_sametype_chans(src, x_out, y_out, src->channels);
  if (!result)
    return NULL;

  accum_row_bytes = sizeof(i_fcolor) * src->xsize;
  if (accum_row_bytes / sizeof(i_fcolor) != src->xsize) {
    i_push_error(0, "integer overflow allocating accumulator row buffer");
    return NULL;
  }

  accum_row  = mymalloc(accum_row_bytes);

#code src->bits <= 8
  IM_COLOR *in_row = NULL;
  IM_COLOR *xscale_row = NULL;
  int in_row_bytes, out_row_bytes;

  in_row_bytes = sizeof(IM_COLOR) * src->xsize;
  if (in_row_bytes / sizeof(IM_COLOR) != src->xsize) {
    i_push_error(0, "integer overflow allocating input row buffer");
    return NULL;
  }
  out_row_bytes = sizeof(IM_COLOR) * x_out;
  if (out_row_bytes / sizeof(IM_COLOR) != x_out) {
    i_push_error(0, "integer overflow allocating output row buffer");
    return NULL;
  }

  in_row     = mymalloc(in_row_bytes);
  xscale_row = mymalloc(out_row_bytes);

  rowsread = 0;
  rowsleft = 0.0;
  for (y = 0; y < y_out; ++y) {
    if (y_out == src->ysize) {
      /* no vertical scaling, just load it */
#ifdef IM_EIGHT_BIT
      int x, ch;
      /* load and convert to doubles */
      IM_GLIN(src, 0, src->xsize, y, in_row);
      for (x = 0; x < src->xsize; ++x) {
        for (ch = 0; ch < src->channels; ++ch) {
          accum_row[x].channel[ch] = in_row[x].channel[ch];
        }
      }
#else
      IM_GLIN(src, 0, src->xsize, y, accum_row);
#endif
      /* alpha adjust if needed */
      if (src->channels == 2 || src->channels == 4) {
	for (x = 0; x < src->xsize; ++x) {
	  for (ch = 0; ch < src->channels-1; ++ch) {
	    accum_row[x].channel[ch] *=
	      accum_row[x].channel[src->channels-1] / IM_SAMPLE_MAX;
	  }
	}
      }
    }
    else {
      fracrowtofill = 1.0;
      zero_row(accum_row, src->xsize, src->channels);
      while (fracrowtofill > 0) {
	if (rowsleft <= 0) {
	  if (rowsread < src->ysize) {
	    IM_GLIN(src, 0, src->xsize, rowsread, in_row);
	    ++rowsread;
	  }
	  /* else just use the last row read */

	  rowsleft = y_scale;
	}
	if (rowsleft < fracrowtofill) {
	  IM_SUFFIX(accum_output_row)(accum_row, rowsleft, in_row, 
                                      src->xsize, src->channels);
	  fracrowtofill -= rowsleft;
	  rowsleft = 0;
	}
	else {
	  IM_SUFFIX(accum_output_row)(accum_row, fracrowtofill, in_row, 
                                      src->xsize, src->channels);
	  rowsleft -= fracrowtofill;
	  fracrowtofill = 0;
	}
      }
    }
    /* we've accumulated a vertically scaled row */
    if (x_out == src->xsize) {
#if IM_EIGHT_BIT
      int x, ch;
      /* no need to scale, but we need to convert it */
      if (result->channels == 2 || result->channels == 4) {
	int alpha_chan = result->channels - 1;
	for (x = 0; x < x_out; ++x) {
	  double alpha = accum_row[x].channel[alpha_chan] / IM_SAMPLE_MAX;
	  if (alpha) {
	    for (ch = 0; ch < alpha_chan; ++ch) {
	      int val = accum_row[x].channel[ch] / alpha + 0.5;
	      xscale_row[x].channel[ch] = IM_LIMIT(val);
	    }
	  }
	  xscale_row[x].channel[alpha_chan] = IM_LIMIT(accum_row[x].channel[alpha_chan]+0.5);
	}
      }
      else {
	for (x = 0; x < x_out; ++x) {
	  for (ch = 0; ch < result->channels; ++ch)
	    xscale_row[x].channel[ch] = IM_LIMIT(accum_row[x].channel[ch]+0.5);
	}
      }
      IM_PLIN(result, 0, x_out, y, xscale_row);
#else
      IM_PLIN(result, 0, x_out, y, accum_row);
#endif
    }
    else {
      IM_SUFFIX(horizontal_scale)(xscale_row, x_out, accum_row, 
                                  src->xsize, src->channels);
      IM_PLIN(result, 0, x_out, y, xscale_row);
    }
  }
  myfree(in_row);
  myfree(xscale_row);
#/code
  myfree(accum_row);

  return result;
}

static void
zero_row(i_fcolor *row, int width, int channels) {
  int x;
  int ch;

  /* with IEEE floats we could just use memset() but that's not
     safe in general under ANSI C.
     memset() is slightly faster.
  */
  for (x = 0; x < width; ++x) {
    for (ch = 0; ch < channels; ++ch)
      row[x].channel[ch] = 0.0;
  }
}

#code

static void
IM_SUFFIX(accum_output_row)(i_fcolor *accum, double fraction, IM_COLOR const *in,
		 int width, int channels) {
  int x, ch;

  /* it's tempting to change this into a pointer iteration loop but
     modern CPUs do the indexing as part of the instruction */
  if (channels == 2 || channels == 4) {
    for (x = 0; x < width; ++x) {
      for (ch = 0; ch < channels-1; ++ch) {
	accum[x].channel[ch] += in[x].channel[ch] * fraction * in[x].channel[channels-1] / IM_SAMPLE_MAX;
      }
      accum[x].channel[channels-1] += in[x].channel[channels-1] * fraction;
    }
  }
  else {
    for (x = 0; x < width; ++x) {
      for (ch = 0; ch < channels; ++ch) {
	accum[x].channel[ch] += in[x].channel[ch] * fraction;
      }
    }
  }
}

static void
IM_SUFFIX(horizontal_scale)(IM_COLOR *out, int out_width, 
		 i_fcolor const *in, int in_width,
		 int channels) {
  double frac_col_to_fill, frac_col_left;
  int in_x;
  int out_x;
  double x_scale = (double)out_width / in_width;
  int ch;
  double accum[MAXCHANNELS] = { 0 };
  
  frac_col_to_fill = 1.0;
  out_x = 0;
  for (in_x = 0; in_x < in_width; ++in_x) {
    frac_col_left = x_scale;
    while (frac_col_left >= frac_col_to_fill) {
      for (ch = 0; ch < channels; ++ch)
	accum[ch] += frac_col_to_fill * in[in_x].channel[ch];

      if (channels == 2 || channels == 4) {
	int alpha_chan = channels - 1;
	double alpha = accum[alpha_chan] / IM_SAMPLE_MAX;
	if (alpha) {
	  for (ch = 0; ch < alpha_chan; ++ch) {
	    IM_WORK_T val = IM_ROUND(accum[ch] / alpha);
	    out[out_x].channel[ch] = IM_LIMIT(val);
	  }
	}
	out[out_x].channel[alpha_chan] = IM_LIMIT(IM_ROUND(accum[alpha_chan]));
      }
      else {
	for (ch = 0; ch < channels; ++ch) {
	  IM_WORK_T val = IM_ROUND(accum[ch]);
	  out[out_x].channel[ch] = IM_LIMIT(val);
	}
      }
      for (ch = 0; ch < channels; ++ch)
	accum[ch] = 0;
      frac_col_left -= frac_col_to_fill;
      frac_col_to_fill = 1.0;
      ++out_x;
    }

    if (frac_col_left > 0) {
      for (ch = 0; ch < channels; ++ch) {
	accum[ch] += frac_col_left * in[in_x].channel[ch];
      }
      frac_col_to_fill -= frac_col_left;
    }
  }

  if (out_x < out_width-1 || out_x > out_width) {
    i_fatal(3, "Internal error: out_x %d out of range (width %d)", out_x, out_width);
  }
  
  if (out_x < out_width) {
    for (ch = 0; ch < channels; ++ch) {
      accum[ch] += frac_col_to_fill * in[in_width-1].channel[ch];
    }
    if (channels == 2 || channels == 4) {
      int alpha_chan = channels - 1;
      double alpha = accum[alpha_chan] / IM_SAMPLE_MAX;
      if (alpha) {
	for (ch = 0; ch < alpha_chan; ++ch) {
	  IM_WORK_T val = IM_ROUND(accum[ch] / alpha);
	  out[out_x].channel[ch] = IM_LIMIT(val);
	}
      }
      out[out_x].channel[alpha_chan] = IM_LIMIT(IM_ROUND(accum[alpha_chan]));
    }
    else {
      for (ch = 0; ch < channels; ++ch) {
	IM_WORK_T val = IM_ROUND(accum[ch]);
	out[out_x].channel[ch] = IM_LIMIT(val);
      }
    }
  }
}

#/code
Commit	Line	Data
658f724e	1	#include "imager.h"
874c55db	2	#include "imageri.h"
658f724e TC	3
	4	/*
	5	* i_scale_mixing() is based on code contained in pnmscale.c, part of
	6	* the netpbm distribution. No code was copied from pnmscale but
	7	* the algorthm was and for this I thank the netpbm crew.
	8	*
	9	* Tony
	10	*/
	11
	12	/* pnmscale.c - read a portable anymap and scale it
	13	**
	14	** Copyright (C) 1989, 1991 by Jef Poskanzer.
	15	**
	16	** Permission to use, copy, modify, and distribute this software and its
	17	** documentation for any purpose and without fee is hereby granted, provided
	18	** that the above copyright notice appear in all copies and that both that
	19	** copyright notice and this permission notice appear in supporting
	20	** documentation. This software is provided "as is" without express or
	21	** implied warranty.
	22	**
	23	*/
	24
	25
	26	static void
	27	zero_row(i_fcolor *row, int width, int channels);
a10945af TC	28
a10945af TC	29	#code
658f724e	30	static void
a10945af	31	IM_SUFFIX(accum_output_row)(i_fcolor accum, double fraction, IM_COLOR const in,
658f724e TC	32	int width, int channels);
658f724e TC	33	static void
a10945af TC	34	IM_SUFFIX(horizontal_scale)(IM_COLOR *out, int out_width,
	35	i_fcolor const *in, int in_width,
	36	int channels);
	37	#/code
658f724e TC	38
	39	/*
	40	=item i_scale_mixing
	41
	42	Returns a new image scaled to the given size.
	43
	44	Unlike i_scale_axis() this does a simple coverage of pixels from
	45	source to target and doesn't resample.
	46
	47	Adapted from pnmscale.
	48
	49	=cut
	50	*/
	51	i_img *
	52	i_scale_mixing(i_img *src, int x_out, int y_out) {
	53	i_img *result;
658f724e	54	i_fcolor *accum_row = NULL;
874c55db	55	int x, y, ch;
a10945af	56	int accum_row_bytes;
658f724e TC	57	double rowsleft, fracrowtofill;
	58	int rowsread;
	59	double y_scale;
	60
	61	mm_log((1, "i_scale_mixing(src %p, x_out %d, y_out %d)\n",
	62	src, x_out, y_out));
	63
	64	i_clear_error();
	65
	66	if (x_out <= 0) {
	67	i_push_errorf(0, "output width %d invalid", x_out);
	68	return NULL;
	69	}
	70	if (y_out <= 0) {
	71	i_push_errorf(0, "output height %d invalid", y_out);
	72	return NULL;
	73	}
	74
658f724e TC	75	if (x_out == src->xsize && y_out == src->ysize) {
	76	return i_copy(src);
	77	}
	78
	79	y_scale = y_out / (double)src->ysize;
	80
	81	result = i_sametype_chans(src, x_out, y_out, src->channels);
	82	if (!result)
	83	return NULL;
	84
a10945af TC	85	accum_row_bytes = sizeof(i_fcolor) * src->xsize;
	86	if (accum_row_bytes / sizeof(i_fcolor) != src->xsize) {
	87	i_push_error(0, "integer overflow allocating accumulator row buffer");
	88	return NULL;
	89	}
	90
	91	accum_row = mymalloc(accum_row_bytes);
	92
	93	#code src->bits <= 8
	94	IM_COLOR *in_row = NULL;
	95	IM_COLOR *xscale_row = NULL;
	96	int in_row_bytes, out_row_bytes;
	97
	98	in_row_bytes = sizeof(IM_COLOR) * src->xsize;
	99	if (in_row_bytes / sizeof(IM_COLOR) != src->xsize) {
	100	i_push_error(0, "integer overflow allocating input row buffer");
	101	return NULL;
	102	}
	103	out_row_bytes = sizeof(IM_COLOR) * x_out;
	104	if (out_row_bytes / sizeof(IM_COLOR) != x_out) {
	105	i_push_error(0, "integer overflow allocating output row buffer");
	106	return NULL;
	107	}
	108
658f724e	109	in_row = mymalloc(in_row_bytes);
658f724e TC	110	xscale_row = mymalloc(out_row_bytes);
	111
	112	rowsread = 0;
	113	rowsleft = 0.0;
	114	for (y = 0; y < y_out; ++y) {
	115	if (y_out == src->ysize) {
a10945af	116	/* no vertical scaling, just load it */
a10945af	117	#ifdef IM_EIGHT_BIT
e4bf9335	118	int x, ch;
a10945af TC	119	/* load and convert to doubles */
	120	IM_GLIN(src, 0, src->xsize, y, in_row);
	121	for (x = 0; x < src->xsize; ++x) {
	122	for (ch = 0; ch < src->channels; ++ch) {
	123	accum_row[x].channel[ch] = in_row[x].channel[ch];
	124	}
	125	}
	126	#else
	127	IM_GLIN(src, 0, src->xsize, y, accum_row);
	128	#endif
874c55db TC	129	/* alpha adjust if needed */
	130	if (src->channels == 2 \|\| src->channels == 4) {
	131	for (x = 0; x < src->xsize; ++x) {
	132	for (ch = 0; ch < src->channels-1; ++ch) {
	133	accum_row[x].channel[ch] *=
	134	accum_row[x].channel[src->channels-1] / IM_SAMPLE_MAX;
	135	}
	136	}
	137	}
658f724e TC	138	}
	139	else {
	140	fracrowtofill = 1.0;
	141	zero_row(accum_row, src->xsize, src->channels);
	142	while (fracrowtofill > 0) {
	143	if (rowsleft <= 0) {
	144	if (rowsread < src->ysize) {
a10945af	145	IM_GLIN(src, 0, src->xsize, rowsread, in_row);
658f724e TC	146	++rowsread;
	147	}
	148	/* else just use the last row read */
	149
	150	rowsleft = y_scale;
	151	}
	152	if (rowsleft < fracrowtofill) {
a10945af TC	153	IM_SUFFIX(accum_output_row)(accum_row, rowsleft, in_row,
a10945af TC	154	src->xsize, src->channels);
658f724e TC	155	fracrowtofill -= rowsleft;
	156	rowsleft = 0;
	157	}
	158	else {
a10945af TC	159	IM_SUFFIX(accum_output_row)(accum_row, fracrowtofill, in_row,
a10945af TC	160	src->xsize, src->channels);
658f724e TC	161	rowsleft -= fracrowtofill;
	162	fracrowtofill = 0;
	163	}
	164	}
a10945af TC	165	}
	166	/* we've accumulated a vertically scaled row */
	167	if (x_out == src->xsize) {
a10945af	168	#if IM_EIGHT_BIT
e4bf9335	169	int x, ch;
a10945af	170	/* no need to scale, but we need to convert it */
874c55db TC	171	if (result->channels == 2 \|\| result->channels == 4) {
	172	int alpha_chan = result->channels - 1;
	173	for (x = 0; x < x_out; ++x) {
	174	double alpha = accum_row[x].channel[alpha_chan] / IM_SAMPLE_MAX;
	175	if (alpha) {
	176	for (ch = 0; ch < alpha_chan; ++ch) {
	177	int val = accum_row[x].channel[ch] / alpha + 0.5;
	178	xscale_row[x].channel[ch] = IM_LIMIT(val);
	179	}
	180	}
	181	xscale_row[x].channel[alpha_chan] = IM_LIMIT(accum_row[x].channel[alpha_chan]+0.5);
	182	}
	183	}
	184	else {
	185	for (x = 0; x < x_out; ++x) {
	186	for (ch = 0; ch < result->channels; ++ch)
	187	xscale_row[x].channel[ch] = IM_LIMIT(accum_row[x].channel[ch]+0.5);
	188	}
658f724e	189	}
a10945af TC	190	IM_PLIN(result, 0, x_out, y, xscale_row);
	191	#else
	192	IM_PLIN(result, 0, x_out, y, accum_row);
	193	#endif
	194	}
	195	else {
	196	IM_SUFFIX(horizontal_scale)(xscale_row, x_out, accum_row,
	197	src->xsize, src->channels);
	198	IM_PLIN(result, 0, x_out, y, xscale_row);
658f724e TC	199	}
658f724e TC	200	}
658f724e	201	myfree(in_row);
658f724e	202	myfree(xscale_row);
a10945af TC	203	#/code
a10945af TC	204	myfree(accum_row);
658f724e TC	205
	206	return result;
	207	}
	208
	209	static void
	210	zero_row(i_fcolor *row, int width, int channels) {
	211	int x;
	212	int ch;
	213
	214	/* with IEEE floats we could just use memset() but that's not
a10945af TC	215	safe in general under ANSI C.
	216	memset() is slightly faster.
	217	*/
658f724e TC	218	for (x = 0; x < width; ++x) {
	219	for (ch = 0; ch < channels; ++ch)
	220	row[x].channel[ch] = 0.0;
	221	}
	222	}
	223
a10945af TC	224	#code
a10945af TC	225
658f724e	226	static void
a10945af	227	IM_SUFFIX(accum_output_row)(i_fcolor accum, double fraction, IM_COLOR const in,
658f724e TC	228	int width, int channels) {
	229	int x, ch;
	230
a10945af TC	231	/* it's tempting to change this into a pointer iteration loop but
a10945af TC	232	modern CPUs do the indexing as part of the instruction */
874c55db TC	233	if (channels == 2 \|\| channels == 4) {
	234	for (x = 0; x < width; ++x) {
	235	for (ch = 0; ch < channels-1; ++ch) {
	236	accum[x].channel[ch] += in[x].channel[ch] * fraction * in[x].channel[channels-1] / IM_SAMPLE_MAX;
	237	}
	238	accum[x].channel[channels-1] += in[x].channel[channels-1] * fraction;
	239	}
	240	}
	241	else {
	242	for (x = 0; x < width; ++x) {
	243	for (ch = 0; ch < channels; ++ch) {
	244	accum[x].channel[ch] += in[x].channel[ch] * fraction;
	245	}
658f724e TC	246	}
	247	}
	248	}
	249
	250	static void
a10945af	251	IM_SUFFIX(horizontal_scale)(IM_COLOR *out, int out_width,
658f724e TC	252	i_fcolor const *in, int in_width,
	253	int channels) {
	254	double frac_col_to_fill, frac_col_left;
	255	int in_x;
	256	int out_x;
	257	double x_scale = (double)out_width / in_width;
	258	int ch;
	259	double accum[MAXCHANNELS] = { 0 };
	260
	261	frac_col_to_fill = 1.0;
	262	out_x = 0;
	263	for (in_x = 0; in_x < in_width; ++in_x) {
	264	frac_col_left = x_scale;
	265	while (frac_col_left >= frac_col_to_fill) {
	266	for (ch = 0; ch < channels; ++ch)
	267	accum[ch] += frac_col_to_fill * in[in_x].channel[ch];
	268
874c55db TC	269	if (channels == 2 \|\| channels == 4) {
	270	int alpha_chan = channels - 1;
	271	double alpha = accum[alpha_chan] / IM_SAMPLE_MAX;
	272	if (alpha) {
	273	for (ch = 0; ch < alpha_chan; ++ch) {
	274	IM_WORK_T val = IM_ROUND(accum[ch] / alpha);
	275	out[out_x].channel[ch] = IM_LIMIT(val);
	276	}
	277	}
	278	out[out_x].channel[alpha_chan] = IM_LIMIT(IM_ROUND(accum[alpha_chan]));
658f724e	279	}
874c55db TC	280	else {
	281	for (ch = 0; ch < channels; ++ch) {
	282	IM_WORK_T val = IM_ROUND(accum[ch]);
	283	out[out_x].channel[ch] = IM_LIMIT(val);
	284	}
	285	}
	286	for (ch = 0; ch < channels; ++ch)
	287	accum[ch] = 0;
658f724e TC	288	frac_col_left -= frac_col_to_fill;
	289	frac_col_to_fill = 1.0;
	290	++out_x;
	291	}
	292
	293	if (frac_col_left > 0) {
	294	for (ch = 0; ch < channels; ++ch) {
	295	accum[ch] += frac_col_left * in[in_x].channel[ch];
	296	}
	297	frac_col_to_fill -= frac_col_left;
	298	}
	299	}
	300
	301	if (out_x < out_width-1 \|\| out_x > out_width) {
	302	i_fatal(3, "Internal error: out_x %d out of range (width %d)", out_x, out_width);
	303	}
	304
	305	if (out_x < out_width) {
	306	for (ch = 0; ch < channels; ++ch) {
	307	accum[ch] += frac_col_to_fill * in[in_width-1].channel[ch];
874c55db TC	308	}
	309	if (channels == 2 \|\| channels == 4) {
	310	int alpha_chan = channels - 1;
	311	double alpha = accum[alpha_chan] / IM_SAMPLE_MAX;
	312	if (alpha) {
	313	for (ch = 0; ch < alpha_chan; ++ch) {
	314	IM_WORK_T val = IM_ROUND(accum[ch] / alpha);
	315	out[out_x].channel[ch] = IM_LIMIT(val);
	316	}
	317	}
	318	out[out_x].channel[alpha_chan] = IM_LIMIT(IM_ROUND(accum[alpha_chan]));
	319	}
	320	else {
	321	for (ch = 0; ch < channels; ++ch) {
	322	IM_WORK_T val = IM_ROUND(accum[ch]);
	323	out[out_x].channel[ch] = IM_LIMIT(val);
	324	}
658f724e TC	325	}
	326	}
	327	}
a10945af TC	328
a10945af TC	329	#/code