/**************************************************************************
 *{@C
 *      Copyright:      1988-2020 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         ImageProcessing
 *      Filename:       IP_Filter.c
 *
 *      Author:         Paul Obermeier
 *
 *      Description:    Various image filtering functions.
 *
 *      Additional documentation:
 *                      None.
 *
 *      Exported functions:
 *                      IP_ChangeGamma
 *                      IP_RemapValues
 *                      IP_Derivates
 *                      IP_Median
 *                      IP_MedianSequence
 *                      IP_Dilatation
 *                      IP_Erosion
 *                      IP_Threshold
 *                      IP_CutOff
 *                      IP_MarkNonZeroPixels
 *                      IP_DifferenceImage
 *
 **************************************************************************/

#include <stdio.h>

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

#define  USE_IMG_MACROS
#include "FF_Image.h"

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

/* Gamma-filter and scale one channel of an image. For an explanation
   of the parameters, see the comments on function IP_ChangeGamma. */

static UT_Bool change_gamma
        (IP_ImageId img, FF_ImgChanType chn,
         Float32 gamma, Float32 scl, Float32 off)
{
    UInt8       ubgtab[256], *ubdata, *ubstop;
    Float32     fgtab[GTABSIZE], *fdata;
    Int32       y, i, tmp;

    switch (img->pimg.channel[chn]) {
        case FF_ImgFmtTypeNone: {
            break;
        }
        case FF_ImgFmtTypeUByte: {
            if (!IP_UByteGammaTable (1.0 / gamma, ubgtab)) {
                return UT_False;
            }
            for (i = 0; i < 256; ++i) {
                tmp = (Int32) (ubgtab[i] * scl + off * 255.0);
                ubgtab[i] = UT_MAX (0, UT_MIN (tmp, 255));
            }
            for (y = 0; y < img->pimg.desc.height; ++y) {
                ubdata = FF_ImgUByteRow (&img->pimg, chn, y);
                ubstop = ubdata + img->pimg.desc.width;
                while (ubdata < ubstop) {
                    *ubdata = ubgtab[*ubdata];
                    ++ubdata;
                }
            }
            break;
        }
        case FF_ImgFmtTypeFloat: {
            if (!IP_FloatGammaTable (1.0 / gamma, fgtab)) {
                return UT_False;
            }
            for (i = 0; i < GTABSIZE; ++i) {
                fgtab[i] = fgtab[i] * scl + off;
            }
            for (y = 0; y < img->pimg.desc.height; ++y) {
                fdata = FF_ImgFloatRow (&img->pimg, chn, y);
                IP_GammaFloat2Float
                    (img->pimg.desc.width, fdata, fgtab, fdata);
            }
            break;
        }
    }
    return UT_True;
}

/* Remap pixel values using a lookup table. For an explanation of the
   parameters, see the comments on function IP_RemapValues. */

static void remap_values
        (IP_ImageId img, FF_ImgChanType chn,
         Int32 lutsize, const Float32 lut[])
{
    UInt8       *ubdata, *ubstop, ublut[256];
    Float32     *fdata, *fstop, lutscale, ftmp, s, t;
    Int32       lutclip, i, y, a, b, c;

    switch (img->pimg.channel[chn]) {
        case FF_ImgFmtTypeNone: {
            break;
        }
        case FF_ImgFmtTypeUByte: {
            lutscale = (lutsize - 1) / 255.0;
            lutclip = lutsize - 2;
            for (i = 0; i < 256; ++i) {
                ftmp = i * lutscale;
                a = (Int32)ftmp;
                a = UT_MAX (0, UT_MIN (a, lutclip));
                b = a + 1;
                s = ftmp - a;
                t = b - ftmp;
                c = (Int32) (255.0 * (t * lut[a] + s * lut[b]));
                ublut[i] = UT_MAX (0, UT_MIN (c, 255));
            }
            for (y = 0; y < img->pimg.desc.height; ++y) {
                ubdata = FF_ImgUByteRow (&img->pimg, chn, y);
                ubstop = ubdata + img->pimg.desc.width;
                while (ubdata < ubstop) {
                    *ubdata = ublut[*ubdata];
                    ++ubdata;
                }
            }
            break;
        }
        case FF_ImgFmtTypeFloat: {
            lutscale = (Float32) (lutsize - 1);
            lutclip = lutsize - 2;
            for (y = 0; y < img->pimg.desc.height; ++y) {
                fdata = FF_ImgFloatRow (&img->pimg, chn, y);
                fstop = fdata + img->pimg.desc.width;
                while (fdata < fstop) {
                    ftmp = *fdata * lutscale;
                    a = (Int32)ftmp;
                    a = UT_MAX (0, UT_MIN (a, lutclip));
                    b = a + 1;
                    s = ftmp - a;
                    t = b - ftmp;
                    *(fdata++) = t * lut[a] + s * lut[b];
                }
            }
            break;
        }
    }
    return;
}

/* Copy channel "chn" of scan line "y" of image "img" into buffer "buf",
   converting the pixel data to FF_ImgFmtTypeFloat format.
   Function "scan_to_buf" observes the "wrap" flag as in IP_Derivates. */

