1/* $Id: tiffmedian.c 276 2010-06-30 12:18:30Z nijtmans $ */
2
3/*
4 * Apply median cut on an image.
5 *
6 * tiffmedian [-c n] [-f] input output
7 *     -C n		- set colortable size.  Default is 256.
8 *     -f		- use Floyd-Steinberg dithering.
9 *     -c lzw		- compress output with LZW
10 *     -c none		- use no compression on output
11 *     -c packbits	- use packbits compression on output
12 *     -r n		- create output with n rows/strip of data
13 * (by default the compression scheme and rows/strip are taken
14 *  from the input file)
15 *
16 * Notes:
17 *
18 * [1] Floyd-Steinberg dither:
19 *  I should point out that the actual fractions we used were, assuming
20 *  you are at X, moving left to right:
21 *
22 *		    X     7/16
23 *	     3/16   5/16  1/16
24 *
25 *  Note that the error goes to four neighbors, not three.  I think this
26 *  will probably do better (at least for black and white) than the
27 *  3/8-3/8-1/4 distribution, at the cost of greater processing.  I have
28 *  seen the 3/8-3/8-1/4 distribution described as "our" algorithm before,
29 *  but I have no idea who the credit really belongs to.
30
31 *  Also, I should add that if you do zig-zag scanning (see my immediately
32 *  previous message), it is sufficient (but not quite as good) to send
33 *  half the error one pixel ahead (e.g. to the right on lines you scan
34 *  left to right), and half one pixel straight down.  Again, this is for
35 *  black and white;  I've not tried it with color.
36 *  --
37 *					    Lou Steinberg
38 *
39 * [2] Color Image Quantization for Frame Buffer Display, Paul Heckbert,
40 *	Siggraph '82 proceedings, pp. 297-307
41 */
42
43#include "tif_config.h"
44
45#include <stdio.h>
46#include <stdlib.h>
47#include <string.h>
48
49#ifdef HAVE_UNISTD_H
50# include <unistd.h>
51#endif
52
53#include "tiffio.h"
54
55#define	MAX_CMAP_SIZE	256
56
57#define	streq(a,b)	(strcmp(a,b) == 0)
58#define	strneq(a,b,n)	(strncmp(a,b,n) == 0)
59
60#define	COLOR_DEPTH	8
61#define	MAX_COLOR	256
62
63#define	B_DEPTH		5		/* # bits/pixel to use */
64#define	B_LEN		(1L<<B_DEPTH)
65
66#define	C_DEPTH		2
67#define	C_LEN		(1L<<C_DEPTH)	/* # cells/color to use */
68
69#define	COLOR_SHIFT	(COLOR_DEPTH-B_DEPTH)
70
71typedef	struct colorbox {
72	struct	colorbox *next, *prev;
73	int	rmin, rmax;
74	int	gmin, gmax;
75	int	bmin, bmax;
76	uint32	total;
77} Colorbox;
78
79typedef struct {
80	int	num_ents;
81	int	entries[MAX_CMAP_SIZE][2];
82} C_cell;
83
84uint16	rm[MAX_CMAP_SIZE], gm[MAX_CMAP_SIZE], bm[MAX_CMAP_SIZE];
85int	num_colors;
86uint32	histogram[B_LEN][B_LEN][B_LEN];
87Colorbox *freeboxes;
88Colorbox *usedboxes;
89C_cell	**ColorCells;
90TIFF	*in, *out;
91uint32	rowsperstrip = (uint32) -1;
92uint16	compression = (uint16) -1;
93uint16	bitspersample = 1;
94uint16	samplesperpixel;
95uint32	imagewidth;
96uint32	imagelength;
97uint16	predictor = 0;
98
99static	void get_histogram(TIFF*, Colorbox*);
100static	void splitbox(Colorbox*);
101static	void shrinkbox(Colorbox*);
102static	void map_colortable(void);
103static	void quant(TIFF*, TIFF*);
104static	void quant_fsdither(TIFF*, TIFF*);
105static	Colorbox* largest_box(void);
106
107static	void usage(void);
108static	int processCompressOptions(char*);
109
110#define	CopyField(tag, v) \
111	if (TIFFGetField(in, tag, &v)) TIFFSetField(out, tag, v)
112
113int
114main(int argc, char* argv[])
115{
116	int i, dither = 0;
117	uint16 shortv, config, photometric;
118	Colorbox *box_list, *ptr;
119	float floatv;
120	uint32 longv;
121	int c;
122	extern int optind;
123	extern char* optarg;
124
125	num_colors = MAX_CMAP_SIZE;
126	while ((c = getopt(argc, argv, "c:C:r:f")) != -1)
127		switch (c) {
128		case 'c':		/* compression scheme */
129			if (!processCompressOptions(optarg))
130				usage();
131			break;
132		case 'C':		/* set colormap size */
133			num_colors = atoi(optarg);
134			if (num_colors > MAX_CMAP_SIZE) {
135				fprintf(stderr,
136				   "-c: colormap too big, max %d\n",
137				   MAX_CMAP_SIZE);
138				usage();
139			}
140			break;
141		case 'f':		/* dither */
142			dither = 1;
143			break;
144		case 'r':		/* rows/strip */
145			rowsperstrip = atoi(optarg);
146			break;
147		case '?':
148			usage();
149			/*NOTREACHED*/
150		}
151	if (argc - optind != 2)
152		usage();
153	in = TIFFOpen(argv[optind], "r");
154	if (in == NULL)
155		return (-1);
156	TIFFGetField(in, TIFFTAG_IMAGEWIDTH, &imagewidth);
157	TIFFGetField(in, TIFFTAG_IMAGELENGTH, &imagelength);
158	TIFFGetField(in, TIFFTAG_BITSPERSAMPLE, &bitspersample);
159	TIFFGetField(in, TIFFTAG_SAMPLESPERPIXEL, &samplesperpixel);
160	if (bitspersample != 8 && bitspersample != 16) {
161		fprintf(stderr, "%s: Image must have at least 8-bits/sample\n",
162		    argv[optind]);
163		return (-3);
164	}
165	if (!TIFFGetField(in, TIFFTAG_PHOTOMETRIC, &photometric) ||
166	    photometric != PHOTOMETRIC_RGB || samplesperpixel < 3) {
167		fprintf(stderr, "%s: Image must have RGB data\n", argv[optind]);
168		return (-4);
169	}
170	TIFFGetField(in, TIFFTAG_PLANARCONFIG, &config);
171	if (config != PLANARCONFIG_CONTIG) {
172		fprintf(stderr, "%s: Can only handle contiguous data packing\n",
173		    argv[optind]);
174		return (-5);
175	}
176
177	/*
178	 * STEP 1:  create empty boxes
179	 */
180	usedboxes = NULL;
181	box_list = freeboxes = (Colorbox *)_TIFFmalloc(num_colors*sizeof (Colorbox));
182	freeboxes[0].next = &freeboxes[1];
183	freeboxes[0].prev = NULL;
184	for (i = 1; i < num_colors-1; ++i) {
185		freeboxes[i].next = &freeboxes[i+1];
186		freeboxes[i].prev = &freeboxes[i-1];
187	}
188	freeboxes[num_colors-1].next = NULL;
189	freeboxes[num_colors-1].prev = &freeboxes[num_colors-2];
190
191	/*
192	 * STEP 2: get histogram, initialize first box
193	 */
194	ptr = freeboxes;
195	freeboxes = ptr->next;
196	if (freeboxes)
197		freeboxes->prev = NULL;
198	ptr->next = usedboxes;
199	usedboxes = ptr;
200	if (ptr->next)
201		ptr->next->prev = ptr;
202	get_histogram(in, ptr);
203
204	/*
205	 * STEP 3: continually subdivide boxes until no more free
206	 * boxes remain or until all colors assigned.
207	 */
208	while (freeboxes != NULL) {
209		ptr = largest_box();
210		if (ptr != NULL)
211			splitbox(ptr);
212		else
213			freeboxes = NULL;
214	}
215
216	/*
217	 * STEP 4: assign colors to all boxes
218	 */
219	for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next) {
220		rm[i] = ((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2;
221		gm[i] = ((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2;
222		bm[i] = ((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2;
223	}
224
225	/* We're done with the boxes now */
226	_TIFFfree(box_list);
227	freeboxes = usedboxes = NULL;
228
229	/*
230	 * STEP 5: scan histogram and map all values to closest color
231	 */
232	/* 5a: create cell list as described in Heckbert[2] */
233	ColorCells = (C_cell **)_TIFFmalloc(C_LEN*C_LEN*C_LEN*sizeof (C_cell*));
234	_TIFFmemset(ColorCells, 0, C_LEN*C_LEN*C_LEN*sizeof (C_cell*));
235	/* 5b: create mapping from truncated pixel space to color
236	   table entries */
237	map_colortable();
238
239	/*
240	 * STEP 6: scan image, match input values to table entries
241	 */
242	out = TIFFOpen(argv[optind+1], "w");
243	if (out == NULL)
244		return (-2);
245
246	CopyField(TIFFTAG_SUBFILETYPE, longv);
247	CopyField(TIFFTAG_IMAGEWIDTH, longv);
248	TIFFSetField(out, TIFFTAG_BITSPERSAMPLE, (short)COLOR_DEPTH);
249	if (compression != (uint16)-1) {
250		TIFFSetField(out, TIFFTAG_COMPRESSION, compression);
251		switch (compression) {
252		case COMPRESSION_LZW:
253		case COMPRESSION_DEFLATE:
254			if (predictor != 0)
255				TIFFSetField(out, TIFFTAG_PREDICTOR, predictor);
256			break;
257		}
258	} else
259		CopyField(TIFFTAG_COMPRESSION, compression);
260	TIFFSetField(out, TIFFTAG_PHOTOMETRIC, (short)PHOTOMETRIC_PALETTE);
261	CopyField(TIFFTAG_ORIENTATION, shortv);
262	TIFFSetField(out, TIFFTAG_SAMPLESPERPIXEL, (short)1);
263	CopyField(TIFFTAG_PLANARCONFIG, shortv);
264	TIFFSetField(out, TIFFTAG_ROWSPERSTRIP,
265	    TIFFDefaultStripSize(out, rowsperstrip));
266	CopyField(TIFFTAG_MINSAMPLEVALUE, shortv);
267	CopyField(TIFFTAG_MAXSAMPLEVALUE, shortv);
268	CopyField(TIFFTAG_RESOLUTIONUNIT, shortv);
269	CopyField(TIFFTAG_XRESOLUTION, floatv);
270	CopyField(TIFFTAG_YRESOLUTION, floatv);
271	CopyField(TIFFTAG_XPOSITION, floatv);
272	CopyField(TIFFTAG_YPOSITION, floatv);
273
274	if (dither)
275		quant_fsdither(in, out);
276	else
277		quant(in, out);
278	/*
279	 * Scale colormap to TIFF-required 16-bit values.
280	 */
281#define	SCALE(x)	(((x)*((1L<<16)-1))/255)
282	for (i = 0; i < MAX_CMAP_SIZE; ++i) {
283		rm[i] = SCALE(rm[i]);
284		gm[i] = SCALE(gm[i]);
285		bm[i] = SCALE(bm[i]);
286	}
287	TIFFSetField(out, TIFFTAG_COLORMAP, rm, gm, bm);
288	(void) TIFFClose(out);
289	return (0);
290}
291
292static int
293processCompressOptions(char* opt)
294{
295	if (streq(opt, "none"))
296		compression = COMPRESSION_NONE;
297	else if (streq(opt, "packbits"))
298		compression = COMPRESSION_PACKBITS;
299	else if (strneq(opt, "lzw", 3)) {
300		char* cp = strchr(opt, ':');
301		if (cp)
302			predictor = atoi(cp+1);
303		compression = COMPRESSION_LZW;
304	} else if (strneq(opt, "zip", 3)) {
305		char* cp = strchr(opt, ':');
306		if (cp)
307			predictor = atoi(cp+1);
308		compression = COMPRESSION_DEFLATE;
309	} else
310		return (0);
311	return (1);
312}
313
314char* stuff[] = {
315"usage: tiffmedian [options] input.tif output.tif",
316"where options are:",
317" -r #		make each strip have no more than # rows",
318" -C #		create a colormap with # entries",
319" -f		use Floyd-Steinberg dithering",
320" -c lzw[:opts]	compress output with Lempel-Ziv & Welch encoding",
321" -c zip[:opts]	compress output with deflate encoding",
322" -c packbits	compress output with packbits encoding",
323" -c none	use no compression algorithm on output",
324"",
325"LZW and deflate options:",
326" #		set predictor value",
327"For example, -c lzw:2 to get LZW-encoded data with horizontal differencing",
328NULL
329};
330
331static void
332usage(void)
333{
334	char buf[BUFSIZ];
335	int i;
336
337	setbuf(stderr, buf);
338        fprintf(stderr, "%s\n\n", TIFFGetVersion());
339	for (i = 0; stuff[i] != NULL; i++)
340		fprintf(stderr, "%s\n", stuff[i]);
341	exit(-1);
342}
343
344static void
345get_histogram(TIFF* in, Colorbox* box)
346{
347	register unsigned char *inptr;
348	register int red, green, blue;
349	register uint32 j, i;
350	unsigned char *inputline;
351
352	inputline = (unsigned char *)_TIFFmalloc(TIFFScanlineSize(in));
353	if (inputline == NULL) {
354		fprintf(stderr, "No space for scanline buffer\n");
355		exit(-1);
356	}
357	box->rmin = box->gmin = box->bmin = 999;
358	box->rmax = box->gmax = box->bmax = -1;
359	box->total = imagewidth * imagelength;
360
361	{ register uint32 *ptr = &histogram[0][0][0];
362	  for (i = B_LEN*B_LEN*B_LEN; i-- > 0;)
363		*ptr++ = 0;
364	}
365	for (i = 0; i < imagelength; i++) {
366		if (TIFFReadScanline(in, inputline, i, 0) <= 0)
367			break;
368		inptr = inputline;
369		for (j = imagewidth; j-- > 0;) {
370			red = *inptr++ >> COLOR_SHIFT;
371			green = *inptr++ >> COLOR_SHIFT;
372			blue = *inptr++ >> COLOR_SHIFT;
373			if (red < box->rmin)
374				box->rmin = red;
375		        if (red > box->rmax)
376				box->rmax = red;
377		        if (green < box->gmin)
378				box->gmin = green;
379		        if (green > box->gmax)
380				box->gmax = green;
381		        if (blue < box->bmin)
382				box->bmin = blue;
383		        if (blue > box->bmax)
384				box->bmax = blue;
385		        histogram[red][green][blue]++;
386		}
387	}
388	_TIFFfree(inputline);
389}
390
391static Colorbox *
392largest_box(void)
393{
394	register Colorbox *p, *b;
395	register uint32 size;
396
397	b = NULL;
398	size = 0;
399	for (p = usedboxes; p != NULL; p = p->next)
400		if ((p->rmax > p->rmin || p->gmax > p->gmin ||
401		    p->bmax > p->bmin) &&  p->total > size)
402		        size = (b = p)->total;
403	return (b);
404}
405
406static void
407splitbox(Colorbox* ptr)
408{
409	uint32		hist2[B_LEN];
410	int		first=0, last=0;
411	register Colorbox	*new;
412	register uint32	*iptr, *histp;
413	register int	i, j;
414	register int	ir,ig,ib;
415	register uint32 sum, sum1, sum2;
416	enum { RED, GREEN, BLUE } axis;
417
418	/*
419	 * See which axis is the largest, do a histogram along that
420	 * axis.  Split at median point.  Contract both new boxes to
421	 * fit points and return
422	 */
423	i = ptr->rmax - ptr->rmin;
424	if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
425		axis = RED;
426	else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
427		axis = GREEN;
428	else
429		axis = BLUE;
430	/* get histogram along longest axis */
431	switch (axis) {
432	case RED:
433		histp = &hist2[ptr->rmin];
434	        for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
435			*histp = 0;
436			for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
437				iptr = &histogram[ir][ig][ptr->bmin];
438				for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
439					*histp += *iptr++;
440			}
441			histp++;
442	        }
443	        first = ptr->rmin;
444		last = ptr->rmax;
445	        break;
446	case GREEN:
447	        histp = &hist2[ptr->gmin];
448	        for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
449			*histp = 0;
450			for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
451				iptr = &histogram[ir][ig][ptr->bmin];
452				for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
453					*histp += *iptr++;
454			}
455			histp++;
456	        }
457	        first = ptr->gmin;
458		last = ptr->gmax;
459	        break;
460	case BLUE:
461	        histp = &hist2[ptr->bmin];
462	        for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) {
463			*histp = 0;
464			for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
465				iptr = &histogram[ir][ptr->gmin][ib];
466				for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
467					*histp += *iptr;
468					iptr += B_LEN;
469				}
470			}
471			histp++;
472	        }
473	        first = ptr->bmin;
474		last = ptr->bmax;
475	        break;
476	}
477	/* find median point */
478	sum2 = ptr->total / 2;
479	histp = &hist2[first];
480	sum = 0;
481	for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
482		;
483	if (i == first)
484		i++;
485
486	/* Create new box, re-allocate points */
487	new = freeboxes;
488	freeboxes = new->next;
489	if (freeboxes)
490		freeboxes->prev = NULL;
491	if (usedboxes)
492		usedboxes->prev = new;
493	new->next = usedboxes;
494	usedboxes = new;
495
496	histp = &hist2[first];
497	for (sum1 = 0, j = first; j < i; j++)
498		sum1 += *histp++;
499	for (sum2 = 0, j = i; j <= last; j++)
500	    sum2 += *histp++;
501	new->total = sum1;
502	ptr->total = sum2;
503
504	new->rmin = ptr->rmin;
505	new->rmax = ptr->rmax;
506	new->gmin = ptr->gmin;
507	new->gmax = ptr->gmax;
508	new->bmin = ptr->bmin;
509	new->bmax = ptr->bmax;
510	switch (axis) {
511	case RED:
512		new->rmax = i-1;
513	        ptr->rmin = i;
514	        break;
515	case GREEN:
516	        new->gmax = i-1;
517	        ptr->gmin = i;
518	        break;
519	case BLUE:
520	        new->bmax = i-1;
521	        ptr->bmin = i;
522	        break;
523	}
524	shrinkbox(new);
525	shrinkbox(ptr);
526}
527
528static void
529shrinkbox(Colorbox* box)
530{
531	register uint32 *histp;
532	register int	ir, ig, ib;
533
534	if (box->rmax > box->rmin) {
535		for (ir = box->rmin; ir <= box->rmax; ++ir)
536			for (ig = box->gmin; ig <= box->gmax; ++ig) {
537				histp = &histogram[ir][ig][box->bmin];
538			        for (ib = box->bmin; ib <= box->bmax; ++ib)
539					if (*histp++ != 0) {
540						box->rmin = ir;
541						goto have_rmin;
542					}
543			}
544	have_rmin:
545		if (box->rmax > box->rmin)
546			for (ir = box->rmax; ir >= box->rmin; --ir)
547				for (ig = box->gmin; ig <= box->gmax; ++ig) {
548					histp = &histogram[ir][ig][box->bmin];
549					ib = box->bmin;
550					for (; ib <= box->bmax; ++ib)
551						if (*histp++ != 0) {
552							box->rmax = ir;
553							goto have_rmax;
554						}
555			        }
556	}
557have_rmax:
558	if (box->gmax > box->gmin) {
559		for (ig = box->gmin; ig <= box->gmax; ++ig)
560			for (ir = box->rmin; ir <= box->rmax; ++ir) {
561				histp = &histogram[ir][ig][box->bmin];
562			        for (ib = box->bmin; ib <= box->bmax; ++ib)
563				if (*histp++ != 0) {
564					box->gmin = ig;
565					goto have_gmin;
566				}
567			}
568	have_gmin:
569		if (box->gmax > box->gmin)
570			for (ig = box->gmax; ig >= box->gmin; --ig)
571				for (ir = box->rmin; ir <= box->rmax; ++ir) {
572					histp = &histogram[ir][ig][box->bmin];
573					ib = box->bmin;
574					for (; ib <= box->bmax; ++ib)
575						if (*histp++ != 0) {
576							box->gmax = ig;
577							goto have_gmax;
578						}
579			        }
580	}
581have_gmax:
582	if (box->bmax > box->bmin) {
583		for (ib = box->bmin; ib <= box->bmax; ++ib)
584			for (ir = box->rmin; ir <= box->rmax; ++ir) {
585				histp = &histogram[ir][box->gmin][ib];
586			        for (ig = box->gmin; ig <= box->gmax; ++ig) {
587					if (*histp != 0) {
588						box->bmin = ib;
589						goto have_bmin;
590					}
591					histp += B_LEN;
592			        }
593		        }
594	have_bmin:
595		if (box->bmax > box->bmin)
596			for (ib = box->bmax; ib >= box->bmin; --ib)
597				for (ir = box->rmin; ir <= box->rmax; ++ir) {
598					histp = &histogram[ir][box->gmin][ib];
599					ig = box->gmin;
600					for (; ig <= box->gmax; ++ig) {
601						if (*histp != 0) {
602							box->bmax = ib;
603							goto have_bmax;
604						}
605						histp += B_LEN;
606					}
607			        }
608	}
609have_bmax:
610	;
611}
612
613static C_cell *
614create_colorcell(int red, int green, int blue)
615{
616	register int ir, ig, ib, i;
617	register C_cell *ptr;
618	int mindist, next_n;
619	register int tmp, dist, n;
620
621	ir = red >> (COLOR_DEPTH-C_DEPTH);
622	ig = green >> (COLOR_DEPTH-C_DEPTH);
623	ib = blue >> (COLOR_DEPTH-C_DEPTH);
624	ptr = (C_cell *)_TIFFmalloc(sizeof (C_cell));
625	*(ColorCells + ir*C_LEN*C_LEN + ig*C_LEN + ib) = ptr;
626	ptr->num_ents = 0;
627
628	/*
629	 * Step 1: find all colors inside this cell, while we're at
630	 *	   it, find distance of centermost point to furthest corner
631	 */
632	mindist = 99999999;
633	for (i = 0; i < num_colors; ++i) {
634		if (rm[i]>>(COLOR_DEPTH-C_DEPTH) != ir  ||
635		    gm[i]>>(COLOR_DEPTH-C_DEPTH) != ig  ||
636		    bm[i]>>(COLOR_DEPTH-C_DEPTH) != ib)
637			continue;
638		ptr->entries[ptr->num_ents][0] = i;
639		ptr->entries[ptr->num_ents][1] = 0;
640		++ptr->num_ents;
641	        tmp = rm[i] - red;
642	        if (tmp < (MAX_COLOR/C_LEN/2))
643			tmp = MAX_COLOR/C_LEN-1 - tmp;
644	        dist = tmp*tmp;
645	        tmp = gm[i] - green;
646	        if (tmp < (MAX_COLOR/C_LEN/2))
647			tmp = MAX_COLOR/C_LEN-1 - tmp;
648	        dist += tmp*tmp;
649	        tmp = bm[i] - blue;
650	        if (tmp < (MAX_COLOR/C_LEN/2))
651			tmp = MAX_COLOR/C_LEN-1 - tmp;
652	        dist += tmp*tmp;
653	        if (dist < mindist)
654			mindist = dist;
655	}
656
657	/*
658	 * Step 3: find all points within that distance to cell.
659	 */
660	for (i = 0; i < num_colors; ++i) {
661		if (rm[i] >> (COLOR_DEPTH-C_DEPTH) == ir  &&
662		    gm[i] >> (COLOR_DEPTH-C_DEPTH) == ig  &&
663		    bm[i] >> (COLOR_DEPTH-C_DEPTH) == ib)
664			continue;
665		dist = 0;
666	        if ((tmp = red - rm[i]) > 0 ||
667		    (tmp = rm[i] - (red + MAX_COLOR/C_LEN-1)) > 0 )
668			dist += tmp*tmp;
669	        if ((tmp = green - gm[i]) > 0 ||
670		    (tmp = gm[i] - (green + MAX_COLOR/C_LEN-1)) > 0 )
671			dist += tmp*tmp;
672	        if ((tmp = blue - bm[i]) > 0 ||
673		    (tmp = bm[i] - (blue + MAX_COLOR/C_LEN-1)) > 0 )
674			dist += tmp*tmp;
675	        if (dist < mindist) {
676			ptr->entries[ptr->num_ents][0] = i;
677			ptr->entries[ptr->num_ents][1] = dist;
678			++ptr->num_ents;
679	        }
680	}
681
682	/*
683	 * Sort color cells by distance, use cheap exchange sort
684	 */
685	for (n = ptr->num_ents - 1; n > 0; n = next_n) {
686		next_n = 0;
687		for (i = 0; i < n; ++i)
688			if (ptr->entries[i][1] > ptr->entries[i+1][1]) {
689				tmp = ptr->entries[i][0];
690				ptr->entries[i][0] = ptr->entries[i+1][0];
691				ptr->entries[i+1][0] = tmp;
692				tmp = ptr->entries[i][1];
693				ptr->entries[i][1] = ptr->entries[i+1][1];
694				ptr->entries[i+1][1] = tmp;
695				next_n = i;
696		        }
697	}
698	return (ptr);
699}
700
701static void
702map_colortable(void)
703{
704	register uint32 *histp = &histogram[0][0][0];
705	register C_cell *cell;
706	register int j, tmp, d2, dist;
707	int ir, ig, ib, i;
708
709	for (ir = 0; ir < B_LEN; ++ir)
710		for (ig = 0; ig < B_LEN; ++ig)
711			for (ib = 0; ib < B_LEN; ++ib, histp++) {
712				if (*histp == 0) {
713					*histp = -1;
714					continue;
715				}
716				cell = *(ColorCells +
717				    (((ir>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2) +
718				    ((ig>>(B_DEPTH-C_DEPTH)) << C_DEPTH) +
719				    (ib>>(B_DEPTH-C_DEPTH))));
720				if (cell == NULL )
721					cell = create_colorcell(
722					    ir << COLOR_SHIFT,
723					    ig << COLOR_SHIFT,
724					    ib << COLOR_SHIFT);
725				dist = 9999999;
726				for (i = 0; i < cell->num_ents &&
727				    dist > cell->entries[i][1]; ++i) {
728					j = cell->entries[i][0];
729					d2 = rm[j] - (ir << COLOR_SHIFT);
730					d2 *= d2;
731					tmp = gm[j] - (ig << COLOR_SHIFT);
732					d2 += tmp*tmp;
733					tmp = bm[j] - (ib << COLOR_SHIFT);
734					d2 += tmp*tmp;
735					if (d2 < dist) {
736						dist = d2;
737						*histp = j;
738					}
739				}
740			}
741}
742
743/*
744 * straight quantization.  Each pixel is mapped to the colors
745 * closest to it.  Color values are rounded to the nearest color
746 * table entry.
747 */
748static void
749quant(TIFF* in, TIFF* out)
750{
751	unsigned char	*outline, *inputline;
752	register unsigned char	*outptr, *inptr;
753	register uint32 i, j;
754	register int red, green, blue;
755
756	inputline = (unsigned char *)_TIFFmalloc(TIFFScanlineSize(in));
757	outline = (unsigned char *)_TIFFmalloc(imagewidth);
758	for (i = 0; i < imagelength; i++) {
759		if (TIFFReadScanline(in, inputline, i, 0) <= 0)
760			break;
761		inptr = inputline;
762		outptr = outline;
763		for (j = 0; j < imagewidth; j++) {
764			red = *inptr++ >> COLOR_SHIFT;
765			green = *inptr++ >> COLOR_SHIFT;
766			blue = *inptr++ >> COLOR_SHIFT;
767			*outptr++ = (unsigned char)histogram[red][green][blue];
768		}
769		if (TIFFWriteScanline(out, outline, i, 0) < 0)
770			break;
771	}
772	_TIFFfree(inputline);
773	_TIFFfree(outline);
774}
775
776#define	SWAP(type,a,b)	{ type p; p = a; a = b; b = p; }
777
778#define	GetInputLine(tif, row, bad)				\
779	if (TIFFReadScanline(tif, inputline, row, 0) <= 0)	\
780		bad;						\
781	inptr = inputline;					\
782	nextptr = nextline;					\
783	for (j = 0; j < imagewidth; ++j) {			\
784		*nextptr++ = *inptr++;				\
785		*nextptr++ = *inptr++;				\
786		*nextptr++ = *inptr++;				\
787	}
788#define	GetComponent(raw, cshift, c)				\
789	cshift = raw;						\
790	if (cshift < 0)						\
791		cshift = 0;					\
792	else if (cshift >= MAX_COLOR)				\
793		cshift = MAX_COLOR-1;				\
794	c = cshift;						\
795	cshift >>= COLOR_SHIFT;
796
797static void
798quant_fsdither(TIFF* in, TIFF* out)
799{
800	unsigned char *outline, *inputline, *inptr;
801	short *thisline, *nextline;
802	register unsigned char	*outptr;
803	register short *thisptr, *nextptr;
804	register uint32 i, j;
805	uint32 imax, jmax;
806	int lastline, lastpixel;
807
808	imax = imagelength - 1;
809	jmax = imagewidth - 1;
810	inputline = (unsigned char *)_TIFFmalloc(TIFFScanlineSize(in));
811	thisline = (short *)_TIFFmalloc(imagewidth * 3 * sizeof (short));
812	nextline = (short *)_TIFFmalloc(imagewidth * 3 * sizeof (short));
813	outline = (unsigned char *) _TIFFmalloc(TIFFScanlineSize(out));
814
815	GetInputLine(in, 0, goto bad);		/* get first line */
816	for (i = 1; i <= imagelength; ++i) {
817		SWAP(short *, thisline, nextline);
818		lastline = (i >= imax);
819		if (i <= imax)
820			GetInputLine(in, i, break);
821		thisptr = thisline;
822		nextptr = nextline;
823		outptr = outline;
824		for (j = 0; j < imagewidth; ++j) {
825			int red, green, blue;
826			register int oval, r2, g2, b2;
827
828			lastpixel = (j == jmax);
829			GetComponent(*thisptr++, r2, red);
830			GetComponent(*thisptr++, g2, green);
831			GetComponent(*thisptr++, b2, blue);
832			oval = histogram[r2][g2][b2];
833			if (oval == -1) {
834				int ci;
835				register int cj, tmp, d2, dist;
836				register C_cell	*cell;
837
838				cell = *(ColorCells +
839				    (((r2>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2) +
840				    ((g2>>(B_DEPTH-C_DEPTH)) << C_DEPTH ) +
841				    (b2>>(B_DEPTH-C_DEPTH))));
842				if (cell == NULL)
843					cell = create_colorcell(red,
844					    green, blue);
845				dist = 9999999;
846				for (ci = 0; ci < cell->num_ents && dist > cell->entries[ci][1]; ++ci) {
847					cj = cell->entries[ci][0];
848					d2 = (rm[cj] >> COLOR_SHIFT) - r2;
849					d2 *= d2;
850					tmp = (gm[cj] >> COLOR_SHIFT) - g2;
851					d2 += tmp*tmp;
852					tmp = (bm[cj] >> COLOR_SHIFT) - b2;
853					d2 += tmp*tmp;
854					if (d2 < dist) {
855						dist = d2;
856						oval = cj;
857					}
858				}
859				histogram[r2][g2][b2] = oval;
860			}
861			*outptr++ = oval;
862			red -= rm[oval];
863			green -= gm[oval];
864			blue -= bm[oval];
865			if (!lastpixel) {
866				thisptr[0] += blue * 7 / 16;
867				thisptr[1] += green * 7 / 16;
868				thisptr[2] += red * 7 / 16;
869			}
870			if (!lastline) {
871				if (j != 0) {
872					nextptr[-3] += blue * 3 / 16;
873					nextptr[-2] += green * 3 / 16;
874					nextptr[-1] += red * 3 / 16;
875				}
876				nextptr[0] += blue * 5 / 16;
877				nextptr[1] += green * 5 / 16;
878				nextptr[2] += red * 5 / 16;
879				if (!lastpixel) {
880					nextptr[3] += blue / 16;
881				        nextptr[4] += green / 16;
882				        nextptr[5] += red / 16;
883				}
884				nextptr += 3;
885			}
886		}
887		if (TIFFWriteScanline(out, outline, i-1, 0) < 0)
888			break;
889	}
890bad:
891	_TIFFfree(inputline);
892	_TIFFfree(thisline);
893	_TIFFfree(nextline);
894	_TIFFfree(outline);
895}
896/*
897 * Local Variables:
898 * mode: c
899 * c-basic-offset: 8
900 * fill-column: 78
901 * End:
902 */
903