#include <math.h>
#include "batypes.h"
#include "batype.h"
#include "balibtyp.h"
#include "barwtyp.h"
#include "baresour.h"
#include "baimage.h"
#include "badebug.h"
#include "palquant.h"

#if (!defined(DOXYGEN))
static const char rcsid[] __RWUNUSED__ =
    "@@(#)$Id: palquant.c,v 1.28 2001/08/20 15:08:26 johns Exp $";
#endif /* (!defined(DOXYGEN)) */

#if defined (__MWERKS__)
#if (defined(RWVERBOSE))
#pragma message (__FILE__ "/" _SKY_EXPAND(__LINE__) ": __MWERKS__ == " _SKY_EXPAND(__MWERKS__))
#endif /* (defined (__MWERKS__)) */
#if (__option (global_optimizer))
#pragma always_inline off
#endif /* (__option (global_optimizer)) */
#endif /*  defined (__MWERKS__) */

#define	RED	1
#define	GREEN	2
#define BLUE	3
#define ALPHA	4

#define MAXDEPTH 5
#define MAXCOLOR 256

#define RwRGBARealLengthSq(a)	((a)->red*(a)->red + (a)->green*(a)->green + (a)->blue*(a)->blue + (a)->alpha*(a)->alpha)

typedef unsigned long OctantMap;

/* local */
static OctantMap    splice[256];

/*************************************************************************/
static LeafNode    *
InitLeaf(LeafNode * Leaf)
{

    RWFUNCTION(RWSTRING("InitLeaf"));
    RWASSERT(Leaf);

    Leaf->palIndex = -1;
    Leaf->weight = 0.0f;
    Leaf->ac.red = 0.0f;
    Leaf->ac.green = 0.0f;
    Leaf->ac.blue = 0.0f;
    Leaf->ac.alpha = 0.0f;
    Leaf->m2 = 0.0f;
    RWRETURN(Leaf);
}

/*************************************************************************/
static BranchNode  *
InitBranch(BranchNode * Branch)
{
    int                 i;

    RWFUNCTION(RWSTRING("InitBranch"));
    RWASSERT(Branch);

    for (i = 0; i < 16; i++)
    {
        Branch->dir[i] = (OctNode *)NULL;
    }

    RWRETURN(Branch);
}

/*************************************************************************/
static OctNode     *
CreateCube(RwFreeList * fl)
{
    OctNode            *cube;

    RWFUNCTION(RWSTRING("CreateCube"));
    RWASSERT(fl);

    cube = (OctNode *) RwFreeListAlloc(fl);
    RWRETURN(cube);
}

/*************************************************************************/
static              OctantMap
GetOctAdr(RwRGBA * c)
{
    int                 cs = 8 - MAXDEPTH;

    RWFUNCTION(RWSTRING("GetOctAdr"));
    RWASSERT(c);

    RWRETURN((splice[c->red >> cs] << 3) | (splice[c->green >> cs] << 2) |
             (splice[c->blue >> cs] << 1) | (splice[c->alpha >> cs] << 0));
}

/*************************************************************************/
static OctNode     *
AllocateToLeaf(rwPalQuant * pq, OctNode * root, OctantMap Octs, int depth)
{

    RWFUNCTION(RWSTRING("AllocateToLeaf"));
    RWASSERT(pq);
    RWASSERT(root);

    /* return leaf */
    if (depth == 0)
    {
        RWRETURN(root);
    }

    /* populate branch */
    if (!root->Branch.dir[Octs & 15])
    {
        OctNode            *node;

        node = CreateCube(pq->cubefreelist);
        root->Branch.dir[Octs & 15] = node;
        if (depth == 1)
        {
            InitLeaf(&node->Leaf);
        }
        else
        {
            InitBranch(&node->Branch);
        }
    }

    RWRETURN(AllocateToLeaf
             (pq, root->Branch.dir[Octs & 15], Octs >> 4, depth - 1));
}