static void scan_to_buf
        (IP_ImageId img, FF_ImgChanType chn, Int32 y, 
         Float32 *buf, UT_Bool wrap)
{
    UInt8       *ubsrc;
    Float32     *fsrc, *fdest, *fstop;

    if (y < 0)
        y = wrap? img->pimg.desc.height - 1: 0;
    fdest = buf + 1;
    fstop = fdest + img->pimg.desc.width;
    switch (img->pimg.channel[chn]) {
        case FF_ImgFmtTypeNone:
            while (fdest < fstop)
                *(fdest++) = 0.0;
            break;
        case FF_ImgFmtTypeUByte:
            ubsrc = FF_ImgUByteRow (&img->pimg, chn, y);
            while (fdest < fstop)
                *(fdest++) = (1.0 / 255.0) * *(ubsrc++);
            break;
        case FF_ImgFmtTypeFloat:
            fsrc = FF_ImgFloatRow (&img->pimg, chn, y);
            while (fdest < fstop)
                *(fdest++) = *(fsrc++);
            break;
    }
    buf[0] = buf[wrap? img->pimg.desc.width: 1];
    return;
}

/* Copy the contents of buffer "buf" into channel "chn" of scan line "y" of
   image "img", converting the format of the pixel data if necessary. */

static void buf_to_scan
        (Float32 *buf, IP_ImageId img, FF_ImgChanType chn, Int32 y, Int32 width)
{
    Int32       tmp;
    UInt8       *ubdest;
    Float32     *fdest, *fsrc, *fstop;

    fsrc = buf;
    fstop = buf + width;
    switch (img->pimg.channel[chn]) {
        case FF_ImgFmtTypeNone: {
            break;
        }
        case FF_ImgFmtTypeUByte: {
            ubdest = FF_ImgUByteRow (&img->pimg, chn, y);
            while (fsrc < fstop) {
                tmp = (Int32) (255.0 * *(fsrc++));
                *(ubdest++) = UT_MAX (0, UT_MIN (tmp, 255));
            }
            break;
        }
        case FF_ImgFmtTypeFloat: {
            fdest = FF_ImgFloatRow (&img->pimg, chn, y);
            while (fsrc < fstop) {
                *(fdest++) = *(fsrc++);
            }
            break;
        }
    }
    return;
}

/* Compute differences between adjacent pixels in horizontal and in
   vertical direction. For an explanation of the parameters, see the
   comments on function IP_Derivates. */

static UT_Bool hvdiff
        (IP_ImageId srcimg, FF_ImgChanType srcchn,
         IP_ImageId destimg, FF_ImgChanType hchn, FF_ImgChanType vchn,
         UT_Bool wrap)
{
    Float32     *scan0, *scan1, *diff, *src0, *src1, *dest, *stop;
    Int32       width, height, y;

    if (!(diff = UT_MemTempArray (srcimg->pimg.desc.width, Float32)))
        return UT_False;
    if (!(scan0 = UT_MemTempArray (srcimg->pimg.desc.width + 1, Float32)) ||
        !(scan1 = UT_MemTempArray (srcimg->pimg.desc.width + 1, Float32)))
        return UT_False;
    width = UT_MIN (srcimg->pimg.desc.width, destimg->pimg.desc.width);
    height = UT_MIN (srcimg->pimg.desc.height, destimg->pimg.desc.height);
    scan_to_buf (srcimg, srcchn, -1, scan0, wrap);
    for (y = 0; y < height; ++y) {
        scan_to_buf (srcimg, srcchn, y, scan1, wrap);
        stop = diff + width;
        dest = diff;
        src0 = scan1;
        src1 = scan1 + 1;
        while (dest < stop)
            *(dest++) = 0.5 + *(src1++) - *(src0++);
        buf_to_scan (diff, destimg, hchn, y, width);
        dest = diff;
        src0 = scan0 + 1;
        src1 = scan1 + 1;
        while (dest < stop)
            *(dest++) = 0.5 + *(src1++) - *(src0++);
        buf_to_scan (diff, destimg, vchn, y, width);
        UT_SWAP (scan0, scan1, Float32 *);
    }
    return UT_True;
}

/* Apply a three pixel wide horizontal median filter to all channels in image
   "img" for which the current drawing mask is UT_True. */

static void median_h (IP_ImageId img)
{
    Int32       x, y, chn;
    UInt8       ub0, ub1, ub2;
    Float32     f0, f1, f2;
    FF_ImgDesc  *pimg;

    pimg = &img->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: {
                for (y = 0; y < pimg->desc.height; ++y) {
                    ub0 = *FF_ImgUBytePixel (pimg, chn, 0, y);
                    ub1 = *FF_ImgUBytePixel (pimg, chn, 1, y);
                    for (x = 2; x < pimg->desc.width; ++x) {
                        ub2 = *FF_ImgUBytePixel (pimg, chn, x, y);
                        *FF_ImgUBytePixel (pimg, chn, x - 1, y) =
                            (ub0 > ub1?
                             (ub0 > ub2?
                              (ub1 > ub2? ub1: ub2)
                              :
                              ub0
                              )
                             :
                             (ub0 > ub2?
                              ub0
                              :
                              (ub1 > ub2? ub2: ub1)
                              )
                            );
                        ub0 = ub1;
                        ub1 = ub2;
                    }
                }
                break;
            }
            case FF_ImgFmtTypeFloat: {
                for (y = 0; y < pimg->desc.height; ++y) {
                    f0 = *FF_ImgFloatPixel (pimg, chn, 0, y);
                    f1 = *FF_ImgFloatPixel (pimg, chn, 1, y);
                    for (x = 2; x < pimg->desc.width; ++x) {
                        f2 = *FF_ImgFloatPixel (pimg, chn, x, y);
                        *FF_ImgFloatPixel (pimg, chn, x - 1, y) =
                            (f0 > f1?
                             (f0 > f2?
                              (f1 > f2? f1: f2)
                              :
                              f0
                              )
                             :
                             (f0 > f2?
                              f0
                              :
                              (f1 > f2? f2: f1)
                              )
                            );
                        f0 = f1;
                        f1 = f2;
                    }
                }
                break;
            }
        }
    }
    return;
}

