colorspace.c

Go to the documentation of this file.

/*====================================================================*
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
 -
 -  Redistribution and use in source and binary forms, with or without
 -  modification, are permitted provided that the following conditions
 -  are met:
 -  1. Redistributions of source code must retain the above copyright
 -     notice, this list of conditions and the following disclaimer.
 -  2. Redistributions in binary form must reproduce the above
 -     copyright notice, this list of conditions and the following
 -     disclaimer in the documentation and/or other materials
 -     provided with the distribution.
 -
 -  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 -  ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 -  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 -  A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL ANY
 -  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 -  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 -  PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 -  PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 -  OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 -  NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 -  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *====================================================================*/

#include <string.h>
#include <math.h>
#include "allheaders.h"

#ifndef  NO_CONSOLE_IO
#define  DEBUG_HISTO       0
#define  SLOW_CUBE_ROOT    0
#endif  /* ~NO_CONSOLE_IO */

    /* Functions used in xyz <--> lab conversions */
static l_float32 lab_forward(l_float32 v);
static l_float32 lab_reverse(l_float32 v);


/*---------------------------------------------------------------------------*
 *                  Colorspace conversion between RGB and HSB                *
 *---------------------------------------------------------------------------*/
PIX *
pixConvertRGBToHSV(PIX  *pixd,
                   PIX  *pixs)
{
l_int32    w, h, d, wpl, i, j, rval, gval, bval, hval, sval, vval;
l_uint32  *line, *data;
PIXCMAP   *cmap;

    PROCNAME("pixConvertRGBToHSV");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixd && pixd != pixs)
        return (PIX *)ERROR_PTR("pixd defined and not inplace", procName, pixd);

    d = pixGetDepth(pixs);
    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("not cmapped or rgb", procName, pixd);

    if (!pixd)
        pixd = pixCopy(NULL, pixs);

    cmap = pixGetColormap(pixd);
    if (cmap) {   /* just convert the colormap */
        pixcmapConvertRGBToHSV(cmap);
        return pixd;
    }

        /* Convert RGB image */
    pixGetDimensions(pixd, &w, &h, NULL);
    wpl = pixGetWpl(pixd);
    data = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        for (j = 0; j < w; j++) {
            extractRGBValues(line[j], &rval, &gval, &bval);
            convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
            line[j] = (hval << 24) | (sval << 16) | (vval << 8);
        }
    }

    return pixd;
}


PIX *
pixConvertHSVToRGB(PIX  *pixd,
                   PIX  *pixs)
{
l_int32    w, h, d, wpl, i, j, rval, gval, bval, hval, sval, vval;
l_uint32   pixel;
l_uint32  *line, *data;
PIXCMAP   *cmap;

    PROCNAME("pixConvertHSVToRGB");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixd && pixd != pixs)
        return (PIX *)ERROR_PTR("pixd defined and not inplace", procName, pixd);

    d = pixGetDepth(pixs);
    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("not cmapped or hsv", procName, pixd);

    if (!pixd)
        pixd = pixCopy(NULL, pixs);

    cmap = pixGetColormap(pixd);
    if (cmap) {   /* just convert the colormap */
        pixcmapConvertHSVToRGB(cmap);
        return pixd;
    }

        /* Convert HSV image */
    pixGetDimensions(pixd, &w, &h, NULL);
    wpl = pixGetWpl(pixd);
    data = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        for (j = 0; j < w; j++) {
            pixel = line[j];
            hval = pixel >> 24;
            sval = (pixel >> 16) & 0xff;
            vval = (pixel >> 8) & 0xff;
            convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
            composeRGBPixel(rval, gval, bval, line + j);
        }
    }

    return pixd;
}


l_ok
convertRGBToHSV(l_int32   rval,
                l_int32   gval,
                l_int32   bval,
                l_int32  *phval,
                l_int32  *psval,
                l_int32  *pvval)
{
l_int32    minrg, maxrg, min, max, delta;
l_float32  h;

    PROCNAME("convertRGBToHSV");

    if (phval) *phval = 0;
    if (psval) *psval = 0;
    if (pvval) *pvval = 0;
    if (!phval || !psval || !pvval)
        return ERROR_INT("&hval, &sval, &vval not all defined", procName, 1);

    minrg = L_MIN(rval, gval);
    min = L_MIN(minrg, bval);
    maxrg = L_MAX(rval, gval);
    max = L_MAX(maxrg, bval);
    delta = max - min;

    *pvval = max;
    if (delta == 0) {  /* gray; no chroma */
        *phval = 0;
        *psval = 0;
    } else {
        *psval = (l_int32)(255. * (l_float32)delta / (l_float32)max + 0.5);
        if (rval == max)  /* between magenta and yellow */
            h = (l_float32)(gval - bval) / (l_float32)delta;
        else if (gval == max)  /* between yellow and cyan */
            h = 2. + (l_float32)(bval - rval) / (l_float32)delta;
        else  /* between cyan and magenta */
            h = 4. + (l_float32)(rval - gval) / (l_float32)delta;
        h *= 40.0;
        if (h < 0.0)
            h += 240.0;
        if (h >= 239.5)
            h = 0.0;
        *phval = (l_int32)(h + 0.5);
    }

    return 0;
}


l_ok
convertHSVToRGB(l_int32   hval,
                l_int32   sval,
                l_int32   vval,
                l_int32  *prval,
                l_int32  *pgval,
                l_int32  *pbval)
{
l_int32   i, x, y, z;
l_float32 h, f, s;

    PROCNAME("convertHSVToRGB");

    if (prval) *prval = 0;
    if (pgval) *pgval = 0;
    if (pbval) *pbval = 0;
    if (!prval || !pgval || !pbval)
        return ERROR_INT("&rval, &gval, &bval not all defined", procName, 1);

    if (sval == 0) {  /* gray */
        *prval = vval;
        *pgval = vval;
        *pbval = vval;
    } else {
        if (hval < 0 || hval > 240)
            return ERROR_INT("invalid hval", procName, 1);
        if (hval == 240)
            hval = 0;
        h = (l_float32)hval / 40.;
        i = (l_int32)h;
        f = h - i;
        s = (l_float32)sval / 255.;
        x = (l_int32)(vval * (1. - s) + 0.5);
        y = (l_int32)(vval * (1. - s * f) + 0.5);
        z = (l_int32)(vval * (1. - s * (1. - f)) + 0.5);
        switch (i)
        {
        case 0:
            *prval = vval;
            *pgval = z;
            *pbval = x;
            break;
        case 1:
            *prval = y;
            *pgval = vval;
            *pbval = x;
            break;
        case 2:
            *prval = x;
            *pgval = vval;
            *pbval = z;
            break;
        case 3:
            *prval = x;
            *pgval = y;
            *pbval = vval;
            break;
        case 4:
            *prval = z;
            *pgval = x;
            *pbval = vval;
            break;
        case 5:
            *prval = vval;
            *pgval = x;
            *pbval = y;
            break;
        default:  /* none possible */
            return 1;
        }
    }

    return 0;
}


