src/rt/colortab.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 * colortab.c - allocate and control dynamic color table.
 *
 *      We start off with a uniform partition of color space.
 *      As pixels are sent to the frame buffer, a histogram is built.
 *      When a new color table is requested, the histogram is used
 *      to make a pseudo-optimal partition, after which the
 *      histogram is cleared.  This algorithm
 *      performs only as well as the next drawing's color
 *      distribution is correlated to the last.
 *
 *  External symbols declared in drvier.h
 */

/* ====================================================================
 * The Radiance Software License, Version 1.0
 *
 * Copyright (c) 1990 - 2002 The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory.   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.
 *
 * 3. The end-user documentation included with the redistribution,
 *           if any, must include the following acknowledgment:
 *             "This product includes Radiance software
 *                 (http://radsite.lbl.gov/)
 *                 developed by the Lawrence Berkeley National Laboratory
 *               (http://www.lbl.gov/)."
 *       Alternately, this acknowledgment may appear in the software itself,
 *       if and wherever such third-party acknowledgments normally appear.
 *
 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
 *       and "The Regents of the University of California" must
 *       not be used to endorse or promote products derived from this
 *       software without prior written permission. For written
 *       permission, please contact [email protected].
 *
 * 5. Products derived from this software may not be called "Radiance",
 *       nor may "Radiance" appear in their name, without prior written
 *       permission of Lawrence Berkeley National Laboratory.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED 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 Lawrence Berkeley National Laboratory OR
 * ITS 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.
 * ====================================================================
 *
 * This software consists of voluntary contributions made by many
 * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
 * information on Lawrence Berkeley National Laboratory, please see
 * <http://www.lbl.gov/>.
 */

#include "standard.h"

#include "color.h"
                                /* histogram resolution */
#define NRED            24
#define NGRN            32
#define NBLU            16
#define HMAX            NGRN
                                /* minimum box count for adaptive partition */
#define MINSAMP         7
                                /* maximum distance^2 before color reassign */
#define MAXDST2         12
                                /* map a color */
#define map_col(c,p)    clrmap[p][ colval(c,p)<1. ? \
                                (int)(colval(c,p)*256.) : 255 ]
                                /* color partition tree */
#define CNODE           short
#define set_branch(p,c) ((c)<<2|(p))
#define set_pval(pv)    ((pv)<<2|3)
#define is_branch(cn)   (((cn)&3)!=3)
#define is_pval(cn)     (((cn)&3)==3)
#define part(cn)        ((cn)>>2)
#define prim(cn)        ((cn)&3)
#define pval(cn)        ((cn)>>2)
                                /* our color table */
static struct tabent {
        long    sum[3];         /* sum of colors using this entry */
        int     n;              /* number of colors */
        BYTE    ent[3];         /* current table value */
}       *clrtab = NULL;
                                /* color cube partition */
static CNODE    *ctree = NULL;
                                /* our color correction map */
static BYTE     clrmap[3][256];
                                /* histogram of colors used */
static unsigned short   histo[NRED][NGRN][NBLU];
                                /* initial color cube boundary */
static int      CLRCUBE[3][2] = {{0,NRED},{0,NGRN},{0,NBLU}};

static int      split();
static void     cut();


int
new_ctab(ncolors)               /* start new color table with max ncolors */
int     ncolors;
{
        int     treesize;

        if (ncolors < 1)
                return(0);
                                /* free old tables */
        if (clrtab != NULL)
                free((void *)clrtab);
        if (ctree != NULL)
                free((void *)ctree);
                                /* get new tables */
        for (treesize = 1; treesize < ncolors; treesize <<= 1)
                ;
        treesize <<= 1;
        clrtab = (struct tabent *)calloc(ncolors, sizeof(struct tabent));
        ctree = (CNODE *)malloc(treesize*sizeof(CNODE));
        if (clrtab == NULL || ctree == NULL)
                return(0);
                                /* partition color space */
        cut(ctree, 0, CLRCUBE, 0, ncolors);
                                /* clear histogram */
        bzero((char *)histo, sizeof(histo));
                                /* return number of colors used */
        return(ncolors);
}