/* Apply a three pixel wide vertical median filter to all channels in image
   "img" for which the current drawing mask is UT_True. */

static void median_v (IP_ImageId img)
{
    Int32       x, y, chn;
    UInt8       ub0, ub1, ub2;
    Float32     f0, f1, f2;
    FF_ImgDesc  *pimg;

    pimg = &img->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: {
                for (x = 0; x < pimg->desc.width; ++x) {
                    ub0 = *FF_ImgUBytePixel (pimg, chn, x, 0);
                    ub1 = *FF_ImgUBytePixel (pimg, chn, x, 1);
                    for (y = 2; y < pimg->desc.height; ++y) {
                        ub2 = *FF_ImgUBytePixel (pimg, chn, x, y);
                        *FF_ImgUBytePixel (pimg, chn, x, y - 1) =
                            (ub0 > ub1?
                             (ub0 > ub2?
                              (ub1 > ub2? ub1: ub2)
                              :
                              ub0
                              )
                             :
                             (ub0 > ub2?
                              ub0
                              :
                              (ub1 > ub2? ub2: ub1)
                              )
                            );
                        ub0 = ub1;
                        ub1 = ub2;
                    }
                }
                break;
            }
            case FF_ImgFmtTypeFloat: {
                for (x = 0; x < pimg->desc.width; ++x) {
                    f0 = *FF_ImgFloatPixel (pimg, chn, x, 0);
                    f1 = *FF_ImgFloatPixel (pimg, chn, x, 1);
                    for (y = 2; y < pimg->desc.height; ++y) {
                        f2 = *FF_ImgFloatPixel (pimg, chn, x, y);
                        *FF_ImgFloatPixel (pimg, chn, x, y - 1) =
                            (f0 > f1?
                             (f0 > f2?
                              (f1 > f2? f1: f2)
                              :
                              f0
                              )
                             :
                             (f0 > f2?
                              f0
                              :
                              (f1 > f2? f2: f1)
                              )
                            );
                        f0 = f1;
                        f1 = f2;
                    }
                }
                break;
            }
        }
    }
    return;
}

/* Apply a median filter to all channels in image sequence "srcimg1",
   "srcimg2", "srcimg3", for which the current drawing mask is UT_True. */

static void median_sequence
        (IP_ImageId srcimg1,
         IP_ImageId srcimg2,
         IP_ImageId srcimg3,
         IP_ImageId destimg)
{
    Int32       x, y, chn, itmp, width, height;
    Float32     src1 = 0.0, src2 = 0.0, src3 = 0.0, dest = 0.0;
    FF_ImgDesc  *psrc1, *psrc2, *psrc3, *pdest;

    pdest = &destimg->pimg;
    psrc1 = &srcimg1->pimg;
    psrc2 = &srcimg2->pimg;
    psrc3 = &srcimg3->pimg;
    width = pdest->desc.width;
    width = UT_MIN (width, psrc1->desc.width);
    width = UT_MIN (width, psrc2->desc.width);
    width = UT_MIN (width, psrc3->desc.width);
    height = pdest->desc.height;
    height = UT_MIN (height, psrc1->desc.height);
    height = UT_MIN (height, psrc2->desc.height);
    height = UT_MIN (height, psrc3->desc.height);
    for (chn = 0; chn < FF_NumImgChanTypes; ++chn) {
        if (!IP_StateDrawMask[chn] || !pdest->channel[chn]) {
            continue;
        }
        for (y = 0; y < height; ++y) {
            for (x = 0; x < width; ++x) {
                switch (psrc1->channel[chn]) {
                    case FF_ImgFmtTypeNone: {
                        src1 = 0.0;
                        break;
                    }
                    case FF_ImgFmtTypeUByte: {
                        src1 = (1.0 / 255.0) *
                            *FF_ImgUBytePixel (psrc1, chn, x, y);
                        break;
                    }
                    case FF_ImgFmtTypeFloat: {
                        src1 = *FF_ImgFloatPixel (psrc1, chn, x, y);
                        break;
                    }
                }
                switch (psrc2->channel[chn]) {
                    case FF_ImgFmtTypeNone: {
                        src2 = 0.0;
                        break;
                    }
                    case FF_ImgFmtTypeUByte: {
                        src2 = (1.0 / 255.0) *
                            *FF_ImgUBytePixel (psrc2, chn, x, y);
                        break;
                    }
                    case FF_ImgFmtTypeFloat: {
                        src2 = *FF_ImgFloatPixel (psrc2, chn, x, y);
                        break;
                    }
                }
                switch (psrc3->channel[chn]) {
                    case FF_ImgFmtTypeNone: {
                        src3 = 0.0;
                        break;
                    }
                    case FF_ImgFmtTypeUByte: {
                        src3 = (1.0 / 255.0) *
                            *FF_ImgUBytePixel (psrc3, chn, x, y);
                        break;
                    }
                    case FF_ImgFmtTypeFloat: {
                        src3 = *FF_ImgFloatPixel (psrc3, chn, x, y);
                        break;
                    }
                }
                dest =
                    (src1 > src2?
                     (src1 > src3?
                      (src2 > src3? src2: src3)
                      :
                      src1
                      )
                     :
                     (src1 > src3?
                      src1
                      :
                      (src2 > src3? src3: src2)
                      )
                    );
                switch (pdest->channel[chn]) {
                    case FF_ImgFmtTypeNone: {
                        break;
                    }
                    case FF_ImgFmtTypeUByte: {
                        itmp = (Int32) (255.0 * dest);
                        *FF_ImgUBytePixel (pdest, chn, x, y) =
                            UT_MAX (0, UT_MIN (itmp, 255));
                        break;
                    }
                    case FF_ImgFmtTypeFloat: {
                        *FF_ImgFloatPixel (pdest, chn, x, y) = dest;
                        break;
                    }
                }
            }
        }
    }
    return;
}

