/**************************************************************************
 *{@C
 *	Copyright:	1988-2020 Paul Obermeier (obermeier@poSoft.de)
 *
 *	Module:		ImageProcessing
 *	Filename:	IP_Blur.c
 *
 *	Author:		Paul Obermeier
 *
 *	Description: 	Functions to blur images using a box filter.
 *
 *	Additional documentation:
 *			None.
 *
 *	Exported functions:
 *			IP_Blur
 *			IP_DepthBlur
 *
 **************************************************************************/

#include <stdio.h>

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

#define  USE_IMG_MACROS
#include "FF_Image.h"

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

/* A struct to pass parameters to the concurrently
invoked image blurring functions */

typedef struct {
	IP_ImageId	img;		/* The image to be blurred */
	Int32		chn;		/* Image channel to be blurred */
	void		*buf;		/* Buffer for intermediate results */
	const UInt32	*mtable;	/* Multiplication table */
	Int32		n;		/* Amount of blurring to be applied */
	Float32		weight;		/* Input pixel weighting factor */
	} blurparms;

/* Amount of memory set aside for the concurrently invoked
functions' auto variable stacks */

#define stacksize	(1 << 12)
#define IS_IN_RANGE(a)	((a) > mindepth && (a) < maxdepth)

/* 
  Blur one UByte channel of an image horizontally.  This routine is an integer
  version of "blur_h_Float", which has been optimized for high speed:

  - All "pixel times weight" multiply operations have been replaced by shifts
    and lookups using a precomputed table.

  - The innermost loop, "for (x = 0; x < width; ++x) {...}", has been broken
    up into three parts.  The first and the last part deal with pixels whose
    colors are influenced by the left and right borders of the image.  The
    middle part of the loop processes all pixels, which are not influenced by
    the image borders Since it doesn't have to care about image boundaries, the
    middle part of the loop can run significantly faster than the first and the
    last part.
*/

static void blur_h_UByte (void *ptr, Int32 n_conc, Int32 i_conc)
{
    blurparms	*parms;
    FF_ImgDesc	*pimg;
    Int32	chn, m, n, x, x1, x2, y, xc, xmax, ymin, ymax, width, height;
    UInt8	*src;
    UInt32	*buf, *dst, average;
    const UInt32 *mtable;

    /* Main loop initialization */

    parms = ptr;
    buf = parms->buf;
    mtable = parms->mtable;
    pimg = &parms->img->pimg;
    chn = parms->chn;
    n = parms->n;
    m = n + 1;
    width = pimg->desc.width;
    xmax = width - 1;
    x1 = UT_MIN (n, width);
    x2 = width - m;
    height = pimg->desc.height;
    n_conc = (n_conc == 0? 1: n_conc);
    ymin = (height * i_conc) / n_conc;
    ymax = (height * (i_conc + 1)) / n_conc;

    /* Main loop */

    for (y = ymin; y < ymax; ++y) {
	/* Process next horizontal row. */

	src = FF_ImgUByteRow (pimg, chn, y);

	/* Per row initialization: compute the color of the leftmost pixel. */

	average = 0;
	for (x = -n; x <= n; ++x) {
	    if (x <= 0)
		xc = 0;
	    else if (x >= xmax)
		xc = xmax;
	    else
		xc = x;
	    average += mtable[src[xc]];
	}

	/* Process pixels influenced by the left end of the row. */

	dst = buf + y * width;
	x = 0;
	while (x < x1) {
	    dst[x] = average >> 16;
	    average -= mtable[src[0]];
	    xc = x + m;
	    if (xc > xmax)
		xc = xmax;
	    average += mtable[src[xc]];
	    ++x;
	}

	/* Process pixels not influenced by either end of the row. */

	while (x < x2) {
	    dst[x] = average >> 16;
	    average -= mtable[src[x - n]];
	    average += mtable[src[x + m]];
	    ++x;
	}

	/* Process pixels influenced by the right end of the row. */

	while (x < width) {
	    dst[x] = average >> 16;
	    average -= mtable[src[x - n]];
	    average += mtable[src[xmax]];
	    ++x;
	}
    }
    return;
}