l_ok
pixcmapConvertRGBToHSV(PIXCMAP  *cmap)
{
l_int32   i, ncolors, rval, gval, bval, hval, sval, vval;

    PROCNAME("pixcmapConvertRGBToHSV");

    if (!cmap)
        return ERROR_INT("cmap not defined", procName, 1);

    ncolors = pixcmapGetCount(cmap);
    for (i = 0; i < ncolors; i++) {
        pixcmapGetColor(cmap, i, &rval, &gval, &bval);
        convertRGBToHSV(rval, gval, bval, &hval, &sval, &vval);
        pixcmapResetColor(cmap, i, hval, sval, vval);
    }
    return 0;
}


l_ok
pixcmapConvertHSVToRGB(PIXCMAP  *cmap)
{
l_int32   i, ncolors, rval, gval, bval, hval, sval, vval;

    PROCNAME("pixcmapConvertHSVToRGB");

    if (!cmap)
        return ERROR_INT("cmap not defined", procName, 1);

    ncolors = pixcmapGetCount(cmap);
    for (i = 0; i < ncolors; i++) {
        pixcmapGetColor(cmap, i, &hval, &sval, &vval);
        convertHSVToRGB(hval, sval, vval, &rval, &gval, &bval);
        pixcmapResetColor(cmap, i, rval, gval, bval);
    }
    return 0;
}


PIX *
pixConvertRGBToHue(PIX  *pixs)
{
l_int32    w, h, d, wplt, wpld;
l_int32    i, j, rval, gval, bval, hval, minrg, min, maxrg, max, delta;
l_float32  fh;
l_uint32   pixel;
l_uint32  *linet, *lined, *datat, *datad;
PIX       *pixt, *pixd;

    PROCNAME("pixConvertRGBToHue");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);

    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32 && !pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("not cmapped or rgb", procName, NULL);
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);

        /* Convert RGB image */
    pixd = pixCreate(w, h, 8);
    pixCopyResolution(pixd, pixs);
    wplt = pixGetWpl(pixt);
    datat = pixGetData(pixt);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            minrg = L_MIN(rval, gval);
            min = L_MIN(minrg, bval);
            maxrg = L_MAX(rval, gval);
            max = L_MAX(maxrg, bval);
            delta = max - min;
            if (delta == 0) {  /* gray; no chroma */
                hval = 0;
            } else {
                if (rval == max)  /* between magenta and yellow */
                    fh = (l_float32)(gval - bval) / (l_float32)delta;
                else if (gval == max)  /* between yellow and cyan */
                    fh = 2. + (l_float32)(bval - rval) / (l_float32)delta;
                else  /* between cyan and magenta */
                    fh = 4. + (l_float32)(rval - gval) / (l_float32)delta;
                fh *= 40.0;
                if (fh < 0.0)
                    fh += 240.0;
                hval = (l_int32)(fh + 0.5);
            }
            SET_DATA_BYTE(lined, j, hval);
        }
    }
    pixDestroy(&pixt);

    return pixd;
}



PIX *
pixConvertRGBToSaturation(PIX  *pixs)
{
l_int32    w, h, d, wplt, wpld;
l_int32    i, j, rval, gval, bval, sval, minrg, min, maxrg, max, delta;
l_uint32   pixel;
l_uint32  *linet, *lined, *datat, *datad;
PIX       *pixt, *pixd;

    PROCNAME("pixConvertRGBToSaturation");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);

    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32 && !pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("not cmapped or rgb", procName, NULL);
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);

        /* Convert RGB image */
    pixd = pixCreate(w, h, 8);
    pixCopyResolution(pixd, pixs);
    wplt = pixGetWpl(pixt);
    datat = pixGetData(pixt);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            minrg = L_MIN(rval, gval);
            min = L_MIN(minrg, bval);
            maxrg = L_MAX(rval, gval);
            max = L_MAX(maxrg, bval);
            delta = max - min;
            if (delta == 0)  /* gray; no chroma */
                sval = 0;
            else
                sval = (l_int32)(255. *
                                 (l_float32)delta / (l_float32)max + 0.5);
            SET_DATA_BYTE(lined, j, sval);
        }
    }

    pixDestroy(&pixt);
    return pixd;
}


PIX *
pixConvertRGBToValue(PIX  *pixs)
{
l_int32    w, h, d, wplt, wpld;
l_int32    i, j, rval, gval, bval, maxrg, max;
l_uint32   pixel;
l_uint32  *linet, *lined, *datat, *datad;
PIX       *pixt, *pixd;

    PROCNAME("pixConvertRGBToValue");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);

    pixGetDimensions(pixs, &w, &h, &d);
    if (d != 32 && !pixGetColormap(pixs))
        return (PIX *)ERROR_PTR("not cmapped or rgb", procName, NULL);
    pixt = pixRemoveColormap(pixs, REMOVE_CMAP_TO_FULL_COLOR);

        /* Convert RGB image */
    pixd = pixCreate(w, h, 8);
    pixCopyResolution(pixd, pixs);
    wplt = pixGetWpl(pixt);
    datat = pixGetData(pixt);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            extractRGBValues(pixel, &rval, &gval, &bval);
            maxrg = L_MAX(rval, gval);
            max = L_MAX(maxrg, bval);
            SET_DATA_BYTE(lined, j, max);
        }
    }

    pixDestroy(&pixt);
    return pixd;
}


/*---------------------------------------------------------------------------*
 *            Selection and display of range of colors in HSV space          *
 *---------------------------------------------------------------------------*/
PIX *
pixMakeRangeMaskHS(PIX     *pixs,
                   l_int32  huecenter,
                   l_int32  huehw,
                   l_int32  satcenter,
                   l_int32  sathw,
                   l_int32  regionflag)
{
l_int32    i, j, w, h, wplt, wpld, hstart, hend, sstart, send, hval, sval;
l_int32   *hlut, *slut;
l_uint32   pixel;
l_uint32  *datat, *datad, *linet, *lined;
PIX       *pixt, *pixd;

    PROCNAME("pixMakeRangeMaskHS");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
    if (regionflag != L_INCLUDE_REGION && regionflag != L_EXCLUDE_REGION)
        return (PIX *)ERROR_PTR("invalid regionflag", procName, NULL);

        /* Set up LUTs for hue and saturation.  These have the value 1
         * within the specified intervals of hue and saturation. */
    hlut = (l_int32 *)LEPT_CALLOC(240, sizeof(l_int32));
    slut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    sstart = L_MAX(0, satcenter - sathw);
    send = L_MIN(255, satcenter + sathw);
    for (i = sstart; i <= send; i++)
        slut[i] = 1;
    hstart = (huecenter - huehw + 240) % 240;
    hend = (huecenter + huehw + 240) % 240;
    if (hstart < hend) {
        for (i = hstart; i <= hend; i++)
            hlut[i] = 1;
    } else {  /* wrap */
        for (i = hstart; i < 240; i++)
            hlut[i] = 1;
        for (i = 0; i <= hend; i++)
            hlut[i] = 1;
    }

        /* Generate the mask */
    pixt = pixConvertRGBToHSV(NULL, pixs);
    pixGetDimensions(pixs, &w, &h, NULL);
    pixd = pixCreateNoInit(w, h, 1);
    if (regionflag == L_INCLUDE_REGION)
        pixClearAll(pixd);
    else  /* L_EXCLUDE_REGION */
        pixSetAll(pixd);
    datat = pixGetData(pixt);
    datad = pixGetData(pixd);
    wplt = pixGetWpl(pixt);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            hval = (pixel >> L_RED_SHIFT) & 0xff;
            sval = (pixel >> L_GREEN_SHIFT) & 0xff;
            if (hlut[hval] == 1 && slut[sval] == 1) {
                if (regionflag == L_INCLUDE_REGION)
                    SET_DATA_BIT(lined, j);
                else  /* L_EXCLUDE_REGION */
                    CLEAR_DATA_BIT(lined, j);
            }
        }
    }

    LEPT_FREE(hlut);
    LEPT_FREE(slut);
    pixDestroy(&pixt);
    return pixd;
}


