/**************************************************************************
 *{@C
 *      Copyright:      1988-2020 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         ImageProcessing
 *      Filename:       IP_Pyramid.c
 *
 *      Author:         Paul Obermeier
 *
 *	Description:	Functions to create and sample pixel pyramids.
 *
 *      Additional documentation:
 *                      None.
 *
 *	Exported functions:
 *			IP_GenPyramid
 *			IP_SamplePyramid
 *
 **************************************************************************/

#include <stdio.h>

#include "UT_Compat.h"
#include "UT_Const.h"
#include "UT_Macros.h"
#include "UT_Error.h"
#include "UT_Memory.h"

#include "FF_Image.h"

#include "IP_Image.h"
#include "IP_ImagePrivate.h"

/* Find the resolution for a pixel pyramid. */

static Int32 find_maxlevel (Int32 width, Int32 height)
{
    Int32 size, maxlevel;

    size = UT_MAX (width, height);
    maxlevel = 1;
    while (size > (1 << maxlevel))
	++maxlevel;
    return maxlevel;
}


/* Allocate and initialize memory for a pixel pyramid. */

static pyramid *new_pyramid
	(Int32 width, Int32 height, const FF_ImgFmtType channel[])
{
    pyramid		*pyr;
    struct pyrlevel	*plv;
    Int32		w, h, i, j, k, level, n;

    if (!(pyr = UT_MemTemp (pyramid)))
	return NULL;
    pyr->numchan = 0;
    for (i = 0; i < FF_NumImgChanTypes; ++i)
	if (channel[i])
	    pyr->channel[pyr->numchan++] = i;
    pyr->maxlevel = find_maxlevel (width, height);
    pyr->xmax = (Float32) (width - 1);
    pyr->ymax = (Float32) (height - 1);
    pyr->xyscale1 = 1.0 / (Float32) ((1 << pyr->maxlevel) - 1);
    w = width;
    h = height;
    level = pyr->maxlevel;
    while (level >= 0)
    {
	plv = pyr->pyrlevel + level;
	plv->rowbytes = pyr->numchan * w;
	n = plv->rowbytes * h;
	if (!(plv->pixels = UT_MemTempArray (n, UInt8)))
	    return NULL;
	plv->xyscale2 = (Float32) ((1 << level) - 1);
	plv->width = w;
	plv->height = h;
	--level;
	if (w & 1)
	    w = (w + 1) >> 1;
	else
	    w >>= 1;
	if (w == 0)
	    w = 1;
	if (h & 1)
	    h = (h + 1) >> 1;
	else
	    h >>= 1;
	if (h == 0)
	    h = 1;
    }
    n = (1 << (pyr->maxlevel - 1)) + 1;
    if (!(pyr->levtable = UT_MemTempArray (n, Int32)))
	return NULL;
    pyr->lvtscale = (Float32) ((1 << (pyr->maxlevel - 1)));
    pyr->levtable[0] = pyr->maxlevel;
    i = 1;
    j = 1;
    k = pyr->maxlevel;
    while (i <= (1 << (pyr->maxlevel - 1)))
    {
	if (i == j)
	{
	    j <<= 1;
	    --k;
	}
	pyr->levtable[i] = k;
	++i;
    }
    if (k != 0)
	UT_ErrFatal ("new_pyramid", "level table inconsistent");
    return pyr;
}


/* Read the pixel data for the channels specified in pixel pyramid "pyr"
   from image "img" into "pyr". */