/* 
  Blur one UByte channel of an image vertically.  This routine is an integer
  version of "blur_v_Float", which has been optimized for high speed:

  - All "pixel times weight" multiply operations have been replaced by shifts
    and lookups using a precomputed table.

  - The innermost loop, "for (y = 0; y < height; ++y) {...}", has been broken
    up into three parts, analogous to "blur_h_UByte".
    
  - Most of the multiply operations resulting from column-by-column processing
    of the image have been eliminated by strength reduction.
*/

static void blur_v_UByte (void *ptr, Int32 n_conc, Int32 i_conc)
{
    blurparms	*parms;
    FF_ImgDesc	*pimg;
    Int32	chn, m, n, x, y, y1, y2, yc, xmin, xmax, ymax, width, height;
    UInt32	average;
    const UInt32 *mtable, *buf, *src, *src1, *src2;

    /* Main loop initialization */

    parms = ptr;
    buf = parms->buf;
    mtable = parms->mtable;
    pimg = &parms->img->pimg;
    chn = parms->chn;
    n = parms->n;
    m = n + 1;
    width = pimg->desc.width;
    n_conc = (n_conc == 0? 1: n_conc);
    xmin = (width * i_conc) / n_conc;
    xmax = (width * (i_conc + 1)) / n_conc;
    height = pimg->desc.height;
    ymax = height - 1;
    y1 = UT_MIN (n, height);
    y2 = height - m;

    /* Main loop */

    for (x = xmin; x < xmax; ++x) {
	/* Process next vertical column. */

	src = buf + x;

	/* Per column initialization: compute the color of the bottom pixel. */

	average = 0;
	for (y = -n; y <= n; ++y) {
	    if (y <= 0)
		yc = 0;
	    else if (y >= ymax)
		yc = ymax;
	    else
		yc = y;
	    average += mtable[src[yc * width]];
	}

	/* Process pixels influenced by the bottom end of the column. */

	y = 0;
	while (y < y1) {
	    *FF_ImgUBytePixel (pimg, chn, x, y) = average >> 24;
	    average -= mtable[src[0]];
	    yc = y + m;
	    if (yc > ymax)
		yc = ymax;
	    average += mtable[src[yc * width]];
	    ++y;
	}

	/* Process pixels not influenced by either end of the column. */

	src1 = src;
	src2 = src + (n + m) * width;
	while (y < y2) {
	    *FF_ImgUBytePixel (pimg, chn, x, y) = average >> 24;
	    average -= mtable[*src1];
	    src1 += width;
	    average += mtable[*src2];
	    src2 += width;
	    ++y;
	}

	/* Process pixels influenced by the top end of the column. */

	src1 = src + (y - n) * width;
	src2 = src +  ymax   * width;
	while (y < height) {
	    *FF_ImgUBytePixel (pimg, chn, x, y) = average >> 24;
	    average -= mtable[*src1];
	    src1 += width;
	    average += mtable[*src2];
	    ++y;
	}
    }
    return;
}

/* Blur one Float channel of an image horizontally. */

static void blur_h_Float (void *ptr, Int32 n_conc, Int32 i_conc)
{
    blurparms	*parms;
    FF_ImgDesc	*pimg;
    Int32	chn, n, x, y, xc, ymin, ymax, width, height;
    Float32	*buf, average, weight;

    parms = ptr;
    buf = parms->buf;
    pimg = &parms->img->pimg;
    chn = parms->chn;
    n = parms->n;
    weight = parms->weight;
    width = pimg->desc.width;
    height = pimg->desc.height;
    n_conc = (n_conc == 0? 1: n_conc);
    ymin = (height * i_conc) / n_conc;
    ymax = (height * (i_conc + 1)) / n_conc;

    for (y = ymin; y < ymax; ++y) {
	average = 0.0;
	for (x = -n; x <= n; ++x) {
	    xc = UT_MAX (0, UT_MIN (x, width - 1));
	    average += weight * *FF_ImgFloatPixel (pimg, chn, xc, y);
	}
	for (x = 0; x < width; ++x) {
	    buf[x + y * width] = average;
	    xc = UT_MAX (0, x - n);
	    average -= weight * *FF_ImgFloatPixel (pimg, chn, xc, y);
	    xc = UT_MIN (x + n + 1, width - 1);
	    average += weight * *FF_ImgFloatPixel (pimg, chn, xc, y);
	}
    }
    return;
}