PIX *
pixMakeRangeMaskHV(PIX     *pixs,
                   l_int32  huecenter,
                   l_int32  huehw,
                   l_int32  valcenter,
                   l_int32  valhw,
                   l_int32  regionflag)
{
l_int32    i, j, w, h, wplt, wpld, hstart, hend, vstart, vend, hval, vval;
l_int32   *hlut, *vlut;
l_uint32   pixel;
l_uint32  *datat, *datad, *linet, *lined;
PIX       *pixt, *pixd;

    PROCNAME("pixMakeRangeMaskHV");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
    if (regionflag != L_INCLUDE_REGION && regionflag != L_EXCLUDE_REGION)
        return (PIX *)ERROR_PTR("invalid regionflag", procName, NULL);

        /* Set up LUTs for hue and maximum intensity (val).  These have
         * the value 1 within the specified intervals of hue and value. */
    hlut = (l_int32 *)LEPT_CALLOC(240, sizeof(l_int32));
    vlut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    vstart = L_MAX(0, valcenter - valhw);
    vend = L_MIN(255, valcenter + valhw);
    for (i = vstart; i <= vend; i++)
        vlut[i] = 1;
    hstart = (huecenter - huehw + 240) % 240;
    hend = (huecenter + huehw + 240) % 240;
    if (hstart < hend) {
        for (i = hstart; i <= hend; i++)
            hlut[i] = 1;
    } else {
        for (i = hstart; i < 240; i++)
            hlut[i] = 1;
        for (i = 0; i <= hend; i++)
            hlut[i] = 1;
    }

        /* Generate the mask */
    pixt = pixConvertRGBToHSV(NULL, pixs);
    pixGetDimensions(pixs, &w, &h, NULL);
    pixd = pixCreateNoInit(w, h, 1);
    if (regionflag == L_INCLUDE_REGION)
        pixClearAll(pixd);
    else  /* L_EXCLUDE_REGION */
        pixSetAll(pixd);
    datat = pixGetData(pixt);
    datad = pixGetData(pixd);
    wplt = pixGetWpl(pixt);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            hval = (pixel >> L_RED_SHIFT) & 0xff;
            vval = (pixel >> L_BLUE_SHIFT) & 0xff;
            if (hlut[hval] == 1 && vlut[vval] == 1) {
                if (regionflag == L_INCLUDE_REGION)
                    SET_DATA_BIT(lined, j);
                else  /* L_EXCLUDE_REGION */
                    CLEAR_DATA_BIT(lined, j);
            }
        }
    }

    LEPT_FREE(hlut);
    LEPT_FREE(vlut);
    pixDestroy(&pixt);
    return pixd;
}


PIX *
pixMakeRangeMaskSV(PIX     *pixs,
                   l_int32  satcenter,
                   l_int32  sathw,
                   l_int32  valcenter,
                   l_int32  valhw,
                   l_int32  regionflag)
{
l_int32    i, j, w, h, wplt, wpld, sval, vval, sstart, send, vstart, vend;
l_int32   *slut, *vlut;
l_uint32   pixel;
l_uint32  *datat, *datad, *linet, *lined;
PIX       *pixt, *pixd;

    PROCNAME("pixMakeRangeMaskSV");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
    if (regionflag != L_INCLUDE_REGION && regionflag != L_EXCLUDE_REGION)
        return (PIX *)ERROR_PTR("invalid regionflag", procName, NULL);

        /* Set up LUTs for saturation and max intensity (val).
         * These have the value 1 within the specified intervals of
         * saturation and max intensity. */
    slut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    vlut = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
    sstart = L_MAX(0, satcenter - sathw);
    send = L_MIN(255, satcenter + sathw);
    vstart = L_MAX(0, valcenter - valhw);
    vend = L_MIN(255, valcenter + valhw);
    for (i = sstart; i <= send; i++)
        slut[i] = 1;
    for (i = vstart; i <= vend; i++)
        vlut[i] = 1;

        /* Generate the mask */
    pixt = pixConvertRGBToHSV(NULL, pixs);
    pixGetDimensions(pixs, &w, &h, NULL);
    pixd = pixCreateNoInit(w, h, 1);
    if (regionflag == L_INCLUDE_REGION)
        pixClearAll(pixd);
    else  /* L_EXCLUDE_REGION */
        pixSetAll(pixd);
    datat = pixGetData(pixt);
    datad = pixGetData(pixd);
    wplt = pixGetWpl(pixt);
    wpld = pixGetWpl(pixd);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            sval = (pixel >> L_GREEN_SHIFT) & 0xff;
            vval = (pixel >> L_BLUE_SHIFT) & 0xff;
            if (slut[sval] == 1 && vlut[vval] == 1) {
                if (regionflag == L_INCLUDE_REGION)
                    SET_DATA_BIT(lined, j);
                else  /* L_EXCLUDE_REGION */
                    CLEAR_DATA_BIT(lined, j);
            }
        }
    }

    LEPT_FREE(slut);
    LEPT_FREE(vlut);
    pixDestroy(&pixt);
    return pixd;
}


PIX *
pixMakeHistoHS(PIX     *pixs,
               l_int32  factor,
               NUMA   **pnahue,
               NUMA   **pnasat)
{
l_int32    i, j, w, h, wplt, hval, sval, nd;
l_uint32   pixel;
l_uint32  *datat, *linet;
void     **lined32;
NUMA      *nahue, *nasat;
PIX       *pixt, *pixd;

    PROCNAME("pixMakeHistoHS");

    if (pnahue) *pnahue = NULL;
    if (pnasat) *pnasat = NULL;
    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);

    if (pnahue) {
        nahue = numaCreate(240);
        numaSetCount(nahue, 240);
        *pnahue = nahue;
    }
    if (pnasat) {
        nasat = numaCreate(256);
        numaSetCount(nasat, 256);
        *pnasat = nasat;
    }

    if (factor <= 1)
        pixt = pixClone(pixs);
    else
        pixt = pixScaleBySampling(pixs, 1.0 / (l_float32)factor,
                                  1.0 / (l_float32)factor);

        /* Create the hue-saturation histogram */
    pixd = pixCreate(256, 240, 32);
    lined32 = pixGetLinePtrs(pixd, NULL);
    pixGetDimensions(pixt, &w, &h, NULL);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            hval = (pixel >> L_RED_SHIFT) & 0xff;

#if  DEBUG_HISTO
            if (hval > 239) {
                fprintf(stderr, "hval = %d for (%d,%d)\n", hval, i, j);
                continue;
            }