/*************************************************************************/
void
_rwPalQuantAddImage(rwPalQuant * pq, RwImage * img, RwReal weight)
{
    RwInt32             ColShift;
    RwInt32             width, height, stride;
    RwUInt8            *pixels;
    RwRGBA             *palette;

    RWFUNCTION(RWSTRING("_rwPalQuantAddImage"));
    RWASSERT(pq);
    RWASSERT(img);

    ColShift = 8 - MAXDEPTH;

    stride = RwImageGetStride(img);
    pixels = RwImageGetPixels(img);
    palette = RwImageGetPalette(img);
    height = RwImageGetHeight(img);

    switch (RwImageGetDepth(img))
    {
        case (4):
        case (8):
            {
                while (height--)
                {
                    RwUInt8            *linePixels = pixels;

                    width = RwImageGetWidth(img);
                    while (width--)
                    {
                        RwRGBA             *color = &palette[*linePixels];
                        RwRGBAReal          rColor;
                        OctNode            *leaf;
                        OctantMap           Octs;

                        /* build down to leaf */
                        Octs = GetOctAdr(color);
                        leaf = AllocateToLeaf(pq, pq->root, Octs, MAXDEPTH);

                        RwRGBARealFromRwRGBA(&rColor, color);
                        leaf->Leaf.weight += weight;
                        RwRGBARealScale(&rColor, &rColor, weight);
                        RwRGBARealAdd(&leaf->Leaf.ac, &leaf->Leaf.ac,
                                      &rColor);
                        leaf->Leaf.m2 += weight*RwRGBARealLengthSq(&rColor);

                        linePixels++;
                    }

                    pixels += stride;
                }
                break;
            }

        case (32):
            {
                while (height--)
                {
                    RwRGBA             *color = (RwRGBA *) pixels;

                    width = RwImageGetWidth(img);
                    while (width--)
                    {
                        RwRGBAReal          rColor;
                        OctNode            *leaf;
                        OctantMap           Octs;

                        /* build down to leaf */
                        Octs = GetOctAdr(color);
                        leaf = AllocateToLeaf(pq, pq->root, Octs, MAXDEPTH);

                        RwRGBARealFromRwRGBA(&rColor, color);
                        leaf->Leaf.weight += weight;
                        RwRGBARealScale(&rColor, &rColor, weight);
                        RwRGBARealAdd(&leaf->Leaf.ac, &leaf->Leaf.ac,
                                      &rColor);
                        leaf->Leaf.m2 += weight*RwRGBARealLengthSq(&rColor);

                        color++;
                    }

                    pixels += stride;
                }
                break;
            }
    }

    RWRETURNVOID();
}

/*************************************************************************/

/*************************************************************************/
static void
assignindex(OctNode * root, RwRGBA * origin, int depth, box * region,
            RwInt32 palIndex)
{
    int                 width, dr, dg, db, da, dR, dG, dB, dA;
    int                 i;

    RWFUNCTION(RWSTRING("assignindex"));
    RWASSERT(origin);
    RWASSERT(region);

    if (!root)
        RWRETURNVOID();

    width = 1 << depth;

    dr = origin->red - region->col1.red;
    dg = origin->green - region->col1.green;
    db = origin->blue - region->col1.blue;
    da = origin->alpha - region->col1.alpha;
    if (dr >= 0 || dg >= 0 || db >= 0 || da >= 0)
    {
        RWRETURNVOID();
    }

    dR = region->col0.red - origin->red;
    dG = region->col0.green - origin->green;
    dB = region->col0.blue - origin->blue;
    dA = region->col0.alpha - origin->alpha;
    if (dR >= width || dG >= width || dB >= width || dA >= width)
    {
        RWRETURNVOID();
    }

    /* wholly inside region and a leaf? */
    if (dr <= -width && dg <= -width && db <= -width && da <= -width)
        if (dR <= 0 && dG <= 0 && dB <= 0 && dA <= 0)
            if (depth == 0)
            {
                root->Leaf.palIndex = (RwUInt8) palIndex;
                RWRETURNVOID();
            }

    /* try children */
    depth--;
    for (i = 0; i < 16; i++)
    {
        RwRGBA              suborigin;

        suborigin.red = origin->red + (((i >> 3) & 1) << depth);
        suborigin.green = origin->green + (((i >> 2) & 1) << depth);
        suborigin.blue = origin->blue + (((i >> 1) & 1) << depth);
        suborigin.alpha = origin->alpha + (((i >> 0) & 1) << depth);

        assignindex(root->Branch.dir[i], &suborigin, depth, region, palIndex);
    }

    RWRETURNVOID();
}

/*************************************************************************/

/* Assign palIndex to leaves */
static void
nAssign(RwInt32 palIndex, OctNode * root, box * cube)
{
    RwRGBA              origin;

    RWFUNCTION(RWSTRING("nAssign"));
    RWASSERT(root);
    RWASSERT(cube);

    origin.red = 0;
    origin.green = 0;
    origin.blue = 0;
    origin.alpha = 0;
    assignindex(root, &origin, MAXDEPTH, cube, palIndex);
    RWRETURNVOID();
}