/* Blur one Float channel of an image vertically. */

static void blur_v_Float (void *ptr, Int32 n_conc, Int32 i_conc)
{
    blurparms	*parms;
    FF_ImgDesc	*pimg;
    Int32	chn, n, x, y, yc, xmin, xmax, width, height;
    Float32	*buf, average, weight;

    parms = ptr;
    buf = parms->buf;
    pimg = &parms->img->pimg;
    chn = parms->chn;
    n = parms->n;
    weight = parms->weight;
    width = pimg->desc.width;
    height = pimg->desc.height;
    n_conc = (n_conc == 0? 1: n_conc);
    xmin = (width * i_conc) / n_conc;
    xmax = (width * (i_conc + 1)) / n_conc;

    for (x = xmin; x < xmax; ++x) {
	average = 0.0;
	for (y = -n; y <= n; ++y) {
	    yc = UT_MAX (0, UT_MIN (y, height - 1));
	    average += weight * buf[x + yc * width];
	}
	for (y = 0; y < height; ++y) {
	    *FF_ImgFloatPixel (pimg, chn, x, y) = average;
	    yc = UT_MAX (0, y - n);
	    average -= weight * buf[x + yc * width];
	    yc = UT_MIN (y + n + 1, height - 1);
	    average += weight * buf[x + yc * width];
	}
    }
    return;
}

/* Blur an image. For an explanation of the parameters, see the comments
   on function IP_Blur. */

static UT_Bool blur (IP_ImageId img, Float32 f, Float32 g)
{
    Float32	weight_h, weight_v;
    Int32	size, n_h, n_v, i;
    blurparms	parms;
    UInt32	uiweight;
    static Int32  mtable_n_h = -1, mtable_n_v = -1;
    static UInt32 mtable_h[1 << 8], mtable_v[1 << 16];

    /* Allocate memory for temporary buffer used by the concurrent
       blurring routines to hold intermediate results. */

    size = img->pimg.desc.width * img->pimg.desc.height *
 	   UT_MAX (sizeof (Float32), sizeof (UInt32));
    if (!(parms.buf = UT_MemTempAlloc (size, "blur temp buf" UT_MemIdString))) {
	return UT_False;
    }

    /* Precompute some data used by the concurrent blurring routines. */

    parms.img = img;
    n_h = (Int32) (f * 0.5);
    n_h = UT_MAX (0, n_h);
    weight_h = 1.0 / (2 * n_h + 1);
    n_v = (Int32) (g * 0.5);
    n_v = UT_MAX (0, n_v);
    weight_v = 1.0 / (2 * n_v + 1);

    /* Set up multiplication tables, so that most of the expensive
       multiply operations in "blur_h_UByte" and "blur_v_UByte" can be
       replaced with cheaper table lookups */

    if (mtable_n_h != n_h) {
	uiweight = (UInt32) (weight_h * (1 << 24) + 0.5);
	for (i = 0; i < (1 << 8); ++i) {
	    mtable_h[i] = i * uiweight;
	}
	mtable_n_h = n_h;
    }
    if (mtable_n_v != n_v) {
	uiweight = (UInt32) (weight_v * (1 << 16) + 0.5);
	for (i = 0; i < (1 << 16); ++i) {
	    mtable_v[i] = i * uiweight;
	}
	mtable_n_v = n_v;
    }

    /* Blur the image, distributing the work onto multiple threads
       running concurrently. */

    for (parms.chn = 0; parms.chn < FF_NumImgChanTypes; ++parms.chn) {
	if (!IP_StateDrawMask[parms.chn]) {
	    continue;
	}
	switch (img->pimg.channel[parms.chn]) {
	    case FF_ImgFmtTypeNone: {
		break;
	    }
	    case FF_ImgFmtTypeUByte: {
		parms.n = n_h;
		parms.mtable = mtable_h;
		if (!UT_THREAD_EXEC (blur_h_UByte, &parms,
				     IP_StateNumThreads, stacksize)) {
		    return UT_False;
		}
		parms.n = n_v;
		parms.mtable = mtable_v;
		if (!UT_THREAD_EXEC (blur_v_UByte, &parms,
				     IP_StateNumThreads, stacksize)) {
		    return UT_False;
		}
		break;
	    }
	    case FF_ImgFmtTypeFloat: {
		parms.n = n_h;
		parms.weight = weight_h;
		if (!UT_THREAD_EXEC (blur_h_Float, &parms,
				    IP_StateNumThreads, stacksize)) {
		    return UT_False;
		}
		parms.n = n_v;
		parms.weight = weight_v;
		if (!UT_THREAD_EXEC (blur_v_Float, &parms,
				     IP_StateNumThreads, stacksize)) {
		    return UT_False;
		}
		break;
	    }
	}
    }
    return UT_True;
}