#endif  /* DEBUG_HISTO */

            sval = (pixel >> L_GREEN_SHIFT) & 0xff;
            if (pnahue)
                numaShiftValue(nahue, hval, 1.0);
            if (pnasat)
                numaShiftValue(nasat, sval, 1.0);
            nd = GET_DATA_FOUR_BYTES(lined32[hval], sval);
            SET_DATA_FOUR_BYTES(lined32[hval], sval, nd + 1);
        }
    }

    LEPT_FREE(lined32);
    pixDestroy(&pixt);
    return pixd;
}


PIX *
pixMakeHistoHV(PIX     *pixs,
               l_int32  factor,
               NUMA   **pnahue,
               NUMA   **pnaval)
{
l_int32    i, j, w, h, wplt, hval, vval, nd;
l_uint32   pixel;
l_uint32  *datat, *linet;
void     **lined32;
NUMA      *nahue, *naval;
PIX       *pixt, *pixd;

    PROCNAME("pixMakeHistoHV");

    if (pnahue) *pnahue = NULL;
    if (pnaval) *pnaval = NULL;
    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);

    if (pnahue) {
        nahue = numaCreate(240);
        numaSetCount(nahue, 240);
        *pnahue = nahue;
    }
    if (pnaval) {
        naval = numaCreate(256);
        numaSetCount(naval, 256);
        *pnaval = naval;
    }

    if (factor <= 1)
        pixt = pixClone(pixs);
    else
        pixt = pixScaleBySampling(pixs, 1.0 / (l_float32)factor,
                                  1.0 / (l_float32)factor);

        /* Create the hue-value histogram */
    pixd = pixCreate(256, 240, 32);
    lined32 = pixGetLinePtrs(pixd, NULL);
    pixGetDimensions(pixt, &w, &h, NULL);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            hval = (pixel >> L_RED_SHIFT) & 0xff;
            vval = (pixel >> L_BLUE_SHIFT) & 0xff;
            if (pnahue)
                numaShiftValue(nahue, hval, 1.0);
            if (pnaval)
                numaShiftValue(naval, vval, 1.0);
            nd = GET_DATA_FOUR_BYTES(lined32[hval], vval);
            SET_DATA_FOUR_BYTES(lined32[hval], vval, nd + 1);
        }
    }

    LEPT_FREE(lined32);
    pixDestroy(&pixt);
    return pixd;
}


PIX *
pixMakeHistoSV(PIX     *pixs,
               l_int32  factor,
               NUMA   **pnasat,
               NUMA   **pnaval)
{
l_int32    i, j, w, h, wplt, sval, vval, nd;
l_uint32   pixel;
l_uint32  *datat, *linet;
void     **lined32;
NUMA      *nasat, *naval;
PIX       *pixt, *pixd;

    PROCNAME("pixMakeHistoSV");

    if (pnasat) *pnasat = NULL;
    if (pnaval) *pnaval = NULL;
    if (!pixs || pixGetDepth(pixs) != 32)
        return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);

    if (pnasat) {
        nasat = numaCreate(256);
        numaSetCount(nasat, 256);
        *pnasat = nasat;
    }
    if (pnaval) {
        naval = numaCreate(256);
        numaSetCount(naval, 256);
        *pnaval = naval;
    }

    if (factor <= 1)
        pixt = pixClone(pixs);
    else
        pixt = pixScaleBySampling(pixs, 1.0 / (l_float32)factor,
                                  1.0 / (l_float32)factor);

        /* Create the hue-value histogram */
    pixd = pixCreate(256, 256, 32);
    lined32 = pixGetLinePtrs(pixd, NULL);
    pixGetDimensions(pixt, &w, &h, NULL);
    datat = pixGetData(pixt);
    wplt = pixGetWpl(pixt);
    for (i = 0; i < h; i++) {
        linet = datat + i * wplt;
        for (j = 0; j < w; j++) {
            pixel = linet[j];
            sval = (pixel >> L_GREEN_SHIFT) & 0xff;
            vval = (pixel >> L_BLUE_SHIFT) & 0xff;
            if (pnasat)
                numaShiftValue(nasat, sval, 1.0);
            if (pnaval)
                numaShiftValue(naval, vval, 1.0);
            nd = GET_DATA_FOUR_BYTES(lined32[sval], vval);
            SET_DATA_FOUR_BYTES(lined32[sval], vval, nd + 1);
        }
    }

    LEPT_FREE(lined32);
    pixDestroy(&pixt);
    return pixd;
}


l_ok
pixFindHistoPeaksHSV(PIX       *pixs,
                     l_int32    type,
                     l_int32    width,
                     l_int32    height,
                     l_int32    npeaks,
                     l_float32  erasefactor,
                     PTA      **ppta,
                     NUMA     **pnatot,
                     PIXA     **ppixa)
{
l_int32   i, xmax, ymax, ewidth, eheight;
l_uint32  maxval;
BOX      *box;
NUMA     *natot;
PIX      *pixh, *pixw, *pix1, *pix2, *pix3;
PTA      *pta;

    PROCNAME("pixFindHistoPeaksHSV");

    if (ppixa) *ppixa = NULL;
    if (ppta) *ppta = NULL;
    if (pnatot) *pnatot = NULL;
    if (!pixs || pixGetDepth(pixs) != 32)
        return ERROR_INT("pixs undefined or not 32 bpp", procName, 1);
    if (!ppta || !pnatot)
        return ERROR_INT("&pta and &natot not both defined", procName, 1);
    if (type != L_HS_HISTO && type != L_HV_HISTO && type != L_SV_HISTO)
        return ERROR_INT("invalid HSV histo type", procName, 1);

    if ((pta = ptaCreate(npeaks)) == NULL)
        return ERROR_INT("pta not made", procName, 1);
    *ppta = pta;
    if ((natot = numaCreate(npeaks)) == NULL)
        return ERROR_INT("natot not made", procName, 1);
    *pnatot = natot;

    *ppta = pta;
    if (type == L_SV_HISTO)
        pixh = pixAddMirroredBorder(pixs, width + 1, width + 1, height + 1,
                                    height + 1);
    else  /* type == L_HS_HISTO or type == L_HV_HISTO */
        pixh = pixAddMixedBorder(pixs, width + 1, width + 1, height + 1,
                                 height + 1);

        /* Get the total count in the sliding window.  If the window
         * fully covers the peak, this will be the integrated
         * volume under the peak. */
    pixw = pixWindowedMean(pixh, width, height, 1, 0);
    pixDestroy(&pixh);

        /* Sequentially identify and erase peaks in the histogram.
         * If requested for debugging, save a pixa of the sequence of
         * false color histograms. */
    if (ppixa)
        *ppixa = pixaCreate(0);
    for (i = 0; i < npeaks; i++) {
        pixGetMaxValueInRect(pixw, NULL, &maxval, &xmax, &ymax);
        if (maxval == 0) break;
        numaAddNumber(natot, maxval);
        ptaAddPt(pta, xmax, ymax);
        ewidth = (l_int32)(width * erasefactor);
        eheight = (l_int32)(height * erasefactor);
        box = boxCreate(xmax - ewidth, ymax - eheight, 2 * ewidth + 1,
                        2 * eheight + 1);

        if (ppixa) {
            pix1 = pixMaxDynamicRange(pixw, L_LINEAR_SCALE);
            pixaAddPix(*ppixa, pix1, L_INSERT);
            pix2 = pixConvertGrayToFalseColor(pix1, 1.0);
            pixaAddPix(*ppixa, pix2, L_INSERT);
            pix1 = pixMaxDynamicRange(pixw, L_LOG_SCALE);
            pix2 = pixConvertGrayToFalseColor(pix1, 1.0);
            pixaAddPix(*ppixa, pix2, L_INSERT);
            pix3 = pixConvertTo32(pix1);
            pixRenderHashBoxArb(pix3, box, 6, 2, L_NEG_SLOPE_LINE,
                                1, 255, 100, 100);
            pixaAddPix(*ppixa, pix3, L_INSERT);
            pixDestroy(&pix1);
        }

        pixClearInRect(pixw, box);
        boxDestroy(&box);
        if (type == L_HS_HISTO || type == L_HV_HISTO) {
                /* clear wraps at bottom and top */
            if (ymax - eheight < 0) {  /* overlap to bottom */
                box = boxCreate(xmax - ewidth, 240 + ymax - eheight,
                                2 * ewidth + 1, eheight - ymax);
            } else if (ymax + eheight > 239) {  /* overlap to top */
                box = boxCreate(xmax - ewidth, 0, 2 * ewidth + 1,
                                ymax + eheight - 239);
            } else {
                box = NULL;
            }
            if (box) {
                pixClearInRect(pixw, box);
                boxDestroy(&box);
            }
        }
    }

    pixDestroy(&pixw);
    return 0;
}


