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