2 * Image Library -- Computation of image signatures
4 * (c) 2006 Pavel Charvat <pchar@ucw.cz>
6 * This software may be freely distributed and used according to the terms
7 * of the GNU Lesser General Public License.
12 #include "sherlock/sherlock.h"
14 #include "lib/fastbuf.h"
15 #include "images/images.h"
16 #include "images/color.h"
17 #include "images/signature.h"
20 static double image_sig_inertia_scale[3] = { 3, 1, 0.3 };
23 u32 area; /* block area in pixels (usually 16) */
24 u32 l, u, v; /* average Luv coefficients */
25 u32 lh, hl, hh; /* energies in Daubechies wavelet bands */
26 u32 x, y; /* block position */
33 u32 sum_l, sum_u, sum_v;
34 u32 sum_lh, sum_hl, sum_hh;
49 dump_segmentation(struct region *regions, uns regions_count, uns cols, uns rows)
51 uns size = (cols + 1) * rows;
54 for (uns i = 0; i < regions_count; i++)
56 byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
57 for (struct block *b = regions[i].blocks; b; b = b->next)
58 buf[b->x + b->y * (cols + 1)] = c;
60 for (uns i = 0; i < rows; i++)
61 log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
67 compute_k_means(struct block *blocks, uns blocks_count, struct region *regions, uns regions_count)
69 ASSERT(regions_count <= blocks_count);
70 struct block *mean[IMAGE_REG_MAX], *b, *blocks_end = blocks + blocks_count;
71 struct region *r, *regions_end = regions + regions_count;
73 /* Select means_count random blocks as initial regions pivots */
74 if (regions_count <= blocks_count - regions_count)
76 for (b = blocks; b != blocks_end; b++)
78 for (uns i = 0; i < regions_count; )
80 uns j = random_max(blocks_count);
83 b->next = mean[i++] = b;
89 for (uns i = regions_count; i; j--)
90 if (random_max(j) <= i)
91 mean[--i] = blocks + j - 1;
94 for (uns i = 0; i < regions_count; i++, r++)
105 /* Convergation cycle */
106 for (uns conv_i = 8; ; conv_i--)
108 for (r = regions; r != regions_end; r++)
110 r->sum_l = r->sum_u = r->sum_v = r->sum_lh = r->sum_hl = r->sum_hh = r->count = 0;
114 /* Find nearest regions and accumulate averages */
115 for (b = blocks; b != blocks_end; b++)
118 struct region *best_r = NULL;
119 for (r = regions; r != regions_end; r++)
134 best_r->sum_l += b->l;
135 best_r->sum_u += b->u;
136 best_r->sum_v += b->v;
137 best_r->sum_lh += b->lh;
138 best_r->sum_hl += b->hl;
139 best_r->sum_hh += b->hh;
141 b->next = best_r->blocks;
145 /* Compute new averages */
146 for (r = regions; r != regions_end; r++)
149 r->l = r->sum_l / r->count;
150 r->u = r->sum_u / r->count;
151 r->v = r->sum_v / r->count;
152 r->lh = r->sum_lh / r->count;
153 r->hl = r->sum_hl / r->count;
154 r->hh = r->sum_hh / r->count;
158 break; // FIXME: convergation criteria
161 /* Remove empty regions */
162 struct region *r2 = regions;
163 for (r = regions; r != regions_end; r++)
170 compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image)
172 bzero(sig, sizeof(*sig));
173 ASSERT((image->flags & IMAGE_PIXEL_FORMAT) == COLOR_SPACE_RGB);
174 uns cols = image->cols;
175 uns rows = image->rows;
176 uns row_size = image->row_size;
178 uns w = (cols + 3) >> 2;
179 uns h = (rows + 3) >> 2;
181 DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h);
183 uns blocks_count = w * h;
184 struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks;
186 /* Every block of 4x4 pixels */
187 byte *row_start = image->pixels;
188 for (uns block_y = 0; block_y < h; block_y++, row_start += row_size * 4)
191 for (uns block_x = 0; block_x < w; block_x++, p += 12, block++)
193 int t[16], s[16], *tp = t;
197 /* Convert pixels to Luv color space and compute average coefficients */
202 if ((!(cols & 3) || block_x < w - 1) && (!(rows & 3) || block_y < h - 1))
204 for (uns y = 0; y < 4; y++, p2 += row_size - 12)
205 for (uns x = 0; x < 4; x++, p2 += 3)
208 srgb_to_luv_pixel(luv, p2);
209 l_sum += *tp++ = luv[0];
214 block->l = (l_sum >> 4);
215 block->u = (u_sum >> 4);
216 block->v = (v_sum >> 4);
218 /* Incomplete square near the edge */
222 uns square_cols = (block_x < w - 1) ? 4 : cols & 3;
223 uns square_rows = (block_y < h - 1) ? 4 : rows & 3;
224 for (y = 0; y < square_rows; y++, p2 += row_size)
227 for (x = 0; x < square_cols; x++, p3 += 3)
230 srgb_to_luv_pixel(luv, p3);
231 l_sum += *tp++ = luv[0];
237 *tp = tp[-square_cols];
242 for (x = 0; x < 4; x++)
244 *tp = tp[-square_rows * 4];
247 block->area = square_cols * square_rows;
248 uns div = 0x10000 / block->area;
249 block->l = (l_sum * div) >> 16;
250 block->u = (u_sum * div) >> 16;
251 block->v = (v_sum * div) >> 16;
254 /* Apply Daubechies wavelet transformation */
256 # define DAUB_0 31651 /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
257 # define DAUB_1 54822 /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
258 # define DAUB_2 14689 /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
259 # define DAUB_3 -8481 /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
261 /* ... to the rows */
263 for (i = 0; i < 16; i += 4)
265 s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
266 s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
267 s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
268 s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
271 /* ... and to the columns... skip LL band */
272 for (i = 0; i < 2; i++)
274 t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
275 t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
279 t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
280 t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
281 t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
282 t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
285 /* Extract energies in LH, HL and HH bands */
286 block->lh = CLAMP((int)(sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]) / 16), 0, 255);
287 block->hl = CLAMP((int)(sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]) / 16), 0, 255);
288 block->hh = CLAMP((int)(sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]) / 16), 0, 255);
292 /* FIXME: simple average is for testing pusposes only */
299 for (uns i = 0; i < blocks_count; i++)
301 l_sum += blocks[i].l;
302 u_sum += blocks[i].u;
303 v_sum += blocks[i].v;
304 lh_sum += blocks[i].lh;
305 hl_sum += blocks[i].hl;
306 hh_sum += blocks[i].hh;
309 sig->vec.f[0] = l_sum / blocks_count;
310 sig->vec.f[1] = u_sum / blocks_count;
311 sig->vec.f[2] = v_sum / blocks_count;
312 sig->vec.f[3] = lh_sum / blocks_count;
313 sig->vec.f[4] = hl_sum / blocks_count;
314 sig->vec.f[5] = hh_sum / blocks_count;
316 if (cols < image_sig_min_width || rows < image_sig_min_height)
319 /* Quantize blocks to image regions */
320 struct region regions[IMAGE_REG_MAX];
321 sig->len = compute_k_means(blocks, blocks_count, regions, MIN(blocks_count, IMAGE_REG_MAX));
323 /* For each region */
325 uns w_border = (MIN(w, h) + 3) / 4;
326 uns w_mul = 127 * 256 / w_border;
327 for (uns i = 0; i < sig->len; i++)
329 struct region *r = regions + i;
330 DBG("Processing region %u: count=%u", i, r->count);
333 /* Copy texture properties */
334 sig->reg[i].f[0] = r->l;
335 sig->reg[i].f[1] = r->u;
336 sig->reg[i].f[2] = r->v;
337 sig->reg[i].f[3] = r->lh;
338 sig->reg[i].f[4] = r->hl;
339 sig->reg[i].f[5] = r->hh;
341 /* Compute coordinates centroid and region weight */
342 u64 x_avg = 0, y_avg = 0, w_sum = 0;
343 for (struct block *b = r->blocks; b; b = b->next)
349 d = MIN(d, w - b->x - 1);
350 d = MIN(d, h - b->y - 1);
354 w_sum += 128 + (d - w_border) * w_mul / 256;
360 DBG(" centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
362 /* Compute normalized inertia */
363 u64 sum1 = 0, sum2 = 0, sum3 = 0;
364 for (struct block *b = r->blocks; b; b = b->next)
366 uns inc2 = dist(x_avg, b->x) + dist(y_avg, b->y);
367 uns inc1 = sqrt(inc2);
372 sig->reg[i].h[0] = CLAMP(image_sig_inertia_scale[0] * sum1 * ((3 * M_PI * M_PI) / 2) * pow(r->count, -1.5), 0, 65535);
373 sig->reg[i].h[1] = CLAMP(image_sig_inertia_scale[1] * sum2 * ((4 * M_PI * M_PI * M_PI) / 2) / ((u64)r->count * r->count), 0, 65535);
374 sig->reg[i].h[2] = CLAMP(image_sig_inertia_scale[2] * sum3 * ((5 * M_PI * M_PI * M_PI * M_PI) / 2) * pow(r->count, -2.5), 0, 65535);
378 /* Compute average differences */
389 for (uns i = 0; i < sig->len; i++)
390 for (uns j = i + 1; j < sig->len; j++)
393 for (uns k = 0; k < IMAGE_REG_F; k++)
394 d += dist(sig->reg[i].f[k], sig->reg[j].f[k]);
397 for (uns k = 0; k < IMAGE_REG_H; k++)
398 d += dist(sig->reg[i].h[k], sig->reg[j].h[k]);
402 sig->df = CLAMP(df / cnt, 1, 255);
403 sig->dh = CLAMP(dh / cnt, 1, 65535);
405 DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
407 /* Compute normalized weights */
408 uns wa = 128, wb = 128;
409 for (uns i = sig->len; --i > 0; )
411 struct region *r = regions + i;
412 wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
413 wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
418 /* Dump regions features */
420 for (uns i = 0; i < sig->len; i++)
422 byte buf[IMAGE_REGION_DUMP_MAX];
423 image_region_dump(buf, sig->reg + i);
424 DBG("region %u: features=%s", i, buf);
426 dump_segmentation(regions, sig->len, w, h);