hide or rename any symbols that are likely to conflict with other
[imager.git] / quant.c
CommitLineData
02d1d628
AMH
1/* quant.c - provides general image quantization
2 currently only used by gif.c, but maybe we'll support producing
3 8-bit (or bigger indexed) png files at some point
4*/
5#include "image.h"
6
7static void makemap_addi(i_quantize *, i_img **imgs, int count);
97c4effc 8static void makemap_mediancut(i_quantize *, i_img **imgs, int count);
02d1d628
AMH
9
10static
11void
12setcol(i_color *cl,unsigned char r,unsigned char g,unsigned char b,unsigned char a) {
13 cl->rgba.r=r;
14 cl->rgba.g=g;
15 cl->rgba.b=b;
16 cl->rgba.a=a;
17}
18
19
20
21/* make a colour map overwrites mc_existing/mc_count in quant Note
22 that i_makemap will be called once for each image if mc_perimage is
23 set and the format support multiple colour maps per image.
24
25 This means we don't need any special processing at this level to
26 handle multiple colour maps.
27*/
28
29void
30quant_makemap(i_quantize *quant, i_img **imgs, int count) {
97c4effc
TC
31
32 if (quant->translate == pt_giflib) {
33 /* giflib does it's own color table generation */
34 /* previously we used giflib's quantizer, but it didn't handle multiple
35 images, which made it hard to build a global color map
36 We've implemented our own median cut code so we can ignore
37 the giflib version */
38 makemap_mediancut(quant, imgs, count);
02d1d628 39 return;
97c4effc
TC
40 }
41
02d1d628
AMH
42 switch (quant->make_colors & mc_mask) {
43 case mc_none:
44 /* use user's specified map */
45 break;
46 case mc_web_map:
47 {
48 int r, g, b;
49 int i = 0;
50 for (r = 0; r < 256; r+=0x33)
51 for (g = 0; g < 256; g+=0x33)
52 for (b = 0; b < 256; b += 0x33)
53 setcol(quant->mc_colors+i++, r, g, b, 0);
54 quant->mc_count = i;
55 }
56 break;
57
97c4effc
TC
58 case mc_median_cut:
59 makemap_mediancut(quant, imgs, count);
60 break;
61
02d1d628
AMH
62 case mc_addi:
63 default:
64 makemap_addi(quant, imgs, count);
65 break;
66 }
67}
68
69#ifdef HAVE_LIBGIF
70static void translate_giflib(i_quantize *, i_img *, i_palidx *);
71#endif
72static void translate_closest(i_quantize *, i_img *, i_palidx *);
73static void translate_errdiff(i_quantize *, i_img *, i_palidx *);
74static void translate_addi(i_quantize *, i_img *, i_palidx *);
75
76/* Quantize the image given the palette in quant.
77
78 The giflib quantizer ignores the palette.
79*/
80i_palidx *quant_translate(i_quantize *quant, i_img *img) {
2ff8ed30 81 i_palidx *result;
a73aeb5f 82 mm_log((1, "quant_translate(quant %p, img %p)\n", quant, img));
2ff8ed30 83
a73aeb5f 84 result = mymalloc(img->xsize * img->ysize);
2ff8ed30 85
02d1d628 86 switch (quant->translate) {
02d1d628 87 case pt_closest:
97c4effc 88 case pt_giflib:
02d1d628
AMH
89 translate_closest(quant, img, result);
90 break;
a73aeb5f 91
02d1d628
AMH
92 case pt_errdiff:
93 translate_errdiff(quant, img, result);
94 break;
a73aeb5f 95
02d1d628
AMH
96 case pt_perturb:
97 default:
98 translate_addi(quant, img, result);
99 break;
100 }
a73aeb5f 101
02d1d628
AMH
102 return result;
103}
104
18accb2a
TC
105#ifdef HAVE_LIBGIF_THIS_NOT_USED
106
02d1d628
AMH
107#include "gif_lib.h"
108
109#define GET_RGB(im, x, y, ri, gi, bi, col) \
110 i_gpix((im),(x),(y),&(col)); (ri)=(col).rgb.r; \
111 if((im)->channels==3) { (bi)=(col).rgb.b; (gi)=(col).rgb.g; }
112
113static int
114quant_replicate(i_img *im, i_palidx *output, i_quantize *quant);
115
18accb2a 116
02d1d628
AMH
117/* Use the gif_lib quantization functions to quantize the image */
118static void translate_giflib(i_quantize *quant, i_img *img, i_palidx *out) {
119 int x,y,ColorMapSize,colours_in;
120 unsigned long Size;
121 int i;
122
123 GifByteType *RedBuffer = NULL, *GreenBuffer = NULL, *BlueBuffer = NULL;
124 GifByteType *RedP, *GreenP, *BlueP;
125 ColorMapObject *OutputColorMap = NULL;
126
127 i_color col;
128
7fd765fe 129 mm_log((1,"translate_giflib(quant %p, img %p, out %p)\n", quant, img, out));
02d1d628
AMH
130
131 /*if (!(im->channels==1 || im->channels==3)) { fprintf(stderr,"Unable to write gif, improper colorspace.\n"); exit(3); }*/
132
133 ColorMapSize = quant->mc_size;
134
135 Size = ((long) img->xsize) * img->ysize * sizeof(GifByteType);
136
137
138 if ((OutputColorMap = MakeMapObject(ColorMapSize, NULL)) == NULL)
139 m_fatal(0,"Failed to allocate memory for Output colormap.");
140 /* if ((OutputBuffer = (GifByteType *) mymalloc(im->xsize * im->ysize * sizeof(GifByteType))) == NULL)
141 m_fatal(0,"Failed to allocate memory for output buffer.");*/
142
143 /* ******************************************************* */
144 /* count the number of colours in the image */
145 colours_in=i_count_colors(img, OutputColorMap->ColorCount);
146
147 if(colours_in != -1) { /* less then the number wanted */
148 /* so we copy them over as-is */
149 mm_log((2,"image has %d colours, which fits in %d. Copying\n",
150 colours_in,ColorMapSize));
151 quant_replicate(img, out, quant);
152 /* saves the colors, so don't fall through */
153 return;
154 } else {
155
156 mm_log((2,"image has %d colours, more then %d. Quantizing\n",colours_in,ColorMapSize));
157
158 if (img->channels >= 3) {
159 if ((RedBuffer = (GifByteType *) mymalloc((unsigned int) Size)) == NULL) {
160 m_fatal(0,"Failed to allocate memory required, aborted.");
161 return;
162 }
163 if ((GreenBuffer = (GifByteType *) mymalloc((unsigned int) Size)) == NULL) {
164 m_fatal(0,"Failed to allocate memory required, aborted.");
862b614c 165 myfree(RedBuffer);
02d1d628
AMH
166 return;
167 }
168
169 if ((BlueBuffer = (GifByteType *) mymalloc((unsigned int) Size)) == NULL) {
170 m_fatal(0,"Failed to allocate memory required, aborted.");
862b614c
AMH
171 myfree(RedBuffer);
172 myfree(GreenBuffer);
02d1d628
AMH
173 return;
174 }
175
176 RedP = RedBuffer;
177 GreenP = GreenBuffer;
178 BlueP = BlueBuffer;
179
180 for (y=0; y< img->ysize; y++) for (x=0; x < img->xsize; x++) {
181 i_gpix(img,x,y,&col);
182 *RedP++ = col.rgb.r;
183 *GreenP++ = col.rgb.g;
184 *BlueP++ = col.rgb.b;
185 }
186
187 } else {
188
189 if ((RedBuffer = (GifByteType *) mymalloc((unsigned int) Size))==NULL) {
190 m_fatal(0,"Failed to allocate memory required, aborted.");
191 return;
192 }
193
194 GreenBuffer=BlueBuffer=RedBuffer;
195 RedP = RedBuffer;
196 for (y=0; y< img->ysize; y++) for (x=0; x < img->xsize; x++) {
197 i_gpix(img,x,y,&col);
198 *RedP++ = col.rgb.r;
199 }
200 }
201
202 if (QuantizeBuffer(img->xsize, img->ysize, &ColorMapSize, RedBuffer, GreenBuffer, BlueBuffer,
203 out, OutputColorMap->Colors) == GIF_ERROR) {
204 mm_log((1,"Error in QuantizeBuffer, unable to write image.\n"));
205 }
206 }
207
862b614c
AMH
208 myfree(RedBuffer);
209 if (img->channels == 3) { myfree(GreenBuffer); myfree(BlueBuffer); }
02d1d628
AMH
210
211 /* copy over the color map */
212 for (i = 0; i < ColorMapSize; ++i) {
213 quant->mc_colors[i].rgb.r = OutputColorMap->Colors[i].Red;
214 quant->mc_colors[i].rgb.g = OutputColorMap->Colors[i].Green;
215 quant->mc_colors[i].rgb.b = OutputColorMap->Colors[i].Blue;
216 }
217 quant->mc_count = ColorMapSize;
218}
219
220static
221int
222quant_replicate(i_img *im, GifByteType *output, i_quantize *quant) {
223 int x, y, alloced, r, g=0, b=0, idx ;
224 i_color col;
225
226 alloced=0;
227 for(y=0; y<im->ysize; y++) {
228 for(x=0; x<im->xsize; x++) {
229
230 GET_RGB(im, x,y, r,g,b, col);
231
232 for(idx=0; idx<alloced; idx++) { /* linear search for an index */
233 if(quant->mc_colors[idx].rgb.r==r &&
234 quant->mc_colors[idx].rgb.g==g &&
235 quant->mc_colors[idx].rgb.b==b) {
236 break;
237 }
238 }
239
240 if(idx >= alloced) { /* if we haven't already, we */
241 idx=alloced++; /* add the colour to the map */
242 if(quant->mc_size < alloced) {
243 mm_log((1,"Tried to allocate more then %d colours.\n",
244 quant->mc_size));
245 return 0;
246 }
247 quant->mc_colors[idx].rgb.r=r;
248 quant->mc_colors[idx].rgb.g=g;
249 quant->mc_colors[idx].rgb.b=b;
250 }
251 *output=idx; /* fill output buffer */
252 output++; /* with colour indexes */
253 }
254 }
255 quant->mc_count = alloced;
256 return 1;
257}
258
259#endif
260
261static void translate_closest(i_quantize *quant, i_img *img, i_palidx *out) {
262 quant->perturb = 0;
263 translate_addi(quant, img, out);
264}
265
266#define PWR2(x) ((x)*(x))
267
268typedef int (*cmpfunc)(const void*, const void*);
269
270typedef struct {
271 unsigned char r,g,b;
36e67d0b
TC
272 char fixed;
273 char used;
02d1d628
AMH
274 int dr,dg,db;
275 int cdist;
276 int mcount;
277} cvec;
278
279typedef struct {
280 int cnt;
281 int vec[256];
282} hashbox;
283
284typedef struct {
285 int boxnum;
286 int pixcnt;
287 int cand;
288 int pdc;
289} pbox;
290
18accb2a 291static void prescan(i_img **im,int count, int cnum, cvec *clr, i_sample_t *line);
02d1d628
AMH
292static void reorder(pbox prescan[512]);
293static int pboxcmp(const pbox *a,const pbox *b);
294static void boxcenter(int box,cvec *cv);
295static float frandn(void);
296static void boxrand(int box,cvec *cv);
297static void bbox(int box,int *r0,int *r1,int *g0,int *g1,int *b0,int *b1);
298static void cr_hashindex(cvec clr[256],int cnum,hashbox hb[512]);
299static int mindist(int boxnum,cvec *cv);
300static int maxdist(int boxnum,cvec *cv);
301
302/* Some of the simpler functions are kept here to aid the compiler -
303 maybe some of them will be inlined. */
304
305static int
306pixbox(i_color *ic) { return ((ic->channel[0] & 224)<<1)+ ((ic->channel[1]&224)>>2) + ((ic->channel[2] &224) >> 5); }
307
18accb2a
TC
308static int
309pixbox_ch(i_sample_t *chans) { return ((chans[0] & 224)<<1)+ ((chans[1]&224)>>2) + ((chans[2] &224) >> 5); }
310
02d1d628
AMH
311static unsigned char
312g_sat(int in) {
313 if (in>255) { return 255; }
314 else if (in>0) return in;
315 return 0;
316}
317
318static
319float
320frand(void) {
321 return rand()/(RAND_MAX+1.0);
322}
323
324static
325int
326eucl_d(cvec* cv,i_color *cl) { return PWR2(cv->r-cl->channel[0])+PWR2(cv->g-cl->channel[1])+PWR2(cv->b-cl->channel[2]); }
327
18accb2a
TC
328static
329int
330eucl_d_ch(cvec* cv,i_sample_t *chans) {
331 return PWR2(cv->r - chans[0]) + PWR2(cv->g - chans[1])
332 + PWR2(cv->b - chans[2]);
333}
334
02d1d628
AMH
335static
336int
337ceucl_d(i_color *c1, i_color *c2) { return PWR2(c1->channel[0]-c2->channel[0])+PWR2(c1->channel[1]-c2->channel[1])+PWR2(c1->channel[2]-c2->channel[2]); }
338
18accb2a
TC
339static const int
340gray_samples[] = { 0, 0, 0 };
341
02d1d628
AMH
342/*
343
344This quantization algorithm and implementation routines are by Arnar
345M. Hrafnkelson. In case any new ideas are here they are mine since
346this was written from scratch.
347
348The algorithm uses local means in the following way:
349
350 For each point in the colormap we find which image points
351 have that point as it's closest point. We calculate the mean
352 of those points and in the next iteration it will be the new
353 entry in the colormap.
354
355In order to speed this process up (i.e. nearest neighbor problem) We
356divied the r,g,b space up in equally large 512 boxes. The boxes are
357numbered from 0 to 511. Their numbering is so that for a given vector
358it is known that it belongs to the box who is formed by concatenating the
3593 most significant bits from each component of the RGB triplet.
360
361For each box we find the list of points from the colormap who might be
362closest to any given point within the box. The exact solution
363involves finding the Voronoi map (or the dual the Delauny
364triangulation) and has many issues including numerical stability.
365
366So we use this approximation:
367
3681. Find which point has the shortest maximum distance to the box.
3692. Find all points that have a shorter minimum distance than that to the box
370
371This is a very simple task and is not computationally heavy if one
372takes into account that the minimum distances from a pixel to a box is
373always found by checking if it's inside the box or is closest to some
374side or a corner. Finding the maximum distance is also either a side
375or a corner.
376
377This approach results 2-3 times more than the actual points needed but
378is still a good gain over the complete space. Usually when one has a
379256 Colorcolor map a search over 30 is often obtained.
380
381A bit of an enhancement to this approach is to keep a seperate list
382for each side of the cube, but this will require even more memory.
383
384 Arnar M. Hrafnkelsson (addi@umich.edu);
385
386*/
387/*
388 Extracted from gifquant.c, removed dependencies on gif_lib,
389 and added support for multiple images.
390 starting from 1nov2000 by TonyC <tony@develop-help.com>.
391
392*/
393
394static void
395makemap_addi(i_quantize *quant, i_img **imgs, int count) {
396 cvec *clr;
36e67d0b 397 int cnum, i, x, y, bst_idx=0, ld, cd, iter, currhb, img_num;
18accb2a 398 i_sample_t *val;
02d1d628 399 float dlt, accerr;
9cfd5724 400 hashbox *hb;
18accb2a
TC
401 i_mempool mp;
402 int maxwidth = 0;
403 i_sample_t *line;
404 const int *sample_indices;
02d1d628 405
18accb2a
TC
406 i_mempool_init(&mp);
407
408 clr = i_mempool_alloc(&mp, sizeof(cvec) * quant->mc_size);
409 hb = i_mempool_alloc(&mp, sizeof(hashbox) * 512);
02d1d628
AMH
410 for (i=0; i < quant->mc_count; ++i) {
411 clr[i].r = quant->mc_colors[i].rgb.r;
412 clr[i].g = quant->mc_colors[i].rgb.g;
413 clr[i].b = quant->mc_colors[i].rgb.b;
36e67d0b
TC
414 clr[i].fixed = 1;
415 clr[i].mcount = 0;
02d1d628
AMH
416 }
417 /* mymalloc doesn't clear memory, so I think we need this */
418 for (; i < quant->mc_size; ++i) {
36e67d0b
TC
419 clr[i].fixed = 0;
420 clr[i].mcount = 0;
02d1d628
AMH
421 }
422 cnum = quant->mc_size;
423 dlt = 1;
424
18accb2a
TC
425 for (img_num = 0; img_num < count; ++img_num) {
426 if (imgs[img_num]->xsize > maxwidth)
427 maxwidth = imgs[img_num]->xsize;
428 }
429 line = i_mempool_alloc(&mp, 3 * maxwidth * sizeof(*line));
430
431 prescan(imgs, count, cnum, clr, line);
02d1d628
AMH
432 cr_hashindex(clr, cnum, hb);
433
434 for(iter=0;iter<3;iter++) {
435 accerr=0.0;
436
36e67d0b
TC
437 for (img_num = 0; img_num < count; ++img_num) {
438 i_img *im = imgs[img_num];
18accb2a
TC
439 sample_indices = im->channels >= 3 ? NULL : gray_samples;
440 for(y=0;y<im->ysize;y++) {
441 i_gsamp(im, 0, im->xsize, y, line, sample_indices, 3);
442 val = line;
443 for(x=0;x<im->xsize;x++) {
444 ld=196608;
445 /*i_gpix(im,x,y,&val);*/
446 currhb=pixbox_ch(val);
447 /* printf("box = %d \n",currhb); */
448 for(i=0;i<hb[currhb].cnt;i++) {
449 /* printf("comparing: pix (%d,%d,%d) vec (%d,%d,%d)\n",val.channel[0],val.channel[1],val.channel[2],clr[hb[currhb].vec[i]].r,clr[hb[currhb].vec[i]].g,clr[hb[currhb].vec[i]].b); */
450
451 cd=eucl_d_ch(&clr[hb[currhb].vec[i]],val);
452 if (cd<ld) {
453 ld=cd; /* shortest distance yet */
454 bst_idx=hb[currhb].vec[i]; /* index of closest vector yet */
455 }
456 }
457
458 clr[bst_idx].mcount++;
459 accerr+=(ld);
460 clr[bst_idx].dr+=val[0];
461 clr[bst_idx].dg+=val[1];
462 clr[bst_idx].db+=val[2];
463
464 val += 3; /* next 3 samples (next pixel) */
465 }
02d1d628
AMH
466 }
467 }
18accb2a
TC
468
469 for(i=0;i<cnum;i++)
470 if (clr[i].mcount) {
471 clr[i].dr/=clr[i].mcount;
472 clr[i].dg/=clr[i].mcount;
473 clr[i].db/=clr[i].mcount;
474 }
475
02d1d628 476 /* for(i=0;i<cnum;i++) printf("vec(%d)=(%d,%d,%d) dest=(%d,%d,%d) matchcount=%d\n",
18accb2a
TC
477 i,clr[i].r,clr[i].g,clr[i].b,clr[i].dr,clr[i].dg,clr[i].db,clr[i].mcount); */
478
02d1d628 479 /* printf("total error: %.2f\n",sqrt(accerr)); */
18accb2a 480
02d1d628 481 for(i=0;i<cnum;i++) {
36e67d0b 482 if (clr[i].fixed) continue; /* skip reserved colors */
18accb2a 483
02d1d628 484 if (clr[i].mcount) {
18accb2a
TC
485 clr[i].used = 1;
486 clr[i].r=clr[i].r*(1-dlt)+dlt*clr[i].dr;
487 clr[i].g=clr[i].g*(1-dlt)+dlt*clr[i].dg;
488 clr[i].b=clr[i].b*(1-dlt)+dlt*clr[i].db;
02d1d628 489 } else {
18accb2a
TC
490 /* let's try something else */
491 clr[i].used = 0;
492 clr[i].r=rand();
493 clr[i].g=rand();
494 clr[i].b=rand();
02d1d628 495 }
18accb2a 496
02d1d628
AMH
497 clr[i].dr=0;
498 clr[i].dg=0;
499 clr[i].db=0;
500 clr[i].mcount=0;
501 }
502 cr_hashindex(clr,cnum,hb);
503 }
504
505
506#ifdef NOTEF
507 for(i=0;i<cnum;i++) {
508 cd=eucl_d(&clr[i],&val);
509 if (cd<ld) {
510 ld=cd;
511 bst_idx=i;
512 }
513 }
514#endif
515
36e67d0b
TC
516 /* if defined, we only include colours with an mcount or that were
517 supplied in the fixed palette, giving us a smaller output palette */
518#define ONLY_USE_USED
519#ifdef ONLY_USE_USED
520 /* transfer the colors back */
521 quant->mc_count = 0;
522 for (i = 0; i < cnum; ++i) {
523 if (clr[i].fixed || clr[i].used) {
524 /*printf("Adding %d (%d,%d,%d)\n", i, clr[i].r, clr[i].g, clr[i].b);*/
525 quant->mc_colors[quant->mc_count].rgb.r = clr[i].r;
526 quant->mc_colors[quant->mc_count].rgb.g = clr[i].g;
527 quant->mc_colors[quant->mc_count].rgb.b = clr[i].b;
528 ++quant->mc_count;
529 }
530 }
531#else
02d1d628
AMH
532 /* transfer the colors back */
533 for (i = 0; i < cnum; ++i) {
534 quant->mc_colors[i].rgb.r = clr[i].r;
535 quant->mc_colors[i].rgb.g = clr[i].g;
536 quant->mc_colors[i].rgb.b = clr[i].b;
537 }
538 quant->mc_count = cnum;
36e67d0b 539#endif
02d1d628 540
18accb2a 541 i_mempool_destroy(&mp);
02d1d628
AMH
542}
543
97c4effc
TC
544typedef struct {
545 i_sample_t rgb[3];
546 int count;
547} quant_color_entry;
548
549#define MEDIAN_CUT_COLORS 32768
550
551#define MED_CUT_INDEX(c) ((((c).rgb.r & 0xF8) << 7) | \
552 (((c).rgb.g & 0xF8) << 2) | (((c).rgb.b & 0xF8) >> 3))
553
18accb2a
TC
554#define MED_CUT_GRAY_INDEX(c) ((((c).rgb.r & 0xF8) << 7) | \
555 (((c).rgb.r & 0xF8) << 2) | (((c).rgb.r & 0xF8) >> 3))
556
97c4effc
TC
557/* scale these to cover the whole range */
558#define MED_CUT_RED(index) ((((index) & 0x7C00) >> 10) * 255 / 31)
559#define MED_CUT_GREEN(index) ((((index) & 0x3E0) >> 5) * 255 / 31)
560#define MED_CUT_BLUE(index) (((index) & 0x1F) * 255 / 31)
561
562typedef struct {
563 i_sample_t min[3]; /* minimum for each channel */
564 i_sample_t max[3]; /* maximum for each channel */
565 i_sample_t width[3]; /* width for each channel */
566 int start, size; /* beginning and size of the partition */
567 int pixels; /* number of pixels represented by this partition */
568} medcut_partition;
569
570/*
571=item calc_part(part, colors)
572
573Calculates the new color limits for the given partition.
574
575Giflib assumes that the limits for the non-split channels stay the
576same, but this strikes me as incorrect, especially if the colors tend
577to be color ramps.
578
579Of course this could be optimized by not recalculating the channel we
580just sorted on, but it's not worth the effort right now.
581
582=cut
583*/
584static void calc_part(medcut_partition *part, quant_color_entry *colors) {
585 int i, ch;
586
587 for (ch = 0; ch < 3; ++ch) {
588 part->min[ch] = 255;
589 part->max[ch] = 0;
590 }
591 for (i = part->start; i < part->start + part->size; ++i) {
592 for (ch = 0; ch < 3; ++ch) {
593 if (part->min[ch] > colors[i].rgb[ch])
594 part->min[ch] = colors[i].rgb[ch];
595 if (part->max[ch] < colors[i].rgb[ch])
596 part->max[ch] = colors[i].rgb[ch];
597 }
598 }
599 for (ch = 0; ch < 3; ++ch) {
600 part->width[ch] = part->max[ch] - part->min[ch];
601 }
602}
603
604/* simple functions to sort by each channel - we could use a global, but
605 that would be bad */
606
607static int
608color_sort_red(void const *left, void const *right) {
609 return ((quant_color_entry *)left)->rgb[0] - ((quant_color_entry *)right)->rgb[0];
610}
611
612static int
613color_sort_green(void const *left, void const *right) {
614 return ((quant_color_entry *)left)->rgb[1] - ((quant_color_entry *)right)->rgb[1];
615}
616
617static int
618color_sort_blue(void const *left, void const *right) {
619 return ((quant_color_entry *)left)->rgb[2] - ((quant_color_entry *)right)->rgb[2];
620}
621
622static int (*sorters[])(void const *, void const *) =
623{
624 color_sort_red,
625 color_sort_green,
626 color_sort_blue,
627};
628
629static void
630makemap_mediancut(i_quantize *quant, i_img **imgs, int count) {
631 quant_color_entry *colors;
632 i_mempool mp;
633 int imgn, x, y, i, ch;
634 int max_width;
635 i_color *line;
636 int color_count;
637 int total_pixels;
638 medcut_partition *parts;
639 int part_num;
640 int in, out;
18accb2a
TC
641 /* number of channels we search for the best channel to partition
642 this isn't terribly efficient, but it should work */
643 int chan_count;
97c4effc
TC
644
645 /*printf("images %d pal size %d\n", count, quant->mc_size);*/
646
647 i_mempool_init(&mp);
648
649 colors = i_mempool_alloc(&mp, sizeof(*colors) * MEDIAN_CUT_COLORS);
650 for (i = 0; i < MEDIAN_CUT_COLORS; ++i) {
651 colors[i].rgb[0] = MED_CUT_RED(i);
652 colors[i].rgb[1] = MED_CUT_GREEN(i);
653 colors[i].rgb[2] = MED_CUT_BLUE(i);
654 colors[i].count = 0;
655 }
656
657 max_width = -1;
658 for (imgn = 0; imgn < count; ++imgn) {
659 if (imgs[imgn]->xsize > max_width)
660 max_width = imgs[imgn]->xsize;
661 }
662 line = i_mempool_alloc(&mp, sizeof(i_color) * max_width);
663
664 /* build the stats */
665 total_pixels = 0;
18accb2a 666 chan_count = 1; /* assume we just have grayscale */
97c4effc
TC
667 for (imgn = 0; imgn < count; ++imgn) {
668 total_pixels += imgs[imgn]->xsize * imgs[imgn]->ysize;
669 for (y = 0; y < imgs[imgn]->ysize; ++y) {
670 i_glin(imgs[imgn], 0, imgs[imgn]->xsize, y, line);
18accb2a
TC
671 if (imgs[imgn]->channels > 2) {
672 chan_count = 3;
673 for (x = 0; x < imgs[imgn]->xsize; ++x) {
674 ++colors[MED_CUT_INDEX(line[x])].count;
675 }
676 }
677 else {
678 /* a gray-scale image, just use the first channel */
679 for (x = 0; x < imgs[imgn]->xsize; ++x) {
680 ++colors[MED_CUT_GRAY_INDEX(line[x])].count;
681 }
97c4effc
TC
682 }
683 }
684 }
685
686 /* eliminate the empty colors */
687 out = 0;
688 for (in = 0; in < MEDIAN_CUT_COLORS; ++in) {
689 if (colors[in].count) {
690 colors[out++] = colors[in];
691 }
692 }
693 /*printf("out %d\n", out);
694
695 for (i = 0; i < out; ++i) {
696 if (colors[i].count) {
697 printf("%d: (%d,%d,%d) -> %d\n", i, colors[i].rgb[0], colors[i].rgb[1],
698 colors[i].rgb[2], colors[i].count);
699 }
700 }*/
701
702 if (out < quant->mc_size) {
703 /* just copy them into the color table */
704 for (i = 0; i < out; ++i) {
705 for (ch = 0; ch < 3; ++ch) {
706 quant->mc_colors[i].channel[ch] = colors[i].rgb[ch];
707 }
708 }
709 quant->mc_count = out;
710 }
711 else {
712 /* build the starting partition */
713 parts = i_mempool_alloc(&mp, sizeof(*parts) * quant->mc_size);
714 parts[0].start = 0;
715 parts[0].size = out;
716 parts[0].pixels = total_pixels;
717 calc_part(parts, colors);
718 color_count = 1;
719
720 while (color_count < quant->mc_size) {
721 int max_index, max_ch; /* index/channel with biggest spread */
722 int max_size;
723 medcut_partition *workpart;
724 int cum_total;
725 int half;
726
727 /* find the partition with the most biggest span with more than
728 one color */
729 max_size = -1;
730 for (i = 0; i < color_count; ++i) {
18accb2a 731 for (ch = 0; ch < chan_count; ++ch) {
97c4effc
TC
732 if (parts[i].width[ch] > max_size
733 && parts[i].size > 1) {
734 max_index = i;
735 max_ch = ch;
736 max_size = parts[i].width[ch];
737 }
738 }
739 }
740
741 /* nothing else we can split */
742 if (max_size == -1)
743 break;
744
745 workpart = parts+max_index;
746 /*printf("splitting partition %d (pixels %ld, start %d, size %d)\n", max_index, workpart->pixels, workpart->start, workpart->size);*/
747 qsort(colors + workpart->start, workpart->size, sizeof(*colors),
748 sorters[max_ch]);
749
750 /* find the median or something like it we need to make sure both
751 sides of the split have at least one color in them, so we don't
752 test at the first or last entry */
753 i = workpart->start;
754 cum_total = colors[i].count;
755 ++i;
756 half = workpart->pixels / 2;
757 while (i < workpart->start + workpart->size - 1
758 && cum_total < half) {
759 cum_total += colors[i++].count;
760 }
761 /*printf("Split at %d to make %d (half %ld, cumtotal %ld)\n", i, color_count, half, cum_total);*/
762
763 /* found the spot to split */
764 parts[color_count].start = i;
765 parts[color_count].size = workpart->start + workpart->size - i;
766 workpart->size = i - workpart->start;
767 parts[color_count].pixels = workpart->pixels - cum_total;
768 workpart->pixels = cum_total;
769
770 /* recalculate the limits */
771 calc_part(workpart, colors);
772 calc_part(parts+color_count, colors);
773 ++color_count;
774 }
775
776 /* fill in the color table - since we could still have partitions
777 that have more than one color, we need to average the colors */
778 for (part_num = 0; part_num < color_count; ++part_num) {
779 long sums[3];
780 medcut_partition *workpart;
781
782 workpart = parts+part_num;
783 for (ch = 0; ch < 3; ++ch)
784 sums[ch] = 0;
785
786 for (i = workpart->start; i < workpart->start + workpart->size; ++i) {
787 for (ch = 0; ch < 3; ++ch) {
788 sums[ch] += colors[i].rgb[ch] * colors[i].count;
789 }
790 }
791 for (ch = 0; ch < 3; ++ch) {
792 quant->mc_colors[part_num].channel[ch] = sums[ch] / workpart->pixels;
793 }
794 }
795 quant->mc_count = color_count;
796 }
797 /*printf("out %d colors\n", quant->mc_count);*/
798 i_mempool_destroy(&mp);
799}
800
02d1d628
AMH
801#define pboxjump 32
802
803/* Define one of the following 4 symbols to choose a colour search method
804 The idea is to try these out, including benchmarking, to see which
805 is fastest in a good spread of circumstances.
806 I'd expect IM_CFLINSEARCH to be fastest for very small palettes, and
807 IM_CFHASHBOX for large images with large palettes.
808
809 Some other possibilities include:
810 - search over entries sorted by luminance
811
812 Initially I was planning on testing using the macros and then
813 integrating the code directly into each function, but this means if
814 we find a bug at a late stage we will need to update N copies of
815 the same code. Also, keeping the code in the macros means that the
816 code in the translation functions is much more to the point,
817 there's no distracting colour search code to remove attention from
818 what makes _this_ translation function different. It may be
819 advisable to move the setup code into functions at some point, but
820 it should be possible to do this fairly transparently.
821
822 If IM_CF_COPTS is defined then CFLAGS must have an appropriate
823 definition.
824
825 Each option needs to define 4 macros:
826 CF_VARS - variables to define in the function
827 CF_SETUP - code to setup for the colour search, eg. allocating and
828 initializing lookup tables
829 CF_FIND - code that looks for the color in val and puts the best
830 matching index in bst_idx
831 CF_CLEANUP - code to clean up, eg. releasing memory
832*/
833#ifndef IM_CF_COPTS
834/*#define IM_CFLINSEARCH*/
835#define IM_CFHASHBOX
836/*#define IM_CFSORTCHAN*/
837/*#define IM_CFRAND2DIST*/
838#endif
839
840#ifdef IM_CFHASHBOX
841
842/* The original version I wrote for this used the sort.
843 If this is defined then we use a sort to extract the indices for
844 the hashbox */
845#define HB_SORT
846
847/* assume i is available */
9cfd5724 848#define CF_VARS hashbox *hb = mymalloc(sizeof(hashbox) * 512); \
02d1d628
AMH
849 int currhb; \
850 long ld, cd
851
852#ifdef HB_SORT
853
854static long *gdists; /* qsort is annoying */
855/* int might be smaller than long, so we need to do a real compare
856 rather than a subtraction*/
857static int distcomp(void const *a, void const *b) {
858 long ra = gdists[*(int const *)a];
859 long rb = gdists[*(int const *)b];
860 if (ra < rb)
861 return -1;
862 else if (ra > rb)
863 return 1;
864 else
865 return 0;
866}
867
868#endif
869
870/* for each hashbox build a list of colours that are in the hb or is closer
871 than other colours
872 This is pretty involved. The original gifquant generated the hashbox
873 as part of it's normal processing, but since the map generation is now
874 separated from the translation we need to do this on the spot.
875 Any optimizations, even if they don't produce perfect results would be
876 welcome.
877 */
878static void hbsetup(i_quantize *quant, hashbox *hb) {
879 long *dists, mind, maxd, cd;
880 int cr, cb, cg, hbnum, i;
881 i_color cenc;
882#ifdef HB_SORT
883 int *indices = mymalloc(quant->mc_count * sizeof(int));
884#endif
885
886 dists = mymalloc(quant->mc_count * sizeof(long));
887 for (cr = 0; cr < 8; ++cr) {
888 for (cg = 0; cg < 8; ++cg) {
889 for (cb = 0; cb < 8; ++cb) {
890 /* centre of the hashbox */
891 cenc.channel[0] = cr*pboxjump+pboxjump/2;
892 cenc.channel[1] = cg*pboxjump+pboxjump/2;
893 cenc.channel[2] = cb*pboxjump+pboxjump/2;
894 hbnum = pixbox(&cenc);
895 hb[hbnum].cnt = 0;
896 /* order indices in the order of distance from the hashbox */
897 for (i = 0; i < quant->mc_count; ++i) {
898#ifdef HB_SORT
899 indices[i] = i;
900#endif
901 dists[i] = ceucl_d(&cenc, quant->mc_colors+i);
902 }
903#ifdef HB_SORT
904 /* it should be possible to do this without a sort
905 but so far I'm too lazy */
906 gdists = dists;
907 qsort(indices, quant->mc_count, sizeof(int), distcomp);
908 /* any colors that can match are within mind+diagonal size of
909 a hashbox */
910 mind = dists[indices[0]];
911 i = 0;
912 maxd = (sqrt(mind)+pboxjump)*(sqrt(mind)+pboxjump);
913 while (i < quant->mc_count && dists[indices[i]] < maxd) {
914 hb[hbnum].vec[hb[hbnum].cnt++] = indices[i++];
915 }
916#else
917 /* work out the minimum */
918 mind = 256*256*3;
919 for (i = 0; i < quant->mc_count; ++i) {
920 if (dists[i] < mind) mind = dists[i];
921 }
922 /* transfer any colours that might be closest to a colour in
923 this hashbox */
924 maxd = (sqrt(mind)+pboxjump)*(sqrt(mind)+pboxjump);
925 for (i = 0; i < quant->mc_count; ++i) {
926 if (dists[i] < maxd)
927 hb[hbnum].vec[hb[hbnum].cnt++] = i;
928 }
929#endif
930 }
931 }
862b614c 932 }
02d1d628
AMH
933#ifdef HB_SORT
934 myfree(indices);
935#endif
936 myfree(dists) ;
937}
938#define CF_SETUP hbsetup(quant, hb)
939
940#define CF_FIND \
941 currhb = pixbox(&val); \
942 ld = 196608; \
943 for (i = 0; i < hb[currhb].cnt; ++i) { \
944 cd = ceucl_d(quant->mc_colors+hb[currhb].vec[i], &val); \
945 if (cd < ld) { ld = cd; bst_idx = hb[currhb].vec[i]; } \
946 }
947
9cfd5724 948#define CF_CLEANUP myfree(hb)
02d1d628
AMH
949
950#endif
951
952#ifdef IM_CFLINSEARCH
953/* as simple as it gets */
954#define CF_VARS long ld, cd
955#define CF_SETUP /* none needed */
956#define CF_FIND \
957 ld = 196608; \
958 for (i = 0; i < quant->mc_count; ++i) { \
959 cd = ceucl_d(quant->mc_colors+i, &val); \
960 if (cd < ld) { ld = cd; bst_idx = i; } \
961 }
962#define CF_CLEANUP
963#endif
964
965#ifdef IM_CFSORTCHAN
966static int gsortchan;
967static i_quantize *gquant;
968static int chansort(void const *a, void const *b) {
969 return gquant->mc_colors[*(int const *)a].channel[gsortchan] -
970 gquant->mc_colors[*(int const *)b].channel[gsortchan];
971}
972#define CF_VARS int *indices, sortchan, diff; \
973 long ld, cd; \
974 int vindex[256] /* where to find value i of chan */
975
976static void chansetup(i_img *img, i_quantize *quant, int *csortchan,
977 int *vindex, int **cindices) {
978 int *indices, sortchan, chan, i, chval;
979 int chanmins[MAXCHANNELS], chanmaxs[MAXCHANNELS], maxrange;
980
981 /* find the channel with the maximum range */
982 /* the maximum stddev would probably be better */
983 for (chan = 0; chan < img->channels; ++chan) {
984 chanmins[chan] = 256; chanmaxs[chan] = 0;
985 for (i = 0; i < quant->mc_count; ++i) {
986 if (quant->mc_colors[i].channel[chan] < chanmins[chan])
987 chanmins[chan] = quant->mc_colors[i].channel[chan];
988 if (quant->mc_colors[i].channel[chan] > chanmaxs[chan])
989 chanmaxs[chan] = quant->mc_colors[i].channel[chan];
990 }
991 }
992 maxrange = -1;
993 for (chan = 0; chan < img->channels; ++chan) {
994 if (chanmaxs[chan]-chanmins[chan] > maxrange) {
995 maxrange = chanmaxs[chan]-chanmins[chan];
996 sortchan = chan;
997 }
998 }
999 indices = mymalloc(quant->mc_count * sizeof(int)) ;
1000 for (i = 0; i < quant->mc_count; ++i) {
1001 indices[i] = i;
1002 }
1003 gsortchan = sortchan;
1004 gquant = quant;
1005 qsort(indices, quant->mc_count, sizeof(int), chansort) ;
1006 /* now a lookup table to find entries faster */
1007 for (chval=0, i=0; i < quant->mc_count; ++i) {
1008 while (chval < 256 &&
1009 chval < quant->mc_colors[indices[i]].channel[sortchan]) {
1010 vindex[chval++] = i;
1011 }
1012 }
1013 while (chval < 256) {
1014 vindex[chval++] = quant->mc_count-1;
1015 }
1016 *csortchan = sortchan;
1017 *cindices = indices;
1018}
1019
1020#define CF_SETUP \
1021 chansetup(img, quant, &sortchan, vindex, &indices)
1022
1023int chanfind(i_color val, i_quantize *quant, int *indices, int *vindex,
1024 int sortchan) {
1025 int i, bst_idx, diff, maxdiff;
1026 long ld, cd;
1027
1028 i = vindex[val.channel[sortchan]];
1029 bst_idx = indices[i];
1030 ld = 196608;
1031 diff = 0;
1032 maxdiff = quant->mc_count;
1033 while (diff < maxdiff) {
1034 if (i+diff < quant->mc_count) {
1035 cd = ceucl_d(&val, quant->mc_colors+indices[i+diff]);
1036 if (cd < ld) {
1037 bst_idx = indices[i+diff];
1038 ld = cd;
1039 maxdiff = sqrt(ld);
1040 }
1041 }
1042 if (i-diff >= 0) {
1043 cd = ceucl_d(&val, quant->mc_colors+indices[i-diff]);
1044 if (cd < ld) {
1045 bst_idx = indices[i-diff];
1046 ld = cd;
1047 maxdiff = sqrt(ld);
1048 }
1049 }
1050 ++diff;
1051 }
1052
1053 return bst_idx;
1054}
1055
1056#define CF_FIND \
1057 bst_idx = chanfind(val, quant, indices, vindex, sortchan)
1058
1059
1060#define CF_CLEANUP myfree(indices)
1061
1062#endif
1063
1064#ifdef IM_CFRAND2DIST
1065
1066/* This is based on a method described by Addi in the #imager channel
1067 on the 28/2/2001. I was about 1am Sydney time at the time, so I
1068 wasn't at my most cogent. Well, that's my excuse :)
1069
1070<TonyC> what I have at the moment is: hashboxes, with optimum hash box
1071filling; simple linear search; and a lookup in the widest channel
1072(currently the channel with the maximum range)
1073<Addi> There is one more way that might be simple to implement.
1074<Addi> You want to hear?
1075<TonyC> what's that?
1076<purl> somebody said that was not true
1077<Addi> For each of the colors in the palette start by creating a
1078sorted list of the form:
1079<Addi> [distance, color]
1080<Addi> Where they are sorted by distance.
1081<TonyC> distance to where?
1082<Addi> Where the elements in the lists are the distances and colors of
1083the other colors in the palette
1084<TonyC> ok
1085<Addi> So if you are at color 0
1086<Addi> ok - now to search for the closest color when you are creating
1087the final image is done like this:
1088<Addi> a) pick a random color from the palette
1089<Addi> b) calculate the distance to it
1090<Addi> c) only check the vectors that are within double the distance
1091in the list of the color you picked from the palette.
1092<Addi> Does that seem logical?
1093<Addi> Lets imagine that we only have grayscale to make an example:
1094<Addi> Our palette has 1 4 10 20 as colors.
1095<Addi> And we want to quantize the color 11
1096<Addi> lets say we picked 10 randomly
1097<Addi> the double distance is 2
1098<Addi> since abs(10-11)*2 is 2
1099<Addi> And the list at vector 10 is this:
1100<Addi> [0, 10], [6 4], [9, 1], [10, 20]
1101<Addi> so we look at the first one (but not the second one since 6 is
1102at a greater distance than 2.
1103<Addi> Any of that make sense?
1104<TonyC> yes, though are you suggesting another random jump to one of
1105the colours with the possible choices? or an exhaustive search?
1106<Addi> TonyC: It's possible to come up with a recursive/iterative
1107enhancement but this is the 'basic' version.
1108<Addi> Which would do an iterative search.
1109<Addi> You can come up with conditions where it pays to switch to a new one.
1110<Addi> And the 'random' start can be switched over to a small tree.
1111<Addi> So you would have a little index at the start.
1112<Addi> to get you into the general direction
1113<Addi> Perhaps just an 8 split.
1114<Addi> that is - split each dimension in half.
1115<TonyC> yep
1116<TonyC> I get the idea
1117<Addi> But this would seem to be a good approach in our case since we
1118usually have few codevectors.
1119<Addi> So we only need 256*256 entries in a table.
1120<Addi> We could even only index some of them that were deemed as good
1121candidates.
1122<TonyC> I was considering adding paletted output support for PNG and
1123TIFF at some point, which support 16-bit palettes
1124<Addi> ohh.
1125<Addi> 'darn' ;)
1126
1127
1128*/
1129
1130
1131typedef struct i_dists {
1132 int index;
1133 long dist;
1134} i_dists;
1135
1136#define CF_VARS \
1137 i_dists *dists;
1138
1139static int dists_sort(void const *a, void const *b) {
1140 return ((i_dists *)a)->dist - ((i_dists *)b)->dist;
1141}
1142
1143static void rand2dist_setup(i_quantize *quant, i_dists **cdists) {
1144 i_dists *dists =
1145 mymalloc(sizeof(i_dists)*quant->mc_count*quant->mc_count);
1146 int i, j;
1147 long cd;
1148 for (i = 0; i < quant->mc_count; ++i) {
1149 i_dists *ldists = dists + quant->mc_count * i;
1150 i_color val = quant->mc_colors[i];
1151 for (j = 0; j < quant->mc_count; ++j) {
1152 ldists[j].index = j;
1153 ldists[j].dist = ceucl_d(&val, quant->mc_colors+j);
1154 }
1155 qsort(ldists, quant->mc_count, sizeof(i_dists), dists_sort);
1156 }
1157 *cdists = dists;
1158}
1159
1160#define CF_SETUP \
1161 bst_idx = rand() % quant->mc_count; \
1162 rand2dist_setup(quant, &dists)
1163
1164static int rand2dist_find(i_color val, i_quantize *quant, i_dists *dists, int index) {
1165 i_dists *cdists;
1166 long cd, ld;
1167 long maxld;
1168 int i;
1169 int bst_idx;
1170
1171 cdists = dists + index * quant->mc_count;
1172 ld = 3 * 256 * 256;
1173 maxld = 8 * ceucl_d(&val, quant->mc_colors+index);
1174 for (i = 0; i < quant->mc_count && cdists[i].dist <= maxld; ++i) {
1175 cd = ceucl_d(&val, quant->mc_colors+cdists[i].index);
1176 if (cd < ld) {
1177 bst_idx = cdists[i].index;
1178 ld = cd;
1179 }
1180 }
1181 return bst_idx;
1182}
1183
1184#define CF_FIND bst_idx = rand2dist_find(val, quant, dists, bst_idx)
1185
1186#define CF_CLEANUP myfree(dists)
1187
1188
1189#endif
1190
1191static void translate_addi(i_quantize *quant, i_img *img, i_palidx *out) {
1192 int x, y, i, k, bst_idx;
1193 i_color val;
1194 int pixdev = quant->perturb;
1195 CF_VARS;
1196
1197 CF_SETUP;
1198
18accb2a
TC
1199 if (img->channels >= 3) {
1200 if (pixdev) {
1201 k=0;
1202 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1203 i_gpix(img,x,y,&val);
1204 val.channel[0]=g_sat(val.channel[0]+(int)(pixdev*frandn()));
1205 val.channel[1]=g_sat(val.channel[1]+(int)(pixdev*frandn()));
1206 val.channel[2]=g_sat(val.channel[2]+(int)(pixdev*frandn()));
1207 CF_FIND;
1208 out[k++]=bst_idx;
1209 }
1210 } else {
1211 k=0;
1212 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1213 i_gpix(img,x,y,&val);
1214 CF_FIND;
1215 out[k++]=bst_idx;
1216 }
02d1d628 1217 }
18accb2a
TC
1218 }
1219 else {
1220 if (pixdev) {
1221 k=0;
1222 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1223 i_gpix(img,x,y,&val);
1224 val.channel[1] = val.channel[2] =
1225 val.channel[0]=g_sat(val.channel[0]+(int)(pixdev*frandn()));
1226 CF_FIND;
1227 out[k++]=bst_idx;
1228 }
1229 } else {
1230 k=0;
1231 for(y=0;y<img->ysize;y++) for(x=0;x<img->xsize;x++) {
1232 i_gpix(img,x,y,&val);
1233 val.channel[1] = val.channel[2] = val.channel[0];
1234 CF_FIND;
1235 out[k++]=bst_idx;
1236 }
02d1d628
AMH
1237 }
1238 }
1239 CF_CLEANUP;
1240}
1241
1242static int floyd_map[] =
1243{
1244 0, 0, 7,
1245 3, 5, 1
1246};
1247
1248static int jarvis_map[] =
1249{
1250 0, 0, 0, 7, 5,
1251 3, 5, 7, 5, 3,
1252 1, 3, 5, 3, 1
1253};
1254
1255static int stucki_map[] =
1256{
1257 0, 0, 0, 8, 4,
1258 2, 4, 8, 4, 2,
1259 1, 2, 4, 2, 1
1260};
1261
1262struct errdiff_map {
1263 int *map;
1264 int width, height, orig;
1265};
1266
1267static struct errdiff_map maps[] =
1268{
1269 { floyd_map, 3, 2, 1 },
1270 { jarvis_map, 5, 3, 2 },
1271 { stucki_map, 5, 3, 2 },
1272};
1273
1274typedef struct errdiff_tag {
1275 int r, g, b;
1276} errdiff_t;
1277
1278/* perform an error diffusion dither */
1279static
1280void
1281translate_errdiff(i_quantize *quant, i_img *img, i_palidx *out) {
1282 int *map;
1283 int mapw, maph, mapo;
1284 int i;
1285 errdiff_t *err;
1286 int errw;
1287 int difftotal;
1288 int x, y, dx, dy;
1289 int minr, maxr, ming, maxg, minb, maxb, cr, cg, cb;
1290 i_color find;
1291 int bst_idx;
1292 CF_VARS;
1293
1294 if ((quant->errdiff & ed_mask) == ed_custom) {
1295 map = quant->ed_map;
1296 mapw = quant->ed_width;
1297 maph = quant->ed_height;
1298 mapo = quant->ed_orig;
1299 }
1300 else {
1301 int index = quant->errdiff & ed_mask;
1302 if (index >= ed_custom) index = ed_floyd;
1303 map = maps[index].map;
1304 mapw = maps[index].width;
1305 maph = maps[index].height;
1306 mapo = maps[index].orig;
1307 }
1308
1309 errw = img->xsize+mapw;
1310 err = mymalloc(sizeof(*err) * maph * errw);
1311 /*errp = err+mapo;*/
1312 memset(err, 0, sizeof(*err) * maph * errw);
1313
1314 difftotal = 0;
1315 for (i = 0; i < maph * mapw; ++i)
1316 difftotal += map[i];
1317 /*printf("map:\n");
1318 for (dy = 0; dy < maph; ++dy) {
1319 for (dx = 0; dx < mapw; ++dx) {
1320 printf("%2d", map[dx+dy*mapw]);
1321 }
1322 putchar('\n');
1323 }*/
1324
1325 CF_SETUP;
1326
1327 for (y = 0; y < img->ysize; ++y) {
1328 for (x = 0; x < img->xsize; ++x) {
1329 i_color val;
1330 long ld, cd;
1331 errdiff_t perr;
1332 i_gpix(img, x, y, &val);
18accb2a
TC
1333 if (img->channels < 3) {
1334 val.channel[1] = val.channel[2] = val.channel[0];
1335 }
02d1d628
AMH
1336 perr = err[x+mapo];
1337 perr.r = perr.r < 0 ? -((-perr.r)/difftotal) : perr.r/difftotal;
1338 perr.g = perr.g < 0 ? -((-perr.g)/difftotal) : perr.g/difftotal;
1339 perr.b = perr.b < 0 ? -((-perr.b)/difftotal) : perr.b/difftotal;
1340 /*printf("x %3d y %3d in(%3d, %3d, %3d) di(%4d,%4d,%4d)\n", x, y, val.channel[0], val.channel[1], val.channel[2], perr.r, perr.g, perr.b);*/
1341 val.channel[0] = g_sat(val.channel[0]-perr.r);
1342 val.channel[1] = g_sat(val.channel[1]-perr.g);
1343 val.channel[2] = g_sat(val.channel[2]-perr.b);
1344 CF_FIND;
1345 /* save error */
1346 perr.r = quant->mc_colors[bst_idx].channel[0] - val.channel[0];
1347 perr.g = quant->mc_colors[bst_idx].channel[1] - val.channel[1];
1348 perr.b = quant->mc_colors[bst_idx].channel[2] - val.channel[2];
1349 /*printf(" out(%3d, %3d, %3d) er(%4d, %4d, %4d)\n", quant->mc_colors[bst_idx].channel[0], quant->mc_colors[bst_idx].channel[1], quant->mc_colors[bst_idx].channel[2], perr.r, perr.g, perr.b);*/
1350 for (dx = 0; dx < mapw; ++dx) {
1351 for (dy = 0; dy < maph; ++dy) {
1352 err[x+dx+dy*errw].r += perr.r * map[dx+mapw*dy];
1353 err[x+dx+dy*errw].g += perr.g * map[dx+mapw*dy];
1354 err[x+dx+dy*errw].b += perr.b * map[dx+mapw*dy];
1355 }
1356 }
1357 *out++ = bst_idx;
1358 }
1359 /* shift up the error matrix */
1360 for (dy = 0; dy < maph-1; ++dy) {
1361 memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1362 }
1363 memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1364 }
1365 CF_CLEANUP;
7fd765fe 1366 myfree(err);
02d1d628
AMH
1367}
1368/* Prescan finds the boxes in the image that have the highest number of colors
1369 and that result is used as the initial value for the vectores */
1370
1371
18accb2a 1372static void prescan(i_img **imgs,int count, int cnum, cvec *clr, i_sample_t *line) {
02d1d628 1373 int i,k,j,x,y;
18accb2a
TC
1374 i_sample_t *val;
1375 const int *chans;
02d1d628
AMH
1376
1377 pbox prebox[512];
1378 for(i=0;i<512;i++) {
1379 prebox[i].boxnum=i;
1380 prebox[i].pixcnt=0;
1381 prebox[i].cand=1;
1382 }
1383
1384 /* process each image */
1385 for (i = 0; i < count; ++i) {
1386 i_img *im = imgs[i];
18accb2a
TC
1387 chans = im->channels >= 3 ? NULL : gray_samples;
1388 for(y=0;y<im->ysize;y++) {
1389 i_gsamp(im, 0, im->xsize, y, line, chans, 3);
1390 val = line;
1391 for(x=0;x<im->xsize;x++) {
1392 prebox[pixbox_ch(val)].pixcnt++;
1393 }
02d1d628
AMH
1394 }
1395 }
1396
1397 for(i=0;i<512;i++) prebox[i].pdc=prebox[i].pixcnt;
1398 qsort(prebox,512,sizeof(pbox),(cmpfunc)pboxcmp);
1399
1400 for(i=0;i<cnum;i++) {
1401 /* printf("Color %d\n",i);
1402 for(k=0;k<10;k++) { printf("box=%03d %04d %d %04d \n",prebox[k].boxnum,prebox[k].pixcnt,prebox[k].cand,prebox[k].pdc); }
1403 printf("\n\n"); */
1404 reorder(prebox);
1405 }
1406
1407 /* for(k=0;k<cnum;k++) { printf("box=%03d %04d %d %04d \n",prebox[k].boxnum,prebox[k].pixcnt,prebox[k].cand,prebox[k].pdc); } */
1408
1409 k=0;
1410 j=1;
1411 i=0;
1412 while(i<cnum) {
1413 /* printf("prebox[%d].cand=%d\n",k,prebox[k].cand); */
36e67d0b 1414 if (clr[i].fixed) { i++; continue; } /* reserved go to next */
02d1d628
AMH
1415 if (j>=prebox[k].cand) { k++; j=1; } else {
1416 if (prebox[k].cand == 2) boxcenter(prebox[k].boxnum,&(clr[i]));
1417 else boxrand(prebox[k].boxnum,&(clr[i]));
1418 /* printf("(%d,%d) %d %d -> (%d,%d,%d)\n",k,j,prebox[k].boxnum,prebox[k].pixcnt,clr[i].r,clr[i].g,clr[i].b); */
1419 j++;
1420 i++;
1421 }
1422 }
1423}
1424
1425
1426static void reorder(pbox prescan[512]) {
1427 int nidx;
1428 pbox c;
1429
1430 nidx=0;
1431 c=prescan[0];
1432
1433 c.cand++;
1434 c.pdc=c.pixcnt/(c.cand*c.cand);
1435 /* c.pdc=c.pixcnt/c.cand; */
1436 while(c.pdc < prescan[nidx+1].pdc && nidx < 511) {
1437 prescan[nidx]=prescan[nidx+1];
1438 nidx++;
1439 }
1440 prescan[nidx]=c;
1441}
1442
1443static int
1444pboxcmp(const pbox *a,const pbox *b) {
1445 if (a->pixcnt > b->pixcnt) return -1;
1446 if (a->pixcnt < b->pixcnt) return 1;
1447 return 0;
1448}
1449
1450static void
1451boxcenter(int box,cvec *cv) {
1452 cv->r=15+((box&448)>>1);
1453 cv->g=15+((box&56)<<2);
1454 cv->b=15+((box&7)<<5);
1455}
1456
1457static void
1458bbox(int box,int *r0,int *r1,int *g0,int *g1,int *b0,int *b1) {
1459 *r0=(box&448)>>1;
1460 *r1=(*r0)|31;
1461 *g0=(box&56)<<2;
1462 *g1=(*g0)|31;
1463 *b0=(box&7)<<5;
1464 *b1=(*b0)|31;
1465}
1466
1467static void
1468boxrand(int box,cvec *cv) {
1469 cv->r=6+(rand()%25)+((box&448)>>1);
1470 cv->g=6+(rand()%25)+((box&56)<<2);
1471 cv->b=6+(rand()%25)+((box&7)<<5);
1472}
1473
1474static float
1475frandn(void) {
1476
1477 float u1,u2,w;
1478
1479 w=1;
1480
1481 while (w >= 1 || w == 0) {
1482 u1 = 2 * frand() - 1;
1483 u2 = 2 * frand() - 1;
1484 w = u1*u1 + u2*u2;
1485 }
1486
1487 w = sqrt((-2*log(w))/w);
1488 return u1*w;
1489}
1490
1491/* Create hash index */
1492static
1493void
1494cr_hashindex(cvec clr[256],int cnum,hashbox hb[512]) {
1495
1496 int bx,mind,cd,cumcnt,bst_idx,i;
1497/* printf("indexing... \n");*/
1498
1499 cumcnt=0;
1500 for(bx=0; bx<512; bx++) {
1501 mind=196608;
1502 for(i=0; i<cnum; i++) {
1503 cd = maxdist(bx,&clr[i]);
1504 if (cd < mind) { mind=cd; bst_idx=i; }
1505 }
1506
1507 hb[bx].cnt=0;
1508 for(i=0;i<cnum;i++) if (mindist(bx,&clr[i])<mind) hb[bx].vec[hb[bx].cnt++]=i;
1509 /*printf("box %d -> approx -> %d\n",bx,hb[bx].cnt); */
1510 /* statbox(bx,cnum,clr); */
1511 cumcnt+=hb[bx].cnt;
1512 }
1513
1514/* printf("Average search space: %d\n",cumcnt/512); */
1515}
1516
1517static int
1518maxdist(int boxnum,cvec *cv) {
1519 int r0,r1,g0,g1,b0,b1;
1520 int r,g,b,mr,mg,mb;
1521
1522 r=cv->r;
1523 g=cv->g;
1524 b=cv->b;
1525
1526 bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1527
b33c08f8
TC
1528 mr=i_max(abs(b-b0),abs(b-b1));
1529 mg=i_max(abs(g-g0),abs(g-g1));
1530 mb=i_max(abs(r-r0),abs(r-r1));
02d1d628
AMH
1531
1532 return PWR2(mr)+PWR2(mg)+PWR2(mb);
1533}
1534
1535static int
1536mindist(int boxnum,cvec *cv) {
1537 int r0,r1,g0,g1,b0,b1;
1538 int r,g,b,mr,mg,mb;
1539
1540 r=cv->r;
1541 g=cv->g;
1542 b=cv->b;
1543
1544 bbox(boxnum,&r0,&r1,&g0,&g1,&b0,&b1);
1545
1546 /* printf("box %d, (%d,%d,%d)-(%d,%d,%d) vec (%d,%d,%d) ",boxnum,r0,g0,b0,r1,g1,b1,r,g,b); */
1547
1548 if (r0<=r && r<=r1 && g0<=g && g<=g1 && b0<=b && b<=b1) return 0;
1549
b33c08f8
TC
1550 mr=i_min(abs(b-b0),abs(b-b1));
1551 mg=i_min(abs(g-g0),abs(g-g1));
1552 mb=i_min(abs(r-r0),abs(r-r1));
02d1d628
AMH
1553
1554 mr=PWR2(mr);
1555 mg=PWR2(mg);
1556 mb=PWR2(mb);
1557
1558 if (r0<=r && r<=r1 && g0<=g && g<=g1) return mb;
1559 if (r0<=r && r<=r1 && b0<=b && b<=b1) return mg;
1560 if (b0<=b && b<=b1 && g0<=g && g<=g1) return mr;
1561
1562 if (r0<=r && r<=r1) return mg+mb;
1563 if (g0<=g && g<=g1) return mr+mb;
1564 if (b0<=b && b<=b1) return mg+mr;
1565
1566 return mr+mg+mb;
1567}
1568
1569static void transparent_threshold(i_quantize *, i_palidx *, i_img *, i_palidx);
1570static void transparent_errdiff(i_quantize *, i_palidx *, i_img *, i_palidx);
1571static void transparent_ordered(i_quantize *, i_palidx *, i_img *, i_palidx);
1572
1573void quant_transparent(i_quantize *quant, i_palidx *data, i_img *img,
1574 i_palidx trans_index)
1575{
1576 switch (quant->transp) {
1577 case tr_none:
1578 break;
1579
1580 default:
1581 quant->tr_threshold = 128;
1582 /* fall through */
1583 case tr_threshold:
1584 transparent_threshold(quant, data, img, trans_index);
1585 break;
1586
1587 case tr_errdiff:
1588 transparent_errdiff(quant, data, img, trans_index);
1589 break;
1590
1591 case tr_ordered:
1592 transparent_ordered(quant, data, img, trans_index);
1593 break;
1594 }
1595}
1596
1597static void
1598transparent_threshold(i_quantize *quant, i_palidx *data, i_img *img,
1599 i_palidx trans_index)
1600{
1601 int x, y;
18accb2a
TC
1602 i_sample_t *line = mymalloc(img->xsize * sizeof(i_sample_t));
1603 int trans_chan = img->channels > 2 ? 3 : 1;
02d1d628
AMH
1604
1605 for (y = 0; y < img->ysize; ++y) {
18accb2a 1606 i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
02d1d628 1607 for (x = 0; x < img->xsize; ++x) {
18accb2a 1608 if (line[x] < quant->tr_threshold)
02d1d628
AMH
1609 data[y*img->xsize+x] = trans_index;
1610 }
1611 }
18accb2a 1612 myfree(line);
02d1d628
AMH
1613}
1614
1615static void
1616transparent_errdiff(i_quantize *quant, i_palidx *data, i_img *img,
1617 i_palidx trans_index)
1618{
1619 int *map;
1620 int index;
1621 int mapw, maph, mapo;
1622 int errw, *err, *errp;
1623 int difftotal, out, error;
1624 int x, y, dx, dy, i;
18accb2a
TC
1625 i_sample_t *line;
1626 int trans_chan = img->channels > 2 ? 3 : 1;
02d1d628
AMH
1627
1628 /* no custom map for transparency (yet) */
1629 index = quant->tr_errdiff & ed_mask;
1630 if (index >= ed_custom) index = ed_floyd;
1631 map = maps[index].map;
1632 mapw = maps[index].width;
1633 maph = maps[index].height;
1634 mapo = maps[index].orig;
1635
1636 errw = img->xsize+mapw-1;
1637 err = mymalloc(sizeof(*err) * maph * errw);
1638 errp = err+mapo;
1639 memset(err, 0, sizeof(*err) * maph * errw);
1640
18accb2a 1641 line = mymalloc(img->xsize * sizeof(i_sample_t));
02d1d628
AMH
1642 difftotal = 0;
1643 for (i = 0; i < maph * mapw; ++i)
1644 difftotal += map[i];
1645 for (y = 0; y < img->ysize; ++y) {
18accb2a 1646 i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
02d1d628 1647 for (x = 0; x < img->xsize; ++x) {
18accb2a
TC
1648 line[x] = g_sat(line[x]-errp[x]/difftotal);
1649 if (line[x] < 128) {
02d1d628
AMH
1650 out = 0;
1651 data[y*img->xsize+x] = trans_index;
1652 }
1653 else {
1654 out = 255;
1655 }
18accb2a 1656 error = out - line[x];
02d1d628
AMH
1657 for (dx = 0; dx < mapw; ++dx) {
1658 for (dy = 0; dy < maph; ++dy) {
1659 errp[x+dx-mapo+dy*errw] += error * map[dx+mapw*dy];
1660 }
1661 }
1662 }
1663 /* shift up the error matrix */
1664 for (dy = 0; dy < maph-1; ++dy)
1665 memcpy(err+dy*errw, err+(dy+1)*errw, sizeof(*err)*errw);
1666 memset(err+(maph-1)*errw, 0, sizeof(*err)*errw);
1667 }
18accb2a
TC
1668 myfree(err);
1669 myfree(line);
02d1d628
AMH
1670}
1671
1672/* builtin ordered dither maps */
b33c08f8
TC
1673static unsigned char
1674orddith_maps[][64] =
02d1d628
AMH
1675{
1676 { /* random
1677 this is purely random - it's pretty awful
1678 */
1679 48, 72, 196, 252, 180, 92, 108, 52,
1680 228, 176, 64, 8, 236, 40, 20, 164,
1681 120, 128, 84, 116, 24, 28, 172, 220,
1682 68, 0, 188, 124, 184, 224, 192, 104,
1683 132, 100, 240, 200, 152, 160, 244, 44,
1684 96, 204, 144, 16, 140, 56, 232, 216,
1685 208, 4, 76, 212, 136, 248, 80, 168,
1686 156, 88, 32, 112, 148, 12, 36, 60,
1687 },
1688 {
1689 /* dot8
1690 perl spot.perl '($x-3.5)*($x-3.5)+($y-3.5)*($y-3.5)'
1691 */
1692 240, 232, 200, 136, 140, 192, 228, 248,
1693 220, 148, 100, 76, 80, 104, 152, 212,
1694 180, 116, 56, 32, 36, 60, 120, 176,
1695 156, 64, 28, 0, 8, 44, 88, 160,
1696 128, 92, 24, 12, 4, 40, 68, 132,
1697 184, 96, 48, 20, 16, 52, 108, 188,
1698 216, 144, 112, 72, 84, 124, 164, 224,
1699 244, 236, 196, 168, 172, 204, 208, 252,
1700 },
1701 { /* dot4
1702 perl spot.perl \
1703 'min(dist(1.5, 1.5),dist(5.5,1.5),dist(1.5,5.5),dist(5.5,5.5))'
1704 */
1705 196, 72, 104, 220, 200, 80, 112, 224,
1706 76, 4, 24, 136, 84, 8, 32, 144,
1707 108, 28, 52, 168, 116, 36, 56, 176,
1708 216, 140, 172, 244, 228, 148, 180, 248,
1709 204, 92, 124, 236, 192, 68, 96, 208,
1710 88, 12, 44, 156, 64, 0, 16, 128,
1711 120, 40, 60, 188, 100, 20, 48, 160,
1712 232, 152, 184, 252, 212, 132, 164, 240,
1713 },
1714 { /* hline
1715 perl spot.perl '$y-3'
1716 */
1717 160, 164, 168, 172, 176, 180, 184, 188,
1718 128, 132, 136, 140, 144, 148, 152, 156,
1719 32, 36, 40, 44, 48, 52, 56, 60,
1720 0, 4, 8, 12, 16, 20, 24, 28,
1721 64, 68, 72, 76, 80, 84, 88, 92,
1722 96, 100, 104, 108, 112, 116, 120, 124,
1723 192, 196, 200, 204, 208, 212, 216, 220,
1724 224, 228, 232, 236, 240, 244, 248, 252,
1725 },
1726 { /* vline
1727 perl spot.perl '$x-3'
1728 */
1729 180, 100, 40, 12, 44, 104, 184, 232,
1730 204, 148, 60, 16, 64, 128, 208, 224,
1731 212, 144, 76, 8, 80, 132, 216, 244,
1732 160, 112, 68, 20, 84, 108, 172, 236,
1733 176, 96, 72, 28, 88, 152, 188, 228,
1734 200, 124, 92, 0, 32, 116, 164, 240,
1735 168, 120, 36, 24, 48, 136, 192, 248,
1736 196, 140, 52, 4, 56, 156, 220, 252,
1737 },
1738 { /* slashline
1739 perl spot.perl '$y+$x-7'
1740 */
1741 248, 232, 224, 192, 140, 92, 52, 28,
1742 240, 220, 196, 144, 108, 60, 12, 64,
1743 216, 180, 148, 116, 76, 20, 80, 128,
1744 204, 152, 104, 44, 16, 72, 100, 160,
1745 164, 96, 68, 24, 56, 112, 168, 176,
1746 124, 40, 8, 36, 88, 136, 184, 212,
1747 84, 4, 32, 120, 156, 188, 228, 236,
1748 0, 48, 132, 172, 200, 208, 244, 252,
1749 },
1750 { /* backline
1751 perl spot.perl '$y-$x'
1752 */
1753 0, 32, 116, 172, 184, 216, 236, 252,
1754 56, 8, 72, 132, 136, 200, 228, 240,
1755 100, 36, 12, 40, 92, 144, 204, 220,
1756 168, 120, 60, 16, 44, 96, 156, 176,
1757 180, 164, 112, 48, 28, 52, 128, 148,
1758 208, 192, 152, 88, 84, 20, 64, 104,
1759 232, 224, 196, 140, 108, 68, 24, 76,
1760 248, 244, 212, 188, 160, 124, 80, 4,
1761 },
11e7329d
TC
1762 {
1763 /* tiny
1764 good for display, bad for print
1765 hand generated
1766 */
1767 0, 128, 32, 192, 8, 136, 40, 200,
1768 224, 64, 160, 112, 232, 72, 168, 120,
1769 48, 144, 16, 208, 56, 152, 24, 216,
1770 176, 96, 240, 80, 184, 104, 248, 88,
1771 12, 140, 44, 204, 4, 132, 36, 196,
1772 236, 76, 172, 124, 228, 68, 164, 116,
1773 60, 156, 28, 220, 52, 148, 20, 212,
1774 188, 108, 252, 92, 180, 100, 244, 84,
1775 },
02d1d628
AMH
1776};
1777
1778static void
1779transparent_ordered(i_quantize *quant, i_palidx *data, i_img *img,
1780 i_palidx trans_index)
1781{
1782 unsigned char *spot;
1783 int x, y;
18accb2a
TC
1784 i_sample_t *line;
1785 int trans_chan = img->channels > 2 ? 3 : 1;
02d1d628
AMH
1786 if (quant->tr_orddith == od_custom)
1787 spot = quant->tr_custom;
1788 else
1789 spot = orddith_maps[quant->tr_orddith];
18accb2a
TC
1790
1791 line = mymalloc(img->xsize * sizeof(i_sample_t));
02d1d628 1792 for (y = 0; y < img->ysize; ++y) {
18accb2a 1793 i_gsamp(img, 0, img->xsize, y, line, &trans_chan, 1);
02d1d628 1794 for (x = 0; x < img->xsize; ++x) {
18accb2a 1795 if (line[x] < spot[(x&7)+(y&7)*8])
02d1d628
AMH
1796 data[x+y*img->xsize] = trans_index;
1797 }
1798 }
18accb2a 1799 myfree(line);
02d1d628 1800}
18accb2a 1801