- /* Compute new averages */
- for (r = regions; r != regions_end; r++)
- if (r->count)
- {
- r->l = r->sum_l / r->count;
- r->u = r->sum_u / r->count;
- r->v = r->sum_v / r->count;
- r->lh = r->sum_lh / r->count;
- r->hl = r->sum_hl / r->count;
- r->hh = r->sum_hh / r->count;
- }
-
- if (!conv_i)
- break; // FIXME: convergation criteria
- }
-
- /* Remove empty regions */
- struct region *r2 = regions;
- for (r = regions; r != regions_end; r++)
- if (r->count)
- *r2++ = *r;
- return r2 - regions;
-}
-
-int
-compute_image_signature(struct image_thread *thread, struct image_signature *sig, struct image *image)
-{
- uns cols = image->cols;
- uns rows = image->rows;
-
- /* FIXME: deal with smaller images */
- if (cols < 4 || rows < 4)
- {
- image_thread_err_format(thread, IMAGE_ERR_INVALID_DIMENSIONS, "Image too small... %ux%u", cols, rows);
- return 0;
- }
-
- uns w = cols >> 2;
- uns h = rows >> 2;
- DBG("Computing signature for image %dx%d... %dx%d blocks", cols, rows, w, h);
- uns blocks_count = w * h;
- struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks; /* FIXME: use mempool */
-
- /* Every block of 4x4 pixels (FIXME: deal with smaller blocks near the edges) */
- byte *p = image->pixels;
- for (uns block_y = 0; block_y < h; block_y++, p += 3 * ((cols & 3) + cols * 3))
- for (uns block_x = 0; block_x < w; block_x++, p -= 3 * (4 * cols - 4), block++)
- {
- block->x = block_x;
- block->y = block_y;
-
- int t[16], s[16], *tp = t;
+ /* Apply Daubechies wavelet transformation */
+
+# define DAUB_0 31651 /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
+# define DAUB_1 54822 /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
+# define DAUB_2 14689 /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
+# define DAUB_3 -8481 /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
+
+ /* ... to the rows */
+ uns i;
+ for (i = 0; i < 16; i += 4)
+ {
+ s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
+ s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
+ s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
+ s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
+ }