/*************************************************************************/
static void
addvolume(OctNode * root, RwRGBA * origin, int depth, box * region,
          LeafNode * volume)
{
    int                 width, dr, dg, db, da, dR, dG, dB, dA;
    int                 i;

    RWFUNCTION(RWSTRING("addvolume"));
    RWASSERT(origin);
    RWASSERT(region);

    if (!root)
        RWRETURNVOID();

    width = 1 << depth;

    dr = origin->red - region->col1.red;
    dg = origin->green - region->col1.green;
    db = origin->blue - region->col1.blue;
    da = origin->alpha - region->col1.alpha;
    if (dr >= 0 || dg >= 0 || db >= 0 || da >= 0)
    {
        RWRETURNVOID();
    }

    dR = region->col0.red - origin->red;
    dG = region->col0.green - origin->green;
    dB = region->col0.blue - origin->blue;
    dA = region->col0.alpha - origin->alpha;
    if (dR >= width || dG >= width || dB >= width || dA >= width)
    {
        RWRETURNVOID();
    }

    /* wholly inside region? */
    if (dr <= -width && dg <= -width && db <= -width && da <= -width)
        if (dR <= 0 && dG <= 0 && dB <= 0 && dA <= 0)
#ifndef CACHEWEIGHTS
            if (depth == 0)    /* we need to visit each leaf */
#endif
            {
                volume->weight += root->Leaf.weight;
                RwRGBARealAdd(&volume->ac, &volume->ac, &root->Leaf.ac);
                volume->m2 += root->Leaf.m2;
                RWRETURNVOID();
            }

    /* try children */
    depth--;
    for (i = 0; i < 16; i++)
    {
        RwRGBA              suborigin;

        suborigin.red = origin->red + (((i >> 3) & 1) << depth);
        suborigin.green = origin->green + (((i >> 2) & 1) << depth);
        suborigin.blue = origin->blue + (((i >> 1) & 1) << depth);
        suborigin.alpha = origin->alpha + (((i >> 0) & 1) << depth);

        addvolume(root->Branch.dir[i], &suborigin, depth, region, volume);
    }
    RWRETURNVOID();
}

/*************************************************************************/

/* Compute sum over a box of any given statistic */
static LeafNode    *
nVol(LeafNode * Vol, OctNode * root, box * cube)
{
    RwRGBA              origin;

    RWFUNCTION(RWSTRING("nVol"));
    RWASSERT(root);
    RWASSERT(cube);

    origin.red = 0;
    origin.green = 0;
    origin.blue = 0;
    origin.alpha = 0;
    InitLeaf(Vol);
    addvolume(root, &origin, MAXDEPTH, cube, Vol);

    RWRETURN(Vol);
}

/*************************************************************************/

/* Compute the weighted variance of a box */

/* NB: as with the raw statistics, this is really the variance * size */
static              RwReal
nVar(OctNode * root, box * cube)
{
    LeafNode            Node;

    RWFUNCTION(RWSTRING("nVar"));
    RWASSERT(root);
    RWASSERT(cube);

    nVol(&Node, root, cube);

    RWRETURN(Node.m2 - (RwRGBARealLengthSq(&Node.ac) / Node.weight));
}

/*************************************************************************/

/* We want to minimize the sum of the variances of two subboxes.
 * The sum(c^2) terms can be ignored since their sum over both subboxes
 * is the same (the sum for the whole box) no matter where we split.
 * The remaining terms have a minus sign in the variance formula,
 * so we drop the minus sign and MAXIMIZE the sum of the two terms.
 */