/* Apply a three by three pixel dilatation/erosion filter to one
   scan line of a UByte channel of an image. */

static void dler_UByte_scan (UInt8 *s1, UInt8 *s2, UInt8 *s3, 
                             Int32 len, Int32 mode, UInt8 *dest)
{
    UInt8 ub;
    UInt8 *stop = s1 + len;

    if (mode == 1) {
        while (s1 < stop) {
            ub = 0;

            ub = UT_MAX (ub, UT_MAX (*(s1-1), UT_MAX (*s1, *(s1+1))));
            ub = UT_MAX (ub, UT_MAX (*(s2-1), UT_MAX (*s2, *(s2+1))));
            ub = UT_MAX (ub, UT_MAX (*(s3-1), UT_MAX (*s3, *(s3+1))));

            *dest++ = ub;
            s1++;
            s2++;
            s3++;
        }
    } else {
        while (s1 < stop) {
            ub = 255;

            ub = UT_MIN (ub, UT_MIN (*(s1-1), UT_MIN (*s1, *(s1+1))));
            ub = UT_MIN (ub, UT_MIN (*(s2-1), UT_MIN (*s2, *(s2+1))));
            ub = UT_MIN (ub, UT_MIN (*(s3-1), UT_MIN (*s3, *(s3+1))));

            *dest++ = ub;
            s1++;
            s2++;
            s3++;
        }
    }
}

/* Apply a three by three pixel dilatation/erosion filter to one
   scan line of a Float32 channel of an image. */

static void dler_Float_scan (Float32 *s1, Float32 *s2, Float32 *s3, 
                             Int32 len, Int32 mode, Float32 *dest)
{
    Float32 f;
    Float32 *stop = s1 + len;

    if (mode == 1) {
        while (s1 < stop) {
            f = 0.0;

            f = UT_MAX (f, UT_MAX (*(s1-1), UT_MAX (*s1, *(s1+1))));
            f = UT_MAX (f, UT_MAX (*(s2-1), UT_MAX (*s2, *(s2+1))));
            f = UT_MAX (f, UT_MAX (*(s3-1), UT_MAX (*s3, *(s3+1))));

            *dest++ = f;
            s1++;
            s2++;
            s3++;
        }
    } else {
        while (s1 < stop) {
            f = 1.0;

            f = UT_MIN (f, UT_MIN (*(s1-1), UT_MIN (*s1, *(s1+1))));
            f = UT_MIN (f, UT_MIN (*(s2-1), UT_MIN (*s2, *(s2+1))));
            f = UT_MIN (f, UT_MIN (*(s3-1), UT_MIN (*s3, *(s3+1))));

            *dest++ = f;
            s1++;
            s2++;
            s3++;
        }
    }
}

/* Apply a three by three pixel dilatation/erosion filter to all channels 
   in image "img" for which the current drawing mask is UT_True. */

