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