int
get_pixel(col, set_pixel)       /* get pixel for color */
COLOR   col;
void    (*set_pixel)();
{
        int     r, g, b;
        int     cv[3];
        register CNODE  *tp;
        register int    h;
                                                /* map color */
        r = map_col(col,RED);
        g = map_col(col,GRN);
        b = map_col(col,BLU);
                                                /* reduce resolution */
        cv[RED] = (r*NRED)>>8;
        cv[GRN] = (g*NGRN)>>8;
        cv[BLU] = (b*NBLU)>>8;
                                                /* add to histogram */
        histo[cv[RED]][cv[GRN]][cv[BLU]]++;
                                                /* find pixel in tree */
        for (tp = ctree, h = 0; is_branch(*tp); h++)
                if (cv[prim(*tp)] < part(*tp))
                        tp += 1<<h;             /* left branch */
                else
                        tp += 1<<(h+1);         /* right branch */
        h = pval(*tp);
                                                /* add to color table */
        clrtab[h].sum[RED] += r;
        clrtab[h].sum[GRN] += g;
        clrtab[h].sum[BLU] += b;
        clrtab[h].n++;
                                        /* recompute average */
        r = clrtab[h].sum[RED] / clrtab[h].n;
        g = clrtab[h].sum[GRN] / clrtab[h].n;
        b = clrtab[h].sum[BLU] / clrtab[h].n;
                                        /* check for movement */
        if (clrtab[h].n == 1 ||
                        (r-clrtab[h].ent[RED])*(r-clrtab[h].ent[RED]) +
                        (g-clrtab[h].ent[GRN])*(g-clrtab[h].ent[GRN]) +
                        (b-clrtab[h].ent[BLU])*(b-clrtab[h].ent[BLU]) > MAXDST2) {
                clrtab[h].ent[RED] = r;
                clrtab[h].ent[GRN] = g; /* reassign pixel */
                clrtab[h].ent[BLU] = b;
#ifdef DEBUG
                sprintf(errmsg, "pixel %d = (%d,%d,%d) (%d refs)\n",
                                h, r, g, b, clrtab[h].n);
                eputs(errmsg);
#endif
                (*set_pixel)(h, r, g, b);
        }
        return(h);                              /* return pixel value */
}


void
make_gmap(gam)                  /* make gamma correction map */
double  gam;
{
        register int    i;
        
        for (i = 0; i < 256; i++)
                clrmap[RED][i] =
                clrmap[GRN][i] =
                clrmap[BLU][i] = 256.0 * pow((i+0.5)/256.0, 1.0/gam);
}


void
set_cmap(rmap, gmap, bmap)      /* set custom color correction map */
BYTE    *rmap, *gmap, *bmap;
{
        bcopy((char *)rmap, (char *)clrmap[RED], 256);
        bcopy((char *)gmap, (char *)clrmap[GRN], 256);
        bcopy((char *)bmap, (char *)clrmap[BLU], 256);
}


void
map_color(rgb, col)             /* map a color to a byte triplet */
BYTE    rgb[3];
COLOR   col;
{
        rgb[RED] = map_col(col,RED);
        rgb[GRN] = map_col(col,GRN);
        rgb[BLU] = map_col(col,BLU);
}


static void
cut(tree, level, box, c0, c1)           /* partition color space */
register CNODE  *tree;
int     level;
register int    box[3][2];
int     c0, c1;
{
        int     kb[3][2];
        
        if (c1-c0 <= 1) {               /* assign pixel */
                *tree = set_pval(c0);
                return;
        }
                                        /* split box */
        *tree = split(box);
        bcopy((char *)box, (char *)kb, sizeof(kb));
                                                /* do left (lesser) branch */
        kb[prim(*tree)][1] = part(*tree);
        cut(tree+(1<<level), level+1, kb, c0, (c0+c1)>>1);
                                                /* do right branch */
        kb[prim(*tree)][0] = part(*tree);
        kb[prim(*tree)][1] = box[prim(*tree)][1];
        cut(tree+(1<<(level+1)), level+1, kb, (c0+c1)>>1, c1);
}


