GEMini Atari/GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso/files/graphics/utility/mgif35s/process.c

C/C++ Source or Header | 1991-06-14 | 37.9 KB | 2,237 lines

#undef DBL_LOOP            /* def for for(x..., for(y... loops */
#undef INL_SORT            /* def to inline piksrt() in median() */

/*
 *    process.c - all the image processing functions.
 */

static char *sccsid = "@(#) process.c 1.3 91/6/9 rosenkra\0\0               ";

#include <stdio.h>
#include "mgif.h"

/*
 *    local functions
 */
int        larger ();        /* frame processes... */
int        zoom ();
int        smaller ();
int        cut ();
int        rotate ();
int        mirror ();
int        convolve ();        /* area processes... */
int        blur ();
int        median ();
int        piksrt ();
int        logscale ();        /* point processes... */
int        Log2x10 ();
int        contrast ();
int        brighten ();
int        invert ();
int        threshold ();
int        histeq ();
int        redistribute ();
int        copyrast ();        /* others... */

/*
 *    external functions
 */
extern int    do_hist ();
extern int    copyrast ();


/*------------------------------*/
/*    larger            */
/*------------------------------*/

#define L_BOTH        0
#define L_HOR        1
#define L_VERT        2

int larger (pras, w, h, opt, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        opt;        /* checked by caller! */
uchar_t           *ptrans;
{

/*
 *    enlarge an image (so far 320x200 -> 640x400). return 1 if
 *    successful, else 0.
 */

    register uchar_t       *pr;
    register uchar_t       *pt;
    register long        x;
    register long        y;
    long            x2;
    long            y2;
    register long        wold;
    long            hold;
    register long        w2;
    long            h2;


    /*
     *   fail if enlarged image will be too big for buffer
     */
#if 0
    if ((w > 320) || (h > 200))
        return (0);
#endif
    switch (opt)
    {
    case L_BOTH:
        if (4L * (long) w * (long) h > MAXIMG)
            return (0);
        break;
    case L_HOR:
        if (2L * (long) w * (long) h > MAXIMG)
            return (0);
        break;
    case L_VERT:
        if (2L * (long) w * (long) h > MAXIMG)
            return (0);
        break;
    }



    /*
     *   here we always do the transform in situ (in ptrans)...
     */
    if (pras == ptrans)
    {
        /*
         *   do in place...
         */
        pr = pras;
        pt = pras;
    }
    else
    {
        /*
         *   copy to new image
         */
        pr = pras;
        pt = ptrans;

        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
        pr = ptrans;
    }


    wold = w;
    hold = h;
    switch (opt)
    {
    case L_BOTH:
        /*
         *   set new sizes
         */
        w2 = wold*2L;
        h2 = hold*2L;


        /*
         *   shift existing data:
         *
         *      0123456789
         *    0|**********|            *=existing pixels
         *    1|**********|            +=new pixels
         *    .
         *    .
         *         |
         *         v
         *                1111111111
         *      01234567890123456789
         *    0|                    |
         *    1| + + + + + + + + + +|
         *    2|                    |
         *    3| + + + + + + + + + +|
         *    .
         *    .
         */
        for (y = hold-1; y >= 0; y--)
        {
            for (x = wold-1; x >= 0; x--)
            {
                x2 = (2*x)+1;
                y2 = (2*y)+1;
                pt[y2*w2 + x2] = pt[y*wold + x];
            }
        }


        /*
         *   now make intermediate points within row (average):
         *
         *                1111111111
         *      01234567890123456789
         *    0|                    |        *=existing pixels
         *    1|+*+*+*+*+*+*+*+*+*+*|        +=new pixels
         *    2|                    |
         *    3|+*+*+*+*+*+*+*+*+*+*|
         *    .
         *    .
         */
        for (y = 1; y < h2; y += 2)
        {
            /* first point is duplicate of second */
            pt[y*w2 + 0] = pt[y*w2 + 1];

            /* others are average */
            for (x = 2; x < w2; x += 2)
            {
                pt[y*w2 + x] = (uchar_t) (
                     ((uint_t) pt[y*w2 + (x-1)]
                    + (uint_t) pt[y*w2 + (x+1)]) / 2);
            }
        }


        /*
         *   first row is duplicate of second row:
         *
         *                1111111111
         *      01234567890123456789
         *    0|++++++++++++++++++++|<---
         *    1|********************|        *=existing pixels
         *    2|                    |        +=new pixels
         *    3|********************|
         *    .
         *    .
         */
        for (x = 0; x < w2; x++)
        {
            pt[x] = pt[w2 + x];
        }


        /*
         *   finally, make new rows (average):
         *
         *                1111111111
         *      01234567890123456789
         *    0|********************|        *=existing pixels
         *    1|********************|        +=new pixels
         *    2|++++++++++++++++++++|<---
         *    3|********************|
         *    .
         *    .
         */
        for (y = 2; y < h2; y += 2)
        {
            for (x = 0; x < w2; x++)
            {
                pt[y*w2 + x] = (uchar_t) (
                     ((uint_t) pt[(y-1)*w2 + x]
                    + (uint_t) pt[(y+1)*w2 + x]) / 2);
            }
        }
        break;


    case L_HOR:
        /*
         *   set new sizes
         */
        w2 = wold*2L;
        h2 = hold;


        /*
         *   shift existing data:
         *
         *      0123456789
         *    0|**********|            *=existing pixels
         *    1|**********|            +=new pixels
         *         .
         *         .
         *         |
         *         v
         *                1111111111
         *      01234567890123456789
         *    0| + + + + + + + + + +|
         *    1| + + + + + + + + + +|
         *    .
         *    .
         */
        for (y = hold-1; y >= 0; y--)
        {
            for (x = wold-1; x >= 0; x--)
            {
                x2 = (2*x)+1;
                pt[y*w2 + x2] = pt[y*wold + x];
            }
        }


        /*
         *   now make intermediate points within row (average):
         *
         *                1111111111
         *      01234567890123456789
         *    0|+*+*+*+*+*+*+*+*+*+*|        *=existing pixels
         *    1|+*+*+*+*+*+*+*+*+*+*|        +=new pixels
         *    .
         *    .
         */
        for (y = 0; y < hold; y++)
        {
            /* first point is duplicate of second */
            pt[y*w2 + 0] = pt[y*w2 + 1];

            /* others are average */
            for (x = 2; x < w2; x += 2)
            {
                pt[y*w2 + x] = (uchar_t) (
                     ((uint_t) pt[y*w2 + (x-1)]
                    + (uint_t) pt[y*w2 + (x+1)]) / 2);
            }
        }
        break;

    case L_VERT:
        /*
         *   set new sizes
         */
        w2 = wold;
        h2 = hold*2L;


        /*
         *   shift existing data:
         *
         *      0123456789
         *    0|**********|            *=existing pixels
         *    1|**********|            +=new pixels
         *    .
         *    .
         *         |
         *         v
         *      0123456789
         *    0|          |
         *    1|++++++++++|
         *    2|          |
         *    3|++++++++++|
         *    .
         *    .
         */
        for (y = hold-1; y >= 0; y--)
        {
            for (x = wold-1; x >= 0; x--)
            {
                y2 = (2*y)+1;
                pt[y2*wold + x] = pt[y*wold + x];
            }
        }


        /*
         *   first row is duplicate of second row:
         *
         *      0123456789
         *    0|++++++++++|<---
         *    1|**********|            *=existing pixels
         *    2|          |            +=new pixels
         *    3|**********|
         *    .
         *    .
         */
        for (x = 0; x < wold; x++)
        {
            pt[x] = pt[wold + x];
        }


        /*
         *   finally, make new rows (average):
         *
         *      0123456789
         *    0|**********|        *=existing pixels
         *    1|**********|        +=new pixels
         *    2|++++++++++|<---
         *    3|**********|
         *    .
         *    .
         */
        for (y = 2; y < h2; y += 2)
        {
            for (x = 0; x < wold; x++)
            {
                pt[y*wold + x] = (uchar_t) (
                     ((uint_t) pt[(y-1)*wold + x]
                    + (uint_t) pt[(y+1)*wold + x]) / 2);
            }
        }
        break;
    }


    return (1);
}




/*------------------------------*/
/*    zoom            */
/*------------------------------*/
int zoom (pras, w, h, xstart, ystart, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        xstart;
int        ystart;
uchar_t           *ptrans;
{

/*
 *    zoom in and enlarge an image (so far 2x) starting at xstart,ystart.
 *    return 1 if successful, else 0.
 */

    register uchar_t       *pr;
    register uchar_t       *pt;
    register long        x;
    register long        y;
    long            x2;
    long            y2;
    register long        wold;
    long            hold;
    register long        w2;
    long            h2;



    /*
     *   enlarged image should ALWAYS fit in buffer!
     *
     *   here we always do the transform in situ (in ptrans)...
     */
    if (pras == ptrans)
    {
        /*
         *   do in place. pr points to "cut" area, pt to start
         */
        pr = &pras[xstart + (ystart-1)*w];
        pt = pras;

        /*
         *   shift smaller image to ul corner of pt
         */
        if ((xstart != 0) && (ystart != 0))
        {
            for (y = 0L; y < h/2; y++)
            {
                for (x = 0L; x < w/2; x++)
                {
                    pt[x + w*(y-1)/2] = pr[x + w*(y-1)];
                }
            }
        }
    }
    else
    {
        /*
         *   copy to new image
         */
        pr = &pras[xstart + (ystart-1)*w];
        pt = ptrans;

        /*
         *   first copy orig image to transform image
         */
        for (y = 0L; y < h/2; y++)
        {
            for (x = 0L; x < w/2; x++)
            {
                pt[x + w*(y-1)/2] = pr[x + w*(y-1)];
            }
        }
    }



    /*
     *   set new sizes
     */
    wold = w/2L;
    hold = h/2L;
    w2   = wold*2L;
    h2   = hold*2L;


    /*
     *   shift existing data:
     *
     *      0123456789
     *    0|**********|            *=existing pixels
     *    1|**********|            +=new pixels
     *    .
     *    .
     *         |
     *         v
     *                1111111111
     *      01234567890123456789
     *    0|                    |
     *    1| + + + + + + + + + +|
     *    2|                    |
     *    3| + + + + + + + + + +|
     *    .
     *    .
     */
    for (y = hold-1; y >= 0; y--)
    {
        for (x = wold-1; x >= 0; x--)
        {
            x2 = (2*x)+1;
            y2 = (2*y)+1;
            pt[y2*w2 + x2] = pt[y*wold + x];
        }
    }


    /*
     *   now make intermediate points within row (average):
     *
     *                1111111111
     *      01234567890123456789
     *    0|                    |        *=existing pixels
     *    1|+*+*+*+*+*+*+*+*+*+*|        +=new pixels
     *    2|                    |
     *    3|+*+*+*+*+*+*+*+*+*+*|
     *    .
     *    .
     */
    for (y = 1; y < h2; y += 2)
    {
        /* first point is duplicate of second */
        pt[y*w2 + 0] = pt[y*w2 + 1];

        /* others are average */
        for (x = 2; x < w2; x += 2)
        {
            pt[y*w2 + x] = (uchar_t) (
                 ((uint_t) pt[y*w2 + (x-1)]
                + (uint_t) pt[y*w2 + (x+1)]) / 2);
        }
    }


    /*
     *   first row is duplicate of second row:
     *
     *                1111111111
     *      01234567890123456789
     *    0|++++++++++++++++++++|<---
     *    1|********************|        *=existing pixels
     *    2|                    |        +=new pixels
     *    3|********************|
     *    .
     *    .
     */
    for (x = 0; x < w2; x++)
    {
        pt[x] = pt[w2 + x];
    }


    /*
     *   finally, make new rows (average):
     *
     *                1111111111
     *      01234567890123456789
     *    0|********************|        *=existing pixels
     *    1|********************|        +=new pixels
     *    2|++++++++++++++++++++|<---
     *    3|********************|
     *    .
     *    .
     */
    for (y = 2; y < h2; y += 2)
    {
        for (x = 0; x < w2; x++)
        {
            pt[y*w2 + x] = (uchar_t) (
                 ((uint_t) pt[(y-1)*w2 + x]
                + (uint_t) pt[(y+1)*w2 + x]) / 2);
        }
    }

    return (1);
}




/*------------------------------*/
/*    smaller            */
/*------------------------------*/

#define S_BOTH        0
#define S_HOR        1
#define S_VERT        2

int smaller (pras, w, h, opt, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        opt;        /* checked by caller! */
uchar_t           *ptrans;
{

/*
 *    reduce an image by half. return 1 if successful, else 0.
 */

    register uchar_t       *pr;
    register uchar_t       *pt;
    register long        x;
    register long        y;
    register long        wold;
    register long        w2;
    long            hold;
    long            h2;


    /*
     *   there is a limit to what is practical...
     */
    if ((w < 10) || (h < 10))
        return (0);


    if (pras == ptrans)
    {
        /*
         *   do in place...
         */
        pr = pras;
        pt = pras;
    }
    else
    {
        /*
         *   copy to new image
         */
        pr = pras;
        pt = ptrans;
    }


    /*
     *   do it. skip odd lines and pixels
     */
    wold = w;
    hold = h;
    switch (opt)
    {
    case S_BOTH:
        /*
         *   set new sizes
         */
        w2 = wold/2L;
        h2 = hold/2L;

        for (y = 0; y < h2; y++)
        {
            for (x = 0; x < w2; x++)
            {
                pt[y*w2 + x] = pr[2*y*wold + x*2];
            }
        }
        break;

    case S_HOR:
        /*
         *   set new sizes
         */
        w2 = wold/2L;
        h2 = hold;

        for (y = 0; y < h2; y++)
        {
            for (x = 0; x < w2; x++)
            {
                pt[y*w2 + x] = pr[y*wold + x*2];
            }
        }
        break;

    case S_VERT:
        /*
         *   set new sizes
         */
        w2 = wold;
        h2 = hold/2L;

        for (y = 0; y < h2; y++)
        {
            for (x = 0; x < w2; x++)
            {
                pt[y*w2 + x] = pr[2*y*wold + x];
            }
        }
        break;
    }

    return (1);
}




/*------------------------------*/
/*    cut            */
/*------------------------------*/
int cut (pras, w, h, x1, y1, x2, y2, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        x1, y1;
int        x2, y2;
uchar_t           *ptrans;
{

/*
 *    cut out a section of the image from pras to ptrans
 */

    register uchar_t       *pr;
    register uchar_t       *pt;
    register long        x;
    register long        y;
    register long        wold;
    register long        w2;
    register long        h2;


    /*
     *   check for bad coords...
     */
    if ((x1 >= x2) || (y1 >= y2))
        return (0);
    if ((x1 < 0) || (x2 < 0) || (y1 < 0) || (y2 < 0))
        return (0);
    if ((x2 >= w) || (y2 >= h))
        return (0);
#if 0
    if ((x2 > 639) || (y2 > 399))
        return (0);
    if ((x2 - x1) > w)        /* don't fail just limit the cut */
        x2 = x1 + w;
    if ((y2 - y1) > h)
        y2 = y1 + h;
#endif


    /*
     *   set up pointers...
     */
    if (pras == ptrans)
    {
        /*
         *   do in place...
         */
        pr = pras;
        pt = pras;
    }
    else
    {
        /*
         *   copy to new image
         */
        pr = pras;
        pt = ptrans;
    }


    /*
     *   do it...
     */
    wold = w;
    w2   = x2 - x1 + 1;
    h2   = y2 - y1 + 1;
    for (y = 0; y < h2; y++)
    {
        for (x = 0; x < w2; x++)
        {
            pt[y*w2 + x] = pr[(y + y1)*wold + (x + x1)];
        }
    }

    return (1);
}




/*------------------------------*/
/*    rotate            */
/*------------------------------*/
int rotate (pras, w, h, angle, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        angle;        /* +90, -90, -180 or +180 */
uchar_t           *ptrans;
{

/*
 *    rotate an image +/-90 or +/-180. does NOT work in place (yet)
 *
 *    based on (rotating about center):
 *
 *        xnew = (x * cos(angle)) + (y * sin(angle))
 *        ynew = (y * cos(angle)) - (x * sin(angle))
 *
 *    note: 45's should be relatively simple:
 *
 *        +45:    xnew = (.707 * x) + (.707 * y)
 *            ynew = (.707 * y) - (.707 * x)
 *
 *        -45:    xnew = (.707 * x) - (.707 * y)
 *            ynew = (.707 * y) + (.707 * x)
 */

    register uchar_t       *pr;
    register uchar_t       *pt;
    register long        xnew;
    register long        ynew;
    register long        wold;
    register long        hold;


    /*
     *   set up pointers...
     */
    if (pras == ptrans)
    {
        /*
         *   can't do in place (yet)...
         */
        return (0);
    }


    /*
     *   copy to new image
     */
    pr = pras;
    pt = ptrans;


    /*
     *   do it...
     */
    wold = (long) w;
    hold = (long) h;
    switch (angle)
    {
    case 90:
        /*
         *   counter clockwise:
         *
         *    xnew = y        ->    y = xnew
         *    ynew = wold - x - 1        x = wold - ynew - 1
         *    wnew = hold
         *    hnew = wold
         */
        for (ynew = 0L; ynew < wold; ynew++)
        {
            for (xnew = 0L; xnew < hold; xnew++)
            {
                pt[ynew*hold + xnew] = pr[xnew*wold + (wold-ynew-1)];
            }
        }
        break;

    case -90:
        /*
         *   clockwise...
         *
         *    xnew = hold - y - 1    ->    y = hold - xnew - 1
         *    ynew = x            x = ynew
         *    wnew = hold
         *    hnew = wold
         */
        for (ynew = 0L; ynew < wold; ynew++)
        {
            for (xnew = 0L; xnew < hold; xnew++)
            {
                pt[ynew*hold + xnew] = pr[(hold-xnew-1)*wold + ynew];
            }
        }
        break;

    case 180:
    case -180:
        /*
         *   clockwise or counter clockwise...
         *
         *    xnew = wold - x - 1    ->    x = wold - xnew - 1
         *    ynew = hold - y - 1        y = hold - ynew - 1
         *    wnew = wold
         *    hnew = hold
         */
        for (ynew = 0L; ynew < hold; ynew++)
        {
            for (xnew = 0L; xnew < wold; xnew++)
            {
                pt[ynew*wold + xnew] = pr[(hold-ynew-1)*wold + (wold-xnew-1)];
            }
        }
        break;

    default:
        /*
         *   undefined angle, fail...
         */
        return (0);
    }

    return (1);
}




/*------------------------------*/
/*    mirror            */
/*------------------------------*/
int mirror (pras, w, h, opt, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        opt;        /* 0=vert mirror, 1=horiz mirror */
uchar_t           *ptrans;
{

/*
 *    mirror image.
 *
 *    vert mirror is:
 *
 *         _______|_______
 *        |    |    |
 *        |    |    |
 *        |    |    |
 *        |_______|_______|
 *            |
 *
 *    horiz mirror is:
 *
 *         _______________
 *        |        |
 *          __|_______________|__
 *        |        |
 *        |_______________|
 */

    register uchar_t       *pr;
    register uchar_t       *pt;
    register long        x;
    register long        y;
    register uchar_t    ctemp;


    /*
     *   set up pointers...
     */
    if (pras == ptrans)
    {
        /*
         *   do in place...
         */
        pr = pras;
        pt = pras;

        if (opt)
        {
            /*
             *   horizontal mirror
             */
            for (y = 0L; y < h/2; y++)
            {
                for (x = 0L; x < w; x++)
                {
                    ctemp           = pt[y*w + x];
                    pt[y*w + x]     = pt[(h-y)*w + x];
                    pt[(h-y)*w + x] = ctemp;
                }
            }
        }
        else
        {
            /*
             *   vertical mirror
             */
            for (y = 0L; y < h; y++)
            {
                for (x = 0L; x < w/2; x++)
                {
                    ctemp           = pt[y*w + x];
                    pt[y*w + x]     = pt[y*w + (w-x)];
                    pt[y*w + (w-x)] = ctemp;
                }
            }
        }
    }
    else
    {
        /*
         *   copy to new image
         */
        pr = pras;
        pt = ptrans;

        if (opt)
        {
            /*
             *   horizontal mirror
             */
            for (y = 0L; y < h; y++)
            {
                for (x = 0L; x < w; x++)
                {
                    pt[y*w + x] = pr[(h-y-1)*w + x];
                }
            }
        }
        else
        {
            /*
             *   vertical mirror
             */
            for (y = 0L; y < h; y++)
            {
                for (x = 0L; x < w; x++)
                {
                    pt[y*w + x] = pr[y*w + (w-x-1)];
                }
            }
        }
    }

    return (1);
}




/*------------------------------*/
/*    convolve        */
/*------------------------------*/

/*
 *    built-in convolution kernels
 *
 *    k[0..8] is kernel (always 3x3)
 *    k[9] is scaling flag: -1=divide, 1=multiply, 0=no scaling
 *    k[10] is scale
 *
 *    example:
 *        for lp1 with flag -1 (divide) and scale 9, use
 *
 *            p = ((p0 * k0) +...+ (p8 * k8)) / 9
 *
 *        for hp2 with flag 0 (no scaling), use
 *
 *            p = (p0 * k0) +...+ (p8 * k8)
 */
int    lp1[11] = { 1, 1, 1,            /* low pass */
            1, 1, 1,
            1, 1, 1,    -1, 9};
int    lp2[11] = { 1, 1, 1,
            1, 2, 1,
            1, 1, 1,    -1, 10};
int    lp3[11] = { 1, 2, 1,
            2, 4, 2,
            1, 2, 1,    -1, 16};

int    hp1[11] = {-1,-1,-1,            /* hi pass */
           -1, 9,-1,
           -1,-1,-1,    0, 1};
int    hp2[11] = { 0,-1, 0,
           -1, 5,-1,
            0,-1, 0,    0, 1};
int    hp3[11] = { 1,-2, 1,
           -2, 5,-2,
            1,-2, 1,    0, 1};

                        /* shift edge */
int    se1[11] = { 0, 0, 0,            /* vertical */
           -1, 1, 0,
            0, 0, 0,    0, 1};
int    se2[11] = { 0,-1, 0,            /* horizontal */
            0, 1, 0,
            0, 0, 0,    0, 1};
int    se3[11] = {-1, 0, 0,            /* hor and vert */
            0, 1, 0,
            0, 0, 0,    0, 1};

                        /* gradient edge */
int    ge1[11] = { 1, 1, 1,            /* north */
            1,-2, 1,
           -1,-1,-1,    0, 1};
int    ge2[11] = { 1, 1, 1,            /* northeast */
           -1,-2, 1,
           -1,-1, 1,    0, 1};
int    ge3[11] = {-1, 1, 1,            /* east */
           -1,-2, 1,
           -1, 1, 1,    0, 1};
int    ge4[11] = {-1,-1, 1,            /* southeast */
           -1,-2, 1,
            1, 1, 1,    0, 1};
int    ge5[11] = {-1,-1,-1,            /* south */
            1,-2, 1,
            1, 1, 1,    0, 1};
int    ge6[11] = { 1,-1,-1,            /* southwest */
            1,-2,-1,
            1, 1, 1,    0, 1};
int    ge7[11] = { 1, 1,-1,            /* west */
            1,-2,-1,
            1, 1,-1,    0, 1};
int    ge8[11] = { 1, 1, 1,            /* northwest */
            1,-2,-1,
            1,-1,-1,    0, 1};

int    le1[11] = { 0, 1, 0,            /* laplace edge */
            1,-4, 1,
            0, 1, 0,    0, 1};
int    le2[11] = {-1,-1,-1,
           -1, 8,-1,
           -1,-1,-1,    0, 1};
int    le3[11] = {-1,-1,-1,
           -1, 9,-1,
           -1,-1,-1,    0, 1};
int    le4[11] = { 1,-2, 1,
           -2, 4,-2,
            1,-2, 1,    0, 1};

int convolve (pras, w, h, opt, kernel, ptrans)
uchar_t           *pras;
register int    w;
int        h;
int        opt;        /* 0=user, 1,2,...=built-in */
int           *kernel;        /* user kernel, if any */
uchar_t           *ptrans;
{

/*
 *    convolute an image (3x3). return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = sum (k  * p ) for i pixels in 3x3 neighborhood
 *                       i    i
 */

    uchar_t               *pr;
    uchar_t               *pt;
    register int        val;
    register long        x;
    register long        y;
    register uchar_t       *ps;
    register int           *pk;


    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   set pointer to user or built-in kernels, depending on opt
     */
    switch (opt)
    {
    case USER_KERN:        pk = kernel;        break;

    case LP1_KERN:        pk = lp1;        break;
    case LP2_KERN:        pk = lp2;        break;
    case LP3_KERN:        pk = lp3;        break;

    case HP1_KERN:        pk = hp1;        break;
    case HP2_KERN:        pk = hp2;        break;
    case HP3_KERN:        pk = hp3;        break;

    case SE1_KERN:        pk = se1;        break;
    case SE2_KERN:        pk = se2;        break;
    case SE3_KERN:        pk = se3;        break;

    case GE1_KERN:        pk = ge1;        break;
    case GE2_KERN:        pk = ge2;        break;
    case GE3_KERN:        pk = ge3;        break;
    case GE4_KERN:        pk = ge4;        break;
    case GE5_KERN:        pk = ge5;        break;
    case GE6_KERN:        pk = ge6;        break;
    case GE7_KERN:        pk = ge7;        break;
    case GE8_KERN:        pk = ge8;        break;

    case LE1_KERN:        pk = le1;        break;
    case LE2_KERN:        pk = le2;        break;
    case LE3_KERN:        pk = le3;        break;
    case LE4_KERN:        pk = le4;        break;

    default:        return (0);
    }


    /*
     *   do it...
     */
    ps = pt;
    for (y = 1; y < h-1; y++)
    {
        for (x = 1; x < w-1; x++)
        {
            val  = ((uint_t) ps[(y-1)*w+(x-1)]) * pk[0];
            val += ((uint_t) ps[(y-1)*w+(x  )]) * pk[1];
            val += ((uint_t) ps[(y-1)*w+(x+1)]) * pk[2];
            val += ((uint_t) ps[(y  )*w+(x-1)]) * pk[3];
            val += ((uint_t) ps[(y  )*w+(x  )]) * pk[4];
            val += ((uint_t) ps[(y  )*w+(x+1)]) * pk[5];
            val += ((uint_t) ps[(y+1)*w+(x-1)]) * pk[6];
            val += ((uint_t) ps[(y+1)*w+(x  )]) * pk[7];
            val += ((uint_t) ps[(y+1)*w+(x+1)]) * pk[8];

            if ((pk[9] < 0) && pk[10])
                val /= pk[10];
            else if (pk[9] > 0)
                val *= pk[10];

            if (val < 0)
                val = 0;
            else if (val > 255)
                val = 255;

            ps[(y*w) + x] = (uchar_t) val;
        }
    }

    return (1);
}




/*------------------------------*/
/*    blur            */
/*------------------------------*/
int blur (pras, w, h, ptrans)
uchar_t           *pras;
register int    w;
int        h;
uchar_t           *ptrans;
{

/*
 *    blur an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = average of neighboring cells (3x3)
 */

    uchar_t               *pr;
    uchar_t               *pt;
    register long        x;
    register long        y;
    register uchar_t       *ps;



    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   do it...
     */
    ps = pt;
    for (y = 1; y < h-1; y++)
    {
        for (x = 1; x < w-1; x++)
        {
            ps[(y*w) + x] = (uchar_t) (
                    ( (uint_t) ps[(y-1)*w + (x-1)]
                    + (uint_t) ps[(y-1)*w + (x  )]
                    + (uint_t) ps[(y-1)*w + (x+1)]
                    + (uint_t) ps[(y  )*w + (x-1)]
                    + (uint_t) ps[(y  )*w + (x  )]
                    + (uint_t) ps[(y  )*w + (x+1)]
                    + (uint_t) ps[(y+1)*w + (x-1)]
                    + (uint_t) ps[(y+1)*w + (x  )]
                    + (uint_t) ps[(y+1)*w + (x+1)]) / 9);
        }
    }

    return (1);
}




/*------------------------------*/
/*    median            */
/*------------------------------*/
int median (pras, w, h, ptrans)
uchar_t           *pras;
int        w;
int        h;
uchar_t           *ptrans;
{

/*
 *    median filter an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = median of neighboring cells
 */

    uchar_t               *pr;
    register uchar_t       *pt;
    register long        x;
    register long        y;
#ifdef INL_SORT
    register int        i;
    register uint_t        a;
#endif
    register uint_t           *pl;
    uint_t            l[9];


    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   do it...
     */
    pl = l;
    for (y = 1; y < h-1; y++)
    {
        for (x = 1; x < w-1; x++)
        {
            pl[0] = (uint_t) pt[(y-1)*w + (x-1)];
            pl[1] = (uint_t) pt[(y-1)*w + (x  )];
            pl[2] = (uint_t) pt[(y-1)*w + (x+1)];
            pl[3] = (uint_t) pt[(y  )*w + (x-1)];
            pl[4] = (uint_t) pt[(y  )*w + (x  )];
            pl[5] = (uint_t) pt[(y  )*w + (x+1)];
            pl[6] = (uint_t) pt[(y+1)*w + (x-1)];
            pl[7] = (uint_t) pt[(y+1)*w + (x  )];
            pl[8] = (uint_t) pt[(y+1)*w + (x+1)];
#ifdef INL_SORT
            /* unroll j loop, too... */
            for (a=pl[1], i=0; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[2], i=1; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[3], i=2; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[4], i=3; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[5], i=4; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[6], i=5; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[7], i=6; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
            for (a=pl[8], i=7; i>=0 && pl[i]>a; i--) pl[i+1]=pl[i];
            pl[i+1] = a;
#else
            piksrt (9, l);
#endif
            pt[(y*w) + x] = pl[4];
        }
    }

    return (1);
}




/*------------------------------*/
/*    piksrt            */
/*------------------------------*/
piksrt (num, arr)
int        num;
uint_t           *arr;
{

/*
 *    straight insertion sort (num. rec. in C, pp 243). is N^2 sort,
 *    probably ok for N=9. shell sort is 3x faster for N=9.
 */

    register uint_t           *parr;
    register int        n;
    register int        i;
    register int        j;
    register uint_t        a;

    n    = num;
    parr = arr;
    for (j = 1; j < n; j++)
    {
        for (a=parr[j], i=j-1; i>=0 && parr[i]>a; i--)
            parr[i+1] = parr[i];

        parr[i+1] = a;
    }
}




/*------------------------------*/
/*    findmedian        */
/*------------------------------*/
uint_t findmedian (num, arr)
int        num;
uint_t           *arr;
{

/*
 *    find median of an array of numbers. this should be faster than sort.
 *    there are pathalogical cases: if all are the same or all but 1 or 2
 *    are the same...
 */

    register uint_t           *parr;
    register int        n;
    register int        i;
    register int        j;
    register int        count;
    register uint_t        a;

    parr = arr;
    n    = num;
    for (j = 0; j < n; j++)
    {
        for (a = parr[j], count = 0, i = 0; i < n; i++)
        {
            if (i == j)
                continue;
            if (a > parr[i])
                count++;
        }
        if (count == (n>>1))
            return (a);
    }

}




/*------------------------------*/
/*    logscale        */
/*------------------------------*/
int logscale (pras, w, h, ptrans)
uchar_t           *pras;
int        w;
int        h;
uchar_t           *ptrans;
{

/*
 *    apply log scaling to an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = maxval * log(old[i]) / log(maxval)
 *    or
 *        new[i] = maxval * log2(old[i]) / log2(maxval)
 */

    uchar_t               *pr;
    uchar_t               *pt;
    register ulong_t    maxval;
    register ulong_t    l2maxval;
    register ulong_t    l2pr;
    uchar_t            uctmp;
    long            x;
    long            y;
    long            val;
    register long        ii;
    register long        lim;
    register uchar_t       *ps;



    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   find maxval. here it is max value a pixel can have, 255.
     */
#if 1
    maxval = 0L;
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            uctmp = pt[(y*w) + x];
            if (uctmp > maxval)
                maxval = (ulong_t) uctmp;
        }
    }
#endif
/*!!!    maxval = 255L;*/


    /*
     *   do it...
     */
#if 1
    l2maxval = (ulong_t) Log2x10 ((uint_t) maxval);
    if (l2maxval == 0)
        return (0);
#endif
/*!!!    l2maxval = 80;*/        /* approx 10 times log2(255) */


#ifdef DBL_LOOP
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            l2pr          = (ulong_t) Log2x10 ((uint_t) pt[(y*w)+x]);
            pt[(y*w) + x] = (uchar_t) ((maxval * l2pr) / l2maxval);
        }
    }
#else
    lim = (long) w * (long) h;
    for (ps = pt, ii = 0L; ii < lim; ii++)
    {
        l2pr  = (ulong_t) Log2x10 ((uint_t) *ps);
        *ps++ = (uchar_t) ((maxval * l2pr) / l2maxval);
    }
#endif

    return (1);
}





/*------------------------------*/
/*    Log2x10            */
/*------------------------------*/
int Log2x10(x)
uint_t    x;
{
/*
 *    VERY approximate log base 2 times 10. basically finds MS bit...
 */
    if (x < 256)        goto lobyte;
    if (x & 0x8000)        return (150);
    if (x & 0x4000)        return (140);
    if (x & 0x2000)        return (130);
    if (x & 0x1000)        return (120);
    if (x & 0x0800)        return (110);
    if (x & 0x0400)        return (100);
    if (x & 0x0200)        return (90);
    if (x & 0x0100)        return (80);
lobyte:
    if (x == 0x00ff)
        return (80);
    if (x == 0x0000)
        return (1);
    if (x & 0x0080)
    {
        if (x & 0x0040)
            return (78);
        else if (x & 0x0020)
            return (75);
        else
            return (71);
    }
    if (x & 0x0040)
    {
        if (x & 0x0020)
            return (68);
        else if (x & 0x0010)
            return (65);
        else
            return (61);
    }
    if (x & 0x0020)
    {
        if (x & 0x0010)
            return (58);
        else if (x & 0x0008)
            return (55);
        else
            return (51);
    }
    if (x & 0x0010)
    {
        if (x & 0x0008)
            return (48);
        else if (x & 0x0004)
            return (45);
        else
            return (41);
    }
    if (x & 0x0008)
    {
        if (x & 0x0004)
            return (38);
        else if (x & 0x0002)
            return (35);
        else
            return (31);
    }
    if (x & 0x0004)
    {
        if (x & 0x0002)
            return (28);
        else if (x & 0x0001)
            return (25);
        else
            return (21);
    }
    if (x & 0x0002)
    {
        if (x & 0x0001)
            return (18);
        else
            return (13);
    }
    if (x & 0x0001)        return (10);
    return (0);
}




/*------------------------------*/
/*    contrast        */
/*------------------------------*/
int contrast (pras, w, h, thresh, hist, ptrans)
uchar_t           *pras;
int        w;
int        h;
long        thresh;
long           *hist;
uchar_t           *ptrans;
{

/*
 *    apply contrast expansion to an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = maxval * (old[i] - lo) / (hi - lo)
 *
 *        where lo to hi is new contrast range
 */

    uchar_t               *pr;
    uchar_t               *pt;
    ulong_t            maxval;
    uchar_t            uctmp;
    long            x;
    long            y;
    long            hi_lo;
    long            newlo;
    long            newhi;
    int            i;
    register long        val;
    register long        rng;
    register long        lim;
    register long        ii;
    register uchar_t       *ps;


    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   get current histogram...
     */
    if (do_hist (pt, w, h, hist) == 0)
        return (0);


    /*
     *    find intensities on either side of hist where pixel count
     *    exceeds thresh
     */
    for (i = 0; i < HISTSIZ; i++)
    {
        if (hist[i] > thresh)
            break;
    }
    newlo = i;
    for (i = HISTSIZ-1; i > newlo; i--)
    {
        if (hist[i] > thresh)
            break;
    }
    newhi = i;


    /*
     *   do it. set pixels below lower threshold 0 and those above to 255.
     *   contrast expand pixels between...
     */
#ifdef DBL_LOOP
    rng = newhi - newlo + 1;
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            val = (ulong_t) pt[(y*w) + x];
            if (val < newlo)
                pt[(y*w) + x] = 0;
            else if (val > newhi)
                pt[(y*w) + x] = 255;
            else
            {
                pt[(y*w) + x] = (uchar_t) ((255L * (val - newlo)) / rng);
            }
        }
    }
#else
    rng = newhi - newlo + 1;
    lim = (long) w * (long) h;
    for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
    {
        val = (ulong_t) *ps;
        if (val < newlo)
            *ps = 0;
        else if (val > newhi)
            *ps = 255;
        else
            *ps = (uchar_t) ((255L * (val - newlo)) / rng);
    }
#endif

    return (1);
}




/*------------------------------*/
/*    brighten        */
/*------------------------------*/
int brighten (pras, w, h, brite, ptrans)
uchar_t           *pras;
int        w;
int        h;
register int    brite;
uchar_t           *ptrans;
{

/*
 *    brighten an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = old[i] + brite
 */

    uchar_t               *pr;
    uchar_t               *pt;
    long            x;
    long            y;
    int            sign;
    register uint_t        val;
    register long        ii;
    register long        lim;
    register uchar_t       *ps;



    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    sign = 0;
    if (brite < 0)
    {
        sign = 1;
        brite = -brite;
    }


    /*
     *   do it...
     */
    if (sign)
    {
#ifdef DBL_LOOP
        for (y = 0; y < h; y++)
        {
            for (x = 0; x < w; x++)
            {
                val = (uint_t) pt[y*w + x];
                if (val < brite)
                    val  = 0;
                else
                    val -= brite;
                pt[(y*w) + x] = (uchar_t) val;
            }
        }
#else
        lim = (long) w * (long) h;
        for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
        {
            if ((uint_t) *ps < brite)
                *ps  = 0;
            else
                *ps -= brite;
        }
#endif
    }
    else
    {
#ifdef DBL_LOOP
        for (y = 0; y < h; y++)
        {
            for (x = 0; x < w; x++)
            {
                val = (uint_t) ((uint_t) pt[y*w + x] + brite);
                if (val > 255)
                    val = 255;
                pt[(y*w) + x] = (uchar_t) val;
            }
        }
#else
        lim = (long) w * (long) h;
        for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
        {
            val = (uint_t) ((uint_t) *ps + brite);
            if (val > 255)
                val = 255;
            *ps = (uchar_t) val;
        }
#endif
    }

    return (1);
}




/*------------------------------*/
/*    invert            */
/*------------------------------*/
int invert (pras, w, h, thresh, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        thresh;
uchar_t           *ptrans;
{

/*
 *    invert (negate) an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] = 255 - old[i]
 */

    uchar_t               *pr;
    uchar_t               *pt;
    long            x;
    long            y;
    int            sign;
    register long        lim;
    register long        ii;
    register int        thr;
    register uchar_t       *ps;


    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    sign = 0;
    if (thresh < 0)
    {
        sign   = 1;
        thresh = -thresh;
    }


    /*
     *   do it...
     */
#ifdef DBL_LOOP
    if (thresh == 0)
    {
        for (y = 0; y < h; y++)
        {
            for (x = 0; x < w; x++)
            {
                pt[(y*w) + x] = 255 - pt[(y*w) + x];
            }
        }
    }
    else if (sign)
    {
        /*
         *   invert only below thresh
         */
        for (y = 0; y < h; y++)
        {
            for (x = 0; x < w; x++)
            {
                if ((uint_t) pt[(y*w) + x] < thresh)
                    pt[(y*w) + x] = 255 - pt[(y*w) + x];
            }
        }
    }
    else
    {
        /*
         *   invert only above thresh
         */
        for (y = 0; y < h; y++)
        {
            for (x = 0; x < w; x++)
            {
                if ((uint_t) pt[(y*w) + x] > thresh)
                    pt[(y*w) + x] = 255 - pt[(y*w) + x];
            }
        }
    }
#else
    lim = (long) h * (long) w;
    thr = thresh;
    if (thr == 0)
    {
        for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
        {
            *ps = 255 - *ps;
        }
    }
    else if (sign)
    {
        /*
         *   invert only below thresh
         */
        for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
        {
            if ((uint_t) *ps < thr)
                *ps = 255 - *ps;
        }
    }
    else
    {
        /*
         *   invert above below thresh
         */
        for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
        {
            if ((uint_t) *ps > thr)
                *ps = 255 - *ps;
        }
    }
#endif

    return (1);
}




/*------------------------------*/
/*    threshold        */
/*------------------------------*/
int threshold (pras, w, h, thresh, ptrans)
uchar_t           *pras;
int        w;
int        h;
int        thresh;
uchar_t           *ptrans;
{

/*
 *    threshold an image. return 1 if successful, else 0.
 *
 *    the basic algorithm is:
 *
 *        new[i] =   0 if old[i] <= thresh
 *        new[i] = 255 if old[i] >  thresh
 */

    uchar_t               *pr;
    uchar_t               *pt;
    long            x;
    long            y;
    uint_t            val;
    register long        ii;
    register long        lim;
    register uchar_t       *ps;


    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   do it...
     */
#ifdef DBL_LOOP
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            val = (uint_t) pt[(y*w) + x];
            if (val > thresh)
                pt[(y*w) + x] = 255;
            else
                pt[(y*w) + x] = 0;
        }
    }
#else
    lim = (long) h * (long) w;
    for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
    {
        if ((uint_t) *ps > thresh)
            *ps = 255;
        else
            *ps = 0;
    }
#endif

    return (1);
}




