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 */
34 compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image)
36 bzero(sig, sizeof(*sig));
37 ASSERT((image->flags & IMAGE_PIXEL_FORMAT) == COLOR_SPACE_RGB);
38 uns cols = image->cols;
39 uns rows = image->rows;
40 uns row_size = image->row_size;
42 uns w = (cols + 3) >> 2;
43 uns h = (rows + 3) >> 2;
45 DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h);
47 uns blocks_count = w * h;
48 struct image_sig_block *blocks = xmalloc(blocks_count * sizeof(struct image_sig_block)), *block = blocks;
50 /* Every block of 4x4 pixels */
51 byte *row_start = image->pixels;
52 for (uns block_y = 0; block_y < h; block_y++, row_start += row_size * 4)
55 for (uns block_x = 0; block_x < w; block_x++, p += 12, block++)
57 int t[16], s[16], *tp = t;
61 /* Convert pixels to Luv color space and compute average coefficients */
66 if ((!(cols & 3) || block_x < w - 1) && (!(rows & 3) || block_y < h - 1))
68 for (uns y = 0; y < 4; y++, p2 += row_size - 12)
69 for (uns x = 0; x < 4; x++, p2 += 3)
72 srgb_to_luv_pixel(luv, p2);
73 l_sum += *tp++ = luv[0];
78 block->v[0] = (l_sum >> 4);
79 block->v[1] = (u_sum >> 4);
80 block->v[2] = (v_sum >> 4);
82 /* Incomplete square near the edge */
86 uns square_cols = (block_x < w - 1 || !(cols & 3)) ? 4 : cols & 3;
87 uns square_rows = (block_y < h - 1 || !(rows & 3)) ? 4 : rows & 3;
88 for (y = 0; y < square_rows; y++, p2 += row_size)
91 for (x = 0; x < square_cols; x++, p3 += 3)
94 srgb_to_luv_pixel(luv, p3);
95 l_sum += *tp++ = luv[0];
101 *tp = tp[-square_cols];
106 for (x = 0; x < 4; x++)
108 *tp = tp[-square_rows * 4];
111 block->area = square_cols * square_rows;
112 uns div = 0x10000 / block->area;
113 block->v[0] = (l_sum * div) >> 16;
114 block->v[1] = (u_sum * div) >> 16;
115 block->v[2] = (v_sum * div) >> 16;
118 /* Apply Daubechies wavelet transformation */
120 # define DAUB_0 31651 /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
121 # define DAUB_1 54822 /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
122 # define DAUB_2 14689 /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
123 # define DAUB_3 -8481 /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
125 /* ... to the rows */
127 for (i = 0; i < 16; i += 4)
129 s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
130 s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
131 s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
132 s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
135 /* ... and to the columns... skip LL band */
136 for (i = 0; i < 2; i++)
138 t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
139 t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
143 t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x2000;
144 t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x2000;
145 t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
146 t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
149 /* Extract energies in LH, HL and HH bands */
150 block->v[3] = fast_sqrt_u16(isqr(t[8]) + isqr(t[9]) + isqr(t[12]) + isqr(t[13]));
151 block->v[4] = fast_sqrt_u16(isqr(t[2]) + isqr(t[3]) + isqr(t[6]) + isqr(t[7]));
152 block->v[5] = fast_sqrt_u16(isqr(t[10]) + isqr(t[11]) + isqr(t[14]) + isqr(t[15]));
156 /* FIXME: simple average is for testing pusposes only */
163 for (uns i = 0; i < blocks_count; i++)
165 l_sum += blocks[i].v[0];
166 u_sum += blocks[i].v[1];
167 v_sum += blocks[i].v[2];
168 lh_sum += blocks[i].v[3];
169 hl_sum += blocks[i].v[4];
170 hh_sum += blocks[i].v[5];
173 sig->vec.f[0] = l_sum / blocks_count;
174 sig->vec.f[1] = u_sum / blocks_count;
175 sig->vec.f[2] = v_sum / blocks_count;
176 sig->vec.f[3] = lh_sum / blocks_count;
177 sig->vec.f[4] = hl_sum / blocks_count;
178 sig->vec.f[5] = hh_sum / blocks_count;
180 if (cols < image_sig_min_width || rows < image_sig_min_height)
186 /* Quantize blocks to image regions */
187 struct image_sig_region regions[IMAGE_REG_MAX];
188 sig->len = image_sig_segmentation(blocks, blocks_count, regions);
190 /* For each region */
192 uns w_border = (MIN(w, h) + 3) / 4;
193 uns w_mul = 127 * 256 / w_border;
194 for (uns i = 0; i < sig->len; i++)
196 struct image_sig_region *r = regions + i;
197 DBG("Processing region %u: count=%u", i, r->count);
200 /* Copy texture properties */
201 sig->reg[i].f[0] = r->a[0];
202 sig->reg[i].f[1] = r->a[1];
203 sig->reg[i].f[2] = r->a[2];
204 sig->reg[i].f[3] = r->a[3];
205 sig->reg[i].f[4] = r->a[4];
206 sig->reg[i].f[5] = r->a[5];
208 /* Compute coordinates centroid and region weight */
209 u64 x_avg = 0, y_avg = 0, w_sum = 0;
210 for (struct image_sig_block *b = r->blocks; b; b = b->next)
216 d = MIN(d, w - b->x - 1);
217 d = MIN(d, h - b->y - 1);
221 w_sum += 128 + (d - w_border) * w_mul / 256;
227 DBG(" centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
229 /* Compute normalized inertia */
230 u64 sum1 = 0, sum2 = 0, sum3 = 0;
231 for (struct image_sig_block *b = r->blocks; b; b = b->next)
233 uns inc2 = isqr(x_avg - b->x) + isqr(y_avg - b->y);
234 uns inc1 = sqrt(inc2);
239 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);
240 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);
241 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);
245 /* Compute average differences */
256 for (uns i = 0; i < sig->len; i++)
257 for (uns j = i + 1; j < sig->len; j++)
260 for (uns k = 0; k < IMAGE_REG_F; k++)
261 d += isqr(sig->reg[i].f[k] - sig->reg[j].f[k]);
264 for (uns k = 0; k < IMAGE_REG_H; k++)
265 d += isqr(sig->reg[i].h[k] - sig->reg[j].h[k]);
269 sig->df = CLAMP(df / cnt, 1, 255);
270 sig->dh = CLAMP(dh / cnt, 1, 65535);
272 DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
274 /* Compute normalized weights */
275 uns wa = 128, wb = 128;
276 for (uns i = sig->len; --i > 0; )
278 struct image_sig_region *r = regions + i;
279 wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
280 wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
285 /* Dump regions features */
287 for (uns i = 0; i < sig->len; i++)
289 byte buf[IMAGE_REGION_DUMP_MAX];
290 image_region_dump(buf, sig->reg + i);
291 DBG("region %u: features=%s", i, buf);