static              RwReal
nMaximize(OctNode * root, box * cube, RwInt32 dir, RwInt32 * cut,
          LeafNode * whole)
{
    box                 infcube;
    LeafNode            left, right;
    RwReal              maxsum, val, lastval;
    RwInt32             i;

    RWFUNCTION(RWSTRING("nMaximize"));
    RWASSERT(root);
    RWASSERT(cube);
    RWASSERT(cut);
    RWASSERT(whole);

    lastval = maxsum = 0.0f;
    *cut = -1;
    infcube = *cube;
    switch (dir)
    {
        case RED:
            for (i = cube->col0.red; i < cube->col1.red; i++)
            {
                infcube.col1.red = (RwUInt8) i;

                nVol(&left, root, &infcube);
                RwRGBARealSub(&right.ac, &whole->ac, &left.ac);
                right.weight = whole->weight - left.weight;
                if ((left.weight > 0.0f) && (right.weight > 0.0f))
                {
                    val = RwRGBARealLengthSq(&left.ac) / left.weight;
                    val += RwRGBARealLengthSq(&right.ac) / right.weight;

                    if (val > maxsum)
                    {
                        maxsum = val;
                        *cut = i;
                    }
                    else if (val < lastval)
                    {
                        /* we've past the peak */
                        break;
                    }

                    lastval = val;
                }
            }
            break;

        case GREEN:
            for (i = cube->col0.green; i < cube->col1.green; i++)
            {
                infcube.col1.green = (RwUInt8) i;

                nVol(&left, root, &infcube);
                RwRGBARealSub(&right.ac, &whole->ac, &left.ac);
                right.weight = whole->weight - left.weight;
                if ((left.weight > 0.0f) && (right.weight > 0.0f))
                {
                    val = RwRGBARealLengthSq(&left.ac) / left.weight;
                    val += RwRGBARealLengthSq(&right.ac) / right.weight;

                    if (val > maxsum)
                    {
                        maxsum = val;
                        *cut = i;
                    }
                    else if (val < lastval)
                    {
                        /* we've past the peak */
                        break;
                    }

                    lastval = val;
                }
            }
            break;

        case BLUE:
            for (i = cube->col0.blue; i < cube->col1.blue; i++)
            {
                infcube.col1.blue = (RwUInt8) i;

                nVol(&left, root, &infcube);
                RwRGBARealSub(&right.ac, &whole->ac, &left.ac);
                right.weight = whole->weight - left.weight;
                if ((left.weight > 0.0f) && (right.weight > 0.0f))
                {
                    val = RwRGBARealLengthSq(&left.ac) / left.weight;
                    val += RwRGBARealLengthSq(&right.ac) / right.weight;

                    if (val > maxsum)
                    {
                        maxsum = val;
                        *cut = i;
                    }
                    else if (val < lastval)
                    {
                        /* we've past the peak */
                        break;
                    }

                    lastval = val;
                }
            }
            break;

        case ALPHA:
            for (i = cube->col0.alpha; i < cube->col1.alpha; i++)
            {
                infcube.col1.alpha = (RwUInt8) i;

                nVol(&left, root, &infcube);
                RwRGBARealSub(&right.ac, &whole->ac, &left.ac);
                right.weight = whole->weight - left.weight;
                if ((left.weight > 0.0f) && (right.weight > 0.0f))
                {
                    val = RwRGBARealLengthSq(&left.ac) / left.weight;
                    val += RwRGBARealLengthSq(&right.ac) / right.weight;

                    if (val > maxsum)
                    {
                        maxsum = val;
                        *cut = i;
                    }
                    else if (val < lastval)
                    {
                        /* we've past the peak */
                        break;
                    }

                    lastval = val;
                }
            }
            break;
    }

    RWRETURN(maxsum);
}

/*************************************************************************/
static              RwBool
nCut(OctNode * root, box * set1, box * set2)
{
    RwInt32             cutr, cutg, cutb, cuta;
    RwReal              maxr, maxg, maxb, maxa;
    LeafNode            whole;

    RWFUNCTION(RWSTRING("nCut"));
    RWASSERT(root);
    RWASSERT(set1);
    RWASSERT(set2);

    nVol(&whole, root, set1);
    maxr = nMaximize(root, set1, RED, &cutr, &whole);
    maxg = nMaximize(root, set1, GREEN, &cutg, &whole);
    maxb = nMaximize(root, set1, BLUE, &cutb, &whole);
    maxa = nMaximize(root, set1, ALPHA, &cuta, &whole);

    /* did we find any splits? */
    if ( (maxr == 0.0f) &&
         (maxg == 0.0f) &&
         (maxb == 0.0f) &&
         (maxa == 0.0f) )
    {
        RWRETURN(FALSE);
    }

    *set2 = *set1;

    /* NEED TO CHECK FOR ALPHA TOO */

    if (maxr >= maxg)
    {
        if (maxr >= maxb)
        {
            if (maxr >= maxa)
            {
                set1->col1.red = set2->col0.red = (RwUInt8) cutr;
            }
            else
            {
                set1->col1.alpha = set2->col0.alpha = (RwUInt8) cuta;
            }
        }
        else
        {
            if (maxb >= maxa)
            {
                set1->col1.blue = set2->col0.blue = (RwUInt8) cutb;
            }
            else
            {
                set1->col1.alpha = set2->col0.alpha = (RwUInt8) cuta;
            }
        }
    }
    else if (maxg >= maxb)
    {
        if (maxg >= maxa)
        {
            set1->col1.green = set2->col0.green = (RwUInt8) cutg;
        }
        else
        {
            set1->col1.alpha = set2->col0.alpha = (RwUInt8) cuta;
        }
    }
    else if (maxb >= maxa)
    {
        set1->col1.blue = set2->col0.blue = (RwUInt8) cutb;
    }
    else
    {
        set1->col1.alpha = set2->col0.alpha = (RwUInt8) cuta;
    }


    RWRETURN(TRUE);
}

