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"
17 #include "images/math.h"
18 #include "images/images.h"
19 #include "images/color.h"
20 #include "images/signature.h"
24 static double image_sig_inertia_scale[3] = { 3, 1, 0.3 };
27 u32 area; /* block area in pixels (usually 16) */
29 u32 x, y; /* block position */
52 dump_segmentation(struct region *regions, uns regions_count, uns cols, uns rows)
54 uns size = (cols + 1) * rows;
57 for (uns i = 0; i < regions_count; i++)
59 byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
60 for (struct block *b = regions[i].blocks; b; b = b->next)
61 buf[b->x + b->y * (cols + 1)] = c;
63 for (uns i = 0; i < rows; i++)
64 log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
68 /* Pre-quantization - recursively split groups of blocks with large error */
71 prequant_init_region(struct region *region)
73 bzero(region, sizeof(*region));
77 prequant_add_block(struct region *region, struct block *block)
79 block->next = region->blocks;
80 region->blocks = block;
82 for (uns i = 0; i < IMAGE_VEC_F; i++)
84 region->b[i] += block->v[i];
85 region->c[i] += isqr(block->v[i]);
90 prequant_finish_region(struct region *region)
92 if (region->count < 2)
93 memcpy(region->c, region->a, sizeof(region->c));
98 for (uns i = 0; i < IMAGE_VEC_F; i++)
100 region->e += region->c[i];
101 a += (u64)region->b[i] * region->b[i];
103 region->e -= a / region->count;
108 prequant_heap_cmp(struct region *a, struct region *b)
113 #define ASORT_PREFIX(x) prequant_##x
114 #define ASORT_KEY_TYPE u32
115 #define ASORT_ELT(i) val[i]
116 #define ASORT_EXTRA_ARGS , u32 *val
117 #include "lib/arraysort.h"
120 prequant(struct block *blocks, uns blocks_count, struct region *regions)
122 DBG("Starting pre-quantization");
124 uns regions_count, heap_count, axis, cov;
125 struct block *blocks_end = blocks + blocks_count, *block, *block2;
126 struct region *heap[IMAGE_REG_MAX + 1], *region, *region2;
128 /* Initialize single region with all blocks */
129 regions_count = heap_count = 1;
131 prequant_init_region(regions);
132 for (block = blocks; block != blocks_end; block++)
133 prequant_add_block(regions, block);
134 prequant_finish_region(regions);
137 while (regions_count < IMAGE_REG_MAX &&
138 regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count)
141 DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u",
142 regions_count, heap_count, region->count, (uns)region->e);
143 if (region->count < 2 ||
144 region->e < image_sig_prequant_thresholds[regions_count - 1] * region->count)
146 HEAP_DELMIN(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
150 /* Select axis to split - the one with maximum covariance */
152 cov = region->count * region->c[0] - region->b[0];
153 for (uns i = 1; i < 6; i++)
155 uns j = region->count * region->c[i] - region->b[i];
162 DBG("Splitting axis %u with covariance %u", axis, cov / region->count);
164 /* Sort values on the split axis */
165 u32 val[region->count];
166 block = region->blocks;
167 for (uns i = 0; i < region->count; i++, block = block->next)
168 val[i] = block->v[axis];
169 prequant_sort(region->count, val);
171 /* Select split value - to minimize error */
172 uns b1 = 0, c1 = 0, b2 = region->b[axis], c2 = region->c[axis];
173 uns i = 0, j = region->count, best_v = 0;
174 u64 best_err = 0xffffffffffffffff;
175 while (i < region->count)
177 u64 err = (u64)i * c1 - (u64)b1 * b1 + (u64)j * c2 - (u64)b2 * b2;
183 uns sqr = isqr(val[i]);
191 uns split_val = best_v;
192 DBG("split_val=%u best_err=%Lu b[axis]=%u c[axis]=%u", split_val, (long long)best_err, region->b[axis], region->c[axis]);
195 block = region->blocks;
196 region2 = regions + regions_count++;
197 prequant_init_region(region);
198 prequant_init_region(region2);
201 block2 = block->next;
202 if (block->v[axis] < split_val)
203 prequant_add_block(region, block);
205 prequant_add_block(region2, block);
208 prequant_finish_region(region);
209 prequant_finish_region(region2);
210 HEAP_INCREASE(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1);
211 heap[++heap_count] = region2;
212 HEAP_INSERT(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
215 DBG("Pre-quantized to %u regions", regions_count);
217 return regions_count;
220 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
223 postquant(struct block *blocks, uns blocks_count, struct region *regions, uns regions_count)
225 DBG("Starting post-quantization");
227 struct block *blocks_end = blocks + blocks_count, *block;
228 struct region *regions_end = regions + regions_count, *region;
229 uns error = 0, last_error;
231 /* Initialize regions and initial segmentation error */
232 for (region = regions; region != regions_end; )
234 uns inv = 0xffffffffU / region->count;
235 for (uns i = 0; i < IMAGE_VEC_F; i++)
237 region->a[i] = ((u64)region->b[i] * inv) >> 32;
238 error += region->c[i] - region->a[i] * region->b[i];
243 /* Convergation cycle */
244 for (uns step = 0; step < image_sig_postquant_max_steps; step++)
249 for (region = regions; region != regions_end; region++)
251 region->blocks = NULL;
253 bzero(region->b, sizeof(region->b));
254 bzero(region->c, sizeof(region->c));
257 /* Assign each block to its nearest pivot and accumulate region variables */
258 for (block = blocks; block != blocks_end; block++)
260 struct region *best_region = NULL;
262 for (region = regions; region != regions_end; region++)
265 isqr(block->v[0] - region->a[0]) +
266 isqr(block->v[1] - region->a[1]) +
267 isqr(block->v[2] - region->a[2]) +
268 isqr(block->v[3] - region->a[3]) +
269 isqr(block->v[4] - region->a[4]) +
270 isqr(block->v[5] - region->a[5]);
271 if (dist <= best_dist)
274 best_region = region;
277 region = best_region;
279 block->next = region->blocks;
280 region->blocks = block;
281 for (uns i = 0; i < IMAGE_VEC_F; i++)
283 region->b[i] += block->v[i];
284 region->c[i] += isqr(block->v[i]);
288 /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
291 for (region = regions; region != regions_end; )
294 uns inv = 0xffffffffU / region->count;
295 for (uns i = 0; i < IMAGE_VEC_F; i++)
297 region->a[i] = ((u64)region->b[i] * inv) >> 32;
298 error += region->c[i] - region->a[i] * region->b[i];
305 *region = *regions_end;
308 DBG("last_error=%u error=%u", last_error, error);
310 /* Convergation criteria */
311 if (step >= image_sig_postquant_min_steps)
313 if (error > last_error)
315 u64 dif = last_error - error;
316 if (dif * image_sig_postquant_threshold < last_error * 100)
321 DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
323 return regions_end - regions;
327 segmentation(struct block *blocks, uns blocks_count, struct region *regions, uns width UNUSED, uns height UNUSED)
330 regions_count = prequant(blocks, blocks_count, regions);
332 dump_segmentation(regions, regions_count, width, height);
334 regions_count = postquant(blocks, blocks_count, regions, regions_count);
336 dump_segmentation(regions, regions_count, width, height);
338 return regions_count;
342 compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image)
344 bzero(sig, sizeof(*sig));
345 ASSERT((image->flags & IMAGE_PIXEL_FORMAT) == COLOR_SPACE_RGB);
346 uns cols = image->cols;
347 uns rows = image->rows;
348 uns row_size = image->row_size;
350 uns w = (cols + 3) >> 2;
351 uns h = (rows + 3) >> 2;
353 DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h);
355 uns blocks_count = w * h;
356 struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks;
358 /* Every block of 4x4 pixels */
359 byte *row_start = image->pixels;
360 for (uns block_y = 0; block_y < h; block_y++, row_start += row_size * 4)
363 for (uns block_x = 0; block_x < w; block_x++, p += 12, block++)
365 int t[16], s[16], *tp = t;
369 /* Convert pixels to Luv color space and compute average coefficients */
374 if ((!(cols & 3) || block_x < w - 1) && (!(rows & 3) || block_y < h - 1))
376 for (uns y = 0; y < 4; y++, p2 += row_size - 12)
377 for (uns x = 0; x < 4; x++, p2 += 3)
380 srgb_to_luv_pixel(luv, p2);
381 l_sum += *tp++ = luv[0];
386 block->v[0] = (l_sum >> 4);
387 block->v[1] = (u_sum >> 4);
388 block->v[2] = (v_sum >> 4);
390 /* Incomplete square near the edge */
394 uns square_cols = (block_x < w - 1 || !(cols & 3)) ? 4 : cols & 3;
395 uns square_rows = (block_y < h - 1 || !(rows & 3)) ? 4 : rows & 3;
396 for (y = 0; y < square_rows; y++, p2 += row_size)
399 for (x = 0; x < square_cols; x++, p3 += 3)
402 srgb_to_luv_pixel(luv, p3);
403 l_sum += *tp++ = luv[0];
409 *tp = tp[-square_cols];
414 for (x = 0; x < 4; x++)
416 *tp = tp[-square_rows * 4];
419 block->area = square_cols * square_rows;
420 uns div = 0x10000 / block->area;
421 block->v[0] = (l_sum * div) >> 16;
422 block->v[1] = (u_sum * div) >> 16;
423 block->v[2] = (v_sum * div) >> 16;
426 /* Apply Daubechies wavelet transformation */
428 # define DAUB_0 31651 /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
429 # define DAUB_1 54822 /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
430 # define DAUB_2 14689 /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
431 # define DAUB_3 -8481 /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
433 /* ... to the rows */
435 for (i = 0; i < 16; i += 4)
437 s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
438 s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
439 s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
440 s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
443 /* ... and to the columns... skip LL band */
444 for (i = 0; i < 2; i++)
446 t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
447 t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
451 t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x2000;
452 t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x2000;
453 t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
454 t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
457 /* Extract energies in LH, HL and HH bands */
458 block->v[3] = fast_sqrt_u16(isqr(t[8]) + isqr(t[9]) + isqr(t[12]) + isqr(t[13]));
459 block->v[4] = fast_sqrt_u16(isqr(t[2]) + isqr(t[3]) + isqr(t[6]) + isqr(t[7]));
460 block->v[5] = fast_sqrt_u16(isqr(t[10]) + isqr(t[11]) + isqr(t[14]) + isqr(t[15]));
464 /* FIXME: simple average is for testing pusposes only */
471 for (uns i = 0; i < blocks_count; i++)
473 l_sum += blocks[i].v[0];
474 u_sum += blocks[i].v[1];
475 v_sum += blocks[i].v[2];
476 lh_sum += blocks[i].v[3];
477 hl_sum += blocks[i].v[4];
478 hh_sum += blocks[i].v[5];
481 sig->vec.f[0] = l_sum / blocks_count;
482 sig->vec.f[1] = u_sum / blocks_count;
483 sig->vec.f[2] = v_sum / blocks_count;
484 sig->vec.f[3] = lh_sum / blocks_count;
485 sig->vec.f[4] = hl_sum / blocks_count;
486 sig->vec.f[5] = hh_sum / blocks_count;
488 if (cols < image_sig_min_width || rows < image_sig_min_height)
494 /* Quantize blocks to image regions */
495 struct region regions[IMAGE_REG_MAX];
496 sig->len = segmentation(blocks, blocks_count, regions, w, h);
498 /* For each region */
500 uns w_border = (MIN(w, h) + 3) / 4;
501 uns w_mul = 127 * 256 / w_border;
502 for (uns i = 0; i < sig->len; i++)
504 struct region *r = regions + i;
505 DBG("Processing region %u: count=%u", i, r->count);
508 /* Copy texture properties */
509 sig->reg[i].f[0] = r->a[0];
510 sig->reg[i].f[1] = r->a[1];
511 sig->reg[i].f[2] = r->a[2];
512 sig->reg[i].f[3] = r->a[3];
513 sig->reg[i].f[4] = r->a[4];
514 sig->reg[i].f[5] = r->a[5];
516 /* Compute coordinates centroid and region weight */
517 u64 x_avg = 0, y_avg = 0, w_sum = 0;
518 for (struct block *b = r->blocks; b; b = b->next)
524 d = MIN(d, w - b->x - 1);
525 d = MIN(d, h - b->y - 1);
529 w_sum += 128 + (d - w_border) * w_mul / 256;
535 DBG(" centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
537 /* Compute normalized inertia */
538 u64 sum1 = 0, sum2 = 0, sum3 = 0;
539 for (struct block *b = r->blocks; b; b = b->next)
541 uns inc2 = dist(x_avg, b->x) + dist(y_avg, b->y);
542 uns inc1 = sqrt(inc2);
547 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);
548 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);
549 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);
553 /* Compute average differences */
564 for (uns i = 0; i < sig->len; i++)
565 for (uns j = i + 1; j < sig->len; j++)
568 for (uns k = 0; k < IMAGE_REG_F; k++)
569 d += dist(sig->reg[i].f[k], sig->reg[j].f[k]);
572 for (uns k = 0; k < IMAGE_REG_H; k++)
573 d += dist(sig->reg[i].h[k], sig->reg[j].h[k]);
577 sig->df = CLAMP(df / cnt, 1, 255);
578 sig->dh = CLAMP(dh / cnt, 1, 65535);
580 DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
582 /* Compute normalized weights */
583 uns wa = 128, wb = 128;
584 for (uns i = sig->len; --i > 0; )
586 struct region *r = regions + i;
587 wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
588 wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
593 /* Dump regions features */
595 for (uns i = 0; i < sig->len; i++)
597 byte buf[IMAGE_REGION_DUMP_MAX];
598 image_region_dump(buf, sig->reg + i);
599 DBG("region %u: features=%s", i, buf);