PIX *
displayHSVColorRange(l_int32  hval,
                     l_int32  sval,
                     l_int32  vval,
                     l_int32  huehw,
                     l_int32  sathw,
                     l_int32  nsamp,
                     l_int32  factor)
{
l_int32  i, j, w, huedelta, satdelta, hue, sat, rval, gval, bval;
PIX     *pixt, *pixd;

    PROCNAME("displayHSVColorRange");

    if (hval < 0 || hval > 240)
        return (PIX *)ERROR_PTR("invalid hval", procName, NULL);
    if (huehw < 5 || huehw > 120)
        return (PIX *)ERROR_PTR("invalid huehw", procName, NULL);
    if (sval - sathw < 0 || sval + sathw > 255)
        return (PIX *)ERROR_PTR("invalid sval/sathw", procName, NULL);
    if (nsamp < 1 || factor < 3)
        return (PIX *)ERROR_PTR("invalid nsamp or rep. factor", procName, NULL);
    if (vval < 0 || vval > 255)
        return (PIX *)ERROR_PTR("invalid vval", procName, NULL);

    w = (2 * nsamp + 1);
    huedelta = (l_int32)((l_float32)huehw / (l_float32)nsamp);
    satdelta = (l_int32)((l_float32)sathw / (l_float32)nsamp);
    pixt = pixCreate(w, w, 32);
    for (i = 0; i < w; i++) {
        hue = hval + huedelta * (i - nsamp);
        if (hue < 0) hue += 240;
        if (hue >= 240) hue -= 240;
        for (j = 0; j < w; j++) {
            sat = sval + satdelta * (j - nsamp);
            convertHSVToRGB(hue, sat, vval, &rval, &gval, &bval);
            pixSetRGBPixel(pixt, j, i, rval, gval, bval);
        }
    }

    pixd = pixExpandReplicate(pixt, factor);
    pixDestroy(&pixt);
    return pixd;
}


/*---------------------------------------------------------------------------*
 *                Colorspace conversion between RGB and YUV                  *
 *---------------------------------------------------------------------------*/
PIX *
pixConvertRGBToYUV(PIX  *pixd,
                   PIX  *pixs)
{
l_int32    w, h, d, wpl, i, j, rval, gval, bval, yval, uval, vval;
l_uint32  *line, *data;
PIXCMAP   *cmap;

    PROCNAME("pixConvertRGBToYUV");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixd && pixd != pixs)
        return (PIX *)ERROR_PTR("pixd defined and not inplace", procName, pixd);

    d = pixGetDepth(pixs);
    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("not cmapped or rgb", procName, pixd);

    if (!pixd)
        pixd = pixCopy(NULL, pixs);

    cmap = pixGetColormap(pixd);
    if (cmap) {   /* just convert the colormap */
        pixcmapConvertRGBToYUV(cmap);
        return pixd;
    }

        /* Convert RGB image */
    pixGetDimensions(pixd, &w, &h, NULL);
    wpl = pixGetWpl(pixd);
    data = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        for (j = 0; j < w; j++) {
            extractRGBValues(line[j], &rval, &gval, &bval);
            convertRGBToYUV(rval, gval, bval, &yval, &uval, &vval);
            line[j] = (yval << 24) | (uval << 16) | (vval << 8);
        }
    }

    return pixd;
}


PIX *
pixConvertYUVToRGB(PIX  *pixd,
                   PIX  *pixs)
{
l_int32    w, h, d, wpl, i, j, rval, gval, bval, yval, uval, vval;
l_uint32   pixel;
l_uint32  *line, *data;
PIXCMAP   *cmap;

    PROCNAME("pixConvertYUVToRGB");

    if (!pixs)
        return (PIX *)ERROR_PTR("pixs not defined", procName, pixd);
    if (pixd && pixd != pixs)
        return (PIX *)ERROR_PTR("pixd defined and not inplace", procName, pixd);

    d = pixGetDepth(pixs);
    cmap = pixGetColormap(pixs);
    if (!cmap && d != 32)
        return (PIX *)ERROR_PTR("not cmapped or hsv", procName, pixd);

    if (!pixd)
        pixd = pixCopy(NULL, pixs);

    cmap = pixGetColormap(pixd);
    if (cmap) {   /* just convert the colormap */
        pixcmapConvertYUVToRGB(cmap);
        return pixd;
    }

        /* Convert YUV image */
    pixGetDimensions(pixd, &w, &h, NULL);
    wpl = pixGetWpl(pixd);
    data = pixGetData(pixd);
    for (i = 0; i < h; i++) {
        line = data + i * wpl;
        for (j = 0; j < w; j++) {
            pixel = line[j];
            yval = pixel >> 24;
            uval = (pixel >> 16) & 0xff;
            vval = (pixel >> 8) & 0xff;
            convertYUVToRGB(yval, uval, vval, &rval, &gval, &bval);
            composeRGBPixel(rval, gval, bval, line + j);
        }
    }

    return pixd;
}


l_ok
convertRGBToYUV(l_int32   rval,
                l_int32   gval,
                l_int32   bval,
                l_int32  *pyval,
                l_int32  *puval,
                l_int32  *pvval)
{
l_float32  norm;

    PROCNAME("convertRGBToYUV");

    if (pyval) *pyval = 0;
    if (puval) *puval = 0;
    if (pvval) *pvval = 0;
    if (!pyval || !puval || !pvval)
        return ERROR_INT("&yval, &uval, &vval not all defined", procName, 1);

    norm = 1.0 / 256.;
    *pyval = (l_int32)(16.0 +
                norm * (65.738 * rval + 129.057 * gval + 25.064 * bval) + 0.5);
    *puval = (l_int32)(128.0 +
                norm * (-37.945 * rval -74.494 * gval + 112.439 * bval) + 0.5);
    *pvval = (l_int32)(128.0 +
                norm * (112.439 * rval - 94.154 * gval - 18.285 * bval) + 0.5);
    return 0;
}


