/**************************************************************************
 *{@C
 *      Copyright:      1988-2023 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;
}
