/*********************************************************************\ Unit to find the best palette of a PICT in RGB555. The palette is given back with color RGB555. References in order of importance: 1) C Implementation of Wu's Color Quantizer (v. 2) (see Graphics Gems vol. II, pp. 126-133) 2) Paul Heckbert's paper "Color Image Quantization for Frame Buffer Display", SIGGRAPH '82 Proceedings, page 297 3) Source of ppmquant.c Copyright (C) 1989, 1991 by Jef Poskanzer. 95/02/15 Adapted for 1) - Laurent \*********************************************************************/ #include #include #include #define MAXCOLOR 256 #define RED 2 #define GREEN 1 #define BLUE 0 typedef float fint32; typedef int int32; typedef unsigned int uint32; typedef short int16; typedef unsigned short uint16; typedef char int8; typedef unsigned char uint8; // clamp the input to the specified range #define CLAMP(v,l,h) ((v)<(l)?(l):(v)>(h)?(h):v) typedef struct TBox { int32 r0; /* min value, exclusive */ int32 r1; /* max value, inclusive */ int32 g0; int32 g1; int32 b0; int32 b1; int32 vol; } TBox; /* At conclusion of the histogram step, we can interpret * wt[r][g][b] = sum over voxel of P(c) * mr[r][g][b] = sum over voxel of r*P(c) , similarly for mg, mb * m2[r][g][b] = sum over voxel of c^2*P(c) * Actually each of these should be divided by 'size' to give the usual * interpretation of P() as ranging from 0 to 1, but we needn't do that here. */ typedef int32 Tint32Table[33][33][33]; typedef fint32 Tfint32Table[33][33][33]; typedef struct STHisto { Tfint32Table m2; // here m2 should be float for huge picture (more than 800x800) Tint32Table wt; // Tint32Table mr; // Tint32Table mg; // Tint32Table mb; // } THisto; /* Histogram is in elements 1..HISTSIZE along each axis, * element 0 is for base or marginal value * NB: these must start out 0! */ static THisto irs_Histo; static THisto *gpHisto = &irs_Histo; static uint8 gTag[33*33*33]; static TBox gcube[MAXCOLOR+1]; static fint32 vv[MAXCOLOR+1]; static int32 gNbCols; /*********************************************************************\ IMPLEMENTATION \*********************************************************************/ /* build 3-D color histogram of counts, r/g/b, c^2 */ // size is the number of pixel in the picture static void BestPal_Hist3d(uint16 *pImage16b, int32 size, uint16 *pMap) { int32 *vwt, *vmr, *vmg, *vmb; fint32 *vm2; int32 ind, r, g, b; int32 i; vm2 = (fint32*)&gpHisto->m2[0]; vwt = (int32*)&gpHisto->wt[0]; vmr = (int32*)&gpHisto->mr[0]; vmg = (int32*)&gpHisto->mg[0]; vmb = (int32*)&gpHisto->mb[0]; memset( (void*)gpHisto, 0, sizeof(THisto) ); if (pMap) { for(i=0; i>10) & 0x1F; g = (r>>5) & 0x1F; r = r & 0x1F; *pMap++=ind=(r+1)*33*33+(g+1)*33+(b+1); vwt[ind]++; vmr[ind] += r; vmg[ind] += g; vmb[ind] += b; vm2[ind] += (fint32)(r*r+g*g+b*b); // could be a lookup table but image is small } else { *pMap++=0; /* transparent ! */ } } } else { for(i=0; i>5) & 0x1F; r = (r>>10) & 0x1F; ind=(r+1)*33*33+(g+1)*33+(b+1); ++vwt[ind]; vmr[ind] += r; vmg[ind] += g; vmb[ind] += b; vm2[ind] += (fint32)(r*r+g*g+b*b); // could be a lookup table but image is small } } } /* We now convert histogram into moments so that we can rapidly calculate * the sums of the above quantities over any desired box. */ /* compute cumulative moments. */ static void BestPal_M3d(void) { int32 *vwt, *vmr, *vmg, *vmb; fint32 *vm2; uint16 ind1, ind2; uint8 i, r, g, b; int32 line, line_r, line_g, line_b; int32 area[33], area_r[33], area_g[33], area_b[33]; fint32 line2, area2[33]; vm2 = (fint32*)&gpHisto->m2[0]; vwt = (int32*)&gpHisto->wt[0]; vmr = (int32*)&gpHisto->mr[0]; vmg = (int32*)&gpHisto->mg[0]; vmb = (int32*)&gpHisto->mb[0]; for(r=1; r<=32; ++r) { for(i=0; i<=32; ++i) { area2[i]=area[i]=area_r[i]=area_g[i]=area_b[i]=0; } for(g=1; g<=32; ++g) { line2 = line = line_r = line_g = line_b = 0; for(b=1; b<=32; ++b) { ind1 = r*33*33+g*33+b; line += vwt[ind1]; line_r += vmr[ind1]; line_g += vmg[ind1]; line_b += vmb[ind1]; line2 += vm2[ind1]; area[b] += line; area_r[b] += line_r; area_g[b] += line_g; area_b[b] += line_b; area2[b] += line2; ind2 = ind1 - 33*33; /* [r-1][g][b] */ vwt[ind1] = vwt[ind2] + area[b]; vmr[ind1] = vmr[ind2] + area_r[b]; vmg[ind1] = vmg[ind2] + area_g[b]; vmb[ind1] = vmb[ind2] + area_b[b]; vm2[ind1] = vm2[ind2] + area2[b]; } } } } /* Compute sum over a box of any given statistic */ static int32 BestPal_Vol(TBox *cube, Tint32Table mmt) { return( mmt[cube->r1][cube->g1][cube->b1] -mmt[cube->r1][cube->g1][cube->b0] -mmt[cube->r1][cube->g0][cube->b1] +mmt[cube->r1][cube->g0][cube->b0] -mmt[cube->r0][cube->g1][cube->b1] +mmt[cube->r0][cube->g1][cube->b0] +mmt[cube->r0][cube->g0][cube->b1] -mmt[cube->r0][cube->g0][cube->b0] ); } /* The next two routines allow a slightly more efficient calculation * of BestPal_Vol() for a proposed subbox of a given box. The sum of Top() * and Bottom() is the BestPal_Vol() of a subbox split in the given direction * and with the specified new upper bound. */ /* Compute part of BestPal_Vol(cube, mmt) that doesn't depend on r1, g1, or b1 */ /* (depending on dir) */ static int32 BestPal_Bottom(TBox *cube, uint8 dir, Tint32Table mmt) { switch(dir) { case RED: return( -mmt[cube->r0][cube->g1][cube->b1] +mmt[cube->r0][cube->g1][cube->b0] +mmt[cube->r0][cube->g0][cube->b1] -mmt[cube->r0][cube->g0][cube->b0] ); break; case GREEN: return( -mmt[cube->r1][cube->g0][cube->b1] +mmt[cube->r1][cube->g0][cube->b0] +mmt[cube->r0][cube->g0][cube->b1] -mmt[cube->r0][cube->g0][cube->b0] ); break; case BLUE: return( -mmt[cube->r1][cube->g1][cube->b0] +mmt[cube->r1][cube->g0][cube->b0] +mmt[cube->r0][cube->g1][cube->b0] -mmt[cube->r0][cube->g0][cube->b0] ); break; } return 0; } /* Compute remainder of BestPal_Vol(cube, mmt), substituting pos for */ /* r1, g1, or b1 (depending on dir) */ static int32 BestPal_Top(TBox *cube, uint8 dir, int32 pos, Tint32Table mmt) { switch(dir) { case RED: return( mmt[pos][cube->g1][cube->b1] -mmt[pos][cube->g1][cube->b0] -mmt[pos][cube->g0][cube->b1] +mmt[pos][cube->g0][cube->b0] ); break; case GREEN: return( mmt[cube->r1][pos][cube->b1] -mmt[cube->r1][pos][cube->b0] -mmt[cube->r0][pos][cube->b1] +mmt[cube->r0][pos][cube->b0] ); break; case BLUE: return( mmt[cube->r1][cube->g1][pos] -mmt[cube->r1][cube->g0][pos] -mmt[cube->r0][cube->g1][pos] +mmt[cube->r0][cube->g0][pos] ); break; } return 0; } /* Compute the weighted variance of a box */ /* NB: as with the raw statistics, this is really the variance * size */ static fint32 BestPal_Var(TBox *cube, Tfint32Table m2) { fint32 dr, dg, db, wt, xx, mm; dr = BestPal_Vol(cube, gpHisto->mr); dg = BestPal_Vol(cube, gpHisto->mg); db = BestPal_Vol(cube, gpHisto->mb); mm = dr*dr+dg*dg+db*db; wt = BestPal_Vol(cube, gpHisto->wt); xx = m2[cube->r1][cube->g1][cube->b1] -m2[cube->r1][cube->g1][cube->b0] -m2[cube->r1][cube->g0][cube->b1] +m2[cube->r1][cube->g0][cube->b0] -m2[cube->r0][cube->g1][cube->b1] +m2[cube->r0][cube->g1][cube->b0] +m2[cube->r0][cube->g0][cube->b1] -m2[cube->r0][cube->g0][cube->b0]; return ( xx - mm / wt ); } /* 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 fint32 BestPal_Maximize(TBox *cube, uint8 dir, int32 first, int32 last, int32 *cut, int32 whole_r, int32 whole_g, int32 whole_b, int32 whole_w) { int32 half_r, half_g, half_b, half_w; int32 base_r, base_g, base_b, base_w; int32 i; fint32 temp, max; base_r = BestPal_Bottom(cube, dir, gpHisto->mr); base_g = BestPal_Bottom(cube, dir, gpHisto->mg); base_b = BestPal_Bottom(cube, dir, gpHisto->mb); base_w = BestPal_Bottom(cube, dir, gpHisto->wt); max = 0.0; *cut = -1; for(i=first; imr); half_g = base_g + BestPal_Top(cube, dir, i, gpHisto->mg); half_b = base_b + BestPal_Top(cube, dir, i, gpHisto->mb); half_w = base_w + BestPal_Top(cube, dir, i, gpHisto->wt); /* now half_x is sum over lower half of box, if split at i */ if (half_w == 0) { /* subbox could be empty of pixels! */ continue; /* never split into an empty box */ } else { temp = ((fint32)half_r*half_r + (fint32)half_g*half_g + (fint32)half_b*half_b)/half_w; } half_r = whole_r - half_r; half_g = whole_g - half_g; half_b = whole_b - half_b; half_w = whole_w - half_w; if (half_w == 0) { /* subbox could be empty of pixels! */ continue; /* never split into an empty box */ } else { temp += ((fint32)half_r*half_r + (fint32)half_g*half_g + (fint32)half_b*half_b)/half_w; } if (temp > max) { *cut = i; max=temp; } else if (temp < max) { return(max); } } return(max); } static int32 BestPal_Cut(TBox *set1, TBox *set2) { uint8 dir; int32 cutr, cutg, cutb; fint32 maxr, maxg, maxb; int32 whole_r, whole_g, whole_b, whole_w; whole_r = BestPal_Vol(set1, gpHisto->mr); whole_g = BestPal_Vol(set1, gpHisto->mg); whole_b = BestPal_Vol(set1, gpHisto->mb); whole_w = BestPal_Vol(set1, gpHisto->wt); maxr = BestPal_Maximize(set1, RED, set1->r0+1, set1->r1, &cutr, whole_r, whole_g, whole_b, whole_w); maxg = BestPal_Maximize(set1, GREEN, set1->g0+1, set1->g1, &cutg, whole_r, whole_g, whole_b, whole_w); maxb = BestPal_Maximize(set1, BLUE, set1->b0+1, set1->b1, &cutb, whole_r, whole_g, whole_b, whole_w); if( (maxr>=maxg)&&(maxr>=maxb) ) { dir = RED; if (cutr < 0) return 0; /* can't split the box */ } else { if ( (maxg>=maxr)&&(maxg>=maxb) ) dir = GREEN; else dir = BLUE; } set2->r1 = set1->r1; set2->g1 = set1->g1; set2->b1 = set1->b1; switch (dir) { case RED: set2->r0 = set1->r1 = cutr; set2->g0 = set1->g0; set2->b0 = set1->b0; break; case GREEN: set2->g0 = set1->g1 = cutg; set2->r0 = set1->r0; set2->b0 = set1->b0; break; case BLUE: set2->b0 = set1->b1 = cutb; set2->r0 = set1->r0; set2->g0 = set1->g0; break; } set1->vol=(set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0); set2->vol=(set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0); return 1; } static void BestPal_Mark(TBox *cube, int8 label) { int32 r, g, b; for(r=cube->r0+1; r<=cube->r1; ++r) for(g=cube->g0+1; g<=cube->g1; ++g) for(b=cube->b0+1; b<=cube->b1; ++b) gTag[r*33*33+g*33+b] = label; } /*********************************************************************\ INPUT PARAMETER: ---------------- pSrcImage16bits Point to the source image in a RAW format RGB555 PixelNumber Is the number of pixel in the source image pPLUT Point to an array of uint16 This array must be at least NbWantedColors wide This array returned the best PLUT NbWantedColors Specify the number of color you request. Maximum color is 256 pDestImage16bits Specify the Destination picture. After execution, each entry contain an index to the PLUT This array must have the same size than the original image because it is used as a temporary array during calculation \*********************************************************************/ uint32 BestPal_BestPalAndRemap( uint16 *pSrcImage16bits, uint32 PixelNumber, uint16 *pPLUT, int32 NbWantedColors, uint16 *pDestImage16bits ) { fint32 temp; int32 lut_r, lut_g, lut_b; int32 next; int32 i, weight; int32 k; if (NbWantedColors > MAXCOLOR) { printf("cannot look for more than %d color but %d asked\n",MAXCOLOR, NbWantedColors); return 0; } gNbCols = NbWantedColors; BestPal_Hist3d(pSrcImage16bits, PixelNumber, pDestImage16bits); BestPal_M3d(); gcube[0].r0 = gcube[0].g0 = gcube[0].b0 = 0; gcube[0].r1 = gcube[0].g1 = gcube[0].b1 = 32; next = 0; for(i=1; i1) ? BestPal_Var(&gcube[next], gpHisto->m2) : 0.0; vv[i] = (gcube[i].vol>1) ? BestPal_Var(&gcube[i], gpHisto->m2) : 0.0; } else { vv[next] = 0.0; /* don't try to split this box again */ i--; /* didn't create box i */ } next = 0; temp = vv[0]; for(k=1; k<=i; k++) { if (vv[k] > temp) { temp = vv[k]; next = k; } } if (temp <= 0.0) { gNbCols = i+2; break; } } for(k=0; kwt); if (weight) { lut_r = BestPal_Vol(&gcube[k], gpHisto->mr) / weight; lut_g = BestPal_Vol(&gcube[k], gpHisto->mg) / weight; lut_b = BestPal_Vol(&gcube[k], gpHisto->mb) / weight; } else { lut_r = lut_g = lut_b = 0; } lut_r = CLAMP(lut_r,0,31); lut_g = CLAMP(lut_g,0,31); lut_b = CLAMP(lut_b,0,31); pPLUT[k+1] = 0x8000 + (lut_b<<10) + (lut_g<<5) + lut_r; } for(i=0; i MAXCOLOR) { printf("cannot look for more than %d color but %d asked\n",MAXCOLOR, NbWantedColors); return 0; } gNbCols = NbWantedColors; BestPal_Hist3d(pSrcImage16bits, PixelNumber, pDestImage16bits); BestPal_M3d(); gcube[0].r0 = gcube[0].g0 = gcube[0].b0 = 0; gcube[0].r1 = gcube[0].g1 = gcube[0].b1 = 32; next = 0; for(i=1; i1) ? BestPal_Var(&gcube[next], gpHisto->m2) : 0.0; vv[i] = (gcube[i].vol>1) ? BestPal_Var(&gcube[i], gpHisto->m2) : 0.0; } else { vv[next] = 0.0; /* don't try to split this box again */ i--; /* didn't create box i */ } next = 0; temp = vv[0]; for(k=1; k<=i; ++k) { if (vv[k] > temp) { temp = vv[k]; next = k; } } if (temp <= 0.0) { gNbCols = i+1; break; } } for(k=0; kwt); if (weight) { lut_r = BestPal_Vol(&gcube[k], gpHisto->mr) / weight; lut_g = BestPal_Vol(&gcube[k], gpHisto->mg) / weight; lut_b = BestPal_Vol(&gcube[k], gpHisto->mb) / weight; } else { lut_r = lut_g = lut_b = 0; } lut_r = CLAMP(lut_r,0,31); lut_g = CLAMP(lut_g,0,31); lut_b = CLAMP(lut_b,0,31); pPLUT[k] = 0x8000 + (lut_b<<10) + (lut_g<<5) + lut_r; } for(i=0; i