/*************************************************************************/
#ifdef CACHEWEIGHTS
static LeafNode    *
CalcNodeWeights(OctNode * root, int depth)
{
    LeafNode           *Leaf;
    int                 i;

    RWFUNCTION(RWSTRING("CalcNodeWeights"));

    Leaf = NULL;
    if (root)
    {
        Leaf = &root->Leaf;

        /* is it a branch? */
        if (depth > 0)
        {
            InitLeaf(Leaf);
            for (i = 0; i < 16; i++)
            {
                LeafNode           *SubNode =
                    CalcNodeWeights(root->Branch.dir[i], depth - 1);

                if (SubNode)
                {
                    Leaf->weight += SubNode->weight;
                    RwRGBARealAdd(&Leaf->ac, &Leaf->ac, &SubNode->ac);
                    Leaf->m2 += SubNode->m2;
                }
            }
        }
    }

    RWRETURN(Leaf);
}
#endif

/*************************************************************************/
static int
CountLeafs(OctNode * root, int depth)
{
    int                 i, n;

    RWFUNCTION(RWSTRING("CountLeafs"));

    n = 0;
    if (root)
    {
        /* is it a branch? */
        if (depth > 0)
        {
            for (i = 0; i < 16; i++)
            {
                n += CountLeafs(root->Branch.dir[i], depth - 1);
            }
        }
        else
        {
            n = 1;
        }
    }

    RWRETURN(n);
}

/*************************************************************************/
static int
CountNodes(OctNode * root, int depth)
{
    int                 i, n;

    RWFUNCTION(RWSTRING("CountNodes"));

    n = 0;
    if (root)
    {
        n = 1;

        /* is it a branch? */
        if (depth > 0)
        {
            for (i = 0; i < 16; i++)
            {
                n += CountNodes(root->Branch.dir[i], depth - 1);
            }
        }
    }

    RWRETURN(n);
}

/*************************************************************************/

/* 
 * Note the use of 255.9999f value when generating the scale variable.
 * This is used instead of 255.0f to prevent rounding errors.
 */

#define node2pal(rgb, node)                                              \
MACRO_START                                                              \
{                                                                        \
    RwInt32             quantize;                                        \
    RwReal              scale = ( ((node)->weight > 0) ?                 \
                                  (255.9999f / (node)->weight) :         \
                                  (RwReal) 0 );                          \
                                                                         \
    quantize = (RwInt32) ((node)->ac.red * scale);                       \
    (rgb)->red = (RwUInt8) quantize;                                     \
    quantize = (RwInt32) ((node)->ac.green * scale);                     \
    (rgb)->green = (RwUInt8) quantize;                                   \
    quantize = (RwInt32) ((node)->ac.blue * scale);                      \
    (rgb)->blue = (RwUInt8) quantize;                                    \
    quantize = (RwInt32) ((node)->ac.alpha * scale);                     \
    (rgb)->alpha = (RwUInt8) quantize;                                   \
}                                                                        \
MACRO_STOP

