1 /*
2 ** © 2009-2018 by Kornel Lesiński.
3 ** © 1989, 1991 by Jef Poskanzer.
4 ** © 1997, 2000, 2002 by Greg Roelofs; based on an idea by Stefan Schneider.
5 **
6 ** See COPYRIGHT file for license.
7 */
8
9 #include <stdlib.h>
10 #include <stddef.h>
11
12 #include "libimagequant.h"
13 #include "pam.h"
14 #include "mediancut.h"
15
16 #define index_of_channel(ch) (offsetof(f_pixel,ch)/sizeof(float))
17
18 static f_pixel averagepixels(unsigned int clrs, const hist_item achv[]);
19
20 struct box {
21 f_pixel color;
22 f_pixel variance;
23 double sum, total_error, max_error;
24 unsigned int ind;
25 unsigned int colors;
26 };
27
28 ALWAYS_INLINE static double variance_diff(double val, const double good_enough);
variance_diff(double val,const double good_enough)29 inline static double variance_diff(double val, const double good_enough)
30 {
31 val *= val;
32 if (val < good_enough*good_enough) return val*0.25;
33 return val;
34 }
35
36 /** Weighted per-channel variance of the box. It's used to decide which channel to split by */
box_variance(const hist_item achv[],const struct box * box)37 static f_pixel box_variance(const hist_item achv[], const struct box *box)
38 {
39 f_pixel mean = box->color;
40 double variancea=0, variancer=0, varianceg=0, varianceb=0;
41
42 for(unsigned int i = 0; i < box->colors; ++i) {
43 const f_pixel px = achv[box->ind + i].acolor;
44 double weight = achv[box->ind + i].adjusted_weight;
45 variancea += variance_diff(mean.a - px.a, 2.0/256.0)*weight;
46 variancer += variance_diff(mean.r - px.r, 1.0/256.0)*weight;
47 varianceg += variance_diff(mean.g - px.g, 1.0/256.0)*weight;
48 varianceb += variance_diff(mean.b - px.b, 1.0/256.0)*weight;
49 }
50
51 return (f_pixel){
52 .a = variancea*(4.0/16.0),
53 .r = variancer*(7.0/16.0),
54 .g = varianceg*(9.0/16.0),
55 .b = varianceb*(5.0/16.0),
56 };
57 }
58
box_max_error(const hist_item achv[],const struct box * box)59 static double box_max_error(const hist_item achv[], const struct box *box)
60 {
61 f_pixel mean = box->color;
62 double max_error = 0;
63
64 for(unsigned int i = 0; i < box->colors; ++i) {
65 const double diff = colordifference(mean, achv[box->ind + i].acolor);
66 if (diff > max_error) {
67 max_error = diff;
68 }
69 }
70 return max_error;
71 }
72
73 ALWAYS_INLINE static double color_weight(f_pixel median, hist_item h);
74
hist_item_swap(hist_item * l,hist_item * r)75 static inline void hist_item_swap(hist_item *l, hist_item *r)
76 {
77 if (l != r) {
78 hist_item t = *l;
79 *l = *r;
80 *r = t;
81 }
82 }
83
84 ALWAYS_INLINE static unsigned int qsort_pivot(const hist_item *const base, const unsigned int len);
qsort_pivot(const hist_item * const base,const unsigned int len)85 inline static unsigned int qsort_pivot(const hist_item *const base, const unsigned int len)
86 {
87 if (len < 32) {
88 return len/2;
89 }
90
91 const unsigned int aidx=8, bidx=len/2, cidx=len-1;
92 const unsigned int a=base[aidx].tmp.sort_value, b=base[bidx].tmp.sort_value, c=base[cidx].tmp.sort_value;
93 return (a < b) ? ((b < c) ? bidx : ((a < c) ? cidx : aidx ))
94 : ((b > c) ? bidx : ((a < c) ? aidx : cidx ));
95 }
96
97 ALWAYS_INLINE static unsigned int qsort_partition(hist_item *const base, const unsigned int len);
qsort_partition(hist_item * const base,const unsigned int len)98 inline static unsigned int qsort_partition(hist_item *const base, const unsigned int len)
99 {
100 unsigned int l = 1, r = len;
101 if (len >= 8) {
102 hist_item_swap(&base[0], &base[qsort_pivot(base,len)]);
103 }
104
105 const unsigned int pivot_value = base[0].tmp.sort_value;
106 while (l < r) {
107 if (base[l].tmp.sort_value >= pivot_value) {
108 l++;
109 } else {
110 while(l < --r && base[r].tmp.sort_value <= pivot_value) {}
111 hist_item_swap(&base[l], &base[r]);
112 }
113 }
114 l--;
115 hist_item_swap(&base[0], &base[l]);
116
117 return l;
118 }
119
120 /** quick select algorithm */
hist_item_sort_range(hist_item base[],unsigned int len,unsigned int sort_start)121 static void hist_item_sort_range(hist_item base[], unsigned int len, unsigned int sort_start)
122 {
123 for(;;) {
124 const unsigned int l = qsort_partition(base, len), r = l+1;
125
126 if (l > 0 && sort_start < l) {
127 len = l;
128 }
129 else if (r < len && sort_start > r) {
130 base += r; len -= r; sort_start -= r;
131 }
132 else break;
133 }
134 }
135
136 /** sorts array to make sum of weights lower than halfvar one side, returns edge between <halfvar and >halfvar parts of the set */
hist_item_sort_halfvar(hist_item base[],unsigned int len,double * const lowervar,const double halfvar)137 static hist_item *hist_item_sort_halfvar(hist_item base[], unsigned int len, double *const lowervar, const double halfvar)
138 {
139 do {
140 const unsigned int l = qsort_partition(base, len), r = l+1;
141
142 // check if sum of left side is smaller than half,
143 // if it is, then it doesn't need to be sorted
144 unsigned int t = 0; double tmpsum = *lowervar;
145 while (t <= l && tmpsum < halfvar) tmpsum += base[t++].color_weight;
146
147 if (tmpsum < halfvar) {
148 *lowervar = tmpsum;
149 } else {
150 if (l > 0) {
151 hist_item *res = hist_item_sort_halfvar(base, l, lowervar, halfvar);
152 if (res) return res;
153 } else {
154 // End of left recursion. This will be executed in order from the first element.
155 *lowervar += base[0].color_weight;
156 if (*lowervar > halfvar) return &base[0];
157 }
158 }
159
160 if (len > r) {
161 base += r; len -= r; // tail-recursive "call"
162 } else {
163 *lowervar += base[r].color_weight;
164 return (*lowervar > halfvar) ? &base[r] : NULL;
165 }
166 } while(1);
167 }
168
169 static f_pixel get_median(const struct box *b, hist_item achv[]);
170
171 typedef struct {
172 unsigned int chan; float variance;
173 } channelvariance;
174
comparevariance(const void * ch1,const void * ch2)175 static int comparevariance(const void *ch1, const void *ch2)
176 {
177 return ((const channelvariance*)ch1)->variance > ((const channelvariance*)ch2)->variance ? -1 :
178 (((const channelvariance*)ch1)->variance < ((const channelvariance*)ch2)->variance ? 1 : 0);
179 }
180
181 /** Finds which channels need to be sorted first and preproceses achv for fast sort */
prepare_sort(struct box * b,hist_item achv[])182 static double prepare_sort(struct box *b, hist_item achv[])
183 {
184 /*
185 ** Sort dimensions by their variance, and then sort colors first by dimension with highest variance
186 */
187 channelvariance channels[4] = {
188 {index_of_channel(a), b->variance.a},
189 {index_of_channel(r), b->variance.r},
190 {index_of_channel(g), b->variance.g},
191 {index_of_channel(b), b->variance.b},
192 };
193
194 qsort(channels, 4, sizeof(channels[0]), comparevariance);
195
196 const unsigned int ind1 = b->ind;
197 const unsigned int colors = b->colors;
198 #if __GNUC__ >= 9 || __clang__
199 #pragma omp parallel for if (colors > 25000) \
200 schedule(static) default(none) shared(achv, channels, colors, ind1)
201 #else
202 #pragma omp parallel for if (colors > 25000) \
203 schedule(static) default(none) shared(achv, channels)
204 #endif
205 for(unsigned int i=0; i < colors; i++) {
206 const float *chans = (const float *)&achv[ind1 + i].acolor;
207 // Only the first channel really matters. When trying median cut many times
208 // with different histogram weights, I don't want sort randomness to influence outcome.
209 achv[ind1 + i].tmp.sort_value = ((unsigned int)(chans[channels[0].chan]*65535.0)<<16) |
210 (unsigned int)((chans[channels[2].chan] + chans[channels[1].chan]/2.0 + chans[channels[3].chan]/4.0)*65535.0);
211 }
212
213 const f_pixel median = get_median(b, achv);
214
215 // box will be split to make color_weight of each side even
216 const unsigned int ind = b->ind, end = ind+b->colors;
217 double totalvar = 0;
218 #pragma omp parallel for if (end - ind > 15000) \
219 schedule(static) default(shared) reduction(+:totalvar)
220 for(unsigned int j=ind; j < end; j++) totalvar += (achv[j].color_weight = color_weight(median, achv[j]));
221 return totalvar / 2.0;
222 }
223
224 /** finds median in unsorted set by sorting only minimum required */
get_median(const struct box * b,hist_item achv[])225 static f_pixel get_median(const struct box *b, hist_item achv[])
226 {
227 const unsigned int median_start = (b->colors-1)/2;
228
229 hist_item_sort_range(&(achv[b->ind]), b->colors,
230 median_start);
231
232 if (b->colors&1) return achv[b->ind + median_start].acolor;
233
234 // technically the second color is not guaranteed to be sorted correctly
235 // but most of the time it is good enough to be useful
236 return averagepixels(2, &achv[b->ind + median_start]);
237 }
238
239 /*
240 ** Find the best splittable box. -1 if no boxes are splittable.
241 */
best_splittable_box(struct box bv[],unsigned int boxes,const double max_mse)242 static int best_splittable_box(struct box bv[], unsigned int boxes, const double max_mse)
243 {
244 int bi=-1; double maxsum=0;
245 for(unsigned int i=0; i < boxes; i++) {
246 if (bv[i].colors < 2) {
247 continue;
248 }
249
250 // looks only at max variance, because it's only going to split by it
251 const double cv = MAX(bv[i].variance.r, MAX(bv[i].variance.g,bv[i].variance.b));
252 double thissum = bv[i].sum * MAX(bv[i].variance.a, cv);
253
254 if (bv[i].max_error > max_mse) {
255 thissum = thissum* bv[i].max_error/max_mse;
256 }
257
258 if (thissum > maxsum) {
259 maxsum = thissum;
260 bi = i;
261 }
262 }
263 return bi;
264 }
265
color_weight(f_pixel median,hist_item h)266 inline static double color_weight(f_pixel median, hist_item h)
267 {
268 float diff = colordifference(median, h.acolor);
269 return sqrt(diff) * (sqrt(1.0+h.adjusted_weight)-1.0);
270 }
271
272 static void set_colormap_from_boxes(colormap *map, struct box bv[], unsigned int boxes, hist_item *achv);
273 static void adjust_histogram(hist_item *achv, const struct box bv[], unsigned int boxes);
274
box_error(const struct box * box,const hist_item achv[])275 static double box_error(const struct box *box, const hist_item achv[])
276 {
277 f_pixel avg = box->color;
278
279 double total_error=0;
280 for (unsigned int i = 0; i < box->colors; ++i) {
281 total_error += colordifference(avg, achv[box->ind + i].acolor) * achv[box->ind + i].perceptual_weight;
282 }
283
284 return total_error;
285 }
286
287
total_box_error_below_target(double target_mse,struct box bv[],unsigned int boxes,const histogram * hist)288 static bool total_box_error_below_target(double target_mse, struct box bv[], unsigned int boxes, const histogram *hist)
289 {
290 target_mse *= hist->total_perceptual_weight;
291 double total_error=0;
292
293 for(unsigned int i=0; i < boxes; i++) {
294 // error is (re)calculated lazily
295 if (bv[i].total_error >= 0) {
296 total_error += bv[i].total_error;
297 }
298 if (total_error > target_mse) return false;
299 }
300
301 for(unsigned int i=0; i < boxes; i++) {
302 if (bv[i].total_error < 0) {
303 bv[i].total_error = box_error(&bv[i], hist->achv);
304 total_error += bv[i].total_error;
305 }
306 if (total_error > target_mse) return false;
307 }
308
309 return true;
310 }
311
box_init(struct box * box,const hist_item * achv,const unsigned int ind,const unsigned int colors,const double sum)312 static void box_init(struct box *box, const hist_item *achv, const unsigned int ind, const unsigned int colors, const double sum) {
313 box->ind = ind;
314 box->colors = colors;
315 box->sum = sum;
316 box->total_error = -1;
317
318 box->color = averagepixels(colors, &achv[ind]);
319 box->variance = box_variance(achv, box);
320 box->max_error = box_max_error(achv, box);
321 }
322
323 /*
324 ** Here is the fun part, the median-cut colormap generator. This is based
325 ** on Paul Heckbert's paper, "Color Image Quantization for Frame Buffer
326 ** Display," SIGGRAPH 1982 Proceedings, page 297.
327 */
mediancut(histogram * hist,unsigned int newcolors,const double target_mse,const double max_mse,void * (* malloc)(size_t),void (* free)(void *))328 LIQ_PRIVATE colormap *mediancut(histogram *hist, unsigned int newcolors, const double target_mse, const double max_mse, void* (*malloc)(size_t), void (*free)(void*))
329 {
330 hist_item *achv = hist->achv;
331 LIQ_ARRAY(struct box, bv, newcolors);
332 unsigned int boxes = 1;
333
334 /*
335 ** Set up the initial box.
336 */
337 {
338 double sum = 0;
339 for(unsigned int i=0; i < hist->size; i++) {
340 sum += achv[i].adjusted_weight;
341 }
342 box_init(&bv[0], achv, 0, hist->size, sum);
343
344
345 /*
346 ** Main loop: split boxes until we have enough.
347 */
348 while (boxes < newcolors) {
349
350 // first splits boxes that exceed quality limit (to have colors for things like odd green pixel),
351 // later raises the limit to allow large smooth areas/gradients get colors.
352 const double current_max_mse = max_mse + (boxes/(double)newcolors)*16.0*max_mse;
353 const int bi = best_splittable_box(bv, boxes, current_max_mse);
354 if (bi < 0) {
355 break; /* ran out of colors! */
356 }
357
358 unsigned int indx = bv[bi].ind;
359 unsigned int clrs = bv[bi].colors;
360
361 /*
362 Classic implementation tries to get even number of colors or pixels in each subdivision.
363
364 Here, instead of popularity I use (sqrt(popularity)*variance) metric.
365 Each subdivision balances number of pixels (popular colors) and low variance -
366 boxes can be large if they have similar colors. Later boxes with high variance
367 will be more likely to be split.
368
369 Median used as expected value gives much better results than mean.
370 */
371
372 const double halfvar = prepare_sort(&bv[bi], achv);
373 double lowervar=0;
374
375 // hist_item_sort_halfvar sorts and sums lowervar at the same time
376 // returns item to break at …minus one, which does smell like an off-by-one error.
377 hist_item *break_p = hist_item_sort_halfvar(&achv[indx], clrs, &lowervar, halfvar);
378 unsigned int break_at = MIN(clrs-1, break_p - &achv[indx] + 1);
379
380 /*
381 ** Split the box.
382 */
383 double sm = bv[bi].sum;
384 double lowersum = 0;
385 for(unsigned int i=0; i < break_at; i++) lowersum += achv[indx + i].adjusted_weight;
386
387 box_init(&bv[bi], achv, indx, break_at, lowersum);
388 box_init(&bv[boxes], achv, indx + break_at, clrs - break_at, sm - lowersum);
389
390 ++boxes;
391
392 if (total_box_error_below_target(target_mse, bv, boxes, hist)) {
393 break;
394 }
395 }
396 }
397
398 colormap *map = pam_colormap(boxes, malloc, free);
399 set_colormap_from_boxes(map, bv, boxes, achv);
400
401 adjust_histogram(achv, bv, boxes);
402
403 return map;
404 }
405
set_colormap_from_boxes(colormap * map,struct box * bv,unsigned int boxes,hist_item * achv)406 static void set_colormap_from_boxes(colormap *map, struct box* bv, unsigned int boxes, hist_item *achv)
407 {
408 /*
409 ** Ok, we've got enough boxes. Now choose a representative color for
410 ** each box. There are a number of possible ways to make this choice.
411 ** One would be to choose the center of the box; this ignores any structure
412 ** within the boxes. Another method would be to average all the colors in
413 ** the box - this is the method specified in Heckbert's paper.
414 */
415
416 for(unsigned int bi = 0; bi < boxes; ++bi) {
417 map->palette[bi].acolor = bv[bi].color;
418
419 /* store total color popularity (perceptual_weight is approximation of it) */
420 map->palette[bi].popularity = 0;
421 for(unsigned int i=bv[bi].ind; i < bv[bi].ind+bv[bi].colors; i++) {
422 map->palette[bi].popularity += achv[i].perceptual_weight;
423 }
424 }
425 }
426
427 /* increase histogram popularity by difference from the final color (this is used as part of feedback loop) */
adjust_histogram(hist_item * achv,const struct box * bv,unsigned int boxes)428 static void adjust_histogram(hist_item *achv, const struct box* bv, unsigned int boxes)
429 {
430 for(unsigned int bi = 0; bi < boxes; ++bi) {
431 for(unsigned int i=bv[bi].ind; i < bv[bi].ind+bv[bi].colors; i++) {
432 achv[i].tmp.likely_colormap_index = bi;
433 }
434 }
435 }
436
averagepixels(unsigned int clrs,const hist_item achv[])437 static f_pixel averagepixels(unsigned int clrs, const hist_item achv[])
438 {
439 double r = 0, g = 0, b = 0, a = 0, sum = 0;
440
441 #pragma omp parallel for if (clrs > 25000) \
442 schedule(static) default(shared) reduction(+:a) reduction(+:r) reduction(+:g) reduction(+:b) reduction(+:sum)
443 for(unsigned int i = 0; i < clrs; i++) {
444 const f_pixel px = achv[i].acolor;
445 const double weight = achv[i].adjusted_weight;
446
447 sum += weight;
448 a += px.a * weight;
449 r += px.r * weight;
450 g += px.g * weight;
451 b += px.b * weight;
452 }
453
454 if (sum) {
455 a /= sum;
456 r /= sum;
457 g /= sum;
458 b /= sum;
459 }
460
461 assert(!isnan(r) && !isnan(g) && !isnan(b) && !isnan(a));
462
463 return (f_pixel){.r=r, .g=g, .b=b, .a=a};
464 }
465