1 /*
2 ** © 2009-2015 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 "libimagequant.h"
10 #include "pam.h"
11 #include "nearest.h"
12 #include "mempool.h"
13 #include <stdlib.h>
14
15 typedef struct vp_sort_tmp {
16 float distance_squared;
17 unsigned int idx;
18 } vp_sort_tmp;
19
20 typedef struct vp_search_tmp {
21 float distance;
22 float distance_squared;
23 unsigned int idx;
24 int exclude;
25 } vp_search_tmp;
26
27 struct leaf {
28 f_pixel color;
29 unsigned int idx;
30 };
31
32 typedef struct vp_node {
33 struct vp_node *near, *far;
34 f_pixel vantage_point;
35 float radius, radius_squared;
36 struct leaf *rest;
37 unsigned short idx;
38 unsigned short restcount;
39 } vp_node;
40
41 struct nearest_map {
42 vp_node *root;
43 const colormap_item *palette;
44 float nearest_other_color_dist[256];
45 mempoolptr mempool;
46 };
47
48 static void vp_search_node(const vp_node *node, const f_pixel *const needle, vp_search_tmp *const best_candidate);
49
vp_compare_distance(const void * ap,const void * bp)50 static int vp_compare_distance(const void *ap, const void *bp) {
51 float a = ((const vp_sort_tmp*)ap)->distance_squared;
52 float b = ((const vp_sort_tmp*)bp)->distance_squared;
53 return a > b ? 1 : -1;
54 }
55
vp_sort_indexes_by_distance(const f_pixel vantage_point,vp_sort_tmp indexes[],int num_indexes,const colormap_item items[])56 static void vp_sort_indexes_by_distance(const f_pixel vantage_point, vp_sort_tmp indexes[], int num_indexes, const colormap_item items[]) {
57 for(int i=0; i < num_indexes; i++) {
58 indexes[i].distance_squared = colordifference(vantage_point, items[indexes[i].idx].acolor);
59 }
60 qsort(indexes, num_indexes, sizeof(indexes[0]), vp_compare_distance);
61 }
62
63 /*
64 * Usually it should pick farthest point, but picking most popular point seems to make search quicker anyway
65 */
vp_find_best_vantage_point_index(vp_sort_tmp indexes[],int num_indexes,const colormap_item items[])66 static int vp_find_best_vantage_point_index(vp_sort_tmp indexes[], int num_indexes, const colormap_item items[]) {
67 int best = 0;
68 float best_popularity = items[indexes[0].idx].popularity;
69 for(int i = 1; i < num_indexes; i++) {
70 if (items[indexes[i].idx].popularity > best_popularity) {
71 best_popularity = items[indexes[i].idx].popularity;
72 best = i;
73 }
74 }
75 return best;
76 }
77
vp_create_node(mempoolptr * m,vp_sort_tmp indexes[],int num_indexes,const colormap_item items[])78 static vp_node *vp_create_node(mempoolptr *m, vp_sort_tmp indexes[], int num_indexes, const colormap_item items[]) {
79 if (num_indexes <= 0) {
80 return NULL;
81 }
82
83 vp_node *node = mempool_alloc(m, sizeof(node[0]), 0);
84
85 if (num_indexes == 1) {
86 *node = (vp_node){
87 .vantage_point = items[indexes[0].idx].acolor,
88 .idx = indexes[0].idx,
89 .radius = MAX_DIFF,
90 .radius_squared = MAX_DIFF,
91 };
92 return node;
93 }
94
95 const int ref = vp_find_best_vantage_point_index(indexes, num_indexes, items);
96 const int ref_idx = indexes[ref].idx;
97
98 // Removes the `ref_idx` item from remaining items, because it's included in the current node
99 num_indexes -= 1;
100 indexes[ref] = indexes[num_indexes];
101
102 vp_sort_indexes_by_distance(items[ref_idx].acolor, indexes, num_indexes, items);
103
104 // Remaining items are split by the median distance
105 const int half_idx = num_indexes/2;
106
107 *node = (vp_node){
108 .vantage_point = items[ref_idx].acolor,
109 .idx = ref_idx,
110 .radius = sqrtf(indexes[half_idx].distance_squared),
111 .radius_squared = indexes[half_idx].distance_squared,
112 };
113 if (num_indexes < 7) {
114 node->rest = mempool_alloc(m, sizeof(node->rest[0]) * num_indexes, 0);
115 node->restcount = num_indexes;
116 for(int i=0; i < num_indexes; i++) {
117 node->rest[i].idx = indexes[i].idx;
118 node->rest[i].color = items[indexes[i].idx].acolor;
119 }
120 } else {
121 node->near = vp_create_node(m, indexes, half_idx, items);
122 node->far = vp_create_node(m, &indexes[half_idx], num_indexes - half_idx, items);
123 }
124
125 return node;
126 }
127
nearest_init(const colormap * map)128 LIQ_PRIVATE struct nearest_map *nearest_init(const colormap *map) {
129 mempoolptr m = NULL;
130 struct nearest_map *handle = mempool_create(&m, sizeof(handle[0]), sizeof(handle[0]) + sizeof(vp_node)*map->colors+16, map->malloc, map->free);
131
132 LIQ_ARRAY(vp_sort_tmp, indexes, map->colors);
133
134 for(unsigned int i=0; i < map->colors; i++) {
135 indexes[i].idx = i;
136 }
137
138 vp_node *root = vp_create_node(&m, indexes, map->colors, map->palette);
139 *handle = (struct nearest_map){
140 .root = root,
141 .palette = map->palette,
142 .mempool = m,
143 };
144
145 for(unsigned int i=0; i < map->colors; i++) {
146 vp_search_tmp best = {
147 .distance = MAX_DIFF,
148 .distance_squared = MAX_DIFF,
149 .exclude = i,
150 };
151 vp_search_node(root, &map->palette[i].acolor, &best);
152 handle->nearest_other_color_dist[i] = best.distance * best.distance / 4.0; // half of squared distance
153 }
154
155 return handle;
156 }
157
vp_search_node(const vp_node * node,const f_pixel * const needle,vp_search_tmp * const best_candidate)158 static void vp_search_node(const vp_node *node, const f_pixel *const needle, vp_search_tmp *const best_candidate) {
159 do {
160 const float distance_squared = colordifference(node->vantage_point, *needle);
161 const float distance = sqrtf(distance_squared);
162
163 if (distance_squared < best_candidate->distance_squared && best_candidate->exclude != node->idx) {
164 best_candidate->distance = distance;
165 best_candidate->distance_squared = distance_squared;
166 best_candidate->idx = node->idx;
167 }
168
169 if (node->restcount) {
170 for(int i=0; i < node->restcount; i++) {
171 const float distance_squared = colordifference(node->rest[i].color, *needle);
172 if (distance_squared < best_candidate->distance_squared && best_candidate->exclude != node->rest[i].idx) {
173 best_candidate->distance = sqrtf(distance_squared);
174 best_candidate->distance_squared = distance_squared;
175 best_candidate->idx = node->rest[i].idx;
176 }
177 }
178 return;
179 }
180
181 // Recurse towards most likely candidate first to narrow best candidate's distance as soon as possible
182 if (distance_squared < node->radius_squared) {
183 if (node->near) {
184 vp_search_node(node->near, needle, best_candidate);
185 }
186 // The best node (final answer) may be just ouside the radius, but not farther than
187 // the best distance we know so far. The vp_search_node above should have narrowed
188 // best_candidate->distance, so this path is rarely taken.
189 if (node->far && distance >= node->radius - best_candidate->distance) {
190 node = node->far; // Fast tail recursion
191 } else {
192 return;
193 }
194 } else {
195 if (node->far) {
196 vp_search_node(node->far, needle, best_candidate);
197 }
198 if (node->near && distance <= node->radius + best_candidate->distance) {
199 node = node->near; // Fast tail recursion
200 } else {
201 return;
202 }
203 }
204 } while(true);
205 }
206
nearest_search(const struct nearest_map * handle,const f_pixel * px,const int likely_colormap_index,float * diff)207 LIQ_PRIVATE unsigned int nearest_search(const struct nearest_map *handle, const f_pixel *px, const int likely_colormap_index, float *diff) {
208 const float guess_diff = colordifference(handle->palette[likely_colormap_index].acolor, *px);
209 if (guess_diff < handle->nearest_other_color_dist[likely_colormap_index]) {
210 if (diff) *diff = guess_diff;
211 return likely_colormap_index;
212 }
213
214 vp_search_tmp best_candidate = {
215 .distance = sqrtf(guess_diff),
216 .distance_squared = guess_diff,
217 .idx = likely_colormap_index,
218 .exclude = -1,
219 };
220 vp_search_node(handle->root, px, &best_candidate);
221 if (diff) {
222 *diff = best_candidate.distance * best_candidate.distance;
223 }
224 return best_candidate.idx;
225 }
226
nearest_free(struct nearest_map * centroids)227 LIQ_PRIVATE void nearest_free(struct nearest_map *centroids)
228 {
229 mempool_destroy(centroids->mempool);
230 }
231