l_ok
convertYUVToRGB(l_int32   yval,
                l_int32   uval,
                l_int32   vval,
                l_int32  *prval,
                l_int32  *pgval,
                l_int32  *pbval)
{
l_int32    rval, gval, bval;
l_float32  norm, ym, um, vm;

    PROCNAME("convertYUVToRGB");

    if (prval) *prval = 0;
    if (pgval) *pgval = 0;
    if (pbval) *pbval = 0;
    if (!prval || !pgval || !pbval)
        return ERROR_INT("&rval, &gval, &bval not all defined", procName, 1);

    norm = 1.0 / 256.;
    ym = yval - 16.0;
    um = uval - 128.0;
    vm = vval - 128.0;
    rval =  (l_int32)(norm * (298.082 * ym + 408.583 * vm) + 0.5);
    gval = (l_int32)(norm * (298.082 * ym - 100.291 * um - 208.120 * vm) +
                       0.5);
    bval = (l_int32)(norm * (298.082 * ym + 516.411 * um) + 0.5);
    *prval = L_MIN(255, L_MAX(0, rval));
    *pgval = L_MIN(255, L_MAX(0, gval));
    *pbval = L_MIN(255, L_MAX(0, bval));

    return 0;
}


l_ok
pixcmapConvertRGBToYUV(PIXCMAP  *cmap)
{
l_int32   i, ncolors, rval, gval, bval, yval, uval, vval;

    PROCNAME("pixcmapConvertRGBToYUV");

    if (!cmap)
        return ERROR_INT("cmap not defined", procName, 1);

    ncolors = pixcmapGetCount(cmap);
    for (i = 0; i < ncolors; i++) {
        pixcmapGetColor(cmap, i, &rval, &gval, &bval);
        convertRGBToYUV(rval, gval, bval, &yval, &uval, &vval);
        pixcmapResetColor(cmap, i, yval, uval, vval);
    }
    return 0;
}


l_ok
pixcmapConvertYUVToRGB(PIXCMAP  *cmap)
{
l_int32   i, ncolors, rval, gval, bval, yval, uval, vval;

    PROCNAME("pixcmapConvertYUVToRGB");

    if (!cmap)
        return ERROR_INT("cmap not defined", procName, 1);

    ncolors = pixcmapGetCount(cmap);
    for (i = 0; i < ncolors; i++) {
        pixcmapGetColor(cmap, i, &yval, &uval, &vval);
        convertYUVToRGB(yval, uval, vval, &rval, &gval, &bval);
        pixcmapResetColor(cmap, i, rval, gval, bval);
    }
    return 0;
}


/*---------------------------------------------------------------------------*
 *                Colorspace conversion between RGB and XYZ                  *
 *---------------------------------------------------------------------------*/
FPIXA *
pixConvertRGBToXYZ(PIX  *pixs)
{
l_int32     w, h, wpls, wpld, i, j, rval, gval, bval;
l_uint32   *lines, *datas;
l_float32   fxval, fyval, fzval;
l_float32  *linex, *liney, *linez, *datax, *datay, *dataz;
FPIX       *fpix;
FPIXA      *fpixa;

    PROCNAME("pixConvertRGBToXYZ");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (FPIXA *)ERROR_PTR("pixs undefined or not rgb", procName, NULL);

        /* Convert RGB image */
    pixGetDimensions(pixs, &w, &h, NULL);
    fpixa = fpixaCreate(3);
    for (i = 0; i < 3; i++) {
        fpix = fpixCreate(w, h);
        fpixaAddFPix(fpixa, fpix, L_INSERT);
    }
    wpls = pixGetWpl(pixs);
    wpld = fpixGetWpl(fpix);
    datas = pixGetData(pixs);
    datax = fpixaGetData(fpixa, 0);
    datay = fpixaGetData(fpixa, 1);
    dataz = fpixaGetData(fpixa, 2);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        linex = datax + i * wpld;
        liney = datay + i * wpld;
        linez = dataz + i * wpld;
        for (j = 0; j < w; j++) {
            extractRGBValues(lines[j], &rval, &gval, &bval);
            convertRGBToXYZ(rval, gval, bval, &fxval, &fyval, &fzval);
            *(linex + j) = fxval;
            *(liney + j) = fyval;
            *(linez + j) = fzval;
        }
    }

    return fpixa;
}


PIX *
fpixaConvertXYZToRGB(FPIXA  *fpixa)
{
l_int32     w, h, wpls, wpld, i, j, rval, gval, bval;
l_float32   fxval, fyval, fzval;
l_float32  *linex, *liney, *linez, *datax, *datay, *dataz;
l_uint32   *lined, *datad;
PIX        *pixd;
FPIX       *fpix;

    PROCNAME("fpixaConvertXYZToRGB");

    if (!fpixa || fpixaGetCount(fpixa) != 3)
        return (PIX *)ERROR_PTR("fpixa undefined or invalid", procName, NULL);

        /* Convert XYZ image */
    if (fpixaGetFPixDimensions(fpixa, 0, &w, &h))
        return (PIX *)ERROR_PTR("fpixa dimensions not found", procName, NULL);
    pixd = pixCreate(w, h, 32);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);
    datax = fpixaGetData(fpixa, 0);
    datay = fpixaGetData(fpixa, 1);
    dataz = fpixaGetData(fpixa, 2);
    fpix = fpixaGetFPix(fpixa, 0, L_CLONE);
    wpls = fpixGetWpl(fpix);
    fpixDestroy(&fpix);
    for (i = 0; i < h; i++) {
        linex = datax + i * wpls;
        liney = datay + i * wpls;
        linez = dataz + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            fxval = linex[j];
            fyval = liney[j];
            fzval = linez[j];
            convertXYZToRGB(fxval, fyval, fzval, 0, &rval, &gval, &bval);
            composeRGBPixel(rval, gval, bval, lined + j);
        }
    }

    return pixd;
}


l_ok
convertRGBToXYZ(l_int32     rval,
                l_int32     gval,
                l_int32     bval,
                l_float32  *pfxval,
                l_float32  *pfyval,
                l_float32  *pfzval)
{
    PROCNAME("convertRGBToXYZ");

    if (pfxval) *pfxval = 0.0;
    if (pfyval) *pfyval = 0.0;
    if (pfzval) *pfzval = 0.0;
    if (!pfxval || !pfyval || !pfzval)
        return ERROR_INT("&xval, &yval, &zval not all defined", procName, 1);

    *pfxval = 0.4125 * rval + 0.3576 * gval + 0.1804 * bval;
    *pfyval = 0.2127 * rval + 0.7152 * gval + 0.0722 * bval;
    *pfzval = 0.0193 * rval + 0.1192 * gval + 0.9502 * bval;
    return 0;
}