static void top_level (IP_ImageId img, pyramid *pyr)
{
    Int32		x, y, i, chn, itmp;
    struct pyrlevel	*plv;
    Float32		*fsrc;
    UInt8		*ubsrc, *ubdest;

    plv = pyr->pyrlevel + pyr->maxlevel;
    for (i = 0; i < pyr->numchan; ++i)
    {
	chn = pyr->channel[i];
	switch (img->pimg.channel[chn])
	{
	case FF_ImgFmtTypeUByte:
	    for (y = 0; y < plv->height; ++y)
	    {
		ubsrc = FF_ImgUByteRow (&img->pimg, chn, y);
		ubdest = plv->pixels + y * plv->rowbytes + i;
		for (x = 0; x < plv->width; ++x)
		{
		    *ubdest = *ubsrc;
		    ubdest += pyr->numchan;
		    ++ubsrc;
		}
	    }
	    break;
	case FF_ImgFmtTypeFloat:
	    for (y = 0; y < plv->height; ++y)
	    {
		fsrc = FF_ImgFloatRow (&img->pimg, chn, y);
		ubdest = plv->pixels + y * plv->rowbytes + i;
		for (x = 0; x < plv->width; ++x)
		{
		    itmp = (Int32) (255.0 * UT_MAX (0.0, UT_MIN (*fsrc, 1.0)));
		    *ubdest = itmp;
		    ubdest += pyr->numchan;
		    ++fsrc;
		}
	    }
	    break;
	default:
	    UT_ErrFatal ("top_level", "invalid pixel format");
	}
    }
    return;
}


/* Compute the pixel data for the lower-resolution levels in pixel pyramid
"pyr" from level number "maxlevel. */

static void lower_levels (pyramid *pyr)
{
    Int32		i, level,
			xsrc1, xsrc2, xsrcmax, xdest,
			ysrc1, ysrc2, ysrcmax, ydest;
    struct pyrlevel	*srcplv, *destplv;
    UInt8		*src1, *src2, *dest;

    for (level = pyr->maxlevel; level > 0; --level)
    {
	srcplv = pyr->pyrlevel + level;
	destplv = srcplv - 1;
	xsrcmax = srcplv->width - 1;
	ysrcmax = srcplv->height - 1;
	for (ydest = 0; ydest < destplv->height; ++ydest)
	{
	    dest = destplv->pixels + ydest * destplv->rowbytes;
	    ysrc1 = UT_MIN (ydest << 1, ysrcmax);
	    ysrc2 = UT_MIN (ysrc1 + 1, ysrcmax);
	    src1 = srcplv->pixels + ysrc1 * srcplv->rowbytes;
	    src2 = srcplv->pixels + ysrc2 * srcplv->rowbytes;
	    for (xdest = 0; xdest < destplv->width; ++xdest)
	    {
		xsrc1 = UT_MIN (xdest << 1, xsrcmax);
		xsrc2 = UT_MIN (xsrc1 + 1, xsrcmax);
		for (i = 0; i < pyr->numchan; ++i)
		    dest[xdest * pyr->numchan + i] =
			(src1[xsrc1 * pyr->numchan + i] +
			 src1[xsrc2 * pyr->numchan + i] +
			 src2[xsrc1 * pyr->numchan + i] +
			 src2[xsrc2 * pyr->numchan + i]) >> 2;
	    }
	}
    }
    return;
}

/***************************************************************************
 *[@e
 *	Name:		IP_GenPyramid
 *
 *	Usage:		Convert an image into a pixel pyramid.
 *
 *	Synopsis:	pyramid *IP_GenPyramid (IP_ImageId img)
 *
 *	Description:	An IP_ImageId is preprocessed into a "pixel pyramid" so
 *			that an approximation of the average color within a
 *			square with arbitrary size and center can be computed
 *			in constant time.
 *
 *	Return value:	A pointer to the pixel pyramid if successful,
 *			else NULL.
 *
 *			Note: IP_GenPyramid allocates only temporary memory
 *			blocks to store the pyramid, so that the pixel pyramid
 *			can be deleted easily using UT_MemRemember and
 *			UT_MemRestore.
 *
 *	See also:
 *
 ***************************************************************************/

pyramid *IP_GenPyramid (IP_ImageId img)
{
    UT_MemState	memstate;
    pyramid	*pyr;

    memstate = UT_MemRemember ();
    if ((pyr = new_pyramid (img->pimg.desc.width,
				img->pimg.desc.height,
				img->pimg.channel)))
    {
	top_level (img, pyr);
	lower_levels (pyr);
    }
    if (!pyr)
	UT_MemRestore (memstate);
    return pyr;
}