/* Copy the pixel of image "isrc" at position (xsrc, ysrc) into image
   "idest" onto position (xdest, ydest). */

static void copy_pixel (FF_ImgDesc *idest, Int32 xdest, Int32 ydest, 
			FF_ImgDesc *isrc,  Int32 xsrc,  Int32 ysrc)
{
    Int32 chn;

    xsrc  = UT_MAX (0, UT_MIN (xsrc,  isrc->desc.width-1));
    ysrc  = UT_MAX (0, UT_MIN (ysrc,  isrc->desc.height-1));
    xdest = UT_MAX (0, UT_MIN (xdest, idest->desc.width-1));
    ydest = UT_MAX (0, UT_MIN (ydest, idest->desc.height-1));

    for (chn=0; chn<FF_NumImgChanTypes; chn++) {
	if (!IP_StateDrawMask[chn]) {
	    continue;
	}
	switch (idest->channel[chn]) {
	    case FF_ImgFmtTypeNone: {
		break;
	    }
	    case FF_ImgFmtTypeUByte: {
		*FF_ImgUBytePixel (idest, chn, xdest, ydest) =
		    *FF_ImgUBytePixel (isrc, chn, xsrc, ysrc);
		break;
	    }
	    case FF_ImgFmtTypeFloat: {
		*FF_ImgFloatPixel (idest, chn, xdest, ydest) =
		    *FF_ImgFloatPixel (isrc, chn, xsrc, ysrc);
		break;
	    }
	}
    }
    return;
}

/* Copy all pixels from image "img" into image "bimg".
   Pixels with depth less than "mindepth" and greater than "maxdepth" are
   copied without change. Pixel with depth between "mindepth" and
   "maxdepth". */

static void prepare_bimg (IP_ImageId img, IP_ImageId bimg, 
			  Float32 mindepth, Float32 maxdepth)
{
    Int32 x, y, width, height;
    Int32 scan1, scan2;
    FF_ImgDesc *pimg, *pbimg;
    UT_Bool in_range;
    Float32 dval;

    pimg = &img->pimg;
    pbimg = &bimg->pimg;
    height = img->pimg.desc.height;
    width  = img->pimg.desc.width;

    for (y=0; y<height; y++)
    {
	scan1 = -1;
	scan2 = -1;
	for (x=0; x<width; x++)
	{
	    dval = *FF_ImgFloatPixel (pimg, FF_ImgChanTypeDepth, x, y);
	    in_range = !IS_IN_RANGE (dval);
	    copy_pixel (pbimg, x, y, pimg, x, y);
	    if (scan1 < 0 && !in_range)
	    {
		;
	    }
	    else if (scan1 < 0 && in_range )
	    {
		scan1 = x;
	    }
	    else if (scan1 >= 0)
	    {
	        if (!in_range)
		    scan2 = x - 1;
		else if (x == width -1)
		    scan2 = width - 1;
	    }

	    if (scan1 >= 0 && scan2 >= 0)
	    {
		if (scan1 == scan2)
		{
		    if (scan1 == 0)
			copy_pixel (pbimg, scan1, y, pbimg, scan1+1, y);
		    else if (scan1 == width - 1)
			copy_pixel (pbimg, scan1, y, pbimg, scan1-1, y);
		    else
			copy_pixel (pbimg, scan1, y, pbimg, scan1+1, y);
		}
		else
		{
		    while (scan1 < scan2)
		    {
			copy_pixel (pbimg, scan1, y, pbimg, scan1-1, y);
			copy_pixel (pbimg, scan2, y, pbimg, scan2+1, y);
			scan1++;
			scan2--;
		    }
		}
		scan1 = -1;
		scan2 = -1;
	    }
	}
    }
    return;
}