static UT_Bool dilat_eros (IP_ImageId img, Int32 n, Int32 mode)
{
    Int32       i, x, y, chn;
    UInt8       *ubuf1 = NULL, *ubuf2 = NULL;
    Float32     *fbuf1 = NULL, *fbuf2 = NULL;
    FF_ImgDesc  *pimg;

    pimg = &img->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   *in_s1, *in_s2, *in_s3, *pix1, *pix2, *pix3, 
                        *buf1_ptr, *buf2_ptr, *out_p1, *out_p2;
                if (!ubuf1 || !ubuf2) {
                    if (!(ubuf1 = UT_MemTempArray (pimg->desc.width, UInt8)) ||
                        !(ubuf2 = UT_MemTempArray (pimg->desc.width, UInt8)))
                        return UT_False;
                }
                for (i = 0; i < n; i++) {
                    in_s1 = FF_ImgUByteRow (pimg, chn, 0);
                    in_s2 = FF_ImgUByteRow (pimg, chn, 1);
                    in_s3 = FF_ImgUByteRow (pimg, chn, 2);
                    pix1 = in_s1 + 1; 
                    pix2 = in_s2 + 1;
                    pix3 = in_s3 + 1;
                    buf1_ptr = ubuf1;
                    buf2_ptr = ubuf2;
                    out_p1 = buf1_ptr + 1;
                    out_p2 = buf2_ptr + 1;
                    dler_UByte_scan (pix1, pix2, pix3, 
                                     pimg->desc.width-2, mode, out_p1);

                    for (y = 2; y < pimg->desc.height -1; ++y) {
                        in_s1 = in_s2;
                        in_s2 = in_s3;
                        pix1 = in_s1 + 1; 
                        pix2 = in_s2 + 1;
                        in_s3 = FF_ImgUByteRow (pimg, chn, y+1);
                        pix3 = in_s3 + 1;
                        dler_UByte_scan (pix1, pix2, pix3, 
                                         pimg->desc.width-2, mode, out_p2);

                        for (x = 1; x < pimg->desc.width -1; ++x)
                            *FF_ImgUBytePixel (pimg, chn, x, y-1) = buf1_ptr[x];
                        UT_SWAP (buf1_ptr, buf2_ptr, UInt8 *);
                        out_p2 = buf2_ptr;
                    }
                    for (x = 1; x < pimg->desc.width -1; ++x)
                        *FF_ImgUBytePixel (pimg, chn, x, y-1) = buf1_ptr[x];
                }
                break;
            }

            case FF_ImgFmtTypeFloat: {
                Float32 *in_s1, *in_s2, *in_s3, *pix1, *pix2, *pix3, 
                        *buf1_ptr, *buf2_ptr, *out_p1, *out_p2;
                if (!fbuf1 || !fbuf2) {
                    if (!(fbuf1 = UT_MemTempArray (pimg->desc.width, Float32)) ||
                        !(fbuf2 = UT_MemTempArray (pimg->desc.width, Float32)))
                        return UT_False;
                }
                for (i = 0; i < n; i++) {
                    in_s1 = FF_ImgFloatRow (pimg, chn, 0);
                    in_s2 = FF_ImgFloatRow (pimg, chn, 1);
                    in_s3 = FF_ImgFloatRow (pimg, chn, 2);
                    pix1 = in_s1 + 1; 
                    pix2 = in_s2 + 1;
                    pix3 = in_s3 + 1;
                    buf1_ptr = fbuf1;
                    buf2_ptr = fbuf2;
                    out_p1 = buf1_ptr + 1;
                    out_p2 = buf2_ptr + 1;
                    dler_Float_scan (pix1, pix2, pix3, 
                                     pimg->desc.width-2, mode, out_p1);

                    for (y = 2; y < pimg->desc.height -1; ++y) {
                        in_s1 = in_s2;
                        in_s2 = in_s3;
                        pix1 = in_s1 + 1; 
                        pix2 = in_s2 + 1;
                        in_s3 = FF_ImgFloatRow (pimg, chn, y+1);
                        pix3 = in_s3 + 1;
                        dler_Float_scan (pix1, pix2, pix3, 
                                         pimg->desc.width-2, mode, out_p2);

                        for (x = 1; x < pimg->desc.width -1; ++x)
                            *FF_ImgFloatPixel (pimg, chn, x, y-1) = buf1_ptr[x];
                        UT_SWAP (buf1_ptr, buf2_ptr, Float32 *);
                        out_p2 = buf2_ptr;
                    }
                    for (x = 1; x < pimg->desc.width -1; ++x)
                        *FF_ImgFloatPixel (pimg, chn, x, y-1) = buf1_ptr[x];
                }
                break;
            }
        }
    }
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_ChangeGamma
 *
 *      Usage:          Change one channel of an image by gamma-filtering
 *                      and linear scaling of the pixel data.
 *
 *      Synopsis:       UT_Bool IP_ChangeGamma
 *                              (IP_ImageId img, FF_ImgChanType chn,
 *                              Float32 gamma, Float32 scl, Float32 off)
 *
 *      Description:    Channel "chn" of every pixel in image "img" is changed:
 *                      The pixel's orignal value, "v", is clipped (values less
 *                      than 0.0 are set to 0.0, values greater than 1.0 are
 *                      set to 1.0). The clipped value is gamma filtered, using
 *                      gamma factor "gamma" and then scaled by "scl". "Off" is
 *                      added to the scaled value. The result replaces the
 *                      pixels original value:
 *
 *                      v = pow (max (0, min (v, 1)), 1/gamma) * scl + off
 *
 *      Return value:   UT_True if successful, else UT_False.
 *
 *      See also:
 *
 ***************************************************************************/

UT_Bool IP_ChangeGamma
        (IP_ImageId img, FF_ImgChanType chn,
         Float32 gamma, Float32 scl, Float32 off)
{
    UT_Bool     success;

    if (chn < 0 || chn >= FF_NumImgChanTypes) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_imgchan, 
                      FF_ImgGetChanTypeName (chn));
        return UT_False;
    }
    if (gamma < MINGAMMA || gamma > MAXGAMMA) {
        UT_ErrSetNum (UT_ErrParamRange, str_range_inclusive,
                     str_gamma, MINGAMMA, MAXGAMMA);
        return UT_False;
    }
    success = change_gamma (img, chn, gamma, scl, off);
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           IP_RemapValues
 *
 *      Usage:          Change one channel of an image by remapping the
 *                      pixel data with a lookup table.
 *
 *      Synopsis:       UT_Bool IP_RemapValues
 *                              (IP_ImageId img, FF_ImgChanType chn,
 *                              Int32 lutsize, const Float32 lut[])
 *
 *      Description:    The original value, "v", of channel "chn" of every
 *                      pixel in image "img" is mapped into a new value which
 *                      replaces the original pixel data:
 *
 *                              v = f (v * (lutsize - 1))
 *
 *                      "F" is a continuous and piecewise linear function so
 *                      that for all x in the range from 0 to (lutsize - 1)
 *
 *                              f (floor(x)) == lut[floor(x)].
 *
 *      Return value:   UT_True if successful, else UT_False.
 *
 *      See also:
 *
 ***************************************************************************/

