Picture Viewer

Development Platform:
Visual C++

MEDIAN.C：Code Content
							/*

** File:        median.c             Copyright (c) Truda Software

** Author:      Anton Kruger         215 Marengo Rd, #2, 

** Date:        March 1992           Oxford, IA 52322-9383

** Revision:    1.0                  March 1992

** 

** Description: Contains an implementation of Heckbert's median-

**              cut color quantization algorithm.

**

** Compilers:  MSC 5.1, 6.0.

**

** Note:       1) Compile in large memory model.

**             2) Delete the "#define FAST_REMAP" statement below

**                in order to deactivate fast remapping.

*/

#include <windows.h>

#define  FAST_REMAP

#include <stdio.h>

#include <stddef.h>              /* for NULL                   */

#include <stdlib.h>              /* for "qsort"                */

#include <float.h>               /* for FLT_MAX, FLT_MIN       */

#define MAXCOLORS 256            /* maximum # of output colors */

#define HSIZE     32768          /* size of image histogram    */

//typedef unsigned  char BYTE;     /* range 0-255                */

//typedef unsigned  short WORD;    /* range 0-65,535             */

//typedef unsigned  long DWORD;    /* range 0-4,294,967,295      */

/* Macros for converting between (r,g,b)-colors and 15-bit     */

/* colors follow.                                              */

#define _RGB(r,g,b) (WORD)(((b)&~7)<<7)|(((g)&~7)<<2)|((r)>>3)

#define RED(x)     (BYTE)(((x)&31)<<3)

#define GREEN(x)   (BYTE)((((x)>>5)&255)<< 3)

#define BLUE(x)    (BYTE)((((x)>>10)&255)<< 3)

typedef  struct {       /* structure for a cube in color space */

   WORD  lower;         /* one corner's index in histogram     */

   WORD  upper;         /* another corner's index in histogram */

   DWORD count;         /* cube's histogram count              */

   int   level;         /* cube's level                        */

   BYTE  rmin,rmax;   

   BYTE  gmin,gmax;   

   BYTE  bmin,bmax;   

} cube_t;

static cube_t list[MAXCOLORS];   /* list of cubes              */

static int longdim;              /* longest dimension of cube  */

static WORD HistPtr[HSIZE];      /* points to colors in "Hist" */

void Shrink(cube_t * Cube);

void InvMap(WORD * Hist, BYTE ColMap[][3],WORD ncubes);

int  compare(const void * a1, const void * a2);

WORD __declspec (dllexport) MedianCut(WORD Hist[],BYTE ColMap[][3], int maxcubes)

{

   /*

   ** Accepts "Hist", a 32,768-element array that contains 15-bit

   ** color counts of the input image. Uses Heckbert's median-cut

   ** algorithm to divide the color space into "maxcubes" cubes,

   ** and returns the centroid (average value) of each cube in

   ** "ColMap". "Hist" is also updated so that it functions as an

   ** inverse color map. MedianCut returns the actual number of

   ** cubes, which may be less that "maxcubes".

   */

   BYTE        lr,lg,lb;

   WORD        i,median,color;

   DWORD       count;

   int         k,level,ncubes,splitpos;

   void        *base;

   size_t      num,width;

   cube_t      Cube,CubeA,CubeB;

   /*

   ** Create the initial cube, which is the whole RGB-cube.

   */

   ncubes = 0;

   Cube.count = 0;

   for (i=0,color=0;i<=HSIZE-1;i++){

      if (Hist[i] != 0){

         HistPtr[color++] = i;

         Cube.count = Cube.count + Hist[i];

      }

   }

   Cube.lower = 0; Cube.upper = color-1;

   Cube.level = 0;

   Shrink(&Cube);

   list[ncubes++] = Cube;

   /*

   ** The main loop follows.  Search the list of cubes for the

   ** next cube to split, which is the lowest level cube.  A

   ** special case is when a cube has only one color, so that it

   ** cannot be split.

   */

   while (ncubes < maxcubes){

      level = 255; splitpos = -1;

      for (k=0;k<=ncubes-1;k++){

         if (list[k].lower == list[k].upper)  

                  ;                            /* single color */

         else if (list[k].level < level){

            level = list[k].level;

            splitpos = k;

         }

      }

      if (splitpos == -1)            /* no more cubes to split */

         break;

      /*

      ** Must split the cube "splitpos" in the list of cubes.

      ** Next find the longest dimension of this cube, and update

      ** the external variable "longdim", which is used by the 

      ** sort routine so that it knows along which axis to sort.

      */

      Cube = list[splitpos];

      lr = Cube.rmax - Cube.rmin;

      lg = Cube.gmax - Cube.gmin;

      lb = Cube.bmax - Cube.bmin;

      if (lr >= lg && lr >= lb) longdim = 0;

      if (lg >= lr && lg >= lb) longdim = 1;

      if (lb >= lr && lb >= lg) longdim = 2;

      /*

      ** Sort along "longdim". This prepares for the next step,

      ** namely, finding the median. Use standard lib's "qsort".

      */

      base = (void *)&HistPtr[Cube.lower];

      num  = (size_t)(Cube.upper - Cube.lower + 1);

      width = (size_t)sizeof(HistPtr[0]);

      qsort(base,num,width,compare);

      /*

      ** Find median by scanning through the cube, and computing

      ** a running sum. When the running sum equals half the

      ** total for the cube, the median has been found.

      */

      count = 0;

      for (i=Cube.lower;i<=Cube.upper-1;i++){

         if (count >= Cube.count/2) break;

         color = HistPtr[i];

         count = count + Hist[color];

      }

      median = i;

      /*

      ** Now split "Cube" at the median. Then add the two new

      ** cubes to the list of cubes.

      */

      CubeA = Cube; CubeA.upper = median-1;

      CubeA.count = count;

      CubeA.level = Cube.level + 1;

      Shrink(&CubeA);

      list[splitpos] = CubeA;               /* add in old slot */

      CubeB = Cube; CubeB.lower = median; 

      CubeB.count = Cube.count - count;

      CubeB.level = Cube.level + 1;

      Shrink(&CubeB);

      list[ncubes++] = CubeB;               /* add in new slot */

      if ((ncubes % 10) == 0)

         fprintf(stderr,".");               /* pacifier        */

   }

   /*

   ** We have enough cubes, or we have split all we can. Now

   ** compute the color map, the inverse color map, and return

   ** the number of colors in the color map.

   */

   InvMap(Hist, ColMap,ncubes);

   return((WORD)ncubes);

}

