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.
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"
26 dump_segmentation(struct image_sig_region *regions, uns regions_count)
28 uns cols = 0, rows = 0;
29 for (uns i = 0; i < regions_count; i++)
30 for (struct image_sig_block *b = regions[i].blocks; b; b = b->next)
32 cols = MAX(cols, b->x);
33 rows = MAX(rows, b->y);
35 uns size = (cols + 1) * rows;
38 for (uns i = 0; i < regions_count; i++)
40 byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
41 for (struct image_sig_block *b = regions[i].blocks; b; b = b->next)
42 buf[b->x + b->y * (cols + 1)] = c;
44 for (uns i = 0; i < rows; i++)
45 log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
49 /* Pre-quantization - recursively split groups of blocks with large error */
52 prequant_init_region(struct image_sig_region *region)
54 bzero(region, sizeof(*region));
58 prequant_add_block(struct image_sig_region *region, struct image_sig_block *block)
60 block->next = region->blocks;
61 region->blocks = block;
63 for (uns i = 0; i < IMAGE_VEC_F; i++)
65 region->b[i] += block->v[i];
66 region->c[i] += isqr(block->v[i]);
71 prequant_finish_region(struct image_sig_region *region)
73 if (region->count < 2)
74 memcpy(region->c, region->a, sizeof(region->c));
79 for (uns i = 0; i < IMAGE_VEC_F; i++)
81 region->e += region->c[i];
82 a += (u64)region->b[i] * region->b[i];
84 region->e -= a / 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 u32
96 #define ASORT_ELT(i) val[i]
97 #define ASORT_EXTRA_ARGS , u32 *val
98 #include "lib/arraysort.h"
101 prequant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions)
103 DBG("Starting pre-quantization");
105 uns regions_count, heap_count, axis, cov;
106 struct image_sig_block *blocks_end = blocks + blocks_count, *block, *block2;
107 struct image_sig_region *heap[IMAGE_REG_MAX + 1], *region, *region2;
109 /* Initialize single region with all blocks */
110 regions_count = heap_count = 1;
112 prequant_init_region(regions);
113 for (block = blocks; block != blocks_end; block++)
114 prequant_add_block(regions, block);
115 prequant_finish_region(regions);
118 while (regions_count < IMAGE_REG_MAX &&
119 regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count)
122 DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u",
123 regions_count, heap_count, region->count, (uns)region->e);
124 if (region->count < 2 ||
125 region->e < image_sig_prequant_thresholds[regions_count - 1] * region->count)
127 HEAP_DELMIN(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
131 /* Select axis to split - the one with maximum covariance */
133 cov = region->count * region->c[0] - region->b[0];
134 for (uns i = 1; i < 6; i++)
136 uns j = region->count * region->c[i] - region->b[i];
143 DBG("Splitting axis %u with covariance %u", axis, cov / region->count);
145 /* Sort values on the split axis */
146 u32 val[region->count];
147 block = region->blocks;
148 for (uns i = 0; i < region->count; i++, block = block->next)
149 val[i] = block->v[axis];
150 prequant_sort(region->count, val);
152 /* Select split value - to minimize error */
153 uns b1 = 0, c1 = 0, b2 = region->b[axis], c2 = region->c[axis];
154 uns i = 0, j = region->count, best_v = 0;
155 u64 best_err = 0xffffffffffffffff;
156 while (i < region->count)
158 u64 err = (u64)i * c1 - (u64)b1 * b1 + (u64)j * c2 - (u64)b2 * b2;
164 uns sqr = isqr(val[i]);
172 uns split_val = best_v;
173 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]);
176 block = region->blocks;
177 region2 = regions + regions_count++;
178 prequant_init_region(region);
179 prequant_init_region(region2);
182 block2 = block->next;
183 if (block->v[axis] < split_val)
184 prequant_add_block(region, block);
186 prequant_add_block(region2, block);
189 prequant_finish_region(region);
190 prequant_finish_region(region2);
191 HEAP_INCREASE(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1);
192 heap[++heap_count] = region2;
193 HEAP_INSERT(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
196 DBG("Pre-quantized to %u regions", regions_count);
198 return regions_count;
201 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
204 postquant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions, uns regions_count)
206 DBG("Starting post-quantization");
208 struct image_sig_block *blocks_end = blocks + blocks_count, *block;
209 struct image_sig_region *regions_end = regions + regions_count, *region;
210 uns error = 0, last_error;
212 /* Initialize regions and initial segmentation error */
213 for (region = regions; region != regions_end; )
215 uns inv = 0xffffffffU / region->count;
216 for (uns i = 0; i < IMAGE_VEC_F; i++)
218 region->a[i] = ((u64)region->b[i] * inv) >> 32;
219 error += region->c[i] - region->a[i] * region->b[i];
224 /* Convergation cycle */
225 for (uns step = 0; step < image_sig_postquant_max_steps; step++)
230 for (region = regions; region != regions_end; region++)
232 region->blocks = NULL;
234 bzero(region->b, sizeof(region->b));
235 bzero(region->c, sizeof(region->c));
238 /* Assign each block to its nearest pivot and accumulate region variables */
239 for (block = blocks; block != blocks_end; block++)
241 struct image_sig_region *best_region = NULL;
243 for (region = regions; region != regions_end; region++)
246 isqr(block->v[0] - region->a[0]) +
247 isqr(block->v[1] - region->a[1]) +
248 isqr(block->v[2] - region->a[2]) +
249 isqr(block->v[3] - region->a[3]) +
250 isqr(block->v[4] - region->a[4]) +
251 isqr(block->v[5] - region->a[5]);
252 if (dist <= best_dist)
255 best_region = region;
258 region = best_region;
260 block->next = region->blocks;
261 region->blocks = block;
262 for (uns i = 0; i < IMAGE_VEC_F; i++)
264 region->b[i] += block->v[i];
265 region->c[i] += isqr(block->v[i]);
269 /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
272 for (region = regions; region != regions_end; )
275 uns inv = 0xffffffffU / region->count;
276 for (uns i = 0; i < IMAGE_VEC_F; i++)
278 region->a[i] = ((u64)region->b[i] * inv) >> 32;
279 error += region->c[i] - region->a[i] * region->b[i];
286 *region = *regions_end;
289 DBG("last_error=%u error=%u", last_error, error);
291 /* Convergation criteria */
292 if (step >= image_sig_postquant_min_steps)
294 if (error > last_error)
296 u64 dif = last_error - error;
297 if (dif * image_sig_postquant_threshold < last_error * 100)
302 DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
304 return regions_end - regions;
308 image_sig_segmentation(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions)
311 regions_count = prequant(blocks, blocks_count, regions);
313 dump_segmentation(regions, regions_count);
315 regions_count = postquant(blocks, blocks_count, regions, regions_count);
317 dump_segmentation(regions, regions_count);
319 return regions_count;