/***************************************************************************
 *[@e
 *	Name:		IP_SamplePyramid
 *
 *	Usage:		Get a color sample from a pixel pyramid.
 *
 *	Synopsis:	void IP_SamplePyramid
 *				(pyramid *pyr, 
 *				 Float32 x, Float32 y,
 *				 Float32 size,
 *				 Float32 cdata[],
 *				 IP_FillModeType oor)
 *
 *	Description:	IP_SamplePyramid returns the average "color" of a
 *			square of size "size" centered around position (x, y)
 *			in pixel pyramid "pyr".
 *			The color data are returned in "cdata". No data are
 *			entered into "cdata" for channels not contained in "p".
 *
 *			The range of values for "x", "y" and size is defined
 *			as follows: The lower left corner of the picture
 *			stored in "pyr" is at position (0.0, 0.0); the upper
 *			left corner is at position (w-1.0, h-1.0), where "w"
 *			and "h" are the width and height of the IP_ImageId from
 *			which "pyr" was constructed.
 *
 *			"Oor" selects how color values for "cdata" are selected
 *			when "x" or "y" are outside the range from 0.0 to 1.0.
 *			The following values for "oor" are acceptable:
 *
 *			IP_FillModeFill	The current drawing color is returned
 *					in "cdata".
 *
 *			IP_FillModeWrap	"X" and "y" are taken modulo "w" and
 *					"h", i.e. the image stored in the
 *					pyramid is repeated in horizontal and
 *					vertical direction.
 *
 *			IP_FillModeClip	Coordinate values less than 0.0 are
 *					set to 0.0; values greater than 1.0
 *					are set to 1.0.
 *
 *	Return value:	None.
 *
 *	See also:
 *
 ***************************************************************************/

