From: Pavel Charvat Date: Fri, 25 Aug 2006 09:48:47 +0000 (+0200) Subject: image segmentation moved to a new file X-Git-Tag: holmes-import~591 X-Git-Url: http://mj.ucw.cz/gitweb/?a=commitdiff_plain;h=956f5794fc8dad3ed98095f578cddbfe342f8db2;p=libucw.git image segmentation moved to a new file --- diff --git a/images/Makefile b/images/Makefile index 47e267d1..3283e44b 100644 --- a/images/Makefile +++ b/images/Makefile @@ -6,7 +6,7 @@ PROGS+=$(addprefix $(o)/images/,image-tool image-dup-test image-sim-test) CONFIGS+=images -LIBIMAGES_MODS=math config image scale color alpha io-main dup-init dup-cmp sig-dump sig-init sig-cmp object +LIBIMAGES_MODS=math config image scale color alpha io-main dup-init dup-cmp sig-dump sig-init sig-seg sig-cmp object LIBIMAGES_LIBS=-lm diff --git a/images/sig-init.c b/images/sig-init.c index 7751a424..6b99e343 100644 --- a/images/sig-init.c +++ b/images/sig-init.c @@ -30,314 +30,6 @@ struct block { struct block *next; }; -struct region { - struct block *blocks; - u32 count; - u32 a[IMAGE_VEC_F]; - u32 b[IMAGE_VEC_F]; - u32 c[IMAGE_VEC_F]; - u64 e; - u64 w_sum; -}; - -static inline uns -dist(uns a, uns b) -{ - int d = a - b; - return d * d; -} - -#ifdef LOCAL_DEBUG -static void -dump_segmentation(struct region *regions, uns regions_count, uns cols, uns rows) -{ - uns size = (cols + 1) * rows; - byte buf[size]; - bzero(buf, size); - for (uns i = 0; i < regions_count; i++) - { - byte c = (i < 10) ? '0' + i : 'A' - 10 + i; - for (struct block *b = regions[i].blocks; b; b = b->next) - buf[b->x + b->y * (cols + 1)] = c; - } - for (uns i = 0; i < rows; i++) - log(L_DEBUG, "%s", &buf[i * (cols + 1)]); -} -#endif - -/* Pre-quantization - recursively split groups of blocks with large error */ - -static inline void -prequant_init_region(struct region *region) -{ - bzero(region, sizeof(*region)); -} - -static inline void -prequant_add_block(struct region *region, struct block *block) -{ - block->next = region->blocks; - region->blocks = block; - region->count++; - for (uns i = 0; i < IMAGE_VEC_F; i++) - { - region->b[i] += block->v[i]; - region->c[i] += isqr(block->v[i]); - } -} - -static void -prequant_finish_region(struct region *region) -{ - if (region->count < 2) - memcpy(region->c, region->a, sizeof(region->c)); - else - { - u64 a = 0; - region->e = 0; - for (uns i = 0; i < IMAGE_VEC_F; i++) - { - region->e += region->c[i]; - a += (u64)region->b[i] * region->b[i]; - } - region->e -= a / region->count; - } -} - -static inline uns -prequant_heap_cmp(struct region *a, struct region *b) -{ - return a->e > b->e; -} - -#define ASORT_PREFIX(x) prequant_##x -#define ASORT_KEY_TYPE u32 -#define ASORT_ELT(i) val[i] -#define ASORT_EXTRA_ARGS , u32 *val -#include "lib/arraysort.h" - -static uns -prequant(struct block *blocks, uns blocks_count, struct region *regions) -{ - DBG("Starting pre-quantization"); - - uns regions_count, heap_count, axis, cov; - struct block *blocks_end = blocks + blocks_count, *block, *block2; - struct region *heap[IMAGE_REG_MAX + 1], *region, *region2; - - /* Initialize single region with all blocks */ - regions_count = heap_count = 1; - heap[1] = regions; - prequant_init_region(regions); - for (block = blocks; block != blocks_end; block++) - prequant_add_block(regions, block); - prequant_finish_region(regions); - - /* Main cycle */ - while (regions_count < IMAGE_REG_MAX && - regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count) - { - region = heap[1]; - DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u", - regions_count, heap_count, region->count, (uns)region->e); - if (region->count < 2 || - region->e < image_sig_prequant_thresholds[regions_count - 1] * region->count) - { - HEAP_DELMIN(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP); - continue; - } - - /* Select axis to split - the one with maximum covariance */ - axis = 0; - cov = region->count * region->c[0] - region->b[0]; - for (uns i = 1; i < 6; i++) - { - uns j = region->count * region->c[i] - region->b[i]; - if (j > cov) - { - axis = i; - cov = j; - } - } - DBG("Splitting axis %u with covariance %u", axis, cov / region->count); - - /* Sort values on the split axis */ - u32 val[region->count]; - block = region->blocks; - for (uns i = 0; i < region->count; i++, block = block->next) - val[i] = block->v[axis]; - prequant_sort(region->count, val); - - /* Select split value - to minimize error */ - uns b1 = 0, c1 = 0, b2 = region->b[axis], c2 = region->c[axis]; - uns i = 0, j = region->count, best_v = 0; - u64 best_err = 0xffffffffffffffff; - while (i < region->count) - { - u64 err = (u64)i * c1 - (u64)b1 * b1 + (u64)j * c2 - (u64)b2 * b2; - if (err < best_err) - { - best_err = err; - best_v = val[i]; - } - uns sqr = isqr(val[i]); - b1 += val[i]; - b2 -= val[i]; - c1 += sqr; - c2 -= sqr; - i++; - j--; - } - uns split_val = best_v; - 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]); - - /* Split region */ - block = region->blocks; - region2 = regions + regions_count++; - prequant_init_region(region); - prequant_init_region(region2); - while (block) - { - block2 = block->next; - if (block->v[axis] < split_val) - prequant_add_block(region, block); - else - prequant_add_block(region2, block); - block = block2; - } - prequant_finish_region(region); - prequant_finish_region(region2); - HEAP_INCREASE(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1); - heap[++heap_count] = region2; - HEAP_INSERT(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP); - } - - DBG("Pre-quantized to %u regions", regions_count); - - return regions_count; -} - -/* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */ - -static uns -postquant(struct block *blocks, uns blocks_count, struct region *regions, uns regions_count) -{ - DBG("Starting post-quantization"); - - struct block *blocks_end = blocks + blocks_count, *block; - struct region *regions_end = regions + regions_count, *region; - uns error = 0, last_error; - - /* Initialize regions and initial segmentation error */ - for (region = regions; region != regions_end; ) - { - uns inv = 0xffffffffU / region->count; - for (uns i = 0; i < IMAGE_VEC_F; i++) - { - region->a[i] = ((u64)region->b[i] * inv) >> 32; - error += region->c[i] - region->a[i] * region->b[i]; - } - region++; - } - - /* Convergation cycle */ - for (uns step = 0; step < image_sig_postquant_max_steps; step++) - { - DBG("Step..."); - - /* Clear regions */ - for (region = regions; region != regions_end; region++) - { - region->blocks = NULL; - region->count = 0; - bzero(region->b, sizeof(region->b)); - bzero(region->c, sizeof(region->c)); - } - - /* Assign each block to its nearest pivot and accumulate region variables */ - for (block = blocks; block != blocks_end; block++) - { - struct region *best_region = NULL; - uns best_dist = ~0U; - for (region = regions; region != regions_end; region++) - { - uns dist = - isqr(block->v[0] - region->a[0]) + - isqr(block->v[1] - region->a[1]) + - isqr(block->v[2] - region->a[2]) + - isqr(block->v[3] - region->a[3]) + - isqr(block->v[4] - region->a[4]) + - isqr(block->v[5] - region->a[5]); - if (dist <= best_dist) - { - best_dist = dist; - best_region = region; - } - } - region = best_region; - region->count++; - block->next = region->blocks; - region->blocks = block; - for (uns i = 0; i < IMAGE_VEC_F; i++) - { - region->b[i] += block->v[i]; - region->c[i] += isqr(block->v[i]); - } - } - - /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */ - last_error = error; - error = 0; - for (region = regions; region != regions_end; ) - if (region->count) - { - uns inv = 0xffffffffU / region->count; - for (uns i = 0; i < IMAGE_VEC_F; i++) - { - region->a[i] = ((u64)region->b[i] * inv) >> 32; - error += region->c[i] - region->a[i] * region->b[i]; - } - region++; - } - else - { - regions_end--; - *region = *regions_end; - } - - DBG("last_error=%u error=%u", last_error, error); - - /* Convergation criteria */ - if (step >= image_sig_postquant_min_steps) - { - if (error > last_error) - break; - u64 dif = last_error - error; - if (dif * image_sig_postquant_threshold < last_error * 100) - break; - } - } - - DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count); - - return regions_end - regions; -} - -static inline uns -segmentation(struct block *blocks, uns blocks_count, struct region *regions, uns width UNUSED, uns height UNUSED) -{ - uns regions_count; - regions_count = prequant(blocks, blocks_count, regions); -#ifdef LOCAL_DEBUG - dump_segmentation(regions, regions_count, width, height); -#endif - regions_count = postquant(blocks, blocks_count, regions, regions_count); -#ifdef LOCAL_DEBUG - dump_segmentation(regions, regions_count, width, height); -#endif - return regions_count; -} - int compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image) { @@ -353,7 +45,7 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h); uns blocks_count = w * h; - struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks; + struct image_sig_block *blocks = xmalloc(blocks_count * sizeof(struct image_sig_block)), *block = blocks; /* Every block of 4x4 pixels */ byte *row_start = image->pixels; @@ -492,8 +184,8 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu } /* Quantize blocks to image regions */ - struct region regions[IMAGE_REG_MAX]; - sig->len = segmentation(blocks, blocks_count, regions, w, h); + struct image_sig_region regions[IMAGE_REG_MAX]; + sig->len = image_sig_segmentation(blocks, blocks_count, regions); /* For each region */ u64 w_total = 0; @@ -501,7 +193,7 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu uns w_mul = 127 * 256 / w_border; for (uns i = 0; i < sig->len; i++) { - struct region *r = regions + i; + struct image_sig_region *r = regions + i; DBG("Processing region %u: count=%u", i, r->count); ASSERT(r->count); @@ -515,7 +207,7 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu /* Compute coordinates centroid and region weight */ u64 x_avg = 0, y_avg = 0, w_sum = 0; - for (struct block *b = r->blocks; b; b = b->next) + for (struct image_sig_block *b = r->blocks; b; b = b->next) { x_avg += b->x; y_avg += b->y; @@ -536,9 +228,9 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu /* Compute normalized inertia */ u64 sum1 = 0, sum2 = 0, sum3 = 0; - for (struct block *b = r->blocks; b; b = b->next) + for (struct image_sig_block *b = r->blocks; b; b = b->next) { - uns inc2 = dist(x_avg, b->x) + dist(y_avg, b->y); + uns inc2 = isqr(x_avg - b->x) + isqr(y_avg - b->y); uns inc1 = sqrt(inc2); sum1 += inc1; sum2 += inc2; @@ -566,11 +258,11 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu { uns d = 0; for (uns k = 0; k < IMAGE_REG_F; k++) - d += dist(sig->reg[i].f[k], sig->reg[j].f[k]); + d += isqr(sig->reg[i].f[k] - sig->reg[j].f[k]); df += sqrt(d); d = 0; for (uns k = 0; k < IMAGE_REG_H; k++) - d += dist(sig->reg[i].h[k], sig->reg[j].h[k]); + d += isqr(sig->reg[i].h[k] - sig->reg[j].h[k]); dh += sqrt(d); cnt++; } @@ -583,7 +275,7 @@ compute_image_signature(struct image_thread *thread UNUSED, struct image_signatu uns wa = 128, wb = 128; for (uns i = sig->len; --i > 0; ) { - struct region *r = regions + i; + struct image_sig_region *r = regions + i; wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i)); wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i)); } diff --git a/images/sig-seg.c b/images/sig-seg.c new file mode 100644 index 00000000..9b5bb268 --- /dev/null +++ b/images/sig-seg.c @@ -0,0 +1,321 @@ +/* + * Image Library -- Image segmentation + * + * (c) 2006 Pavel Charvat + * + * This software may be freely distributed and used according to the terms + * of the GNU Lesser General Public License. + */ + +#undef LOCAL_DEBUG + +#include "sherlock/sherlock.h" +#include "lib/math.h" +#include "lib/fastbuf.h" +#include "lib/conf.h" +#include "lib/heap.h" +#include "images/math.h" +#include "images/images.h" +#include "images/color.h" +#include "images/signature.h" + +#include + +#ifdef LOCAL_DEBUG +static void +dump_segmentation(struct image_sig_region *regions, uns regions_count) +{ + uns cols = 0, rows = 0; + for (uns i = 0; i < regions_count; i++) + for (struct image_sig_block *b = regions[i].blocks; b; b = b->next) + { + cols = MAX(cols, b->x); + rows = MAX(rows, b->y); + } + uns size = (cols + 1) * rows; + byte buf[size]; + bzero(buf, size); + for (uns i = 0; i < regions_count; i++) + { + byte c = (i < 10) ? '0' + i : 'A' - 10 + i; + for (struct image_sig_block *b = regions[i].blocks; b; b = b->next) + buf[b->x + b->y * (cols + 1)] = c; + } + for (uns i = 0; i < rows; i++) + log(L_DEBUG, "%s", &buf[i * (cols + 1)]); +} +#endif + +/* Pre-quantization - recursively split groups of blocks with large error */ + +static inline void +prequant_init_region(struct image_sig_region *region) +{ + bzero(region, sizeof(*region)); +} + +static inline void +prequant_add_block(struct image_sig_region *region, struct image_sig_block *block) +{ + block->next = region->blocks; + region->blocks = block; + region->count++; + for (uns i = 0; i < IMAGE_VEC_F; i++) + { + region->b[i] += block->v[i]; + region->c[i] += isqr(block->v[i]); + } +} + +static void +prequant_finish_region(struct image_sig_region *region) +{ + if (region->count < 2) + memcpy(region->c, region->a, sizeof(region->c)); + else + { + u64 a = 0; + region->e = 0; + for (uns i = 0; i < IMAGE_VEC_F; i++) + { + region->e += region->c[i]; + a += (u64)region->b[i] * region->b[i]; + } + region->e -= a / region->count; + } +} + +static inline uns +prequant_heap_cmp(struct image_sig_region *a, struct image_sig_region *b) +{ + return a->e > b->e; +} + +#define ASORT_PREFIX(x) prequant_##x +#define ASORT_KEY_TYPE u32 +#define ASORT_ELT(i) val[i] +#define ASORT_EXTRA_ARGS , u32 *val +#include "lib/arraysort.h" + +static uns +prequant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions) +{ + DBG("Starting pre-quantization"); + + uns regions_count, heap_count, axis, cov; + struct image_sig_block *blocks_end = blocks + blocks_count, *block, *block2; + struct image_sig_region *heap[IMAGE_REG_MAX + 1], *region, *region2; + + /* Initialize single region with all blocks */ + regions_count = heap_count = 1; + heap[1] = regions; + prequant_init_region(regions); + for (block = blocks; block != blocks_end; block++) + prequant_add_block(regions, block); + prequant_finish_region(regions); + + /* Main cycle */ + while (regions_count < IMAGE_REG_MAX && + regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count) + { + region = heap[1]; + DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u", + regions_count, heap_count, region->count, (uns)region->e); + if (region->count < 2 || + region->e < image_sig_prequant_thresholds[regions_count - 1] * region->count) + { + HEAP_DELMIN(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP); + continue; + } + + /* Select axis to split - the one with maximum covariance */ + axis = 0; + cov = region->count * region->c[0] - region->b[0]; + for (uns i = 1; i < 6; i++) + { + uns j = region->count * region->c[i] - region->b[i]; + if (j > cov) + { + axis = i; + cov = j; + } + } + DBG("Splitting axis %u with covariance %u", axis, cov / region->count); + + /* Sort values on the split axis */ + u32 val[region->count]; + block = region->blocks; + for (uns i = 0; i < region->count; i++, block = block->next) + val[i] = block->v[axis]; + prequant_sort(region->count, val); + + /* Select split value - to minimize error */ + uns b1 = 0, c1 = 0, b2 = region->b[axis], c2 = region->c[axis]; + uns i = 0, j = region->count, best_v = 0; + u64 best_err = 0xffffffffffffffff; + while (i < region->count) + { + u64 err = (u64)i * c1 - (u64)b1 * b1 + (u64)j * c2 - (u64)b2 * b2; + if (err < best_err) + { + best_err = err; + best_v = val[i]; + } + uns sqr = isqr(val[i]); + b1 += val[i]; + b2 -= val[i]; + c1 += sqr; + c2 -= sqr; + i++; + j--; + } + uns split_val = best_v; + 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]); + + /* Split region */ + block = region->blocks; + region2 = regions + regions_count++; + prequant_init_region(region); + prequant_init_region(region2); + while (block) + { + block2 = block->next; + if (block->v[axis] < split_val) + prequant_add_block(region, block); + else + prequant_add_block(region2, block); + block = block2; + } + prequant_finish_region(region); + prequant_finish_region(region2); + HEAP_INCREASE(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1); + heap[++heap_count] = region2; + HEAP_INSERT(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP); + } + + DBG("Pre-quantized to %u regions", regions_count); + + return regions_count; +} + +/* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */ + +static uns +postquant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions, uns regions_count) +{ + DBG("Starting post-quantization"); + + struct image_sig_block *blocks_end = blocks + blocks_count, *block; + struct image_sig_region *regions_end = regions + regions_count, *region; + uns error = 0, last_error; + + /* Initialize regions and initial segmentation error */ + for (region = regions; region != regions_end; ) + { + uns inv = 0xffffffffU / region->count; + for (uns i = 0; i < IMAGE_VEC_F; i++) + { + region->a[i] = ((u64)region->b[i] * inv) >> 32; + error += region->c[i] - region->a[i] * region->b[i]; + } + region++; + } + + /* Convergation cycle */ + for (uns step = 0; step < image_sig_postquant_max_steps; step++) + { + DBG("Step..."); + + /* Clear regions */ + for (region = regions; region != regions_end; region++) + { + region->blocks = NULL; + region->count = 0; + bzero(region->b, sizeof(region->b)); + bzero(region->c, sizeof(region->c)); + } + + /* Assign each block to its nearest pivot and accumulate region variables */ + for (block = blocks; block != blocks_end; block++) + { + struct image_sig_region *best_region = NULL; + uns best_dist = ~0U; + for (region = regions; region != regions_end; region++) + { + uns dist = + isqr(block->v[0] - region->a[0]) + + isqr(block->v[1] - region->a[1]) + + isqr(block->v[2] - region->a[2]) + + isqr(block->v[3] - region->a[3]) + + isqr(block->v[4] - region->a[4]) + + isqr(block->v[5] - region->a[5]); + if (dist <= best_dist) + { + best_dist = dist; + best_region = region; + } + } + region = best_region; + region->count++; + block->next = region->blocks; + region->blocks = block; + for (uns i = 0; i < IMAGE_VEC_F; i++) + { + region->b[i] += block->v[i]; + region->c[i] += isqr(block->v[i]); + } + } + + /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */ + last_error = error; + error = 0; + for (region = regions; region != regions_end; ) + if (region->count) + { + uns inv = 0xffffffffU / region->count; + for (uns i = 0; i < IMAGE_VEC_F; i++) + { + region->a[i] = ((u64)region->b[i] * inv) >> 32; + error += region->c[i] - region->a[i] * region->b[i]; + } + region++; + } + else + { + regions_end--; + *region = *regions_end; + } + + DBG("last_error=%u error=%u", last_error, error); + + /* Convergation criteria */ + if (step >= image_sig_postquant_min_steps) + { + if (error > last_error) + break; + u64 dif = last_error - error; + if (dif * image_sig_postquant_threshold < last_error * 100) + break; + } + } + + DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count); + + return regions_end - regions; +} + +uns +image_sig_segmentation(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions) +{ + uns regions_count; + regions_count = prequant(blocks, blocks_count, regions); +#ifdef LOCAL_DEBUG + dump_segmentation(regions, regions_count); +#endif + regions_count = postquant(blocks, blocks_count, regions, regions_count); +#ifdef LOCAL_DEBUG + dump_segmentation(regions, regions_count); +#endif + return regions_count; +} + diff --git a/images/signature.h b/images/signature.h index fd736f74..54452638 100644 --- a/images/signature.h +++ b/images/signature.h @@ -47,6 +47,27 @@ image_signature_size(uns len) byte *image_vector_dump(byte *buf, struct image_vector *vec); byte *image_region_dump(byte *buf, struct image_region *reg); +struct image_sig_block { + struct image_sig_block *next; + u32 area; /* block area in pixels (usually 16) */ + u32 v[IMAGE_VEC_F]; + u32 x, y; /* block position */ +}; + +struct image_sig_region { + struct image_sig_block *blocks; + u32 count; + u32 a[IMAGE_VEC_F]; + u32 b[IMAGE_VEC_F]; + u32 c[IMAGE_VEC_F]; + u64 e; + u64 w_sum; +}; + +/* sig-seg.c */ + +uns image_sig_segmentation(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions); + /* sig-init.c */ int compute_image_signature(struct image_thread *thread, struct image_signature *sig, struct image *image);