/* Mix the original and the blurred image according to the pixel's depth.
   If the pixel's depth is less than "mindepth", the new pixel value is taken
   from "img". If the depth is greater than "maxdepth" the new pixel value
   is taken from "bimg". A linear interpolation between the pixel values of
   "img" and "bimg" is applied otherwise. */

static void blurmix (IP_ImageId img, IP_ImageId bimg, 
		     Float32 mindepth, Float32 maxdepth, 
		     Float32 imin, Float32 imax)
{
    Int32 chn, y, width, height;
    Float32 *dsrc, fact, tmp;
    FF_ImgDesc *pimg, *pbimg;

    width  = img->pimg.desc.width;
    height = img->pimg.desc.height;
    fact   = 1.0 / (imax - imin);

    pimg  = &img->pimg;
    pbimg = &bimg->pimg;

    for (chn=0; chn<FF_NumImgChanTypes; chn++) {
	if (!IP_StateDrawMask[chn]) {
	    continue;
	}

	switch (img->pimg.channel[chn]) {
	    case FF_ImgFmtTypeNone: {
		break;
	    }
	    case FF_ImgFmtTypeUByte: {
		UInt8 *src, *stop, *bsrc;
		for (y=0; y<height; y++) {
		    src  = FF_ImgUByteRow (pimg, chn, y);
		    bsrc = FF_ImgUByteRow (pbimg, chn, y);
		    dsrc = FF_ImgFloatRow (pimg, FF_ImgChanTypeDepth, y);
		    stop = src + width;
		    while (src < stop) {
			if (IS_IN_RANGE (*dsrc)) {
			    tmp = fact * (*dsrc - imin);
			    *src = (UInt8) (tmp * *src + (1.0 - tmp) * *bsrc);
			}
			src++;
			bsrc++;
			dsrc++;
		    }
		}
		break;
	    }
	    case FF_ImgFmtTypeFloat: {
		Float32 *src, *stop, *bsrc;
		for (y=0; y<height; y++) {
		    src  = FF_ImgFloatRow (pimg, chn, y);
		    bsrc = FF_ImgFloatRow (pbimg, chn, y);
		    dsrc = FF_ImgFloatRow (pimg, FF_ImgChanTypeDepth, y);
		    stop = src + width;
		    while (src < stop) {
			if (IS_IN_RANGE (*dsrc)) {
			    tmp = fact * (*dsrc - imin);
			    *src = tmp * *src + (1.0 - tmp) * *bsrc;
			}
			src++;
			bsrc++;
			dsrc++;
		    }
		}
		break;
	    }
	}
    }
    return;
}


/* Blur an image according to depth information. 
   For an explanation of the parameters, see the 
   comments on function IP_DepthBlur. */

#define restore {		\
    IP_ImageDelete     (bimg);	\
    IP_SetDrawMask  (tmpmsk);	\
    IP_SetFormat (tmpfmt);  }