/*************************************************************************/
RwInt32
_rwPalQuantResolvePalette(RwRGBA * palette, RwInt32 maxcols, rwPalQuant * pq)
{
    RwInt32             numcols, uniquecols, i, k, next;

    RWFUNCTION(RWSTRING("_rwPalQuantResolvePalette"));

#if (defined(RWDEBUG))
    if (256 < maxcols)
    {
        RWMESSAGE(("256 < %d == maxcols", (int) maxcols));
    }
#endif /* (defined(RWDEBUG)) */


    RWASSERT(palette);
    RWASSERT(maxcols <= 256);
    RWASSERT(pq);

    numcols = maxcols;
    uniquecols = CountLeafs(pq->root, MAXDEPTH);
    if (numcols > uniquecols)
    {
        numcols = uniquecols;
    }

#ifdef CACHEWEIGHTS
    /* cache weightings at every node */
    CalcNodeWeights(pq->root, MAXDEPTH);
#endif

    /* divide and conquer */
    pq->Mcube[0].col0.red = 0;
    pq->Mcube[0].col0.green = 0;
    pq->Mcube[0].col0.blue = 0;
    pq->Mcube[0].col0.alpha = 0;
    pq->Mcube[0].col1.red = 1 << MAXDEPTH;
    pq->Mcube[0].col1.green = 1 << MAXDEPTH;
    pq->Mcube[0].col1.blue = 1 << MAXDEPTH;
    pq->Mcube[0].col1.alpha = 1 << MAXDEPTH;
    pq->Mvv[0] = nVar(pq->root, &pq->Mcube[0]);
    for (i = 1; i < numcols; i++)
    {
        float               maxvar;

        /* find best box to split */
        next = -1;
        maxvar = 0.0f;
        for (k = 0; k < i; k++)
        {
            if (pq->Mvv[k] > maxvar)
            {
                maxvar = pq->Mvv[k];
                next = k;
            }
        }

        /* stop if we couldn't find a box to split */
        if (next == -1)
        {
            break;
        }

        /* split box */
        if (nCut(pq->root, &pq->Mcube[next], &pq->Mcube[i]))
        {
            /* volume test ensures we won't try to cut one-cell box */
            pq->Mvv[next] = nVar(pq->root, &pq->Mcube[next]);
            pq->Mvv[i] = nVar(pq->root, &pq->Mcube[i]);
        }
        else
        {
            /* don't try to split this box again */
            pq->Mvv[next] = 0.0f;
            i--;
        }
    }

    /* extract the new palette */
    for (k = 0; k < maxcols; k++)
    {
        if (k < numcols)
        {
            LeafNode            Node;

            nAssign(k, pq->root, &pq->Mcube[k]);
            nVol(&Node, pq->root, &pq->Mcube[k]);
            node2pal(&palette[k], &Node);
        }
        else
        {
            palette[k].red = 0;
            palette[k].green = 0;
            palette[k].blue = 0;
            palette[k].alpha = 0;
        }
    }

    RWRETURN(numcols);
}

/*************************************************************************/
static              RwUInt8
GetIndex(OctNode * root, OctantMap Octs, int depth)
{
    RwUInt8             result;

    RWFUNCTION(RWSTRING("GetIndex"));
    RWASSERT(root);

    if (depth == 0)
    {
        result = root->Leaf.palIndex;
    }
    else
    {
        result = GetIndex(root->Branch.dir[Octs & 15], Octs >> 4, depth - 1);
    }

    RWRETURN(result);
}