UT_Bool IP_RemapValues
        (IP_ImageId img, FF_ImgChanType chn,
         Int32 lutsize, const Float32 lut[])
{
    UT_Bool success;

    if (chn < 0 || chn >= FF_NumImgChanTypes) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_imgchan, 
                      FF_ImgGetChanTypeName (chn));
        return UT_False;
    }
    if (lutsize < 2) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_inotless, str_lutsize, 2);
        return UT_False;
    }
    remap_values (img, chn, lutsize, lut);
    success = UT_True;
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           IP_Derivates
 *
 *      Usage:          Computes differences between adjacent pixels.
 *
 *      Synopsis:       UT_Bool IP_Derivates
 *                              (IP_ImageId srcimg,
 *                              FF_ImgChanType srcchn,
 *                              IP_ImageId destimg,
 *                              FF_ImgChanType hchn, FF_ImgChanType vchn,
 *                              UT_Bool wrap)
 *                      
 *      Description:    For each pixel in channel "srcchn" of image "srcimg",
 *                      IP_Derivates computes the differences between the pixel
 *                      and its left and bottom neighbors. The result is stored
 *                      in the "hchn" and "vchn" channels of image "destimg":
 *
 *                          hchn (x, y) = 0.5 + srcchn (x, y) - srcchn (x-1, y)
 *                          vchn (x, y) = 0.5 + srcchn (x, y) - srcchn (x, y-1)
 *
 *                      The "wrap" flag determines the differences computed
 *                      at the left and bottom borders of the source image.
 *                      If "wrap" is set to UT_True, IP_Derivates assumes that
 *
 *                          srcchn (-1, y) = srcchn (w-1, y) and
 *                          srcchn (x, -1) = srcchn (x, h-1),
 *
 *                      where "w" and "h" are the width and height of "srcimg",
 *                      in pixels.
 *                      If "wrap" is set to UT_False, IP_Derivates assumes that
 *
 *                          srcchn (-1, y) = srcchn (0, y) and
 *                          srcchn (x, -1) = srcchn (x, 0).
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:
 *
 ***************************************************************************/

UT_Bool IP_Derivates (IP_ImageId srcimg, FF_ImgChanType srcchn,
                 IP_ImageId destimg, FF_ImgChanType hchn, FF_ImgChanType vchn,
                 UT_Bool wrap)
{
    UT_MemState memstate;
    UT_Bool     success;

    if (srcchn < 0 || srcchn >= FF_NumImgChanTypes) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_imgchan, 
                      FF_ImgGetChanTypeName (srcchn));
        return UT_False;
    }
    if (hchn < 0 || hchn >= FF_NumImgChanTypes) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_imgchan, 
                      FF_ImgGetChanTypeName (hchn));
        return UT_False;
    }
    if (vchn < 0 || vchn >= FF_NumImgChanTypes) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_imgchan, 
                      FF_ImgGetChanTypeName (vchn));
        return UT_False;
    }
    memstate = UT_MemRemember ();
    success = hvdiff (srcimg, srcchn, destimg, hchn, vchn, wrap);
    UT_MemRestore (memstate);
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           IP_Median
 *
 *      Usage:          Apply a median filter to an image.
 *
 *      Synopsis:       UT_Bool IP_Median (IP_ImageId img)
 *
 *      Description:    A three pixel wide median filter is applied,
 *                      first horizontally, then vertically, to all
 *                      channels in image "img" for which the current
 *                      drawing mask is "UT_True".
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:
 *
 ***************************************************************************/

UT_Bool IP_Median (IP_ImageId img)
{
    UT_Bool success;

    median_h (img);
    median_v (img);
    success = UT_True;
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           IP_MedianSequence
 *
 *      Usage:          Applay a median filter to an image sequence.
 *
 *      Synopsis:       UT_Bool IP_MedianSequence
 *                              (IP_ImageId srcimg1,
 *                              IP_ImageId srcimg2,
 *                              IP_ImageId srcimg3,
 *                              IP_ImageId destimg)
 *                      
 *      Description:    A median filter is applied to all channels in image
 *                      sequence "srcimg1", "srcimg2", "srcimg3", for which
 *                      the current drawing mask is "UT_True":
 *                      Corresponding pixel values in "srcimg1", "srcimg2" and
 *                      "srcimg3" are compared and the second-largest value is
 *                      stored in "destimg".
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:       IP_Median
 *
 ***************************************************************************/

UT_Bool IP_MedianSequence
        (IP_ImageId srcimg1,
         IP_ImageId srcimg2,
         IP_ImageId srcimg3,
         IP_ImageId destimg)
{
    UT_Bool success;

    median_sequence (srcimg1, srcimg2, srcimg3, destimg);
    success = UT_True;
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           IP_Dilatation
 *
 *      Usage:          Dilatation (expansion) of bright areas of an image.
 *
 *      Synopsis:       UT_Bool IP_Dilatation (IP_ImageId img, Int32 n)
 *
 *      Description:    A three by three pixel wide dilatation filter is 
 *                      applied "n" times to all channels in image "img" 
 *                      for which the current drawing mask is "UT_True".
 *
 *                      The dilatation filter expands bright areas in 
 *                      the image and eliminates dark pixels.
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:       IP_Erosion
 *
 ***************************************************************************/

UT_Bool IP_Dilatation (IP_ImageId img, Int32 n)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = dilat_eros (img, n, 1);
    UT_MemRestore (memstate);
    return success;
}

/***************************************************************************
 *[@e
 *      Name:           IP_Erosion
 *
 *      Usage:          Erosion (contraction) of bright areas of an image.
 *
 *      Synopsis:       UT_Bool IP_Erosion (IP_ImageId img, Int32 n)
 *
 *      Description:    A three by three pixel wide erosion filter is 
 *                      applied "n" times to all channels in image "img" 
 *                      for which the current drawing mask is "UT_True".
 *
 *                      The erosion filter contracts bright areas in 
 *                      the image and eliminates bright pixels.
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:       IP_Dilatation
 *
 ***************************************************************************/

UT_Bool IP_Erosion (IP_ImageId img, Int32 n)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = dilat_eros (img, n, -1);
    UT_MemRestore (memstate);
    return success;
}