void IP_SamplePyramid
	(pyramid *pyr,
	 Float32 x, Float32 y, Float32 size,
	 Float32 cdata[],
	 IP_FillModeType oor)
{
    Int32			i, n0, n1, col0, col1, row0, row1;
    Float32			a, b, tmp, fcol, frow, clipx, clipy,
				u0, u1, v0, v1, size2,
				in_horz, in_vert, in, out;
    const UInt8			*LL, *LR, *UL, *UR;
    const struct pyrlevel	*plv;

    /* Compute the indices, "n0" and "n1" of the two pyramid levels we will
    sample, and the weighting factors, "a" and "b" for the samples taken
    from the levels. */

    if (size < 0.0)
	size = 0.0;
    tmp = size * pyr->xyscale1;
    tmp = UT_MIN (tmp, 1.0);
    n0 = pyr->levtable[(int) (tmp * pyr->lvtscale)];
    n1 = n0 - 1;
    if (n1 < 0)
	n1 = 0;
    a = (1 << n0) * tmp - 1.0;
    if (a < 0.0)
	a = 0.0;
    b = 1.0 - a;
    a *= 1.0f / 255.0f;
    b *= 1.0f / 255.0f;

    /* Clip the x and y coordinates of the center of the sampled area. */

    if (oor == IP_FillModeWrap)
    {
	clipx = x / pyr->xmax;
	clipx = (clipx >= 0.0)?
	    (clipx - (int) clipx):
	    (1.0 + clipx + (int) (-clipx));
	clipx *= pyr->xmax;
	clipy = y / pyr->ymax;
	clipy = (clipy >= 0.0)?
	    (clipy - (int) clipy):
	    (1.0 + clipy + (int) (-clipy));
	clipy *= pyr->ymax;
    }
    else
    {
	clipx = (x < 0.0)? 0.0: ((x > pyr->xmax)? pyr->xmax: x);
	clipy = (y < 0.0)? 0.0: ((y > pyr->ymax)? pyr->ymax: y);
    }
    clipx *= pyr->xyscale1;
    clipy *= pyr->xyscale1;

    /* Extract the first sample from level "n0" by bilinear interpolation
    from the four pixels closest to the center fo the sampled area. Scale
    the sampled channel data by "b". */

    plv = pyr->pyrlevel + n0;

    fcol = clipx * plv->xyscale2;
    col0 = (Int32)fcol;
    col1 = UT_MIN (col0 + 1, plv->width - 1);
    if (oor == IP_FillModeWrap)
    {
	if (col0 >= plv->width - 1)
	    col0 = 0;
	if (col1 >= plv->width - 1)
	    col1 = 0;
    }
    u0 = fcol - col0;
    u1 = 1.0 - u0;

    frow = clipy * plv->xyscale2;
    row0 = (Int32)frow;
    row1 = UT_MIN (row0 + 1, plv->height - 1);
    if (oor == IP_FillModeWrap)
    {
	if (row0 >= plv->height - 1)
	    row0 = 0;
	if (row1 >= plv->height - 1)
	    row1 = 0;
    }
    v0 = frow - row0;
    v1 = 1.0 - v0;

    LL = plv->pixels + pyr->numchan * col0 + plv->rowbytes * row0;
    LR = plv->pixels + pyr->numchan * col1 + plv->rowbytes * row0;
    UL = plv->pixels + pyr->numchan * col0 + plv->rowbytes * row1;
    UR = plv->pixels + pyr->numchan * col1 + plv->rowbytes * row1;

    for (i = 0; i < pyr->numchan; ++i)
	cdata[pyr->channel[i]] = b *
	    (u1 * (v1 * LL[i] + v0 * UL[i]) +
	     u0 * (v1 * LR[i] + v0 * UR[i]));

    /* Extract a second sample from level "n1". Scan the sampled channel data
    by "a" and add the result to the sample taken from level "n0". */

    plv = pyr->pyrlevel + n1;

    fcol = clipx * plv->xyscale2;
    col0 = (Int32)fcol;
    col1 = UT_MIN (col0 + 1, plv->width - 1);
    if (oor == IP_FillModeWrap)
    {
	if (col0 >= plv->width - 1)
	    col0 = 0;
	if (col1 >= plv->width - 1)
	    col1 = 0;
    }
    u0 = fcol - col0;
    u1 = 1.0 - u0;

    frow = clipy * plv->xyscale2;
    row0 = (Int32)frow;
    row1 = UT_MIN (row0 + 1, plv->height - 1);
    if (oor == IP_FillModeWrap)
    {
	if (row0 >= plv->height - 1)
	    row0 = 0;
	if (row1 >= plv->height - 1)
	    row1 = 0;
    }
    v0 = frow - row0;
    v1 = 1.0 - v0;

    LL = plv->pixels + pyr->numchan * col0 + plv->rowbytes * row0;
    LR = plv->pixels + pyr->numchan * col1 + plv->rowbytes * row0;
    UL = plv->pixels + pyr->numchan * col0 + plv->rowbytes * row1;
    UR = plv->pixels + pyr->numchan * col1 + plv->rowbytes * row1;

    for (i = 0; i < pyr->numchan; ++i)
	cdata[pyr->channel[i]] += a *
	    (u1 * (v1 * LL[i] + v0 * UL[i]) +
	     u0 * (v1 * LR[i] + v0 * UR[i]));

    if (oor == IP_FillModeFill)
    {
	/* Find out how much of the sampled area is outside the region
	covered by the pyramid, and mix an appropriate amount of the
	current drawing color into the sampled color. */

	if (size > MinNrmFloat32)
	{
	    size2 = size * 0.5;
	    in_horz = size;
	    tmp = x + size2 - pyr->xmax + 0.5;
	    if (tmp > 0.0)
		in_horz -= tmp;
	    tmp = x - size2 + 0.5;
	    if (tmp < 0.0)
		in_horz += tmp;
	    if (in_horz < 0.0)
		in_horz = 0.0;
	    in_vert = size;
	    tmp = y + size2 - pyr->ymax + 0.5;
	    if (tmp > 0.0)
		in_vert -= tmp;
	    tmp = y - size2 + 0.5;
	    if (tmp < 0.0)
		in_vert += tmp;
	    if (in_vert < 0.0)
		in_vert = 0.0;
	    in = (in_horz / size) * (in_vert / size);
	    out = 1.0 - in;
	    for (i = 0; i < pyr->numchan; ++i)
		cdata[pyr->channel[i]] =
		    cdata[pyr->channel[i]] * in +
		    IP_StateDrawColorFloat[pyr->channel[i]] * out;
	}
	else if (x < 0.0 || x > pyr->xmax || y < 0.0 || y > pyr->ymax)
	{
	    for (i = 0; i < pyr->numchan; ++i)
		cdata[pyr->channel[i]] = IP_StateDrawColorFloat[pyr->channel[i]];
	}
    }
    return;
}