void Shrink(cube_t * Cube)

{

   /*

   ** Encloses "Cube" with a tight-fitting cube by updating the

   ** (rmin,gmin,bmin) and (rmax,gmax,bmax) members of "Cube".

   */

   BYTE        r,g,b;

   WORD        i,color;

   Cube->rmin = 255; Cube->rmax = 0;

   Cube->gmin = 255; Cube->gmax = 0;

   Cube->bmin = 255; Cube->bmax = 0;

   for (i=Cube->lower;i<=Cube->upper;i++){

      color = HistPtr[i];

      r = RED(color);

      if (r > Cube->rmax) Cube->rmax = r;

      if (r < Cube->rmin) Cube->rmin = r;

      g = GREEN(color);

      if (g > Cube->gmax) Cube->gmax = g;

      if (g < Cube->gmin) Cube->gmin = g;

      b = BLUE(color);

      if (b > Cube->bmax) Cube->bmax = b;

      if (b < Cube->bmin) Cube->bmin = b;

   }

}

void InvMap(WORD * Hist, BYTE ColMap[][3],WORD ncubes)

{

   /*

   ** For each cube in the list of cubes, computes the centroid

   ** (average value) of the colors enclosed by that cube, and

   ** then loads the centroids in the color map. Next loads the

   ** histogram with indices into the color map. A preprocessor

   ** directive "#define FAST_REMAP" controls whether the cube

   ** centroids become the output color for all the colors in a

   ** cube, or whether a "best remap" method is followed.

   */

   BYTE        r,g,b;

   WORD        i,k,color;

   float       rsum,gsum,bsum;

   cube_t      Cube;

   for (k=0;k<=ncubes-1;k++){

      Cube = list[k];

      rsum = gsum = bsum = (float)0.0;

      for (i=Cube.lower;i<=Cube.upper;i++){

         color = HistPtr[i];

         r = RED(color);  

         rsum += (float)r*(float)Hist[color];

         g = GREEN(color);

         gsum += (float)g*(float)Hist[color];

         b = BLUE(color);

         bsum += (float)b*(float)Hist[color];

      }

      /* Update the color map */  

      ColMap[k][0] = (BYTE)(rsum/(float)Cube.count);

      ColMap[k][1] = (BYTE)(gsum/(float)Cube.count);

      ColMap[k][2] = (BYTE)(bsum/(float)Cube.count);

   }

#ifdef FAST_REMAP

   /*

   ** Fast remap: for each color in each cube, load the corre- 

   ** sponding slot in "Hist" with the centroid of the cube.

   */

   for (k=0;k<=ncubes-1;k++){

      Cube = list[k];

      for (i=Cube.lower;i<=Cube.upper;i++){

         color = HistPtr[i];

         Hist[color] = k;

      }

      if ((k % 10) == 0) fprintf(stderr,".");   /* pacifier    */

   }

#else

   /*

   ** Best remap: for each color in each cube, find the entry in

   ** "ColMap" that has the smallest Euclidian distance from the

   ** color. Record this information in "Hist".

   */

   for (k=0;k<=ncubes-1;k++){

      Cube = list[k];

      for (i=Cube.lower;i<=Cube.upper;i++){

         color = HistPtr[i];

         r = RED(color);  g = GREEN(color); b = BLUE(color);

         /* Search for closest entry in "ColMap" */

         dmin = (float)FLT_MAX;

         for (j=0;j<=ncubes-1;j++){

            dr = (float)ColMap[j][0] - (float)r;

            dg = (float)ColMap[j][1] - (float)g;

            db = (float)ColMap[j][2] - (float)b;

            d = dr*dr + dg*dg + db*db;

            if (d == (float)0.0){

               index = j; break;

            }

            else if (d < dmin){

               dmin = d; index = j;

            }

         }

         Hist[color] = index;

      }

      if ((k % 10) == 0) fprintf(stderr,".");   /* pacifier    */

   }

#endif

   return;

}

int compare(const void * a1, const void * a2)

{

   /*

   ** Called by the sort routine in "MedianCut". Compares two

   ** colors based on the external variable "longdim".

   */

   WORD        color1,color2;

   BYTE        c1,c2;

   color1 = (WORD)*(WORD *)a1;

   color2 = (WORD)*(WORD *)a2;

   switch (longdim){

      case 0:  

         c1 = RED(color1),  c2 = RED(color2);

         break;

      case 1:

         c1 = GREEN(color1), c2 = GREEN(color2);

         break;

      case 2:

         c1 = BLUE(color2), c2 = BLUE(color2);

         break;

   }

   return ((int)(c1-c2));

}