static UT_Bool threshold_filter (IP_ImageId src, IP_ImageId dest,
                                 Float32 t, UT_Bool overwrite)
{
    Int32 x, y, xmax, ymax, ithresh,
          vcenter, vbottom, vright;
    UInt8 *rbottom, *rcenter,
          *gbottom, *gcenter,
          *bbottom, *bcenter;
    FF_ImgDesc *psrc, *pdest;

    /* if (src->pimg.channel[chn]) != FF_ImgChanTypeRed ... */

    psrc  = &src->pimg;
    pdest = &dest->pimg;

    ithresh = (Int32) ((3.0 * 255.0) * t);
    xmax = src->pimg.desc.width  - 2;
    ymax = src->pimg.desc.height - 1;
    for (y = ymax; y >= 1; --y) {
        rbottom = FF_ImgUByteRow (psrc, FF_ImgChanTypeRed,   y-1);
        rcenter = FF_ImgUByteRow (psrc, FF_ImgChanTypeRed,   y);
        gbottom = FF_ImgUByteRow (psrc, FF_ImgChanTypeGreen, y-1);
        gcenter = FF_ImgUByteRow (psrc, FF_ImgChanTypeGreen, y);
        bbottom = FF_ImgUByteRow (psrc, FF_ImgChanTypeBlue,  y-1);
        bcenter = FF_ImgUByteRow (psrc, FF_ImgChanTypeBlue,  y);
        for (x = 1; x <= xmax; ++x) {
            vcenter = rcenter[x]   + gcenter[x]   + bcenter[x];
            vright  = rcenter[x+1] + gcenter[x+1] + bcenter[x+1] - vcenter;
            vbottom = rbottom[x]   + gbottom[x]   + bbottom[x]   - vcenter;
            if (ithresh < UT_ABS (vright) || ithresh < UT_ABS (vbottom)) {
                *FF_ImgUBytePixel (pdest, FF_ImgChanTypeRed,   x, y) = 255;
                *FF_ImgUBytePixel (pdest, FF_ImgChanTypeGreen, x, y) = 255;
                *FF_ImgUBytePixel (pdest, FF_ImgChanTypeBlue,  x, y) = 255;
            } else {
                if (overwrite) {
                    *FF_ImgUBytePixel (pdest, FF_ImgChanTypeRed,   x, y) = 0;
                    *FF_ImgUBytePixel (pdest, FF_ImgChanTypeGreen, x, y) = 0;
                    *FF_ImgUBytePixel (pdest, FF_ImgChanTypeBlue,  x, y) = 0;
                }
            }
        }
    }
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_Threshold
 *
 *      Usage:          Apply a threshold filter to an image.
 *
 *      Synopsis:       UT_Bool IP_Threshold 
 *                              (IP_ImageId src, IP_ImageId dest,
 *                               Float32 t, UT_Bool overwrite)
 *
 *      Description:    TBD
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:       IP_CutOff
 *
 ***************************************************************************/

UT_Bool IP_Threshold (IP_ImageId src, IP_ImageId dest, 
                      Float32 t, UT_Bool overwrite)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = threshold_filter (src, dest, t, overwrite);
    UT_MemRestore (memstate);
    return success;
}

static UT_Bool cutoff_filter (IP_ImageId img, Float32 tmin, Float32 tmax)
{
    Int32 x, y, xmax, ymax, chn;
    UInt8 t1, t2;
    UInt8 *upix;
    Float32 *fpix;
    FF_ImgDesc *pimg;

    pimg = &img->pimg;
    xmax = pimg->desc.width;
    ymax = pimg->desc.height;
    for (chn = 0; chn < FF_NumImgChanTypes; ++chn) {
        if (!IP_StateDrawMask[chn]) {
            continue;
        }
        switch (img->pimg.channel[chn]) {
            case FF_ImgFmtTypeNone:
                break;
            case FF_ImgFmtTypeUByte: {
                t1 = (UInt8) (tmin * 255.0);
                t2 = (UInt8) (tmax * 255.0);
                for (y = 0; y < ymax; ++y) {
                    for (x = 0; x < xmax; ++x) {
                        upix = FF_ImgUBytePixel (pimg, chn, x, y);
                        if (*upix <= t1) {
                            *upix = 0;
                        } else if (*upix >= t2) {
                            *upix = 255;
                        }
                    }
                }
                break;
            }
            case FF_ImgFmtTypeFloat: {
                for (y = 0; y < ymax; ++y) {
                    for (x = 0; x < xmax; ++x) {
                        fpix = FF_ImgFloatPixel (pimg, chn, x, y);
                        if (*fpix <= tmin) {
                            *fpix = 0.0;
                        } else if (*fpix >= tmax) {
                            *fpix = 1.0;
                        }
                    }
                }
                break;
            }
            default:
                UT_ErrFatal ("cutoff_filter", "unknown pixel data format");
        }
    }
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_CutOff
 *
 *      Usage:          Apply a cut-off filter to an image.
 *
 *      Synopsis:       UT_Bool IP_CutOff 
 *                              (IP_ImageId img, Float32 tmin,
 *                               Float32 tmax)
 *
 *      Description:    TBD
 *
 *      Return value:   UT_True if succesfull, else UT_False.
 *
 *      See also:       IP_Threshold
 *
 ***************************************************************************/

UT_Bool IP_CutOff (IP_ImageId img, Float32 tmin, Float32 tmax)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = cutoff_filter (img, tmin, tmax);
    UT_MemRestore (memstate);
    return success;
}