static int
split(box)                              /* find median cut for box */
register int    box[3][2];
{
#define c0      r
        register int    r, g, b;
        int     pri;
        long    t[HMAX], med;
                                        /* find dominant axis */
        pri = RED;
        if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
                pri = GRN;
        if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
                pri = BLU;
                                        /* sum histogram over box */
        med = 0;
        switch (pri) {
        case RED:
                for (r = box[RED][0]; r < box[RED][1]; r++) {
                        t[r] = 0;
                        for (g = box[GRN][0]; g < box[GRN][1]; g++)
                                for (b = box[BLU][0]; b < box[BLU][1]; b++)
                                        t[r] += histo[r][g][b];
                        med += t[r];
                }
                break;
        case GRN:
                for (g = box[GRN][0]; g < box[GRN][1]; g++) {
                        t[g] = 0;
                        for (b = box[BLU][0]; b < box[BLU][1]; b++)
                                for (r = box[RED][0]; r < box[RED][1]; r++)
                                        t[g] += histo[r][g][b];
                        med += t[g];
                }
                break;
        case BLU:
                for (b = box[BLU][0]; b < box[BLU][1]; b++) {
                        t[b] = 0;
                        for (r = box[RED][0]; r < box[RED][1]; r++)
                                for (g = box[GRN][0]; g < box[GRN][1]; g++)
                                        t[b] += histo[r][g][b];
                        med += t[b];
                }
                break;
        }
        if (med < MINSAMP)              /* if too sparse, split at midpoint */
                return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
                                        /* find median position */
        med >>= 1;
        for (c0 = box[pri][0]; med > 0; c0++)
                med -= t[c0];
        if (c0 > (box[pri][0]+box[pri][1])>>1)  /* if past the midpoint */
                c0--;                           /* part left of median */
        return(set_branch(pri,c0));
#undef c0
}
Revision:	2.5
Committed:	Sat Feb 22 02:07:28 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.4:	+69 -10 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.5	static const char RCSid[] = "$Id$";
3	greg	1.1	#endif
4			/*
5			* colortab.c - allocate and control dynamic color table.
6			*
7			* We start off with a uniform partition of color space.
8			* As pixels are sent to the frame buffer, a histogram is built.
9			* When a new color table is requested, the histogram is used
10			* to make a pseudo-optimal partition, after which the
11			* histogram is cleared. This algorithm
12			* performs only as well as the next drawing's color
13			* distribution is correlated to the last.
14	greg	2.5	*
15			* External symbols declared in drvier.h
16			*/
17
18			/* ====================================================================
19			* The Radiance Software License, Version 1.0
20			*
21			* Copyright (c) 1990 - 2002 The Regents of the University of California,
22			* through Lawrence Berkeley National Laboratory. All rights reserved.
23			*
24			* Redistribution and use in source and binary forms, with or without
25			* modification, are permitted provided that the following conditions
26			* are met:
27			*
28			* 1. Redistributions of source code must retain the above copyright
29			* notice, this list of conditions and the following disclaimer.
30			*
31			* 2. Redistributions in binary form must reproduce the above copyright
32			* notice, this list of conditions and the following disclaimer in
33			* the documentation and/or other materials provided with the
34			* distribution.
35			*
36			* 3. The end-user documentation included with the redistribution,
37			* if any, must include the following acknowledgment:
38			* "This product includes Radiance software
39			* (http://radsite.lbl.gov/)
40			* developed by the Lawrence Berkeley National Laboratory
41			* (http://www.lbl.gov/)."
42			* Alternately, this acknowledgment may appear in the software itself,
43			* if and wherever such third-party acknowledgments normally appear.
44			*
45			* 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
46			* and "The Regents of the University of California" must
47			* not be used to endorse or promote products derived from this
48			* software without prior written permission. For written
49			* permission, please contact [email protected].
50			*
51			* 5. Products derived from this software may not be called "Radiance",
52			* nor may "Radiance" appear in their name, without prior written
53			* permission of Lawrence Berkeley National Laboratory.
54			*
55			* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
56			* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
57			* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
58			* DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
59			* ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
60			* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
61			* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
62			* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
63			* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
64			* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
65			* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
66			* SUCH DAMAGE.
67			* ====================================================================
68			*
69			* This software consists of voluntary contributions made by many
70			* individuals on behalf of Lawrence Berkeley National Laboratory. For more
71			* information on Lawrence Berkeley National Laboratory, please see
72			* <http://www.lbl.gov/>.
73	greg	1.1	*/
74
75	greg	1.15	#include "standard.h"
76
77	greg	1.1	#include "color.h"
78			/* histogram resolution */
79	greg	1.3	#define NRED 24
80	greg	1.1	#define NGRN 32
81	greg	1.3	#define NBLU 16
82	greg	1.1	#define HMAX NGRN
83			/* minimum box count for adaptive partition */
84			#define MINSAMP 7
85	greg	1.3	/* maximum distance^2 before color reassign */
86	greg	1.5	#define MAXDST2 12
87	greg	1.1	/* map a color */
88	greg	1.4	#define map_col(c,p) clrmap[p][ colval(c,p)<1. ? \
89			(int)(colval(c,p)*256.) : 255 ]
90	greg	1.1	/* color partition tree */
91			#define CNODE short
92	greg	2.2	#define set_branch(p,c) ((c)<<2\|(p))
93	greg	1.1	#define set_pval(pv) ((pv)<<2\|3)
94	greg	1.7	#define is_branch(cn) (((cn)&3)!=3)
95	greg	1.1	#define is_pval(cn) (((cn)&3)==3)
96			#define part(cn) ((cn)>>2)
97			#define prim(cn) ((cn)&3)
98			#define pval(cn) ((cn)>>2)
99			/* our color table */
100	greg	1.6	static struct tabent {
101	greg	1.3	long sum[3]; /* sum of colors using this entry */
102	greg	1.15	int n; /* number of colors */
103			BYTE ent[3]; /* current table value */
104	greg	1.8	} *clrtab = NULL;
105			/* color cube partition */
106			static CNODE *ctree = NULL;
107	greg	1.1	/* our color correction map */
108	greg	1.4	static BYTE clrmap[3][256];
109	greg	1.1	/* histogram of colors used */
110	greg	1.9	static unsigned short histo[NRED][NGRN][NBLU];
111	greg	1.8	/* initial color cube boundary */
112	greg	2.3	static int CLRCUBE[3][2] = {{0,NRED},{0,NGRN},{0,NBLU}};
113
114	greg	2.5	static int split();
115			static void cut();
116	greg	1.1
117
118	greg	1.3	int
119	greg	1.4	new_ctab(ncolors) /* start new color table with max ncolors */
120	greg	1.1	int ncolors;
121			{
122	greg	1.8	int treesize;
123
124	greg	1.6	if (ncolors < 1)
125	greg	1.3	return(0);
126	greg	1.8	/* free old tables */
127			if (clrtab != NULL)
128	greg	2.5	free((void *)clrtab);
129	greg	1.8	if (ctree != NULL)
130	greg	2.5	free((void *)ctree);
131	greg	1.8	/* get new tables */
132			for (treesize = 1; treesize < ncolors; treesize <<= 1)
133			;
134			treesize <<= 1;
135			clrtab = (struct tabent *)calloc(ncolors, sizeof(struct tabent));
136			ctree = (CNODE )malloc(treesizesizeof(CNODE));
137			if (clrtab == NULL \|\| ctree == NULL)
138			return(0);
139	greg	1.3	/* partition color space */
140	greg	1.7	cut(ctree, 0, CLRCUBE, 0, ncolors);
141	greg	1.1	/* clear histogram */
142	greg	1.11	bzero((char *)histo, sizeof(histo));
143	greg	1.3	/* return number of colors used */
144			return(ncolors);
145	greg	1.1	}
146
147
148			int
149	greg	1.4	get_pixel(col, set_pixel) /* get pixel for color */
150	greg	1.1	COLOR col;
151	greg	2.5	void (*set_pixel)();
152	greg	1.1	{
153	greg	1.3	int r, g, b;
154	greg	1.1	int cv[3];
155	greg	1.10	register CNODE *tp;
156	greg	1.1	register int h;
157	greg	1.3	/* map color */
158			r = map_col(col,RED);
159			g = map_col(col,GRN);
160			b = map_col(col,BLU);
161			/* reduce resolution */
162			cv[RED] = (r*NRED)>>8;
163			cv[GRN] = (g*NGRN)>>8;
164			cv[BLU] = (b*NBLU)>>8;
165			/* add to histogram */
166	greg	1.1	histo[cv[RED]][cv[GRN]][cv[BLU]]++;
167	greg	1.3	/* find pixel in tree */
168	greg	1.10	for (tp = ctree, h = 0; is_branch(*tp); h++)
169			if (cv[prim(tp)] < part(tp))
170			tp += 1<<h; /* left branch */
171	greg	1.1	else
172	greg	1.15	tp += 1<<(h+1); /* right branch */
173	greg	1.10	h = pval(*tp);
174	greg	1.3	/* add to color table */
175	greg	1.10	clrtab[h].sum[RED] += r;
176			clrtab[h].sum[GRN] += g;
177			clrtab[h].sum[BLU] += b;
178			clrtab[h].n++;
179	greg	1.3	/* recompute average */
180	greg	1.10	r = clrtab[h].sum[RED] / clrtab[h].n;
181			g = clrtab[h].sum[GRN] / clrtab[h].n;
182			b = clrtab[h].sum[BLU] / clrtab[h].n;
183	greg	1.3	/* check for movement */
184	greg	1.10	if (clrtab[h].n == 1 \|\|
185			(r-clrtab[h].ent[RED])*(r-clrtab[h].ent[RED]) +
186			(g-clrtab[h].ent[GRN])*(g-clrtab[h].ent[GRN]) +
187			(b-clrtab[h].ent[BLU])*(b-clrtab[h].ent[BLU]) > MAXDST2) {
188			clrtab[h].ent[RED] = r;
189	greg	2.2	clrtab[h].ent[GRN] = g; /* reassign pixel */
190	greg	1.10	clrtab[h].ent[BLU] = b;
191	greg	1.12	#ifdef DEBUG
192			sprintf(errmsg, "pixel %d = (%d,%d,%d) (%d refs)\n",
193	greg	1.10	h, r, g, b, clrtab[h].n);
194	greg	1.12	eputs(errmsg);
195	greg	1.1	#endif
196	greg	1.3	(*set_pixel)(h, r, g, b);
197			}
198			return(h); /* return pixel value */
199	greg	1.1	}
200
201
202	greg	2.5	void
203	greg	1.4	make_gmap(gam) /* make gamma correction map */
204	greg	2.2	double gam;
205	greg	1.1	{
206			register int i;
207
208	greg	1.4	for (i = 0; i < 256; i++)
209			clrmap[RED][i] =
210			clrmap[GRN][i] =
211	greg	1.15	clrmap[BLU][i] = 256.0 * pow((i+0.5)/256.0, 1.0/gam);
212	greg	1.4	}
213
214
215	greg	2.5	void
216	greg	1.4	set_cmap(rmap, gmap, bmap) /* set custom color correction map */
217			BYTE rmap, gmap, *bmap;
218			{
219	greg	1.11	bcopy((char )rmap, (char )clrmap[RED], 256);
220			bcopy((char )gmap, (char )clrmap[GRN], 256);
221			bcopy((char )bmap, (char )clrmap[BLU], 256);
222	greg	1.1	}
223
224
225	greg	2.5	void
226	greg	1.13	map_color(rgb, col) /* map a color to a byte triplet */
227			BYTE rgb[3];
228			COLOR col;
229			{
230			rgb[RED] = map_col(col,RED);
231			rgb[GRN] = map_col(col,GRN);
232			rgb[BLU] = map_col(col,BLU);
233			}
234
235
236	greg	2.5	static void
237	greg	1.7	cut(tree, level, box, c0, c1) /* partition color space */
238	greg	1.1	register CNODE *tree;
239	greg	1.7	int level;
240	greg	1.1	register int box[3][2];
241			int c0, c1;
242			{
243			int kb[3][2];
244
245	greg	1.3	if (c1-c0 <= 1) { /* assign pixel */
246	greg	1.1	*tree = set_pval(c0);
247			return;
248			}
249			/* split box */
250			*tree = split(box);
251	greg	1.11	bcopy((char )box, (char )kb, sizeof(kb));
252	greg	1.3	/* do left (lesser) branch */
253	greg	1.1	kb[prim(tree)][1] = part(tree);
254	greg	1.7	cut(tree+(1<<level), level+1, kb, c0, (c0+c1)>>1);
255	greg	1.3	/* do right branch */
256	greg	1.1	kb[prim(tree)][0] = part(tree);
257			kb[prim(tree)][1] = box[prim(tree)][1];
258	greg	1.7	cut(tree+(1<<(level+1)), level+1, kb, (c0+c1)>>1, c1);
259	greg	1.1	}
260
261
262			static int
263			split(box) /* find median cut for box */
264			register int box[3][2];
265			{
266			#define c0 r
267			register int r, g, b;
268			int pri;
269	greg	2.4	long t[HMAX], med;
270	greg	1.1	/* find dominant axis */
271			pri = RED;
272			if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
273			pri = GRN;
274			if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
275			pri = BLU;
276			/* sum histogram over box */
277			med = 0;
278			switch (pri) {
279			case RED:
280			for (r = box[RED][0]; r < box[RED][1]; r++) {
281			t[r] = 0;
282			for (g = box[GRN][0]; g < box[GRN][1]; g++)
283			for (b = box[BLU][0]; b < box[BLU][1]; b++)
284			t[r] += histo[r][g][b];
285			med += t[r];
286			}
287			break;
288			case GRN:
289			for (g = box[GRN][0]; g < box[GRN][1]; g++) {
290			t[g] = 0;
291			for (b = box[BLU][0]; b < box[BLU][1]; b++)
292			for (r = box[RED][0]; r < box[RED][1]; r++)
293			t[g] += histo[r][g][b];
294			med += t[g];
295			}
296			break;
297			case BLU:
298			for (b = box[BLU][0]; b < box[BLU][1]; b++) {
299			t[b] = 0;
300			for (r = box[RED][0]; r < box[RED][1]; r++)
301			for (g = box[GRN][0]; g < box[GRN][1]; g++)
302			t[b] += histo[r][g][b];
303			med += t[b];
304			}
305			break;
306			}
307			if (med < MINSAMP) /* if too sparse, split at midpoint */
308			return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
309			/* find median position */
310			med >>= 1;
311			for (c0 = box[pri][0]; med > 0; c0++)
312			med -= t[c0];
313			if (c0 > (box[pri][0]+box[pri][1])>>1) /* if past the midpoint */
314			c0--; /* part left of median */
315			return(set_branch(pri,c0));
316			#undef c0
317			}