void
_rwPalQuantMatchImage(RwUInt8 * dstpixels, RwInt32 dststride, RwInt32 dstdepth,
                     RwBool dstPacked, rwPalQuant * pq, RwImage * img)
{
    RwUInt32            width, x, height, stride, maxcol;
    RwUInt8            *pixels, *dstLinePixels;
    RwUInt8             index;
    OctantMap           Octs;

    RWFUNCTION(RWSTRING("_rwPalQuantMatchImage"));
    RWASSERT(dstpixels);
    RWASSERT(pq);
    RWASSERT(img);

    stride = RwImageGetStride(img);
    pixels = RwImageGetPixels(img);
    maxcol = 1 << dstdepth;

    /* store two pixels per byte if the destination is packed
     * (i.e is a raster not an image) */
    if ((dstdepth == 4) && (dstPacked != FALSE))
    {
        switch (RwImageGetDepth(img))
        {
        case 4:
        case 8:
            {
                RwRGBA             *palette = RwImageGetPalette(img);
                height = RwImageGetHeight(img);
                while (height--)
                {
                    RwUInt8            *srcLinePixels = pixels;

                    dstLinePixels = dstpixels;

                    width = RwImageGetWidth(img);
                    for (x = 0; x < width; x++)
                    {
                        RwRGBA *color = &palette[*srcLinePixels++];

                        /* map to new index */
                        Octs = GetOctAdr(color);
                        index = GetIndex(pq->root, Octs, MAXDEPTH);
                        RWASSERT(index < maxcol);

                        if ((x/*+raster->nOffsetX*/) & 1)
                        {
                           *dstLinePixels &= 0x0F;
                           *dstLinePixels |= (index & 0x0F) << 4;
                            dstLinePixels++;
                        }
                        else
                        {
                           *dstLinePixels &= 0xF0;
                           *dstLinePixels |= index & 0x0F;
                        }
                    }

                    pixels += stride;
                    dstpixels += dststride;
                }
                break;
            }
        case 32:
            {
                height = RwImageGetHeight(img);
                while (height--)
                {
                    RwRGBA             *srcLinePixels = (RwRGBA *) pixels;

                    dstLinePixels = dstpixels;

                    width = RwImageGetWidth(img);
                    for (x = 0; x < width; x++)
                    {
                        RwRGBA *color = srcLinePixels++;

                        /* map to new index */
                        Octs = GetOctAdr(color);
                        index = GetIndex(pq->root, Octs, MAXDEPTH);
                        RWASSERT(index < maxcol);

                        if ((x/*+raster->nOffsetX*/) & 1)
                        {
                           *dstLinePixels &= 0x0F;
                           *dstLinePixels |= (index & 0x0F) << 4;
                            dstLinePixels++;
                        }
                        else
                        {
                           *dstLinePixels &= 0xF0;
                           *dstLinePixels |= index & 0x0F;
                        }
                    }

                    pixels += stride;
                    dstpixels += dststride;
                }
                break;
            }
        }
    }
    /* store one pixel per byte */
    else
    {
        switch (RwImageGetDepth(img))
        {
        case 4:
        case 8:
            {
                RwRGBA             *palette = RwImageGetPalette(img);

                height = RwImageGetHeight(img);
                while (height--)
                {
                    RwUInt8            *srcLinePixels = pixels;

                    dstLinePixels = dstpixels;

                    width = RwImageGetWidth(img);
                    while (width--)
                    {
                        RwRGBA *color = &palette[*srcLinePixels++];

                        /* map to new index */
                        Octs = GetOctAdr(color);
                        index = GetIndex(pq->root, Octs, MAXDEPTH);
                        RWASSERT(index < maxcol);
                        *dstLinePixels++ = index;
                    }

                    pixels += stride;
                    dstpixels += dststride;
                }
                break;
            }
        case 32:
            {
                height = RwImageGetHeight(img);
                while (height--)
                {
                    RwRGBA             *srcLinePixels = (RwRGBA *) pixels;

                    dstLinePixels = dstpixels;

                    width = RwImageGetWidth(img);
                    while (width--)
                    {
                        RwRGBA             *color = srcLinePixels++;

                        /* map to new index */
                        Octs = GetOctAdr(color);
                        index = GetIndex(pq->root, Octs, MAXDEPTH);
                        RWASSERT(index < maxcol);
                        *dstLinePixels++ = index;
                    }

                    pixels += stride;
                    dstpixels += dststride;
                }
                break;
            }
        }
    }

    RWRETURNVOID();
}

/*************************************************************************/
RwBool
_rwPalQuantInit(rwPalQuant * pq)
{
    int                 i, j, maxval;

    RWFUNCTION(RWSTRING("_rwPalQuantInit"));
    RWASSERT(pq);

    /* lookup mapping (8) bit-patterns to every 4th bit b31->b00 (least to most) */
    maxval = 1 << MAXDEPTH;
    for (i = 0; i < maxval; i++)
    {
        OctantMap           mask = 0;

        for (j = 0; j < MAXDEPTH; j++)
        {
            mask |= (i & (1 << j)) ? (1 << ((MAXDEPTH - 1 - j) * 4)) : 0;
        }
        splice[i] = mask;
    }

    pq->Mcube = (box *) RwCalloc(sizeof(box), MAXCOLOR);
    pq->Mvv = (float *) RwCalloc(sizeof(float), MAXCOLOR);

    pq->cubefreelist = RwFreeListCreate(sizeof(OctNode), 64, 0);
    pq->root = CreateCube(pq->cubefreelist);
    InitBranch(&pq->root->Branch);

    RWRETURN(TRUE);
}

/*************************************************************************/
static void
DeleteOctTree(rwPalQuant * pq, OctNode * root, int depth)
{
    int                 i;

    RWFUNCTION(RWSTRING("DeleteOctTree"));
    RWASSERT(pq);

    if (root)
    {
        /* is it a branch? */
        if (depth > 0)
        {
            for (i = 0; i < 16; i++)
            {
                DeleteOctTree(pq, root->Branch.dir[i], depth - 1);
            }
        }

        RwFreeListFree(pq->cubefreelist, root);
    }

    RWRETURNVOID();
}