/*------------------------------*/
/*    histeq            */
/*------------------------------*/
int histeq (pras, w, h, hist, ptrans)
uchar_t           *pras;
int        w;
int        h;
long           *hist;
uchar_t           *ptrans;
{

/*
 *    histogram equalization of image. return 1 if successful, else 0.
 */

    long            p_r[HISTSIZ];
    long            s_k[HISTSIZ];
    long            s_kapprx[HISTSIZ];
    uchar_t               *pr;
    uchar_t               *pt;
    int            i;
    long            count;



    pr = pras;
    pt = ptrans;

    if (pras != ptrans)
    {
        /*
         *   first copy orig image to transform image
         */
        copyrast (pr, w, h, pt);
    }


    /*
     *   get current histogram...
     */
    if (do_hist (pt, w, h, hist) == 0)
        return (0);


    /*
     *   find prob density function (p_r), transformation function (s_k),
     *   and new distribution. note the 1000 to avoid floats...
     */
    count = (long) w * (long) h;
    for (i = 0; i < HISTSIZ; i++)
    {
        p_r[i] = (hist[i] * 1000L) / count;
    }
    for (s_k[0] = p_r[0], i = 1; i < HISTSIZ; i++)
    {
        s_k[i] = s_k[i-1] + p_r[i];
    }
    for (i = 0; i < HISTSIZ; i++)
    {
        /* the 500 is for rounding to nearest int... */
        s_kapprx[i] = ((s_k[i] * 255L) + 500L) / 1000L;
    }

#if 0
    printf ("\n    i      hist       p_r       s_k  s_kapprx\n");
    for (count = 0, i = 0; i < HISTSIZ; i++)
    {
        printf ("%5d%10ld%10ld%10ld%10ld\n",
            i,hist[i],p_r[i],s_k[i],s_kapprx[i]);
    }
#endif

    /*
     *   now redistribute. s_kapprx will contain the new value for
     *   a given pixel intensity
     */
    redistribute (pt, w, h, s_kapprx);

    return (1);
}




/*------------------------------*/
/*    redistribute        */
/*------------------------------*/
int redistribute (pras, w, h, s_kapprx)
uchar_t           *pras;
int        w;
int        h;
register long  *s_kapprx;
{

/*
 *    redistribute pixels. each pixel set to value given by s_kapprx
 */

    register uchar_t       *pr;
    register long        ii;
    register long        lim;


    lim = (long) w * (long) h;
    for (pr = pras, ii = 0L; ii < lim; ii++, pr++)
    {
        *pr = (uchar_t) s_kapprx[(uint_t) *pr];
    }
}




/*------------------------------*/
/*    copyrast        */
/*------------------------------*/
int copyrast (ps, w, h, pd)
uchar_t           *ps;        /* source */
int        w;
int        h;
uchar_t           *pd;        /* dest */
{

/*
 *    copy an image. return 1 if successful, else 0.
 */

#ifdef DBL_LOOP
    long        x;
    long        y;

    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            pd[(y*w) + x] = ps[(y*w) + x];
        }
    }
#else
    register long    ii;
    register long    lim;

    lim = (long) w * (long) h;
    for (ii = 0L; ii < lim; ii++)
        *pd++ = *ps++;
#endif
    return (1);
}