1 /* reduce.c
2  *
3  * reduce an image's colormap usage to a set number of colors.  this also
4  * translates a true color image to a TLA-style image of `n' colors.
5  *
6  * this function cannot reduce to any number of colors larger than 32768.
7  *
8  * Based on the file lib/colorquant.c in the Utah Raster Toolkit version 3.0
9  *
10  * > Perform variance-based color quantization on a "full color" image.
11  * > Author:    Craig Kolb
12  * >    Department of Mathematics
13  * >    Yale University
14  * >    kolb@yale.edu
15  * > Date:      Tue Aug 22 1989
16  * > Copyright (C) 1989 Craig E. Kolb
17  *
18  * Adaption for xli, faster colour dithering
19  * and other improvements by Graeme Gill.
20  *
21  */
22 
23 #include "xli.h"
24 
25 #define MAXCOLORS	32768
26 #define INPUTQUANT	5
27 
28 #define REDUCE_GAMMA 2.2	/* Work in gamma compressed space */
29 				/* since this is closer to a linear */
30 				/* perceptual space */
31 
32 /* Ansi uses HUGE_VAL. */
33 #ifndef HUGE
34 #define HUGE HUGE_VAL
35 #endif
36 
37 typedef struct {
38 	float weightedvar,	/* weighted variance */
39 				/* (set to 0 if un-cutable) */
40 	mean[3];		/* centroid */
41 	unsigned long weight,	/* # of pixels in box */
42 	freq[3][256];		/* Projected frequencies */
43 	int low[3], high[3];	/* Box extent  - within low <= col < high */
44 } Box;
45 
46 static void QuantHistogram(Image *image, Box *box);
47 static void BoxStats(Box *box);
48 static void UpdateFrequencies(Box *box1, Box *box2);
49 static void ComputeRGBMap(Box *boxes, int colors, short unsigned int *rgbmap,
50 	int ditherf);
51 static void SetRGBmap(int boxnum, Box *box, short unsigned int *rgbmap);
52 static boolean CutBoxes(Box *boxes, int colors);
53 static int CutBox(Box *box, Box *newbox);
54 static int GreatestVariance(Box *boxes, int n);
55 static boolean FindCutpoint(Box *box, int color, Box *newbox1, Box *newbox2);
56 static void CopyToNewImage(Image *inimage, Image *outimage,
57 	unsigned short *rgbmap, int ditherf, int colors, float gamma,
58 	int verbose);
59 
60 
61 /* this converts a TLA-style pixel into a 15-bit true color pixel
62  */
63 
64 #define TLA_TO_15BIT(TABLE,PIXEL)           \
65   ((((unsigned short)((TABLE).red[PIXEL] & 0xf800)) >> 1) |   \
66    (((unsigned short)((TABLE).green[PIXEL] & 0xf800)) >> 6) | \
67    (((unsigned short)((TABLE).blue[PIXEL] & 0xf800)) >> 11))
68 
69 /* this converts a 24-bit true color pixel into a 15-bit true color pixel
70  */
71 
72 #define TRUE_TO_15BIT(PIXEL)     \
73   ((((unsigned long)((PIXEL) & 0xf80000)) >> 9) | \
74    (((unsigned long)((PIXEL) & 0x00f800)) >> 6) | \
75    (((unsigned long)((PIXEL) & 0x0000f8)) >> 3))
76 
77 /* these macros extract color intensities from a 15-bit true color pixel
78  */
79 
80 #define RED_INTENSITY(P)   (((unsigned short)((P) & 0x7c00)) >> 10)
81 #define GREEN_INTENSITY(P) (((unsigned short)((P) & 0x03e0)) >> 5)
82 #define BLUE_INTENSITY(P)   ((unsigned short)((P) & 0x001f))
83 
84 /*
85  * Value corrersponding to full intensity in colormap.  The values placed
86  * in the colormap are scaled to be between zero and this number.
87  */
88 #define MAX(x,y)	((x) > (y) ? (x) : (y))
89 
90 /*
91  * Readability constants.
92  */
93 #define REDI		0
94 #define GREENI		1
95 #define BLUEI		2
96 
97 static unsigned long *Histogram,	/* image histogram */
98  NPixels;			/* total # of pixels */
99 #define Bits INPUTQUANT
100 #define cBits (8-Bits)
101 #define ColormaxI (1 << Bits)	/* 2 ^ Bits */
102 #define Colormax (ColormaxI - 1)	/* quantized full-intensity */
103 #define ColormaxcI (1 << cBits)	/* 2 ^ (8-Bits) */
104 #define Colormaxc (ColormaxcI - 1)	/* quantized bits lost */
105 
106 /*
107  * if "ditherf" is True, apply color disthering
108  * if "Gamma" != 0.0, compensate for gamma post processing
109  */
110 unsigned int verbose;
reduce(Image * image,unsigned colors,int ditherf,float gamma,int verbose)111 Image *reduce(Image *image, unsigned colors, int ditherf, float gamma,
112 	int verbose)
113 {
114 	unsigned short *rgbmap;
115 	Box *Boxes;		/* Array of color boxes. */
116 	int i;			/* Counter */
117 	int OutColors;		/* # of entries computed */
118 	int depth;
119 	Image *new_image;
120 	char buf[BUFSIZ];
121 
122 	CURRFUNC("reduce");
123 
124 	if (colors > MAXCOLORS)	/* max # of colors we can handle */
125 		colors = MAXCOLORS;
126 
127 	if (BITMAPP(image) || (RGBP(image) && image->rgb.used <= colors))
128 		return (image);
129 
130 	if (TRUEP(image) && image->pixlen != 3) {
131 		fprintf(stderr, "reduce: true color image has strange pixel length?\n");
132 		return (image);
133 	}
134 	/* do color reduction and dithering in linear color space */
135 	if (GAMMA_NOT_EQUAL(image->gamma, REDUCE_GAMMA))
136 		gammacorrect(image, REDUCE_GAMMA, verbose);
137 
138 	NPixels = image->width * image->height;
139 
140 	Histogram = (unsigned long *) lcalloc(ColormaxI * ColormaxI * ColormaxI * sizeof(long));
141 	Boxes = (Box *) lmalloc(colors * sizeof(Box));
142 	rgbmap = (unsigned short *) lmalloc(ColormaxI * ColormaxI * ColormaxI * sizeof(unsigned short));
143 
144 	switch (image->type) {
145 	case IRGB:
146 		if (verbose) {
147 			printf("  Reducing RGB image color usage to %d colors...", colors);
148 			fflush(stdout);
149 		}
150 		QuantHistogram(image, &Boxes[0]);
151 		break;
152 
153 	case ITRUE:
154 		if (verbose) {
155 			printf("  Converting true color image to RGB image with %d colors...",
156 			       colors);
157 			fflush(stdout);
158 		}
159 		QuantHistogram(image, &Boxes[0]);
160 		break;
161 
162 	default:
163 		{
164 			lfree((char *) Histogram);
165 			lfree((char *) Boxes);
166 			lfree((char *) rgbmap);
167 			return (image);		/* not something we can reduce, thank you anyway */
168 		}
169 	}
170 
171 	OutColors = CutBoxes(Boxes, colors);
172 
173 	/*
174 	 * We now know the set of representative colors.  We now
175 	 * must fill in the colormap.
176 	 */
177 
178 	/* get destination image */
179 	depth = colorsToDepth(OutColors);
180 	new_image = newRGBImage(image->width, image->height, depth);
181 	snprintf(buf, sizeof(buf), "%s (%d colors)", image->title, OutColors);
182 	new_image->title = dupString(buf);
183 	new_image->gamma = image->gamma;
184 
185 	for (i = 0; i < OutColors; i++) {
186 		new_image->rgb.red[i] =
187 		    (((unsigned short) (Boxes[i].mean[REDI] + 0.5)) << (cBits + 8));
188 		new_image->rgb.green[i] =
189 		    (((unsigned short) (Boxes[i].mean[GREENI] + 0.5)) << (cBits + 8));
190 		new_image->rgb.blue[i] =
191 		    (((unsigned short) (Boxes[i].mean[BLUEI] + 0.5)) << (cBits + 8));
192 	}
193 	new_image->rgb.used = OutColors;
194 	new_image->rgb.compressed = TRUE;
195 
196 	ComputeRGBMap(Boxes, OutColors, rgbmap, ditherf);
197 	lfree((char *) Histogram);
198 	lfree((char *) Boxes);
199 
200 	/* copy old image into new image */
201 
202 	CopyToNewImage(image, new_image, rgbmap, ditherf, OutColors, gamma, verbose);
203 
204 	lfree((char *) rgbmap);
205 	if (verbose)
206 		printf("done\n");
207 	return (new_image);
208 }
209 
210 /*
211  * Compute the histogram of the image as well as the projected frequency
212  * arrays for the first world-encompassing box.
213  */
QuantHistogram(Image * image,Box * box)214 static void QuantHistogram(Image *image, Box *box)
215 {
216 	register byte *pixel;
217 	register unsigned long *rf, *gf, *bf;
218 	int x, y;
219 
220 	rf = box->freq[0];
221 	gf = box->freq[1];
222 	bf = box->freq[2];
223 
224 	/*
225 	 * Zero-out the projected frequency arrays of the largest box.
226 	 * We compute both the histogram and the proj. frequencies of
227 	 * the first box at the same time to save a pass through the
228 	 * entire image.
229 	 */
230 	bzero(rf, ColormaxI * sizeof(unsigned long));
231 	bzero(gf, ColormaxI * sizeof(unsigned long));
232 	bzero(bf, ColormaxI * sizeof(unsigned long));
233 
234 	if (image->type == IRGB) {
235 		Intensity *red = image->rgb.red, *green = image->rgb.green,
236 		*blue = image->rgb.blue;
237 		pixel = image->data;
238 		if (image->pixlen == 1)		/* special case most common this for speed */
239 			for (y = image->height; y > 0; y--)
240 				for (x = image->width; x > 0; x--) {
241 					register unsigned char r, g,
242 					 b;
243 					r = red[memToVal(pixel, 1)] >> (8 + cBits);
244 					g = green[memToVal(pixel, 1)] >> (8 + cBits);
245 					b = blue[memToVal(pixel, 1)] >> (8 + cBits);
246 					rf[r]++;
247 					gf[g]++;
248 					bf[b]++;
249 					Histogram[(((r << Bits) | g) << Bits) | b]++;
250 					pixel += 1;
251 		} else
252 			for (y = image->height; y > 0; y--)
253 				for (x = image->width; x > 0; x--) {
254 					register unsigned char r, g,
255 					 b;
256 					r = red[memToVal(pixel, image->pixlen)] >> (8 + cBits);
257 					g = green[memToVal(pixel, image->pixlen)] >> (8 + cBits);
258 					b = blue[memToVal(pixel, image->pixlen)] >> (8 + cBits);
259 					rf[r]++;
260 					gf[g]++;
261 					bf[b]++;
262 					Histogram[(((r << Bits) | g) << Bits) | b]++;
263 					pixel += image->pixlen;
264 				}
265 	} else {		/* assume ITRUE */
266 		pixel = image->data;
267 		if (image->pixlen == 3)		/* most common */
268 			for (y = image->height; y > 0; y--)
269 				for (x = image->width; x > 0; x--) {
270 					register unsigned long a, b;
271 					while (x > 4) {		/* Unroll 4 times */
272 						a = (*pixel++) >> cBits;	/* Red */
273 						rf[a]++;
274 						b = a << Bits;
275 						a = (*pixel++) >> cBits;	/* Green */
276 						gf[a]++;
277 						b = (b | a) << Bits;
278 						a = (*pixel++) >> cBits;	/* Blue */
279 						bf[a]++;
280 						b |= a;
281 						Histogram[b]++;
282 						a = (*pixel++) >> cBits;	/* Red */
283 						rf[a]++;
284 						b = a << Bits;
285 						a = (*pixel++) >> cBits;	/* Green */
286 						gf[a]++;
287 						b = (b | a) << Bits;
288 						a = (*pixel++) >> cBits;	/* Blue */
289 						bf[a]++;
290 						b |= a;
291 						Histogram[b]++;
292 						a = (*pixel++) >> cBits;	/* Red */
293 						rf[a]++;
294 						b = a << Bits;
295 						a = (*pixel++) >> cBits;	/* Green */
296 						gf[a]++;
297 						b = (b | a) << Bits;
298 						a = (*pixel++) >> cBits;	/* Blue */
299 						bf[a]++;
300 						b |= a;
301 						Histogram[b]++;
302 						a = (*pixel++) >> cBits;	/* Red */
303 						rf[a]++;
304 						b = a << Bits;
305 						a = (*pixel++) >> cBits;	/* Green */
306 						gf[a]++;
307 						b = (b | a) << Bits;
308 						a = (*pixel++) >> cBits;	/* Blue */
309 						bf[a]++;
310 						b |= a;
311 						Histogram[b]++;
312 						x -= 4;
313 					}
314 					a = (*pixel++) >> cBits;	/* Red */
315 					rf[a]++;
316 					b = a << Bits;
317 					a = (*pixel++) >> cBits;	/* Green */
318 					gf[a]++;
319 					b = (b | a) << Bits;
320 					a = (*pixel++) >> cBits;	/* Blue */
321 					bf[a]++;
322 					b |= a;
323 					Histogram[b]++;
324 		} else		/* less common */
325 			for (y = image->height; y > 0; y--)
326 				for (x = image->width; x > 0; x--) {
327 					register unsigned char r, g,
328 					 b;
329 					r = TRUE_RED(memToVal(pixel, image->pixlen)) >> cBits;
330 					rf[r]++;
331 					g = TRUE_GREEN(memToVal(pixel, image->pixlen)) >> cBits;
332 					gf[g]++;
333 					b = TRUE_BLUE(memToVal(pixel, image->pixlen)) >> cBits;
334 					bf[b]++;
335 					Histogram[(((r << Bits) | g) << Bits) | b]++;
336 					pixel += image->pixlen;
337 				}
338 	}
339 }
340 
341 /*
342    * Interatively cut the boxes.
343  */
CutBoxes(Box * boxes,int colors)344 static int CutBoxes(Box *boxes, int colors)
345 {
346 	int curbox;
347 
348 	boxes[0].low[REDI] = boxes[0].low[GREENI] = boxes[0].low[BLUEI] = 0;
349 	boxes[0].high[REDI] = boxes[0].high[GREENI] =
350 	    boxes[0].high[BLUEI] = ColormaxI;
351 	boxes[0].weight = NPixels;
352 
353 	BoxStats(&boxes[0]);
354 
355 	for (curbox = 1; curbox < colors;) {
356 		int n;
357 		n = GreatestVariance(boxes, curbox);
358 		if (n == curbox)
359 			break;	/* all are un-cutable */
360 		if (CutBox(&boxes[n], &boxes[curbox]))
361 			curbox++;	/* cut successfully */
362 	}
363 
364 	return curbox;
365 }
366 
367 /*
368  * Return the number of the box in 'boxes' with the greatest variance.
369  * Restrict the search to those boxes with indices between 0 and n-1.
370  * Return n if none of the boxes are cutable.
371  */
GreatestVariance(Box * boxes,int n)372 static int GreatestVariance(Box *boxes, int n)
373 {
374 	register int i, whichbox = n;
375 	float max = 0.0;
376 
377 	for (i = 0; i < n; i++) {
378 		if (boxes[i].weightedvar > max) {
379 			max = boxes[i].weightedvar;
380 			whichbox = i;
381 		}
382 	}
383 	return whichbox;
384 }
385 
386 /* Compute mean and weighted variance of the given box. */
BoxStats(Box * box)387 static void BoxStats(Box *box)
388 {
389 	register int i, color;
390 	unsigned long *freq;
391 	float mean, var;
392 
393 	if (box->weight == 0) {
394 		box->weightedvar = 0.0;
395 		return;
396 	}
397 	box->weightedvar = 0.;
398 	for (color = 0; color < 3; color++) {
399 		var = mean = 0;
400 		i = box->low[color];
401 		freq = &box->freq[color][i];
402 		for (; i < box->high[color]; i++, freq++) {
403 			mean += i * *freq;
404 			var += i * i * *freq;
405 		}
406 		box->mean[color] = mean / (float) box->weight;
407 		box->weightedvar += var - box->mean[color] * box->mean[color] *
408 		    (float) box->weight;
409 	}
410 	box->weightedvar /= NPixels;
411 }
412 
413 /*
414  * Cut the given box.  Returns TRUE if the box could be cut,
415  * FALSE (and weightedvar == 0.0) otherwise.
416  */
CutBox(Box * box,Box * newbox)417 static boolean CutBox(Box *box, Box *newbox)
418 {
419 	int i;
420 	double totalvar[3];
421 	Box newboxes[3][2];
422 
423 	if (box->weightedvar == 0.0 || box->weight == 0) {
424 		/*
425 		 * Can't cut this box.
426 		 */
427 		box->weightedvar = 0.0;		/* Don't try cutting it again */
428 		return FALSE;
429 	}
430 	/*
431 	 * Find 'optimal' cutpoint along each of the red, green and blue
432 	 * axes.  Sum the variances of the two boxes which would result
433 	 * by making each cut and store the resultant boxes for
434 	 * (possible) later use.
435 	 */
436 	for (i = 0; i < 3; i++) {
437 		if (FindCutpoint(box, i, &newboxes[i][0], &newboxes[i][1]))
438 			totalvar[i] = newboxes[i][0].weightedvar +
439 			    newboxes[i][1].weightedvar;
440 		else
441 			totalvar[i] = HUGE;	/* Not cutable on that axis */
442 	}
443 
444 	/*
445 	 * Find which of the three cuts minimized the total variance
446 	 * and make that the 'real' cut.
447 	 */
448 	if (totalvar[REDI] < HUGE &&
449 	    totalvar[REDI] <= totalvar[GREENI] &&
450 	    totalvar[REDI] <= totalvar[BLUEI]) {
451 		*box = newboxes[REDI][0];
452 		*newbox = newboxes[REDI][1];
453 		return TRUE;
454 	} else if (totalvar[GREENI] < HUGE &&
455 		   totalvar[GREENI] <= totalvar[REDI] &&
456 		   totalvar[GREENI] <= totalvar[BLUEI]) {
457 		*box = newboxes[GREENI][0];
458 		*newbox = newboxes[GREENI][1];
459 		return TRUE;
460 	} else if (totalvar[BLUEI] < HUGE) {
461 		*box = newboxes[BLUEI][0];
462 		*newbox = newboxes[BLUEI][1];
463 		return TRUE;
464 	}
465 	/* Unable to cut box on any axis - don't try again */
466 	box->weightedvar = 0.0;
467 	return FALSE;
468 }
469 
470 /*
471  * Compute the 'optimal' cutpoint for the given box along the axis
472  * indcated by 'color'.  Store the boxes which result from the cut
473  * in newbox1 and newbox2.
474  * If it is not possible to
475  */
FindCutpoint(Box * box,int color,Box * newbox1,Box * newbox2)476 static boolean FindCutpoint(Box *box, int color, Box *newbox1, Box *newbox2)
477 {
478 	float u, v, max;
479 	int i, maxindex, minindex, cutpoint;
480 	unsigned long optweight, curweight;
481 
482 	if (box->low[color] + 1 == box->high[color])
483 		return FALSE;	/* Cannot be cut. */
484 	minindex = (int) (((float) box->low[color] + box->mean[color]) * 0.5);
485 	maxindex = (int) ((box->mean[color] + (float) box->high[color]) * 0.5);
486 
487 	cutpoint = minindex;
488 	optweight = box->weight;
489 
490 	curweight = 0;
491 	for (i = box->low[color]; i < minindex; i++)
492 		curweight += box->freq[color][i];
493 	u = 0.0;
494 	max = -1.0;
495 	for (i = minindex; i <= maxindex; i++) {
496 		curweight += box->freq[color][i];
497 		if (curweight == box->weight)
498 			break;
499 		u += (float) (i * box->freq[color][i]) /
500 		    (float) box->weight;
501 		v = ((float) curweight / (float) (box->weight - curweight)) *
502 		    (box->mean[color] - u) * (box->mean[color] - u);
503 		if (v > max) {
504 			max = v;
505 			cutpoint = i;
506 			optweight = curweight;
507 		}
508 	}
509 	cutpoint++;
510 	*newbox1 = *newbox2 = *box;
511 	newbox1->weight = optweight;
512 	newbox2->weight -= optweight;
513 	if (newbox1->weight == 0 || newbox2->weight == 0)
514 		return FALSE;	/* Unable to cut on this axis */
515 	newbox1->high[color] = cutpoint;
516 	newbox2->low[color] = cutpoint;
517 	UpdateFrequencies(newbox1, newbox2);
518 	BoxStats(newbox1);
519 	BoxStats(newbox2);
520 
521 	return TRUE;		/* Found cutpoint. */
522 }
523 
524 /*
525  * Update projected frequency arrays for two boxes which used to be
526  * a single box. Also shrink the box sizes to fit the points.
527  */
UpdateFrequencies(Box * box1,Box * box2)528 static void UpdateFrequencies(Box *box1, Box *box2)
529 {
530 	register unsigned long myfreq, *h;
531 	register int b, g, r;
532 	int roff;
533 	register unsigned long *b1f0, *b1f1, *b1f2, *b2f0, *b2f1, *b2f2;
534 
535 	b1f0 = &box1->freq[0][0];
536 	b1f1 = &box1->freq[1][0];
537 	b1f2 = &box1->freq[2][0];
538 	b2f0 = &box2->freq[0][0];
539 	b2f1 = &box2->freq[1][0];
540 	b2f2 = &box2->freq[2][0];
541 
542 	bzero(b1f0, ColormaxI * sizeof(unsigned long));
543 	bzero(b1f1, ColormaxI * sizeof(unsigned long));
544 	bzero(b1f2, ColormaxI * sizeof(unsigned long));
545 
546 	for (r = box1->low[0]; r < box1->high[0]; r++) {
547 		roff = r << Bits;
548 		for (g = box1->low[1]; g < box1->high[1]; g++) {
549 			b = box1->low[2];
550 			h = Histogram + (((roff | g) << Bits) | b);
551 			for (; b < box1->high[2]; b++) {
552 				myfreq = *h++;
553 				if (myfreq == 0)
554 					continue;
555 				b1f0[r] += myfreq;
556 				b1f1[g] += myfreq;
557 				b1f2[b] += myfreq;
558 				b2f0[r] -= myfreq;
559 				b2f1[g] -= myfreq;
560 				b2f2[b] -= myfreq;
561 			}
562 		}
563 	}
564 
565 	/* shrink the boxes to fit the new points */
566 	for (r = 0; r < 3; r++) {
567 		register int l1, l2, h1, h2;
568 		l1 = l2 = ColormaxI;
569 		h1 = h2 = 0;
570 
571 		for (g = 0; g < ColormaxI; g++) {
572 			if (box1->freq[r][g] != 0) {
573 				if (g < l1)
574 					l1 = g;
575 				if (g > h1)
576 					h1 = g;
577 			}
578 			if (box2->freq[r][g] != 0) {
579 				if (g < l2)
580 					l2 = g;
581 				if (g > h2)
582 					h2 = g;
583 			}
584 		}
585 		box1->low[r] = l1;
586 		box2->low[r] = l2;
587 		box1->high[r] = h1 + 1;
588 		box2->high[r] = h2 + 1;
589 	}
590 }
591 
592 /* Compute RGB to colormap index map. */
ComputeRGBMap(Box * boxes,int colors,unsigned short * rgbmap,int ditherf)593 static void ComputeRGBMap(Box *boxes, int colors, unsigned short *rgbmap,
594 	int ditherf)
595 {
596 	register int i;
597 
598 	if (!ditherf) {
599 		/*
600 		 * The centroid of each box serves as the representative
601 		 * for each color in the box.
602 		 */
603 		for (i = 0; i < colors; i++)
604 			SetRGBmap(i, &boxes[i], rgbmap);
605 	} else {
606 		/* Initialise the pixel value to illegal value.
607 		 * Will fill in entry on demand.
608 		 */
609 		bfill((char *) rgbmap, ColormaxI * ColormaxI * ColormaxI
610 			* sizeof(*rgbmap), 0xff);
611 		/* speed things up by pre-filling boxes with input colors */
612 		for (i = 0; i < colors; i++)
613 			SetRGBmap(i, &boxes[i], rgbmap);
614 	}
615 }
616 
617 /*
618  * Make the centroid of "boxnum" serve as the representative for
619  * each color in the box.
620  */
SetRGBmap(int boxnum,Box * box,unsigned short * rgbmap)621 static void SetRGBmap(int boxnum, Box *box, unsigned short *rgbmap)
622 {
623 	register int r, g, b;
624 
625 	for (r = box->low[REDI]; r < box->high[REDI]; r++) {
626 		for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
627 			for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
628 				rgbmap[(((r << Bits) | g) << Bits) | b] = (unsigned short) boxnum;
629 			}
630 		}
631 	}
632 }
633 
634 
635 /* nearest neighbor structure */
636 typedef struct {
637 	int length;		/* number of colors */
638 	unsigned char *red, *green, *blue;
639 	unsigned short *pixel;
640 } NN;
641 
642 #define NNBits 3
643 #define NNcBits (8-NNBits)
644 #define NNmaxI (1 << NNBits)	/* 2 ^ Bits */
645 #define NNmax (NNmaxI - 1)	/* quantized full-intensity */
646 #define NNmaxcI (1 << NNcBits)	/* 2 ^ (8 - Bits) */
647 #define NNmaxc (NNmaxcI - 1)	/* quantized bits lost */
648 
649 /* Initialise the appropriate nn cell */
find_nnearest(int rval,int gval,int bval,NN * nnp,Image * outimage)650 static void find_nnearest(int rval, int gval, int bval, NN *nnp,
651 	Image *outimage)
652 {
653 	long bdist, mdist;
654 	Intensity *red = outimage->rgb.red, *green = outimage->rgb.green,
655 	*blue = outimage->rgb.blue;
656 	long *ds, *de, dists[MAXCOLORS];
657 	int *is, indexes[MAXCOLORS];
658 	int i, length = 0;
659 	double tdouble;
660 
661 	rval = ((rval & ~NNmaxc) << 1) + NNmaxc;	/* closest color to center of cell */
662 	gval = ((gval & ~NNmaxc) << 1) + NNmaxc;	/* (2 * scaled values) */
663 	bval = ((bval & ~NNmaxc) << 1) + NNmaxc;
664 
665 #define SQUARE(aa) ((aa) * ( aa))
666 #define QERR(r1,g1,b1,r2,g2,b2) /* sum of squares distance (in 2 * scale) */ \
667 	(SQUARE(r1 - (r2 >>7)) +  \
668 	 SQUARE(g1 - (g2 >>7)) +  \
669 	 SQUARE(b1 - (b2 >>7)))
670 
671 	/* Calculate the distances of all the colors in the colormap */
672 	/* from the center of this nn cell. */
673 	/* (unroll this loop 4 times) */
674 	for (ds = dists, de = ds + outimage->rgb.used; ds < (de - 3);) {
675 		*ds = QERR(rval, gval, bval, *red, *green, *blue);
676 		ds++;
677 		red++;
678 		green++;
679 		blue++;
680 		*ds = QERR(rval, gval, bval, *red, *green, *blue);
681 		ds++;
682 		red++;
683 		green++;
684 		blue++;
685 		*ds = QERR(rval, gval, bval, *red, *green, *blue);
686 		ds++;
687 		red++;
688 		green++;
689 		blue++;
690 		*ds = QERR(rval, gval, bval, *red, *green, *blue);
691 		ds++;
692 		red++;
693 		green++;
694 		blue++;
695 	}
696 	/* finish off remaining */
697 	while (ds < de) {
698 		*ds = QERR(rval, gval, bval, *red, *green, *blue);
699 		ds++;
700 		red++;
701 		green++;
702 		blue++;
703 	}
704 
705 	/* find nearest colormap colour to the center */
706 	/* (unroll this loop 8 times) */
707 	for (bdist = 0x7fffffff, ds = dists; ds < (de - 7); ds++) {
708 		if (*ds < bdist)
709 			goto found_small;
710 		ds++;
711 		if (*ds < bdist)
712 			goto found_small;
713 		ds++;
714 		if (*ds < bdist)
715 			goto found_small;
716 		ds++;
717 		if (*ds < bdist)
718 			goto found_small;
719 		ds++;
720 		if (*ds < bdist)
721 			goto found_small;
722 		ds++;
723 		if (*ds < bdist)
724 			goto found_small;
725 		ds++;
726 		if (*ds < bdist)
727 			goto found_small;
728 		ds++;
729 		if (*ds >= bdist)
730 			continue;
731 	      found_small:
732 		bdist = *ds;
733 	}
734 	/* finish off remaining */
735 	while (ds < de) {
736 		if (*ds < bdist)
737 			bdist = *ds;
738 		ds++;
739 	}
740 
741 	/* figure out the maximum distance from center */
742 	tdouble = sqrt((double) bdist) + sqrt((double) (NNmaxc * NNmaxc * 12));
743 	mdist = (int) ((tdouble * tdouble) + 0.5);	/* back to square distance */
744 
745 	/* now select the neighborhood groups */
746 	/* (unroll this loop 8 times) */
747 	for (is = indexes, ds = dists; ds < (de - 7); ds++) {
748 		if (*ds <= mdist)
749 			goto found_neighbor;
750 		ds++;
751 		if (*ds <= mdist)
752 			goto found_neighbor;
753 		ds++;
754 		if (*ds <= mdist)
755 			goto found_neighbor;
756 		ds++;
757 		if (*ds <= mdist)
758 			goto found_neighbor;
759 		ds++;
760 		if (*ds <= mdist)
761 			goto found_neighbor;
762 		ds++;
763 		if (*ds <= mdist)
764 			goto found_neighbor;
765 		ds++;
766 		if (*ds <= mdist)
767 			goto found_neighbor;
768 		ds++;
769 		if (*ds > mdist)
770 			continue;
771 	      found_neighbor:
772 		*is++ = ds - dists;
773 	}
774 	/* finish off remaining */
775 	for (; ds < de; ds++) {
776 		if (*ds <= mdist)
777 			*is++ = ds - dists;
778 	}
779 
780 	/* now make up the entries in the neighborhood structure entry */
781 	length = is - indexes;
782 	nnp->length = length;
783 	nnp->red = (unsigned char *) lmalloc(length * 3);
784 	nnp->green = nnp->red + length;
785 	nnp->blue = nnp->green + length;
786 	nnp->pixel = (unsigned short *) lmalloc(length * sizeof(unsigned short));
787 	red = outimage->rgb.red;
788 	green = outimage->rgb.green;
789 	blue = outimage->rgb.blue;
790 	for (i = 0, is = indexes; i < length; i++, is++) {
791 		nnp->red[i] = red[*is] >> 8;
792 		nnp->green[i] = green[*is] >> 8;
793 		nnp->blue[i] = blue[*is] >> 8;
794 		nnp->pixel[i] = *is;
795 	}
796 
797 #undef SQUARE
798 #undef QERR
799 }
800 
801 /*
802  * Find the neareset chosen pixel value to the given color.
803  * Uses nn accelerator.
804  */
find_nearest(int rval,int gval,int bval,NN * nna,Image * outimage)805 static unsigned short find_nearest(int rval, int gval, int bval, NN *nna,
806 	Image *outimage)
807 {
808 	long bdist, ndist;
809 	int i, bindex = 0;
810 	int lastc;
811 	unsigned char *red, *green, *blue;
812 	int nnindex;
813 	NN *nnp;
814 
815 	nnindex = ((((rval >> NNcBits) << NNBits) | (gval >> NNcBits)) << NNBits) | (bval >> NNcBits);
816 
817 	nnp = nna + nnindex;
818 	if (nnp->length == 0)	/* Create nearest neighbor color mapping on demand */
819 		find_nnearest(rval, gval, bval, nnp, outimage);
820 	lastc = nnp->length - 1;
821 	red = &nnp->red[lastc],
822 	    green = &nnp->green[lastc],
823 	    blue = &nnp->blue[lastc];
824 
825 	rval = (rval & ~Colormaxc) + (ColormaxcI >> 1);		/* closest color to center of 15 bit color */
826 	gval = (gval & ~Colormaxc) + (ColormaxcI >> 1);
827 	bval = (bval & ~Colormaxc) + (ColormaxcI >> 1);
828 
829 #define SQUARE(aa) ((aa) * ( aa))
830 #define QERR(r1,g1,b1,r2,g2,b2) \
831 	(SQUARE(r1 - r2) +  \
832 	 SQUARE(g1 - g2) +  \
833 	 SQUARE(b1 - b2))
834 
835 	/* now find the the nearest colour */
836 	/* (try and unroll this loop 8 times) */
837 	for (bdist = 0x7fffffff, i = lastc; i >= 7; i--, red--, green--, blue--) {
838 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
839 			goto found_small;
840 		i--;
841 		red--;
842 		green--;
843 		blue--;
844 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
845 			goto found_small;
846 		i--;
847 		red--;
848 		green--;
849 		blue--;
850 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
851 			goto found_small;
852 		i--;
853 		red--;
854 		green--;
855 		blue--;
856 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
857 			goto found_small;
858 		i--;
859 		red--;
860 		green--;
861 		blue--;
862 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
863 			goto found_small;
864 		i--;
865 		red--;
866 		green--;
867 		blue--;
868 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
869 			goto found_small;
870 		i--;
871 		red--;
872 		green--;
873 		blue--;
874 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist)
875 			goto found_small;
876 		i--;
877 		red--;
878 		green--;
879 		blue--;
880 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) >= bdist)
881 			continue;
882 	      found_small:
883 		bdist = ndist;
884 		bindex = i;
885 	}
886 	/* finish off remaining */
887 	for (; i >= 0; i--, red--, green--, blue--) {
888 		if ((ndist = QERR(rval, gval, bval, *red, *green, *blue)) < bdist) {
889 			bdist = ndist;
890 			bindex = i;
891 		}
892 	}
893 	return nnp->pixel[bindex];
894 #undef SQUARE
895 #undef QERR
896 }
897 
898 /* "rgbmap" is pixel value lookup map */
CopyToNewImage(Image * inimage,Image * outimage,unsigned short * rgbmap,int ditherf,int colors,float gamma,int verbose)899 static void CopyToNewImage(Image *inimage, Image *outimage,
900 	unsigned short *rgbmap, int ditherf, int colors, float gamma,
901 	int verbose)
902 {
903 	byte *pixel, *dpixel;
904 	int x, y;
905 
906 	pixel = inimage->data;
907 	dpixel = outimage->data;
908 
909 	if (!ditherf) {
910 		switch (inimage->type) {
911 		case IRGB:
912 			if (outimage->pixlen == 1 && inimage->pixlen == 1)	/* most common */
913 				for (y = inimage->height; y > 0; y--)
914 					for (x = inimage->width; x > 0; x--) {
915 						register unsigned long temp;
916 						temp = TLA_TO_15BIT(inimage->rgb, *pixel);
917 						temp = rgbmap[temp];
918 						valToMem(temp, dpixel, 1);
919 						pixel++;
920 						dpixel++;
921 			} else	/* less common */
922 				for (y = inimage->height; y > 0; y--)
923 					for (x = inimage->width; x > 0; x--) {
924 						register unsigned long temp;
925 						temp = memToVal(pixel, inimage->pixlen);
926 						temp = TLA_TO_15BIT(inimage->rgb, temp);
927 						temp = rgbmap[temp];
928 						valToMem(temp, dpixel, outimage->pixlen);
929 						pixel += inimage->pixlen;
930 						dpixel += outimage->pixlen;
931 					}
932 			break;
933 
934 		case ITRUE:
935 			if (inimage->pixlen == 3 && outimage->pixlen == 1) {	/* common */
936 				for (y = inimage->height; y > 0; y--) {
937 					register unsigned long temp;
938 					x = inimage->width;
939 					while (x >= 4) {	/* Unroll 4 times */
940 						temp = (((unsigned long) *pixel) & 0xf8) << (Bits + Bits - cBits);
941 						pixel++;
942 						temp |= (((unsigned long) *pixel) & 0xf8) << (Bits - cBits);
943 						pixel++;
944 						temp |= ((unsigned long) *pixel) >> cBits;
945 						*dpixel = rgbmap[temp];
946 						pixel++;
947 						dpixel++;
948 						temp = (((unsigned long) *pixel) & 0xf8) << (Bits + Bits - cBits);
949 						pixel++;
950 						temp |= (((unsigned long) *pixel) & 0xf8) << (Bits - cBits);
951 						pixel++;
952 						temp |= ((unsigned long) *pixel) >> cBits;
953 						*dpixel = rgbmap[temp];
954 						pixel++;
955 						dpixel++;
956 						temp = (((unsigned long) *pixel) & 0xf8) << (Bits + Bits - cBits);
957 						pixel++;
958 						temp |= (((unsigned long) *pixel) & 0xf8) << (Bits - cBits);
959 						pixel++;
960 						temp |= ((unsigned long) *pixel) >> cBits;
961 						*dpixel = rgbmap[temp];
962 						pixel++;
963 						dpixel++;
964 						temp = (((unsigned long) *pixel) & 0xf8) << (Bits + Bits - cBits);
965 						pixel++;
966 						temp |= (((unsigned long) *pixel) & 0xf8) << (Bits - cBits);
967 						pixel++;
968 						temp |= ((unsigned long) *pixel) >> cBits;
969 						*dpixel = rgbmap[temp];
970 						pixel++;
971 						dpixel++;
972 						x -= 4;
973 					}
974 					while (x >= 1) {	/* Finish remaining */
975 						temp = (((unsigned long) *pixel) & 0xf8) << (Bits + Bits - cBits);
976 						pixel++;
977 						temp |= (((unsigned long) *pixel) & 0xf8) << (Bits - cBits);
978 						pixel++;
979 						temp |= ((unsigned long) *pixel) >> cBits;
980 						*dpixel = rgbmap[temp];
981 						pixel++;
982 						dpixel++;
983 						x--;
984 					}
985 				}
986 			} else	/* less common */
987 				for (y = inimage->height; y > 0; y--)
988 					for (x = inimage->width; x > 0; x--) {
989 						register unsigned long temp;
990 						temp = memToVal(pixel, inimage->pixlen);
991 						temp = TRUE_TO_15BIT(temp);
992 						temp = rgbmap[temp];
993 						valToMem(temp, dpixel, outimage->pixlen);
994 						pixel += inimage->pixlen;
995 						dpixel += outimage->pixlen;
996 					}
997 			break;
998 		}
999 	} else {		/* else use dithering */
1000 		NN *nna, *nnp;
1001 		Intensity *ored, *ogreen, *oblue;
1002 		int satary[3 * 256 + 20], *sat = &satary[256 + 10];
1003 		int *epixel, *ip, nextr, nextg, nextb;
1004 		int err1ary[512], *err1 = &err1ary[256];
1005 		int err3ary[512], *err3 = &err3ary[256];
1006 		int err5ary[512], *err5 = &err5ary[256];
1007 		int err7ary[512], *err7 = &err7ary[256];
1008 
1009 		if (verbose) {
1010 			printf("\n  (Using color dithering to reduce the impact of restricted colors)...");
1011 		}
1012 		nna = (NN *) lcalloc((NNmaxI * NNmaxI * NNmaxI) * sizeof(NN));
1013 
1014 		/* Init saturation table */
1015 		for (x = -256 - 10; x < (512 + 10); x++) {
1016 			y = x;
1017 			if (y < 0)
1018 				y = 0;
1019 			else if (y > 255)
1020 				y = 255;
1021 			sat[x] = y;
1022 		}
1023 		{
1024 			int maxerr;	/* Fudge factor */
1025 			maxerr = (int) (256.0 / pow(((double) colors), 0.45));
1026 			/* Init error lookup array */
1027 			for (x = -255; x < 256; x++) {
1028 				int err, ro;
1029 				err = x;
1030 				if (err > maxerr)
1031 					err = maxerr;
1032 				else if (err < -maxerr)
1033 					err = -maxerr;
1034 				err7[x] = (err * 7) / 16;
1035 				err5[x] = (err * 5) / 16;
1036 				err3[x] = (err * 3) / 16;
1037 				err1[x] = (err * 1) / 16;
1038 				ro = err - err1[x] - err3[x] - err5[x] - err7[x];
1039 				switch (err & 3) {	/* Distribute round off error */
1040 				case 0:
1041 					err1[x] += ro;
1042 					break;
1043 				case 1:
1044 					err5[x] += ro;
1045 					break;
1046 				case 2:
1047 					err3[x] += ro;
1048 					break;
1049 				case 3:
1050 					err7[x] += ro;
1051 					break;
1052 				}
1053 			}
1054 		}
1055 		epixel = ((int *) lcalloc((inimage->width + 2) * 3 * sizeof(int))) + 3;
1056 
1057 		/* set up index mapping tables */
1058 		ored = (Intensity *) lmalloc(sizeof(Intensity) * outimage->rgb.size);
1059 		ogreen = (Intensity *) lmalloc(sizeof(Intensity) * outimage->rgb.size);
1060 		oblue = (Intensity *) lmalloc(sizeof(Intensity) * outimage->rgb.size);
1061 		if (UNSET_GAMMA == gamma || (gamma / inimage->gamma) == 1.0)
1062 			for (x = 0; x < outimage->rgb.size; x++) {
1063 				ored[x] = outimage->rgb.red[x] >> 8;
1064 				ogreen[x] = outimage->rgb.green[x] >> 8;
1065 				oblue[x] = outimage->rgb.blue[x] >> 8;
1066 		} else
1067 			/* Create map that accounts for quantizing effect of output */
1068 			/* gamma mapping applied to small values. */
1069 			for (x = 0; x < outimage->rgb.size; x++) {
1070 				int val;
1071 				val = outimage->rgb.red[x] >> 8;
1072 				val = (int) (0.5 + 255 * pow(val / 255.0, inimage->gamma / gamma));
1073 				val = (int) (0.5 + 255 * pow(val / 255.0, gamma / inimage->gamma));
1074 				ored[x] = val;
1075 				val = outimage->rgb.green[x] >> 8;
1076 				val = (int) (0.5 + 255 * pow(val / 255.0, inimage->gamma / gamma));
1077 				val = (int) (0.5 + 255 * pow(val / 255.0, gamma / inimage->gamma));
1078 				ogreen[x] = val;
1079 				val = outimage->rgb.blue[x] >> 8;
1080 				val = (int) (0.5 + 255 * pow(val / 255.0, inimage->gamma / gamma));
1081 				val = (int) (0.5 + 255 * pow(val / 255.0, gamma / inimage->gamma));
1082 				oblue[x] = val;
1083 			}
1084 
1085 		switch (inimage->type) {
1086 		case IRGB:
1087 			{
1088 				Intensity *red = inimage->rgb.red, *green = inimage->rgb.green,
1089 				*blue = inimage->rgb.blue;
1090 				if (outimage->pixlen == 1 && inimage->pixlen == 1) {	/* most common */
1091 					for (x = inimage->width, ip = epixel; x > 0; x--, pixel++) {
1092 						*ip++ = red[*pixel] >> 8;;
1093 						*ip++ = green[*pixel] >> 8;;
1094 						*ip++ = blue[*pixel] >> 8;;
1095 					}
1096 					for (y = inimage->height; y > 0; y--) {
1097 						if (y == 1)	/* protect last line */
1098 							pixel = inimage->data;
1099 						for (x = inimage->width, ip = epixel, nextr = *ip, *ip++ = 0,
1100 						     nextg = *ip, *ip++ = 0, nextb = *ip, *ip++ = 0; x > 0; x--) {
1101 							register unsigned long rgbindex;
1102 							register int rval,
1103 							 gval, bval,
1104 							 color;
1105 
1106 							rgbindex = ((rval = sat[nextr]) & 0xf8) << (Bits + Bits - cBits);
1107 							rgbindex |= ((gval = sat[nextg]) & 0xf8) << (Bits - cBits);
1108 							rgbindex |= (bval = sat[nextb]) >> cBits;
1109 							*dpixel = color = rgbmap[rgbindex];
1110 							if (color == 0xffff)	/* Create nearest color mapping on demand */
1111 								*dpixel = color = rgbmap[rgbindex] = find_nearest(rval, gval, bval, nna, outimage);
1112 							dpixel++;
1113 							rval -= ored[color];
1114 							nextr = *ip++ + err7[rval];	/* this line, next pixel */
1115 							ip[-1] = err1[rval];	/* next line, next pixel */
1116 							ip[-4] += (red[*pixel] >> 8) + err5[rval];	/* next line, this pixel */
1117 							ip[-7] += err3[rval];	/* next line, last pixel */
1118 							gval -= ogreen[color];
1119 							nextg = *ip++ + err7[gval];	/* this line, next pixel */
1120 							ip[-1] = err1[gval];	/* next line, next pixel */
1121 							ip[-4] += (green[*pixel] >> 8) + err5[gval];	/* next line, this pixel */
1122 							ip[-7] += err3[gval];	/* next line, last pixel */
1123 							bval -= oblue[color];
1124 							nextb = *ip++ + err7[bval];	/* this line, next pixel */
1125 							ip[-1] = err1[bval];	/* next line, next pixel */
1126 							ip[-4] += (blue[*pixel++] >> 8) + err5[bval];	/* next line, this pixel */
1127 							ip[-7] += err3[bval];	/* next line, last pixel */
1128 						}
1129 					}
1130 				} else {	/* less common */
1131 					unsigned long temp;
1132 					for (x = inimage->width, ip = epixel; x > 0; x--, pixel += inimage->pixlen) {
1133 						temp = memToVal(pixel, inimage->pixlen);
1134 						*ip++ = red[temp] >> 8;;
1135 						*ip++ = green[temp] >> 8;;
1136 						*ip++ = blue[temp] >> 8;;
1137 					}
1138 					for (y = inimage->height; y > 0; y--) {
1139 						if (y == 1)	/* protect last line */
1140 							pixel = inimage->data;
1141 						for (x = inimage->width, ip = epixel, nextr = *ip, *ip++ = 0,
1142 						     nextg = *ip, *ip++ = 0, nextb = *ip, *ip++ = 0; x > 0; x--) {
1143 							register unsigned long rgbindex;
1144 							register int rval,
1145 							 gval, bval,
1146 							 color;
1147 
1148 							rgbindex = ((rval = sat[nextr]) & 0xf8) << (Bits + Bits - cBits);
1149 							rgbindex |= ((gval = sat[nextg]) & 0xf8) << (Bits - cBits);
1150 							rgbindex |= (bval = sat[nextb]) >> cBits;
1151 							color = rgbmap[rgbindex];
1152 							if (color == 0xffff)	/* Create nearest color mapping on demand */
1153 								color = rgbmap[rgbindex] = find_nearest(rval, gval, bval, nna, outimage);
1154 							valToMem(color, dpixel, outimage->pixlen);
1155 							dpixel += outimage->pixlen;
1156 							temp = memToVal(pixel, inimage->pixlen);
1157 							pixel += inimage->pixlen;
1158 							rval -= ored[color];
1159 							nextr = *ip++ + err7[rval];	/* this line, next pixel */
1160 							ip[-1] = err1[rval];	/* next line, next pixel */
1161 							ip[-4] += (red[temp] >> 8) + err5[rval];	/* next line, this pixel */
1162 							ip[-7] += err3[rval];	/* next line, last pixel */
1163 							gval -= ogreen[color];
1164 							nextg = *ip++ + err7[gval];	/* this line, next pixel */
1165 							ip[-1] = err1[gval];	/* next line, next pixel */
1166 							ip[-4] += (green[temp] >> 8) + err5[gval];	/* next line, this pixel */
1167 							ip[-7] += err3[gval];	/* next line, last pixel */
1168 							bval -= oblue[color];
1169 							nextb = *ip++ + err7[bval];	/* this line, next pixel */
1170 							ip[-1] = err1[bval];	/* next line, next pixel */
1171 							ip[-4] += (blue[temp] >> 8) + err5[bval];	/* next line, this pixel */
1172 							ip[-7] += err3[bval];	/* next line, last pixel */
1173 						}
1174 					}
1175 				}
1176 			}
1177 			break;
1178 
1179 		case ITRUE:
1180 			if (inimage->pixlen == 3 && outimage->pixlen == 1) {	/* common */
1181 				for (x = inimage->width * 3, ip = epixel; x > 0; x--, ip++, pixel++) {
1182 					*ip = *pixel;
1183 				}
1184 				for (y = inimage->height; y > 0; y--) {
1185 					if (y == 1)	/* protect last line */
1186 						pixel = inimage->data;
1187 					for (x = inimage->width, ip = epixel, nextr = *ip, *ip++ = 0,
1188 					     nextg = *ip, *ip++ = 0, nextb = *ip, *ip++ = 0; x > 0; x--) {
1189 						register unsigned long rgbindex;
1190 						register int rval, gval,
1191 						 bval, color;
1192 
1193 						rgbindex = ((rval = sat[nextr]) & 0xf8) << (Bits + Bits - cBits);
1194 						rgbindex |= ((gval = sat[nextg]) & 0xf8) << (Bits - cBits);
1195 						rgbindex |= (bval = sat[nextb]) >> cBits;
1196 						*dpixel = color = rgbmap[rgbindex];
1197 						if (color == 0xffff)	/* Create nearest color mapping on demand */
1198 							*dpixel = color = rgbmap[rgbindex] = find_nearest(rval, gval, bval, nna, outimage);
1199 						dpixel++;
1200 						rval -= ored[color];
1201 						nextr = *ip++ + err7[rval];	/* this line, next pixel */
1202 						ip[-1] = err1[rval];	/* next line, next pixel */
1203 						ip[-4] += *pixel++ + err5[rval];	/* next line, this pixel */
1204 						ip[-7] += err3[rval];	/* next line, last pixel */
1205 						gval -= ogreen[color];
1206 						nextg = *ip++ + err7[gval];	/* this line, next pixel */
1207 						ip[-1] = err1[gval];	/* next line, next pixel */
1208 						ip[-4] += *pixel++ + err5[gval];	/* next line, this pixel */
1209 						ip[-7] += err3[gval];	/* next line, last pixel */
1210 						bval -= oblue[color];
1211 						nextb = *ip++ + err7[bval];	/* this line, next pixel */
1212 						ip[-1] = err1[bval];	/* next line, next pixel */
1213 						ip[-4] += *pixel++ + err5[bval];	/* next line, this pixel */
1214 						ip[-7] += err3[bval];	/* next line, last pixel */
1215 					}
1216 				}
1217 			} else {	/* less common */
1218 				unsigned long temp;
1219 				for (x = inimage->width, ip = epixel; x > 0; x--, pixel += inimage->pixlen) {
1220 					temp = memToVal(pixel, inimage->pixlen);
1221 					*ip++ = TRUE_RED(temp);
1222 					*ip++ = TRUE_GREEN(temp);
1223 					*ip++ = TRUE_BLUE(temp);
1224 				}
1225 				for (y = inimage->height; y > 0; y--) {
1226 					if (y == 1)	/* protect last line */
1227 						pixel = inimage->data;
1228 					for (x = inimage->width, ip = epixel, nextr = *ip, *ip++ = 0,
1229 					     nextg = *ip, *ip++ = 0, nextb = *ip, *ip++ = 0; x > 0; x--) {
1230 						register unsigned long rgbindex;
1231 						register int rval, gval,
1232 						 bval, color;
1233 
1234 						rgbindex = ((rval = sat[nextr]) & 0xf8) << (Bits + Bits - cBits);
1235 						rgbindex |= ((gval = sat[nextg]) & 0xf8) << (Bits - cBits);
1236 						rgbindex |= (bval = sat[nextb]) >> cBits;
1237 						color = rgbmap[rgbindex];
1238 						if (color == 0xffff)	/* Create nearest color mapping on demand */
1239 							color = rgbmap[rgbindex] = find_nearest(rval, gval, bval, nna, outimage);
1240 						valToMem(color, dpixel, outimage->pixlen);
1241 						dpixel += outimage->pixlen;
1242 						temp = memToVal(pixel, inimage->pixlen);
1243 						pixel += inimage->pixlen;
1244 						rval -= ored[color];
1245 						nextr = *ip++ + err7[rval];	/* this line, next pixel */
1246 						ip[-1] = err1[rval];	/* next line, next pixel */
1247 						ip[-4] += err5[rval]	/* next line, this pixel */
1248 						    +TRUE_RED(temp);
1249 						ip[-7] += err3[rval];	/* next line, last pixel */
1250 						gval -= ogreen[color];
1251 						nextg = *ip++ + err7[gval];	/* this line, next pixel */
1252 						ip[-1] = err1[gval];	/* next line, next pixel */
1253 						ip[-4] += err5[gval]	/* next line, this pixel */
1254 						    +TRUE_GREEN(temp);
1255 						ip[-7] += err3[gval];	/* next line, last pixel */
1256 						bval -= oblue[color];
1257 						nextb = *ip++ + err7[bval];	/* this line, next pixel */
1258 						ip[-1] = err1[bval];	/* next line, next pixel */
1259 						ip[-4] += err5[bval]	/* next line, this pixel */
1260 						    +TRUE_BLUE(temp);
1261 						ip[-7] += err3[bval];	/* next line, last pixel */
1262 					}
1263 				}
1264 			}
1265 			break;
1266 		}
1267 
1268 		lfree((byte *) ored);
1269 		lfree((byte *) ogreen);
1270 		lfree((byte *) oblue);
1271 		lfree((byte *) (epixel - 3));
1272 
1273 		/* Clean up nn stuff */
1274 		for (x = NNmaxI * NNmaxI * NNmaxI, nnp = nna; x > 0; x--, nnp++) {
1275 			if (nnp->length != 0) {
1276 				/* free red, green and blue */
1277 				lfree((byte *) nnp->red);
1278 				lfree((byte *) nnp->pixel);
1279 			}
1280 		}
1281 		lfree((char *) nna);
1282 	}
1283 }
1284 
1285 /* expand an image into a true color image */
expandtotrue(Image * image)1286 Image *expandtotrue(Image *image)
1287 {
1288 	Pixel fg, bg;
1289 	Image *new_image;
1290 	int x, y;
1291 	byte *spixel, *dpixel, *line;
1292 	unsigned int linelen;
1293 	byte mask;
1294 
1295 	CURRFUNC("expandtotrue");
1296 	if TRUEP
1297 		(image)
1298 		    return (image);
1299 
1300 	new_image = newTrueImage(image->width, image->height);
1301 	new_image->title = dupString(image->title);
1302 	new_image->gamma = image->gamma;
1303 
1304 	switch (image->type) {
1305 	case IBITMAP:
1306 		fg = RGB_TO_TRUE(image->rgb.red[1], image->rgb.green[1], image->rgb.blue[1]);
1307 		bg = RGB_TO_TRUE(image->rgb.red[0], image->rgb.green[0], image->rgb.blue[0]);
1308 		line = image->data;
1309 		dpixel = new_image->data;
1310 		linelen = (image->width / 8) + (image->width % 8 ? 1 : 0);
1311 		for (y = 0; y < image->height; y++) {
1312 			spixel = line;
1313 			mask = 0x80;
1314 			if (new_image->pixlen == 3)	/* most common */
1315 				for (x = 0; x < image->width; x++) {
1316 					valToMem((mask & *spixel ? fg : bg), dpixel, 3);
1317 					mask >>= 1;
1318 					if (!mask) {
1319 						mask = 0x80;
1320 						spixel++;
1321 					}
1322 					dpixel += 3;
1323 			} else	/* less common */
1324 				for (x = 0; x < image->width; x++) {
1325 					valToMem((mask & *spixel ? fg : bg), dpixel, new_image->pixlen);
1326 					mask >>= 1;
1327 					if (!mask) {
1328 						mask = 0x80;
1329 						spixel++;
1330 					}
1331 					dpixel += new_image->pixlen;
1332 				}
1333 			line += linelen;
1334 		}
1335 		new_image->gamma = 1.0;		/* will be linear now */
1336 		break;
1337 	case IRGB:
1338 		spixel = image->data;
1339 		dpixel = new_image->data;
1340 		if (image->pixlen == 1 && new_image->pixlen == 3)	/* most common */
1341 			for (y = 0; y < image->height; y++)
1342 				for (x = 0; x < image->width; x++) {
1343 					register unsigned long temp;
1344 					temp = memToVal(spixel, 1);
1345 					temp = RGB_TO_TRUE(image->rgb.red[temp],
1346 						  image->rgb.green[temp],
1347 						  image->rgb.blue[temp]);
1348 					valToMem(temp, dpixel, 3);
1349 					spixel += 1;
1350 					dpixel += 3;
1351 		} else		/* less common */
1352 			for (y = 0; y < image->height; y++)
1353 				for (x = 0; x < image->width; x++) {
1354 					register unsigned long temp;
1355 					temp = memToVal(spixel, image->pixlen);
1356 					temp = RGB_TO_TRUE(image->rgb.red[temp],
1357 						  image->rgb.green[temp],
1358 						  image->rgb.blue[temp]);
1359 					valToMem(temp, dpixel, new_image->pixlen);
1360 					spixel += image->pixlen;
1361 					dpixel += new_image->pixlen;
1362 				}
1363 		break;
1364 	}
1365 	return (new_image);
1366 }
1367 
1368 /* expand a bitmap image into an indexed RGB pixmap */
expandbittoirgb(Image * image,int depth)1369 Image *expandbittoirgb(Image *image, int depth)
1370 {
1371 	Image *new_image;
1372 	int x, y;
1373 	byte *spixel, *dpixel, *line;
1374 	unsigned int linelen;
1375 	byte mask;
1376 
1377 	CURRFUNC("expandbittoirgb");
1378 	if (!BITMAPP(image))
1379 		return (image);
1380 
1381 	new_image = newRGBImage(image->width, image->height, depth);
1382 	new_image->title = dupString(image->title);
1383 	new_image->gamma = image->gamma;
1384 
1385 	new_image->rgb.red[0] = image->rgb.red[0];
1386 	new_image->rgb.green[0] = image->rgb.green[0];
1387 	new_image->rgb.blue[0] = image->rgb.blue[0];
1388 	new_image->rgb.red[1] = image->rgb.red[1];
1389 	new_image->rgb.green[1] = image->rgb.green[1];
1390 	new_image->rgb.blue[1] = image->rgb.blue[1];
1391 	new_image->rgb.used = 2;
1392 	line = image->data;
1393 	dpixel = new_image->data;
1394 	linelen = (image->width / 8) + (image->width % 8 ? 1 : 0);
1395 	for (y = 0; y < image->height; y++) {
1396 		spixel = line;
1397 		mask = 0x80;
1398 		if (new_image->pixlen == 1)	/* most common */
1399 			for (x = 0; x < image->width; x++) {
1400 				valToMem((mask & *spixel ? 1 : 0), dpixel, 1);
1401 				mask >>= 1;
1402 				if (!mask) {
1403 					mask = 0x80;
1404 					spixel++;
1405 				}
1406 				dpixel += 1;
1407 		} else		/* less common */
1408 			for (x = 0; x < image->width; x++) {
1409 				valToMem((mask & *spixel ? 1 : 0), dpixel, new_image->pixlen);
1410 				mask >>= 1;
1411 				if (!mask) {
1412 					mask = 0x80;
1413 					spixel++;
1414 				}
1415 				dpixel += new_image->pixlen;
1416 			}
1417 		line += linelen;
1418 	}
1419 	return (new_image);
1420 }
1421 
1422 /* expand an indexed RGB pixmaps depth */
expandirgbdepth(Image * image,int depth)1423 Image *expandirgbdepth(Image *image, int depth)
1424 {
1425 	Image *new_image;
1426 	int i, x, y;
1427 	byte *spixel, *dpixel;
1428 
1429 	CURRFUNC("expandirgb");
1430 	if (TRUEP(image))
1431 		return (image);
1432 	if (BITMAPP(image))
1433 		return expandbittoirgb(image, depth);
1434 
1435 	if (((depth + 7) / 8) == image->pixlen) {	/* no change in storage */
1436 		image->depth = depth;
1437 		resizeRGBMapData(&(image->rgb), depthToColors(depth));
1438 	}
1439 	/* change in storage size */
1440 	new_image = newRGBImage(image->width, image->height, depth);
1441 	new_image->title = dupString(image->title);
1442 	new_image->gamma = image->gamma;
1443 
1444 	/* copy the color map */
1445 	for (i = 0; i < image->rgb.used; i++) {
1446 		new_image->rgb.red[i] = image->rgb.red[i];
1447 		new_image->rgb.green[i] = image->rgb.green[i];
1448 		new_image->rgb.blue[i] = image->rgb.blue[i];
1449 	}
1450 	new_image->rgb.used = image->rgb.used;
1451 
1452 	spixel = image->data;
1453 	dpixel = new_image->data;
1454 	if (image->pixlen == 1 && new_image->pixlen == 2)	/* most common */
1455 		for (y = 0; y < image->height; y++)
1456 			for (x = 0; x < image->width; x++) {
1457 				*dpixel = 0;	/* sign extend bigendian 8bit to 16 bit */
1458 				*(dpixel + 1) = *spixel;
1459 				spixel += 1;
1460 				dpixel += 2;
1461 	} else			/* less common */
1462 		for (y = 0; y < image->height; y++)
1463 			for (x = 0; x < image->width; x++) {
1464 				register unsigned long temp;
1465 				temp = memToVal(spixel, image->pixlen);
1466 				valToMem(temp, dpixel, new_image->pixlen);
1467 				spixel += image->pixlen;
1468 				dpixel += new_image->pixlen;
1469 			}
1470 	return (new_image);
1471 }
1472