static UT_Bool markNonZeroPixels (IP_ImageId src, IP_ImageId dest,
                                  Int32 threshold, Int32 *numMarkedPixels)
{
    Int32 x, y, xmax, ymax;
    Int32 numMarked = 0;
    UInt8 *redSrc,  *greenSrc,  *blueSrc;
    UInt8 *redDest, *greenDest, *blueDest;
    FF_ImgDesc *psrc, *pdest;

    /* TODO Check if (src->pimg.channel[chn]) != FF_ImgChanTypeRed ... */

    psrc  = &src->pimg;
    pdest = &dest->pimg;

    xmax = src->pimg.desc.width  - 1;
    ymax = src->pimg.desc.height - 1;
    for (y = ymax; y >= 0; --y) {
        redSrc    = FF_ImgUByteRow (psrc,  FF_ImgChanTypeRed,   y);
        greenSrc  = FF_ImgUByteRow (psrc,  FF_ImgChanTypeGreen, y);
        blueSrc   = FF_ImgUByteRow (psrc,  FF_ImgChanTypeBlue,  y);
        redDest   = FF_ImgUByteRow (pdest, FF_ImgChanTypeRed,   y);
        greenDest = FF_ImgUByteRow (pdest, FF_ImgChanTypeGreen, y);
        blueDest  = FF_ImgUByteRow (pdest, FF_ImgChanTypeBlue,  y);
        for (x = 0; x <= xmax; ++x) {
            if (threshold < *redSrc || threshold < *greenSrc || threshold < *blueSrc) {
                *redDest   = IP_StateDrawColorUByte[FF_ImgChanTypeRed];
                *greenDest = IP_StateDrawColorUByte[FF_ImgChanTypeGreen];
                *blueDest  = IP_StateDrawColorUByte[FF_ImgChanTypeBlue];
                numMarked++;
            } else {
                *redDest   = 0;
                *greenDest = 0;
                *blueDest  = 0;
            }
            redSrc++;
            greenSrc++;
            blueSrc++;
            redDest++;
            greenDest++;
            blueDest++;
        }
    }
    *numMarkedPixels = numMarked;
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_MarkNonZeroPixels
 *
 *      Usage:          Mark all non-zero pixels of an image.
 *
 *      Synopsis:       UT_Bool IP_MarkNonZeroPixels 
 *                              (IP_ImageId src, IP_ImageId dest,
 *                               Int32 threshold, Int32 *numMarkedPixels)
 *
 *      Description:    TBD
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:       IP_CutOff
 *                      IP_Threshold
 *
 ***************************************************************************/

UT_Bool IP_MarkNonZeroPixels (IP_ImageId src, IP_ImageId dest, 
                              Int32 threshold, Int32 *numMarkedPixels)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = markNonZeroPixels (src, dest, threshold, numMarkedPixels);
    UT_MemRestore (memstate);
    return success;
}

static UT_Bool diffImgs (IP_ImageId src1, IP_ImageId src2, IP_ImageId dest,
                         UT_Bool useAbsValue)
{
    Int32 y, xmin, ymin, xmax, ymax;
    Int32 chn;
    FF_ImgDesc *psrc1, *psrc2, *pdest;

    /* TODO Check if (src->pimg.channel[chn]) != FF_ImgChanTypeRed ... */
    /* TODO Add alignment type: Left, Right, Bottom, Top */

    psrc1 = &src1->pimg;
    psrc2 = &src2->pimg;
    pdest = &dest->pimg;

    xmin = UT_MIN (src1->pimg.desc.width  - 1, src2->pimg.desc.width  - 1);
    ymin = UT_MIN (src1->pimg.desc.height - 1, src2->pimg.desc.height - 1);
    xmax = UT_MAX (src1->pimg.desc.width  - 1, src2->pimg.desc.width  - 1);
    ymax = UT_MAX (src1->pimg.desc.height - 1, src2->pimg.desc.height - 1);

    IP_BlankRect (dest, 0, ymin, xmin, ymax);
    IP_BlankRect (dest, xmin+1, 0, xmax, ymax);
    for (chn = 0; chn < FF_NumImgChanTypes; ++chn) {
        if (!IP_StateDrawMask[chn]) {
            continue;
        }
        switch (pdest->channel[chn]) {
            case FF_ImgFmtTypeNone: {
                break;
            }
            case FF_ImgFmtTypeUByte: {
                UInt8 *src1data, *src2data;
                UInt8 *destdata, *deststop;
                for (y = ymin; y >= 0; --y) {
                    src1data = FF_ImgUByteRow (psrc1, chn, y);
                    src2data = FF_ImgUByteRow (psrc2, chn, y);
                    destdata = FF_ImgUByteRow (pdest, chn, y);
                    deststop = destdata + xmin + 1;
                    while (destdata < deststop) {
                        *destdata = UT_ABS (*src1data - *src2data);
                        src1data++;
                        src2data++;
                        destdata++;
                    }
                }
                break;
            }
            case FF_ImgFmtTypeFloat: {
                UT_ErrSetNum (UT_ErrParamInvalid, str_chn_ubyte);
                return UT_False;
                break;
            }
            default:
                UT_ErrFatal ("diffImgs", "unknown pixel data format");
        }
    }
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_DifferenceImage
 *
 *      Usage:          Generate the difference image of two images.
 *
 *      Synopsis:       UT_Bool IP_DifferenceImage
 *                              (IP_ImageId src1, IP_ImageId src2, IP_ImageId dest,
 *                               UT_Bool useAbsValue)
 *
 *      Description:    TBD
 *
 *      Return value:   UT_True if succesful, else UT_False.
 *
 *      See also:       IP_MarkNonZeroPixels
 *                      IP_CopyRect
 *
 ***************************************************************************/

UT_Bool IP_DifferenceImage (IP_ImageId src1, IP_ImageId src2, IP_ImageId dest,
                            UT_Bool useAbsValue)
{
    UT_MemState memstate;
    UT_Bool     success;

    memstate = UT_MemRemember ();
    success = diffImgs (src1, src2, dest, useAbsValue);
    UT_MemRestore (memstate);
    return success;
}