l_ok
convertXYZToRGB(l_float32  fxval,
                l_float32  fyval,
                l_float32  fzval,
                l_int32    blackout,
                l_int32   *prval,
                l_int32   *pgval,
                l_int32   *pbval)
{
l_int32  rval, gval, bval;

    PROCNAME("convertXYZToRGB");

    if (prval) *prval = 0;
    if (pgval) *pgval = 0;
    if (pbval) *pbval = 0;
    if (!prval || !pgval ||!pbval)
        return ERROR_INT("&rval, &gval, &bval not all defined", procName, 1);
    *prval = *pgval = *pbval = 0;

    rval = (l_int32)(3.2405 * fxval - 1.5372 * fyval - 0.4985 * fzval + 0.5);
    gval = (l_int32)(-0.9693 * fxval + 1.8760 * fyval + 0.0416 * fzval + 0.5);
    bval = (l_int32)(0.0556 * fxval - 0.2040 * fyval + 1.0573 * fzval + 0.5);
    if (blackout == 0) {  /* the usual situation; use nearest rgb color */
        *prval = L_MAX(0, L_MIN(rval, 255));
        *pgval = L_MAX(0, L_MIN(gval, 255));
        *pbval = L_MAX(0, L_MIN(bval, 255));
    } else {  /* use black for out of gamut */
        if (rval >= 0 && rval < 256 && gval >= 0 && gval < 256 &&
            bval >= 0 && bval < 256) {  /* in gamut */
            *prval = rval;
            *pgval = gval;
            *pbval = bval;
        }
    }
    return 0;
}


/*---------------------------------------------------------------------------*
 *               Colorspace conversion between XYZ and LAB                   *
 *---------------------------------------------------------------------------*/
FPIXA *
fpixaConvertXYZToLAB(FPIXA  *fpixas)
{
l_int32     w, h, wpl, i, j;
l_float32   fxval, fyval, fzval, flval, faval, fbval;
l_float32  *linex, *liney, *linez, *datax, *datay, *dataz;
l_float32  *linel, *linea, *lineb, *datal, *dataa, *datab;
FPIX       *fpix;
FPIXA      *fpixad;

    PROCNAME("fpixaConvertXYZToLAB");

    if (!fpixas || fpixaGetCount(fpixas) != 3)
        return (FPIXA *)ERROR_PTR("fpixas undefined/invalid", procName, NULL);

        /* Convert XYZ image */
    if (fpixaGetFPixDimensions(fpixas, 0, &w, &h))
        return (FPIXA *)ERROR_PTR("fpixas sizes not found", procName, NULL);
    fpixad = fpixaCreate(3);
    for (i = 0; i < 3; i++) {
        fpix = fpixCreate(w, h);
        fpixaAddFPix(fpixad, fpix, L_INSERT);
    }
    wpl = fpixGetWpl(fpix);
    datax = fpixaGetData(fpixas, 0);
    datay = fpixaGetData(fpixas, 1);
    dataz = fpixaGetData(fpixas, 2);
    datal = fpixaGetData(fpixad, 0);
    dataa = fpixaGetData(fpixad, 1);
    datab = fpixaGetData(fpixad, 2);

        /* Convert XYZ image */
    for (i = 0; i < h; i++) {
        linex = datax + i * wpl;
        liney = datay + i * wpl;
        linez = dataz + i * wpl;
        linel = datal + i * wpl;
        linea = dataa + i * wpl;
        lineb = datab + i * wpl;
        for (j = 0; j < w; j++) {
            fxval = *(linex + j);
            fyval = *(liney + j);
            fzval = *(linez + j);
            convertXYZToLAB(fxval, fyval, fzval, &flval, &faval, &fbval);
            *(linel + j) = flval;
            *(linea + j) = faval;
            *(lineb + j) = fbval;
        }
    }

    return fpixad;
}


FPIXA *
fpixaConvertLABToXYZ(FPIXA  *fpixas)
{
l_int32     w, h, wpl, i, j;
l_float32   fxval, fyval, fzval, flval, faval, fbval;
l_float32  *linel, *linea, *lineb, *datal, *dataa, *datab;
l_float32  *linex, *liney, *linez, *datax, *datay, *dataz;
FPIX       *fpix;
FPIXA      *fpixad;

    PROCNAME("fpixaConvertLABToXYZ");

    if (!fpixas || fpixaGetCount(fpixas) != 3)
        return (FPIXA *)ERROR_PTR("fpixas undefined/invalid", procName, NULL);

        /* Convert LAB image */
    if (fpixaGetFPixDimensions(fpixas, 0, &w, &h))
        return (FPIXA *)ERROR_PTR("fpixas sizes not found", procName, NULL);
    fpixad = fpixaCreate(3);
    for (i = 0; i < 3; i++) {
        fpix = fpixCreate(w, h);
        fpixaAddFPix(fpixad, fpix, L_INSERT);
    }
    wpl = fpixGetWpl(fpix);
    datal = fpixaGetData(fpixas, 0);
    dataa = fpixaGetData(fpixas, 1);
    datab = fpixaGetData(fpixas, 2);
    datax = fpixaGetData(fpixad, 0);
    datay = fpixaGetData(fpixad, 1);
    dataz = fpixaGetData(fpixad, 2);

        /* Convert XYZ image */
    for (i = 0; i < h; i++) {
        linel = datal + i * wpl;
        linea = dataa + i * wpl;
        lineb = datab + i * wpl;
        linex = datax + i * wpl;
        liney = datay + i * wpl;
        linez = dataz + i * wpl;
        for (j = 0; j < w; j++) {
            flval = *(linel + j);
            faval = *(linea + j);
            fbval = *(lineb + j);
            convertLABToXYZ(flval, faval, fbval, &fxval, &fyval, &fzval);
            *(linex + j) = fxval;
            *(liney + j) = fyval;
            *(linez + j) = fzval;
        }
    }

    return fpixad;
}


l_ok
convertXYZToLAB(l_float32   xval,
                l_float32   yval,
                l_float32   zval,
                l_float32  *plval,
                l_float32  *paval,
                l_float32  *pbval)
{
l_float32  xn, yn, zn, fx, fy, fz;

    PROCNAME("convertXYZToLAB");

    if (plval) *plval = 0.0;
    if (paval) *paval = 0.0;
    if (pbval) *pbval = 0.0;
    if (!plval || !paval || !pbval)
        return ERROR_INT("&lval, &aval, &bval not all defined", procName, 1);

        /* First normalize to the corresponding white values */
    xn = 0.0041259 * xval;
    yn = 0.0039216 * yval;
    zn = 0.0036012 * zval;
        /* Then apply the lab_forward function */
    fx = lab_forward(xn);
    fy = lab_forward(yn);
    fz = lab_forward(zn);
    *plval = 116.0 * fy - 16.0;
    *paval = 500.0 * (fx - fy);
    *pbval = 200.0 * (fy - fz);
    return 0;
}


l_ok
convertLABToXYZ(l_float32   lval,
                l_float32   aval,
                l_float32   bval,
                l_float32  *pxval,
                l_float32  *pyval,
                l_float32  *pzval)
{
l_float32  fx, fy, fz;
l_float32  xw = 242.37;  /* x component corresponding to rgb white */
l_float32  yw = 255.0;  /* y component corresponding to rgb white */
l_float32  zw = 277.69;  /* z component corresponding to rgb white */

    PROCNAME("convertLABToXYZ");

    if (pxval) *pxval = 0.0;
    if (pyval) *pyval = 0.0;
    if (pzval) *pzval = 0.0;
    if (!pxval || !pyval || !pzval)
        return ERROR_INT("&xval, &yval, &zval not all defined", procName, 1);

    fy = 0.0086207 * (16.0 + lval);
    fx = fy + 0.002 * aval;
    fz = fy - 0.005 * bval;
    *pxval = xw * lab_reverse(fx);
    *pyval = yw * lab_reverse(fy);
    *pzval = zw * lab_reverse(fz);
    return 0;
}


