images/sig-cmp.c

   1 /*
   2  *      Image Library -- Comparitions of image signatures
   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 #define LOCAL_DEBUG
  11
  12 #include "lib/lib.h"
  13 #include "lib/math.h"
  14 #include "images/math.h"
  15 #include "images/images.h"
  16 #include "images/signature.h"
  17
  18 #include <stdio.h>
  19
  20 static uns
  21 image_signatures_dist_1(struct image_signature *sig1, struct image_signature *sig2)
  22 {
  23   DBG("image_signatures_dist_1()");
  24
  25   uns cnt1 = sig1->len;
  26   uns cnt2 = sig2->len;
  27   struct image_region *reg1 = sig1->reg;
  28   struct image_region *reg2 = sig2->reg;
  29   uns mf[IMAGE_REG_MAX][IMAGE_REG_MAX], mh[IMAGE_REG_MAX][IMAGE_REG_MAX];
  30   uns lf[IMAGE_REG_MAX * 2], lh[IMAGE_REG_MAX * 2];
  31   uns df = sig1->df + sig2->df, dh = sig1->dh + sig2->dh;
  32
  33   /* Compute distance matrix */
  34   for (uns i = 0; i < cnt1; i++)
  35     for (uns j = 0; j < cnt2; j++)
  36       {
  37         uns d = 0;
  38         for (uns k = 0; k < IMAGE_REG_F; k++)
  39           {
  40             int dif = reg1[i].f[k] - reg2[j].f[k];
  41             d += dif * dif;
  42           }
  43         mf[i][j] = d;
  44         d = 0;
  45         for (uns k = 0; k < IMAGE_REG_H; k++)
  46           {
  47             int dif = reg1[i].h[k] - reg2[j].h[k];
  48             d += dif * dif;
  49           }
  50         mh[i][j] = d;
  51       }
  52
  53   uns lfs = 0, lhs = 0;
  54   for (uns i = 0; i < cnt1; i++)
  55     {
  56       uns f = mf[i][0], h = mh[i][0];
  57       for (uns j = 1; j < cnt2; j++)
  58         {
  59           f = MIN(f, mf[i][j]);
  60           h = MIN(h, mh[i][j]);
  61         }
  62       lf[i] = (df * 0x10000) / (df + fast_sqrt_u32(f));
  63       lh[i] = (dh * 0x10000) / (dh + fast_sqrt_u32(h));
  64       lfs += lf[i] * (6 * reg1[i].wa + 2 * reg1[i].wb);
  65       lhs += lh[i] * reg1[i].wa;
  66     }
  67   for (uns i = 0; i < cnt2; i++)
  68     {
  69       uns f = mf[0][i], h = mh[0][i];
  70       for (uns j = 1; j < cnt1; j++)
  71         {
  72           f = MIN(f, mf[j][i]);
  73           h = MIN(h, mh[j][i]);
  74         }
  75       lf[i + cnt1] = (df * 0x10000) / (df + fast_sqrt_u32(f));
  76       lh[i + cnt1] = (dh * 0x10000) / (dh + fast_sqrt_u32(h));
  77       lfs += lf[i] * (6 * reg2[i].wa + 2 * reg2[i].wb);
  78       lhs += lh[i] * reg2[i].wa;
  79     }
  80
  81   uns measure = lfs * 6 + lhs * 2 * 8;
  82
  83 #ifdef LOCAL_DEBUG
  84   /* Display similarity vectors */
  85   byte buf[2 * IMAGE_REG_MAX * 16 + 3], *b = buf;
  86   for (uns i = 0; i < cnt1 + cnt2; i++)
  87     {
  88       if (i)
  89         *b++ = ' ';
  90       if (i == cnt1)
  91         *b++ = '~', *b++ = ' ';
  92       b += sprintf(b, "%.4f", (double)lf[i] / 0x10000);
  93     }
  94   *b = 0;
  95   DBG("Lf=(%s)", buf);
  96   b = buf;
  97   for (uns i = 0; i < cnt1 + cnt2; i++)
  98     {
  99       if (i)
 100         *b++ = ' ';
 101       if (i == cnt1)
 102         *b++ = '~', *b++ = ' ';
 103       b += sprintf(b, "%.4f", (double)lh[i] / 0x10000);
 104     }
 105   *b = 0;
 106   DBG("Lh=(%s)", buf);
 107   DBG("Lfm=%.4f", lfs / (double)(1 << (3 + 8 + 16)));
 108   DBG("Lhm=%.4f", lhs / (double)(1 << (8 + 16)));
 109   DBG("measure=%.4f", measure / (double)(1 << (3 + 3 + 8 + 16)));
 110 #endif
 111
 112   return (1 << (3 + 3 + 8 + 16)) - measure;
 113 }
 114
 115 #define ASORT_PREFIX(x) image_signatures_dist_2_##x
 116 #define ASORT_KEY_TYPE uns
 117 #define ASORT_ELT(i) items[i]
 118 #define ASORT_EXTRA_ARGS , uns *items
 119 #include "lib/arraysort.h"
 120
 121 static uns
 122 image_signatures_dist_2(struct image_signature *sig1, struct image_signature *sig2)
 123 {
 124   DBG("image_signatures_dist_2()");
 125
 126   uns dist[IMAGE_REG_MAX * IMAGE_REG_MAX], p[IMAGE_REG_MAX], q[IMAGE_REG_MAX];
 127   uns n, i, j, k, l, s, d;
 128   struct image_region *reg1, *reg2;
 129
 130   /* Compute distance matrix */
 131   n = 0;
 132   /* ... for non-textured images */
 133   if (!((sig1->flags | sig2->flags) & IMAGE_SIG_TEXTURED))
 134     for (j = 0, reg2 = sig2->reg; j < sig2->len; j++, reg2++)
 135       for (i = 0, reg1 = sig1->reg; i < sig1->len; i++, reg1++)
 136         {
 137           uns ds =
 138             isqr((int)reg1->h[0] - (int)reg2->h[0]) +
 139             isqr((int)reg1->h[1] - (int)reg2->h[1]) +
 140             isqr((int)reg1->h[2] - (int)reg2->h[2]);
 141           uns dt =
 142             isqr((int)reg1->f[0] - (int)reg2->f[0]) +
 143             isqr((int)reg1->f[1] - (int)reg2->f[1]) +
 144             isqr((int)reg1->f[2] - (int)reg2->f[2]) +
 145             isqr((int)reg1->f[3] - (int)reg2->f[3]) +
 146             isqr((int)reg1->f[4] - (int)reg2->f[4]) +
 147             isqr((int)reg1->f[5] - (int)reg2->f[5]);
 148           if (ds < 1000)
 149             dt *= 8;
 150           else if (ds < 10000)
 151             dt *= 12;
 152           else
 153             dt *= 16;
 154           DBG("[%u][%u] ... ds=%u dt=%u", i, j, ds, dt);
 155           dist[n++] = (dt << 8) + i + (j << 4) ;
 156         }
 157   /* ... for textured images (ignore shape properties) */
 158   else
 159     for (j = 0, reg2 = sig2->reg; j < sig2->len; j++, reg2++)
 160       for (i = 0, reg1 = sig1->reg; i < sig1->len; i++, reg1++)
 161         {
 162           uns dt =
 163             isqr((int)reg1->f[0] - (int)reg2->f[0]) +
 164             isqr((int)reg1->f[1] - (int)reg2->f[1]) +
 165             isqr((int)reg1->f[2] - (int)reg2->f[2]) +
 166             isqr((int)reg1->f[3] - (int)reg2->f[3]) +
 167             isqr((int)reg1->f[4] - (int)reg2->f[4]) +
 168             isqr((int)reg1->f[5] - (int)reg2->f[5]);
 169           dist[n++] = (dt << 12) + i + (j << 4) ;
 170         }
 171
 172   /* One or both signatures have no regions */
 173   if (!n)
 174     return 1 << IMAGE_SIG_DIST_SCALE;
 175
 176   /* Get percentages */
 177   for (i = 0, reg1 = sig1->reg; i < sig1->len; i++, reg1++)
 178     p[i] = reg1->wa;
 179   for (i = 0, reg2 = sig2->reg; i < sig2->len; i++, reg2++)
 180     q[i] = reg2->wa;
 181
 182   /* Sort entries in distance matrix */
 183   image_signatures_dist_2_sort(n, dist);
 184
 185   /* Compute significance matrix and resulting distance */
 186   uns sum = 0;
 187   for (k = 0, l = 128; l; k++)
 188     {
 189       i = dist[k] & 15;
 190       j = (dist[k] >> 4) & 15;
 191       d = dist[k] >> 8;
 192       if (p[i] <= q[j])
 193         {
 194           s = p[i];
 195           q[j] -= p[i];
 196           p[i] = 0;
 197         }
 198       else
 199         {
 200           s = q[j];
 201           p[i] -= q[j];
 202           q[j] = 0;
 203         }
 204       l -= s;
 205       sum += s * d;
 206       DBG("s[%u][%u]=%u d=%u", i, j, s, d);
 207     }
 208
 209   return sum;
 210 }
 211
 212 uns
 213 image_signatures_dist(struct image_signature *sig1, struct image_signature *sig2)
 214 {
 215   switch (image_sig_compare_method)
 216     {
 217       case 1:
 218         return image_signatures_dist_1(sig1, sig2);
 219       case 2:
 220         return image_signatures_dist_2(sig1, sig2);
 221       default:
 222         die("Invalid image signatures compare method.");
 223     }
 224 }
 225