]> mj.ucw.cz Git - libucw.git/blob - images/sig-seg.c
Merge with git+ssh://git.ucw.cz/projects/sherlock/GIT/sherlock.git
[libucw.git] / images / sig-seg.c
1 /*
2  *      Image Library -- Image segmentation
3  *
4  *      (c) 2006 Pavel Charvat <pchar@ucw.cz>
5  *
6  *      This software may be freely distributed and used according to the terms
7  *      of the GNU Lesser General Public License.
8  */
9
10 #undef LOCAL_DEBUG
11
12 #include "lib/lib.h"
13 #include "lib/conf.h"
14 #include "lib/heap.h"
15 #include "images/images.h"
16 #include "images/signature.h"
17 #include "images/math.h"
18
19 #include <string.h>
20
21 #ifdef LOCAL_DEBUG
22 static void
23 dump_segmentation(struct image_sig_region *regions, uns regions_count)
24 {
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)
28       {
29         cols = MAX(cols, b->x + 1);
30         rows = MAX(rows, b->y + 1);
31       }
32   uns size = (cols + 1) * rows;
33   byte buf[size];
34   bzero(buf, size);
35   for (uns i = 0; i < regions_count; i++)
36     {
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;
40     }
41   for (uns i = 0; i < rows; i++)
42     log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
43 }
44 #endif
45
46 /* Pre-quantization - recursively split groups of blocks with large error */
47
48 static inline void
49 prequant_init_region(struct image_sig_region *region)
50 {
51   bzero(region, sizeof(*region));
52 }
53
54 static inline void
55 prequant_add_block(struct image_sig_region *region, struct image_sig_block *block)
56 {
57   block->next = region->blocks;
58   region->blocks = block;
59   region->count++;
60   for (uns i = 0; i < IMAGE_VEC_F; i++)
61     {
62       region->b[i] += block->v[i];
63       region->c[i] += isqr(block->v[i]);
64     }
65 }
66
67 static void
68 prequant_finish_region(struct image_sig_region *region)
69 {
70   if (region->count < 2)
71     {
72       region->e = 0;
73     }
74   else
75     {
76       u64 a = 0;
77       region->e = 0;
78       for (uns i = 0; i < IMAGE_VEC_F; i++)
79         {
80           region->e += region->c[i];
81           a += (u64)region->b[i] * region->b[i];
82         }
83       region->e -= a / region->count;
84       DBG("Finished region %u", (uns)region->e / region->count);
85     }
86 }
87
88 static inline uns
89 prequant_heap_cmp(struct image_sig_region *a, struct image_sig_region *b)
90 {
91   return a->e > b->e;
92 }
93
94 #define ASORT_PREFIX(x) prequant_##x
95 #define ASORT_KEY_TYPE uns
96 #define ASORT_ELT(i) val[i]
97 #define ASORT_EXTRA_ARGS , uns *val
98 #include "lib/arraysort.h"
99
100 static uns
101 prequant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions)
102 {
103   DBG("Starting pre-quantization");
104
105   uns regions_count, heap_count, axis;
106   struct image_sig_block *blocks_end = blocks + blocks_count, *block, *block2;
107   struct image_sig_region *heap[IMAGE_REG_MAX + 1], *region, *region2;
108
109   /* Initialize single region with all blocks */
110   regions_count = heap_count = 1;
111   heap[1] = regions;
112   prequant_init_region(regions);
113   for (block = blocks; block != blocks_end; block++)
114     prequant_add_block(regions, block);
115   prequant_finish_region(regions);
116
117   /* Main cycle */
118   while (regions_count < IMAGE_REG_MAX &&
119       regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count)
120     {
121       region = heap[1];
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] * blocks_count)
126         {
127           HEAP_DELMIN(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
128           continue;
129         }
130
131       /* Select axis to split - the one with maximum average quadratic error */
132       axis = 0;
133       u64 cov = (u64)region->count * region->c[0] - (u64)region->b[0] * region->b[0];
134       for (uns i = 1; i < 6; i++)
135         {
136           uns j = (u64)region->count * region->c[i] - (u64)region->b[i] * region->b[i];
137           if (j > cov)
138             {
139               axis = i;
140               cov = j;
141             }
142         }
143       DBG("Splitting axis %u with average quadratic error %u", axis, (uns)(cov / (region->count * region->count)));
144
145       /* Sort values on the split axis */
146       uns val[256], cnt[256], cval;
147       if (region->count > 64)
148         {
149           bzero(cnt, sizeof(cnt));
150           for (block = region->blocks; block; block = block->next)
151             cnt[block->v[axis]]++;
152           cval = 0;
153           for (uns i = 0; i < 256; i++)
154             if (cnt[i])
155               {
156                 val[cval] = i;
157                 cnt[cval] = cnt[i];
158                 cval++;
159               }
160         }
161       else
162         {
163           block = region->blocks;
164           for (uns i = 0; i < region->count; i++, block = block->next)
165             val[i] = block->v[axis];
166           prequant_sort(region->count, val);
167           cval = 1;
168           cnt[0] = 1;
169           for (uns i = 1; i < region->count; i++)
170             if (val[i] == val[cval - 1])
171               cnt[cval - 1]++;
172             else
173               {
174                 val[cval] = val[i];
175                 cnt[cval] = 1;
176                 cval++;
177               }
178         }
179
180       /* Select split value - to minimize error */
181       uns b1 = val[0] * cnt[0];
182       uns c1 = isqr(val[0]) * cnt[0];
183       uns b2 = region->b[axis] - b1;
184       uns c2 = region->c[axis] - c1;
185       uns i = cnt[0], j = region->count - cnt[0];
186       u64 best_err = c1 - (u64)b1 * b1 / i + c2 - (u64)b2 * b2 / j;
187       uns split_val = val[0];
188       for (uns k = 1; k < cval - 1; k++)
189         {
190           uns b0 = val[k] * cnt[k];
191           uns c0 = isqr(val[k]) * cnt[k];
192           b1 += b0;
193           b2 -= b0;
194           c1 += c0;
195           c2 -= c0;
196           i += cnt[k];
197           j -= cnt[k];
198           u64 err = (u64)c1 - (u64)b1 * b1 / i + (u64)c2 - (u64)b2 * b2 / j;
199           if (err < best_err)
200             {
201               best_err = err;
202               split_val = val[k];
203             }
204         }
205       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
207       /* Split region */
208       block = region->blocks;
209       region2 = regions + regions_count++;
210       prequant_init_region(region);
211       prequant_init_region(region2);
212       while (block)
213         {
214           block2 = block->next;
215           if (block->v[axis] <= split_val)
216             prequant_add_block(region, block);
217           else
218             prequant_add_block(region2, block);
219           block = block2;
220         }
221       prequant_finish_region(region);
222       prequant_finish_region(region2);
223       HEAP_INCREASE(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1);
224       heap[++heap_count] = region2;
225       HEAP_INSERT(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
226     }
227
228   DBG("Pre-quantized to %u regions", regions_count);
229
230   return regions_count;
231 }
232
233 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
234
235 static uns
236 postquant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions, uns regions_count)
237 {
238   DBG("Starting post-quantization");
239
240   struct image_sig_block *blocks_end = blocks + blocks_count, *block;
241   struct image_sig_region *regions_end = regions + regions_count, *region;
242   uns error = 0, last_error;
243
244   /* Initialize regions and initial segmentation error */
245   for (region = regions; region != regions_end; )
246     {
247       uns inv = 0xffffffffU / region->count;
248       for (uns i = 0; i < IMAGE_VEC_F; i++)
249         {
250           region->a[i] = ((u64)region->b[i] * inv) >> 32;
251           error += region->c[i] - region->a[i] * region->b[i];
252         }
253       region++;
254     }
255
256   /* Convergation cycle */
257   for (uns step = 0; step < image_sig_postquant_max_steps; step++)
258     {
259       DBG("Step...");
260
261       /* Clear regions */
262       for (region = regions; region != regions_end; region++)
263         {
264           region->blocks = NULL;
265           region->count = 0;
266           bzero(region->b, sizeof(region->b));
267           bzero(region->c, sizeof(region->c));
268         }
269
270       /* Assign each block to its nearest pivot and accumulate region variables */
271       for (block = blocks; block != blocks_end; block++)
272         {
273           struct image_sig_region *best_region = NULL;
274           uns best_dist = ~0U;
275           for (region = regions; region != regions_end; region++)
276             {
277               uns dist =
278                 isqr(block->v[0] - region->a[0]) +
279                 isqr(block->v[1] - region->a[1]) +
280                 isqr(block->v[2] - region->a[2]) +
281                 isqr(block->v[3] - region->a[3]) +
282                 isqr(block->v[4] - region->a[4]) +
283                 isqr(block->v[5] - region->a[5]);
284               if (dist <= best_dist)
285                 {
286                   best_dist = dist;
287                   best_region = region;
288                 }
289             }
290           region = best_region;
291           region->count++;
292           block->next = region->blocks;
293           region->blocks = block;
294           for (uns i = 0; i < IMAGE_VEC_F; i++)
295             {
296               region->b[i] += block->v[i];
297               region->c[i] += isqr(block->v[i]);
298             }
299         }
300
301       /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
302       last_error = error;
303       error = 0;
304       for (region = regions; region != regions_end; )
305         if (region->count)
306           {
307             uns inv = 0xffffffffU / region->count;
308             for (uns i = 0; i < IMAGE_VEC_F; i++)
309               {
310                 region->a[i] = ((u64)region->b[i] * inv) >> 32;
311                 error += region->c[i] - region->a[i] * region->b[i];
312               }
313             region++;
314           }
315         else
316           {
317             regions_end--;
318             *region = *regions_end;
319           }
320
321       DBG("last_error=%u error=%u", last_error, error);
322
323       /* Convergation criteria */
324       if (step >= image_sig_postquant_min_steps)
325         {
326           if (error > last_error)
327             break;
328           u64 dif = last_error - error;
329           if (dif * image_sig_postquant_threshold < (u64)last_error * 100)
330             break;
331         }
332     }
333
334   DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
335
336   return regions_end - regions;
337 }
338
339 void
340 image_sig_segmentation(struct image_sig_data *data)
341 {
342   data->regions_count = prequant(data->blocks, data->blocks_count, data->regions);
343 #ifdef LOCAL_DEBUG
344   dump_segmentation(data->regions, data->regions_count);
345 #endif
346   data->regions_count = postquant(data->blocks, data->blocks_count, data->regions, data->regions_count);
347 #ifdef LOCAL_DEBUG
348   dump_segmentation(data->regions, data->regions_count);
349 #endif
350 }
351