static UT_Bool depthblur (IP_ImageId img, Float32 f, Float32 g, 
		       Float32 mindepth, Float32 maxdepth)
{
    Int32       chn;
    Float32     imin, imax;
    IP_ImageId bimg;
    UT_Bool      tmpmsk[FF_NumImgChanTypes];
    FF_ImgFmtType tmpfmt[FF_NumImgChanTypes];

    if (!img->pimg.channel[FF_ImgChanTypeDepth])
    {
	UT_ErrSetNum (UT_ErrParamInvalid, str_fg_depth);
	return UT_False;
    }

    IP_GetDrawMask  (tmpmsk);
    IP_GetFormat (tmpfmt);

    for (chn=0; chn<FF_NumImgChanTypes; chn++)
	IP_StateFormat[chn] = img->pimg.channel[chn];
    IP_StateFormat[FF_ImgChanTypeDepth] = FF_ImgFmtTypeNone;
    IP_StateDrawMask[FF_ImgChanTypeDepth] = UT_False;

    /* Allocate memory for the blurred image and copy over all channels,
       for which the current draw mask is on. */
    if (!(bimg = IP_ImageNew (img->pimg.desc.width, 
				 img->pimg.desc.height,
				 img->pimg.desc.aspect)))
    {
	restore;
	return UT_False;
    }

    if (!IP_GetDepthRange (img, &imin, &imax))
    {
	restore;
	return UT_False;
    }

    /* If depth of all pixels is the same, do nothing. */
    if (UT_ABS (imax - imin) < MinNrmFloat32)
	return UT_True;

    if (mindepth < 0.0)
	mindepth = imin;	
    if (maxdepth < 0.0)
	maxdepth = imax;	

    prepare_bimg (img, bimg, mindepth, maxdepth);

    if (!blur (bimg, f, g))
    {
	restore;
	return UT_False;
    }

    blurmix (img, bimg, mindepth, maxdepth, imin, imax);
    restore;
    return UT_True;
}

/***************************************************************************
 *[@e
 *	Name:		IP_Blur
 *
 *	Usage:		Blurs an image.
 *
 *	Synopsis:	UT_Bool IP_Blur (IP_ImageId img, Float32 f, Float32 g)
 *
 *	Description:	All channels in image "img", for which the current
 *			drawing mask is "UT_True", are blurred:
 *			The color of every pixel in "img" is distributed over
 *			a rectangular area of "f" by "g" pixels using a box
 *			filter.
 *
 *			Note that IP_Blur can be run concurrently on
 *			multiple processors; see IP_SetNumThreads.
 *
 *	Return value:	UT_True if successful, else UT_False.
 *
 *	See also:	IP_DepthBlur
 *
 ***************************************************************************/

UT_Bool IP_Blur (IP_ImageId img, Float32 f, Float32 g)
{
    UT_MemState	memstate;
    UT_Bool	success;

    memstate = UT_MemRemember ();
    success = blur (img, f, g);
    UT_MemRestore (memstate);
    return success;
}

/***************************************************************************
 *[@e
 *	Name:		IP_DepthBlur
 *
 *	Usage:		Blurs an image according to a depth channel.
 *
 *	Synopsis:	UT_Bool IP_DepthBlur 
 *				(IP_ImageId img, Float32 f, Float32 g,
 *			         Float32 mindepth, Float32 maxdepth)
 *
 *	Description:	All channels in image "img", for which the current
 *			drawing mask is "UT_True", are blurred:
 *			The color of every pixel in "img" is distributed over
 *			a rectangular area of "f" by "g" pixels using a box
 *			filter.
 *			Only pixels whose depth is between "mindepth" and
 *			"maxdepth" are blurred.
 *
 *	Return value:	UT_True if successful, else UT_False.
 *
 *	See also:	IP_Blur
 *
 ***************************************************************************/

UT_Bool IP_DepthBlur (IP_ImageId img, Float32 f, Float32 g, 
		    Float32 mindepth, Float32 maxdepth)
{
    UT_MemState	memstate;
    UT_Bool	success;

    memstate = UT_MemRemember ();
    success = depthblur (img, f, g, 
                         UT_ABS (maxdepth) < MinNrmFloat32? 0.0: 1.0/maxdepth,
                         UT_ABS (mindepth) < MinNrmFloat32? 0.0: 1.0/mindepth);
    UT_MemRestore (memstate);
    return success;
}
