2 * Image Library -- Image segmentation
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.
15 #include <images/images.h>
16 #include <images/signature.h>
17 #include <images/math.h>
23 dump_segmentation(struct image_sig_region *regions, uns regions_count)
25 uns cols = 0, rows = 0;
26 for (uns i = 0; i < regions_count; i++)
27 for (struct image_sig_block *b = regions[i].blocks; b; b = b->next)
29 cols = MAX(cols, b->x + 1);
30 rows = MAX(rows, b->y + 1);
32 uns size = (cols + 1) * rows;
35 for (uns i = 0; i < regions_count; i++)
37 byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
38 for (struct image_sig_block *b = regions[i].blocks; b; b = b->next)
39 buf[b->x + b->y * (cols + 1)] = c;
41 for (uns i = 0; i < rows; i++)
42 log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
46 /* Pre-quantization - recursively split groups of blocks with large error */
49 prequant_init_region(struct image_sig_region *region)
51 bzero(region, sizeof(*region));
55 prequant_add_block(struct image_sig_region *region, struct image_sig_block *block)
57 block->next = region->blocks;
58 region->blocks = block;
60 for (uns i = 0; i < IMAGE_VEC_F; i++)
62 region->b[i] += block->v[i];
63 region->c[i] += isqr(block->v[i]);
68 prequant_finish_region(struct image_sig_region *region)
70 if (region->count < 2)
78 for (uns i = 0; i < IMAGE_VEC_F; i++)
80 region->e += region->c[i];
81 a += (u64)region->b[i] * region->b[i];
83 region->e -= a / region->count;
84 DBG("Finished region %u", (uns)region->e / region->count);
89 prequant_heap_cmp(struct image_sig_region *a, struct image_sig_region *b)
94 #define ASORT_PREFIX(x) prequant_##x
95 #define ASORT_KEY_TYPE uns
96 #include <ucw/sorter/array-simple.h>
99 prequant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions)
101 DBG("Starting pre-quantization");
103 uns regions_count, heap_count, axis;
104 struct image_sig_block *blocks_end = blocks + blocks_count, *block, *block2;
105 struct image_sig_region *heap[IMAGE_REG_MAX + 1], *region, *region2;
107 /* Initialize single region with all blocks */
108 regions_count = heap_count = 1;
110 prequant_init_region(regions);
111 for (block = blocks; block != blocks_end; block++)
112 prequant_add_block(regions, block);
113 prequant_finish_region(regions);
116 while (regions_count < IMAGE_REG_MAX &&
117 regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count)
120 DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u",
121 regions_count, heap_count, region->count, (uns)region->e);
122 if (region->count < 2 ||
123 region->e < image_sig_prequant_thresholds[regions_count - 1] * blocks_count)
125 HEAP_DELETE_MIN(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
129 /* Select axis to split - the one with maximum average quadratic error */
131 u64 cov = (u64)region->count * region->c[0] - (u64)region->b[0] * region->b[0];
132 for (uns i = 1; i < 6; i++)
134 uns j = (u64)region->count * region->c[i] - (u64)region->b[i] * region->b[i];
141 DBG("Splitting axis %u with average quadratic error %u", axis, (uns)(cov / (region->count * region->count)));
143 /* Sort values on the split axis */
144 uns val[256], cnt[256], cval;
145 if (region->count > 64)
147 bzero(cnt, sizeof(cnt));
148 for (block = region->blocks; block; block = block->next)
149 cnt[block->v[axis]]++;
151 for (uns i = 0; i < 256; i++)
161 block = region->blocks;
162 for (uns i = 0; i < region->count; i++, block = block->next)
163 val[i] = block->v[axis];
164 prequant_sort(val, region->count);
167 for (uns i = 1; i < region->count; i++)
168 if (val[i] == val[cval - 1])
178 /* Select split value - to minimize error */
179 uns b1 = val[0] * cnt[0];
180 uns c1 = isqr(val[0]) * cnt[0];
181 uns b2 = region->b[axis] - b1;
182 uns c2 = region->c[axis] - c1;
183 uns i = cnt[0], j = region->count - cnt[0];
184 u64 best_err = c1 - (u64)b1 * b1 / i + c2 - (u64)b2 * b2 / j;
185 uns split_val = val[0];
186 for (uns k = 1; k < cval - 1; k++)
188 uns b0 = val[k] * cnt[k];
189 uns c0 = isqr(val[k]) * cnt[k];
196 u64 err = (u64)c1 - (u64)b1 * b1 / i + (u64)c2 - (u64)b2 * b2 / j;
203 DBG("split_val=%u best_err=%llu b[axis]=%u c[axis]=%u", split_val, (long long)best_err, region->b[axis], region->c[axis]);
206 block = region->blocks;
207 region2 = regions + regions_count++;
208 prequant_init_region(region);
209 prequant_init_region(region2);
212 block2 = block->next;
213 if (block->v[axis] <= split_val)
214 prequant_add_block(region, block);
216 prequant_add_block(region2, block);
219 prequant_finish_region(region);
220 prequant_finish_region(region2);
221 HEAP_INCREASE(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1, region);
222 HEAP_INSERT(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, region2);
225 DBG("Pre-quantized to %u regions", regions_count);
227 return regions_count;
230 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
233 postquant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions, uns regions_count)
235 DBG("Starting post-quantization");
237 struct image_sig_block *blocks_end = blocks + blocks_count, *block;
238 struct image_sig_region *regions_end = regions + regions_count, *region;
239 uns error = 0, last_error;
241 /* Initialize regions and initial segmentation error */
242 for (region = regions; region != regions_end; )
244 uns inv = 0xffffffffU / region->count;
245 for (uns i = 0; i < IMAGE_VEC_F; i++)
247 region->a[i] = ((u64)region->b[i] * inv) >> 32;
248 error += region->c[i] - region->a[i] * region->b[i];
253 /* Convergation cycle */
254 for (uns step = 0; step < image_sig_postquant_max_steps; step++)
259 for (region = regions; region != regions_end; region++)
261 region->blocks = NULL;
263 bzero(region->b, sizeof(region->b));
264 bzero(region->c, sizeof(region->c));
267 /* Assign each block to its nearest pivot and accumulate region variables */
268 for (block = blocks; block != blocks_end; block++)
270 struct image_sig_region *best_region = NULL;
272 for (region = regions; region != regions_end; region++)
275 isqr(block->v[0] - region->a[0]) +
276 isqr(block->v[1] - region->a[1]) +
277 isqr(block->v[2] - region->a[2]) +
278 isqr(block->v[3] - region->a[3]) +
279 isqr(block->v[4] - region->a[4]) +
280 isqr(block->v[5] - region->a[5]);
281 if (dist <= best_dist)
284 best_region = region;
287 region = best_region;
289 block->next = region->blocks;
290 region->blocks = block;
291 for (uns i = 0; i < IMAGE_VEC_F; i++)
293 region->b[i] += block->v[i];
294 region->c[i] += isqr(block->v[i]);
298 /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
301 for (region = regions; region != regions_end; )
304 uns inv = 0xffffffffU / region->count;
305 for (uns i = 0; i < IMAGE_VEC_F; i++)
307 region->a[i] = ((u64)region->b[i] * inv) >> 32;
308 error += region->c[i] - region->a[i] * region->b[i];
315 *region = *regions_end;
318 DBG("last_error=%u error=%u", last_error, error);
320 /* Convergation criteria */
321 if (step >= image_sig_postquant_min_steps)
323 if (error > last_error)
325 u64 dif = last_error - error;
326 if (dif * image_sig_postquant_threshold < (u64)last_error * 100)
331 DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
333 return regions_end - regions;
337 image_sig_segmentation(struct image_sig_data *data)
339 data->regions_count = prequant(data->blocks, data->blocks_count, data->regions);
341 dump_segmentation(data->regions, data->regions_count);
343 data->regions_count = postquant(data->blocks, data->blocks_count, data->regions, data->regions_count);
345 dump_segmentation(data->regions, data->regions_count);