/*
 * See http://en.wikipedia.org/wiki/Lab_color_space for formulas.
 * This is the forward function: from xyz to lab.  It includes a rational
 * function approximation over [0.008856 ... 1] to the cube root, from
 * "Fast Color Space Transformations Using Minimax Approximations",
 * M. Celebi et al, http://arxiv.org/pdf/1009.0854v1.pdf.
 */
static l_float32
lab_forward(l_float32  v)
{
const l_float32  f_thresh = 0.008856;  /* (6/29)^3  */
const l_float32  f_factor = 7.787;  /* (1/3) * (29/6)^2)  */
const l_float32  f_offset = 0.13793;  /* 4/29 */

    if (v > f_thresh) {
#if  SLOW_CUBE_ROOT
        return powf(v, 0.333333);
#else
        l_float32  num, den;
        num = 4.37089e-04 + v * (9.52695e-02 + v * (1.25201 + v * 1.30273));
        den = 3.91236e-03 + v * (2.95408e-01 + v * (1.71714 + v * 6.34341e-01));
        return num / den;
#endif
    } else {
        return f_factor * v + f_offset;
    }
}


/*
 * See http://en.wikipedia.org/wiki/Lab_color_space for formulas.
 * This is the reverse (inverse) function: from lab to xyz.
 */
static l_float32
lab_reverse(l_float32  v)
{
const l_float32  r_thresh = 0.20690;  /* 6/29  */
const l_float32  r_factor = 0.12842;  /* 3 * (6/29)^2   */
const l_float32  r_offset = 0.13793;  /* 4/29 */

    if (v > r_thresh) {
        return v * v * v;
    } else {
        return r_factor * (v - r_offset);
    }
}


/*---------------------------------------------------------------------------*
 *               Colorspace conversion between RGB and LAB                   *
 *---------------------------------------------------------------------------*/
FPIXA *
pixConvertRGBToLAB(PIX  *pixs)
{
l_int32     w, h, wpls, wpld, i, j, rval, gval, bval;
l_uint32   *lines, *datas;
l_float32   flval, faval, fbval;
l_float32  *linel, *linea, *lineb, *datal, *dataa, *datab;
FPIX       *fpix;
FPIXA      *fpixa;

    PROCNAME("pixConvertRGBToLAB");

    if (!pixs || pixGetDepth(pixs) != 32)
        return (FPIXA *)ERROR_PTR("pixs undefined or not rgb", procName, NULL);

        /* Convert RGB image */
    pixGetDimensions(pixs, &w, &h, NULL);
    fpixa = fpixaCreate(3);
    for (i = 0; i < 3; i++) {
        fpix = fpixCreate(w, h);
        fpixaAddFPix(fpixa, fpix, L_INSERT);
    }
    wpls = pixGetWpl(pixs);
    wpld = fpixGetWpl(fpix);
    datas = pixGetData(pixs);
    datal = fpixaGetData(fpixa, 0);
    dataa = fpixaGetData(fpixa, 1);
    datab = fpixaGetData(fpixa, 2);
    for (i = 0; i < h; i++) {
        lines = datas + i * wpls;
        linel = datal + i * wpld;
        linea = dataa + i * wpld;
        lineb = datab + i * wpld;
        for (j = 0; j < w; j++) {
            extractRGBValues(lines[j], &rval, &gval, &bval);
            convertRGBToLAB(rval, gval, bval, &flval, &faval, &fbval);
            *(linel + j) = flval;
            *(linea + j) = faval;
            *(lineb + j) = fbval;
        }
    }

    return fpixa;
}


PIX *
fpixaConvertLABToRGB(FPIXA  *fpixa)
{
l_int32     w, h, wpls, wpld, i, j, rval, gval, bval;
l_float32   flval, faval, fbval;
l_float32  *linel, *linea, *lineb, *datal, *dataa, *datab;
l_uint32   *lined, *datad;
PIX        *pixd;
FPIX       *fpix;

    PROCNAME("fpixaConvertLABToRGB");

    if (!fpixa || fpixaGetCount(fpixa) != 3)
        return (PIX *)ERROR_PTR("fpixa undefined or invalid", procName, NULL);

        /* Convert LAB image */
    if (fpixaGetFPixDimensions(fpixa, 0, &w, &h))
        return (PIX *)ERROR_PTR("fpixa dimensions not found", procName, NULL);
    pixd = pixCreate(w, h, 32);
    wpld = pixGetWpl(pixd);
    datad = pixGetData(pixd);
    datal = fpixaGetData(fpixa, 0);
    dataa = fpixaGetData(fpixa, 1);
    datab = fpixaGetData(fpixa, 2);
    fpix = fpixaGetFPix(fpixa, 0, L_CLONE);
    wpls = fpixGetWpl(fpix);
    fpixDestroy(&fpix);
    for (i = 0; i < h; i++) {
        linel = datal + i * wpls;
        linea = dataa + i * wpls;
        lineb = datab + i * wpls;
        lined = datad + i * wpld;
        for (j = 0; j < w; j++) {
            flval = linel[j];
            faval = linea[j];
            fbval = lineb[j];
            convertLABToRGB(flval, faval, fbval, &rval, &gval, &bval);
            composeRGBPixel(rval, gval, bval, lined + j);
        }
    }

    return pixd;
}


l_ok
convertRGBToLAB(l_int32     rval,
                l_int32     gval,
                l_int32     bval,
                l_float32  *pflval,
                l_float32  *pfaval,
                l_float32  *pfbval)
{
l_float32  fxval, fyval, fzval;

    PROCNAME("convertRGBToLAB");

    if (pflval) *pflval = 0.0;
    if (pfaval) *pfaval = 0.0;
    if (pfbval) *pfbval = 0.0;
    if (!pflval || !pfaval || !pfbval)
        return ERROR_INT("&flval, &faval, &fbval not all defined", procName, 1);

    convertRGBToXYZ(rval, gval, bval, &fxval, &fyval, &fzval);
    convertXYZToLAB(fxval, fyval, fzval, pflval, pfaval, pfbval);
    return 0;
}


l_ok
convertLABToRGB(l_float32  flval,
                l_float32  faval,
                l_float32  fbval,
                l_int32   *prval,
                l_int32   *pgval,
                l_int32   *pbval)
{
l_float32  fxval, fyval, fzval;

    PROCNAME("convertLABToRGB");

    if (prval) *prval = 0;
    if (pgval) *pgval = 0;
    if (pbval) *pbval = 0;
    if (!prval || !pgval || !pbval)
        return ERROR_INT("&rval, &gval, &bval not all defined", procName, 1);

    convertLABToXYZ(flval, faval, fbval, &fxval, &fyval, &fzval);
    convertXYZToRGB(fxval, fyval, fzval, 0, prval, pgval, pbval);
    return 0;
}