/*************************************************************************/
void
_rwPalQuantTerm(rwPalQuant * pq)
{

    RWFUNCTION(RWSTRING("_rwPalQuantTerm"));
    RWASSERT(pq);

    DeleteOctTree(pq, pq->root, MAXDEPTH);
    pq->root = (OctNode *)NULL;

    RwFreeListDestroy(pq->cubefreelist);
    RwFree(pq->Mvv);
    RwFree(pq->Mcube);

    RWRETURNVOID();
}

/*************************************************************************/
#ifdef AMB_SPECIFIC
#include <stdio.h>



void
_rwPaletteCopy(RwImage * img, RwRGBA * newpal)
{
    RwRGBA             *palette;
    RwInt32             i, maxcol;

    RWFUNCTION(RWSTRING("_rwPaletteCopy"));

    palette = RwImageGetPalette(img);
    maxcol = 1 << RwImageGetDepth(img);
    for (i = 0; i < maxcol; i++)
    {
        palette[i] = newpal[i];
    }

    RWRETURNVOID();
}

void
_rwtest(void)
{
    RwImage            *im, *img[12], *img8[12];
    RwRGBA              palette[256];
    RwInt32             i, w, h, nummips;
    rwPalQuant          pq;

    RWFUNCTION(RWSTRING("_rwtest"));

    _rwPalQuantInit(&pq);

    RwImageSetPath("dffs\\");
    im = RwImageRead("Dhabih.bmp");
    img[0] = RwImageCreate(RwImageGetWidth(im), RwImageGetHeight(im), 32);
    RwImageAllocatePixels(img[0]);
    RwImageCopy(img[0], im);
    RwImageDestroy(im);

    i = 1;
    w = RwImageGetWidth(img[0]);
    h = RwImageGetHeight(img[0]);
    while ((w > 8) && (h > 8))
    {
        w >>= 1;
        h >>= 1;
        img[i] = RwImageCreate(w, h, 32);
        RwImageAllocatePixels(img[i]);
        RwImageResample(img[i], img[i - 1]);
        i++;
    }
    nummips = i;

    for (i = 0; i < nummips; i++)
    {
        _rwPalQuantAddImage(&pq, img[i], (RwReal)(1 << (i << 1)));
    }

    {
        char                buff[128];

        sprintf(buff, "%d unique colours (%d nodes using %d bytes)\n",
                CountLeafs(pq.root, MAXDEPTH),
                CountNodes(pq.root, MAXDEPTH), sizeof(OctNode));
        //OutputDebugString(buff);
    }
    _rwPalQuantResolvePalette(palette, 256, &pq);

    for (i = 0; i < nummips; i++)
    {
        img8[i] =
            RwImageCreate(RwImageGetWidth(img[i]), RwImageGetHeight(img[i]),
                          8);
        RwImageAllocatePixels(img8[i]);
        _rwPaletteCopy(img8[i], palette);
        _rwPalQuantMatchImage(RwImageGetPixels(img8[i]),
                             RwImageGetStride(img8[i]),
                             RwImageGetDepth(img8[i]), &pq, img[i]);

        RwImageDestroy(img[i]);
    }

    for (i = 0; i < nummips; i++)
    {
        RwImage            *img32;
        char                filename[128];

        /* write it */
        img32 =
            RwImageCreate(RwImageGetWidth(img8[i]), RwImageGetHeight(img8[i]),
                          32);
        RwImageAllocatePixels(img32);
        RwImageCopy(img32, img8[i]);
        sprintf(filename, "test%03d.bmp", i);
        RwImageWriteBMP(img32, filename);
        RwImageDestroy(img32);

        RwImageDestroy(img8[i]);
    }

    _rwPalQuantTerm(&pq);

    RWRETURNVOID();
}
#endif

#if defined (__MWERKS__)
#if (defined(RWVERBOSE))
#pragma message (__FILE__ "/" _SKY_EXPAND(__LINE__) ": __MWERKS__ == " _SKY_EXPAND(__MWERKS__))
#endif /* (defined (__MWERKS__)) */
#if (__option (global_optimizer))
#pragma always_inline on
#endif /* (__option (global_optimizer)) */
#endif /*  defined (__MWERKS__) */
