matching.c

   1 /*########################################################################
   2 *
   3 *  Copyright(C) 2002-2007. All Rights Reserved.
   4 *
   5 *  Authors: Shivang Patel
   6 *           Jaap de Haan(jdh)
   7 *
   8 *  Changes: jdh -> Added support for ImageMagick that enables
   9 *                  to export files to more than 40 formats.
  10 *   This is free software; you can redistribute it and/or modify it under
  11 *   the terms of the GNU General Public License as published by the Free
  12 *   Software Foundation; either version 2, or (at your option) any later
  13 *   version.
  14 *
  15 *   This is distributed in the hope that it will be useful, but WITHOUT
  16 *   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  17 *   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  18 *   for more details.
  19 *
  20 *   You should have received a copy of the GNU General Public License with
  21 *   the fvs source package as the
  22 *   file COPYING. If not, write to the Free Software Foundation, Inc.,
  23 *   59 Temple Place - Suite 330, Boston, MA
  24 *   02111-1307, USA.
  25 ########################################################################*/
  26 // this code by pallav gupta (pallav@ieee.org)
  27
  28 // this algorithm is described in two papers.
  29
  30 // xudong jiang and wei-yun yau, "fingerprint minutiae matching based on the
  31 //local and global structure
  32
  33 // shenglin yang and ingrid m. verbauwhede, "a secure fingerprint matching
  34 // technique"
  35
  36 #include <math.h>
  37 #include <string.h>
  38 #include <stdlib.h>
  39 #include <stdio.h>
  40 #include <time.h>
  41
  42 #include "matching.h"
  43
  44 static FvsInt_t CompareLocalStructures(Fvs_LocalStructure_t *ils,
  45                                        Fvs_LocalStructure_t *tls,
  46                                        FvsInt_t ilen, FvsInt_t tlen);
  47
  48 static void FindMaxDistance(FvsMinutia_t *src, FvsInt_t *neigh, FvsInt_t ref,
  49                             FvsFloat_t *d, FvsInt_t *idx);
  50 static void FindNeighbors(FvsMinutia_t *src, FvsInt_t length, FvsInt_t ineigh[][NEIGHBORS]);
  51 static void ComputeLocalStructure(FvsMinutia_t *src, FvsInt_t length,
  52                                   FvsInt_t ineigh[][NEIGHBORS],
  53                                   Fvs_LocalStructure_t *ls);
  54 static inline FvsFloat_t anglediff(FvsFloat_t a, FvsFloat_t b);
  55
  56
  57 FvsError_t MatchingCompareMinutiaSets
  58 (
  59     const FvsMinutiaSet_t set1,
  60     const FvsMinutiaSet_t set2,
  61     FvsFloat_t* pgoodness
  62 )
  63 {
  64   FvsMinutia_t *iminu = MinutiaSetGetBuffer(set1);
  65   FvsInt_t nb_iminu = MinutiaSetGetCount(set1);
  66
  67   FvsMinutia_t *tminu = MinutiaSetGetBuffer(set2);
  68   FvsInt_t nb_tminu = MinutiaSetGetCount(set2);
  69
  70   FvsInt_t ineigh[MAXMINUTIA][NEIGHBORS] = {};
  71   FvsInt_t tneigh[MAXMINUTIA][NEIGHBORS] = {};
  72
  73   Fvs_LocalStructure_t ls_iminu[MAXMINUTIA];
  74   Fvs_LocalStructure_t ls_tminu[MAXMINUTIA];
  75
  76   FvsInt_t matched, max;
  77
  78   // find the neigherest neighbors
  79   FindNeighbors(iminu, nb_iminu, ineigh);
  80   FindNeighbors(tminu, nb_tminu, tneigh);
  81   // compute the local structures of the minutia
  82   ComputeLocalStructure(iminu, nb_iminu, ineigh, ls_iminu);
  83   ComputeLocalStructure(tminu, nb_tminu, tneigh, ls_tminu);
  84   matched = CompareLocalStructures(ls_iminu, ls_tminu, nb_iminu, nb_tminu);
  85   max = (nb_iminu > nb_tminu) ? nb_iminu : nb_tminu;
  86   (*pgoodness) = (FvsFloat_t)matched/max;
  87   return FvsOK;
  88 }
  89
  90
  91 static FvsInt_t CompareLocalStructures(Fvs_LocalStructure_t *ils,
  92                                        Fvs_LocalStructure_t *tls,
  93                                        FvsInt_t ilen, FvsInt_t tlen)
  94 {
  95   Fvs_LocalStructure_t s1, s2;
  96   FvsInt_t marked, matched;
  97   FvsInt_t i, j, k, l;
  98
  99   matched = 0;
 100   marked = 0;
 101
 102   // one minutia from template
 103   for (i = 0; i < ilen; i++) {
 104     s1 = ils[i];
 105
 106     // one minutia from input
 107     for (j = 0; j < tlen; j++) {
 108
 109       s2 = tls[j];
 110
 111       // need to go further only if the local ridge direction are similiar
 112       if (fabs(s1.angle - s2.angle) > DELTA_ANGLE)
 113                                 continue;
 114
 115       marked = 0;
 116
 117       // each neighbor of input minutia
 118       for (k = 0; k < NEIGHBORS; k++) {
 119                                 // each neighbor of template minutia
 120                                 for (l = 0; l < NEIGHBORS; l++) {
 121
 122                                   if (fabs(s1.distance[k] - s2.distance[l]) > DELTA_DISTANCE)
 123                                     continue;
 124                                   if (fabs(s1.varphi[k] - s2.varphi[l]) > DELTA_VARPHI)
 125                                     continue;
 126                                   if (fabs(s1.vartheta[k] - s2.vartheta[l]) > DELTA_VARTHETA)
 127                                     continue;
 128                                   // all conditions matched, consider k'th neighboer and l'th neighbor marked
 129                                   marked++;
 130                                 }
 131       }
 132
 133       if (marked > THRESHOLD_MARKED)
 134                                 matched++;
 135     }
 136   }
 137
 138   return matched;
 139 }
 140
 141
 142 static void ComputeLocalStructure(FvsMinutia_t *src, FvsInt_t length,
 143                                   FvsInt_t ineigh[][NEIGHBORS], Fvs_LocalStructure_t *ls)
 144 {
 145   FvsFloat_t deltax, deltay, distance;
 146   FvsInt_t i, j, p;
 147
 148   for (i = 0; i < length; i++) {
 149     for (j = 0; j < NEIGHBORS; j++) {
 150       // get index of jth neighbor
 151       p = ineigh[i][j];
 152
 153       // compute distance
 154       deltax = src[p].x - src[i].x;
 155       deltay = src[p].y - src[i].y;
 156       distance = sqrt(deltax*deltax + deltay*deltay);
 157       ls[i].distance[j] = distance;
 158
 159       // compute phi
 160       ls[i].varphi[j] = anglediff(src[p].angle, src[i].angle);
 161
 162       // compute vartheta
 163       ls[i].vartheta[j] = anglediff(atan2(deltay, deltax), src[i].angle);
 164
 165       // assign angle
 166       ls[i].angle = src[i].angle;
 167     }
 168   }
 169 }
 170
 171 static inline FvsFloat_t anglediff(FvsFloat_t a, FvsFloat_t b)
 172 {
 173   FvsFloat_t dif;
 174
 175   dif = a-b;
 176
 177   if (dif > -M_PI && dif <= M_PI)
 178     return dif;
 179   else if (dif <= -M_PI)
 180     return (2*M_PI + dif);
 181   else
 182     return (2*M_PI - dif);
 183 }
 184
 185
 186
 187 static void FindNeighbors(FvsMinutia_t *src, FvsInt_t length, FvsInt_t ineigh[][NEIGHBORS])
 188 {
 189
 190   FvsInt_t i, j, k, idx, end, wrap = 0;
 191   FvsFloat_t deltax, deltay, distance, maxdis;
 192
 193   // find the neigherest neighbors for input minutia
 194   for (i = 0; i < length; i++) {
 195     // pick the first 6 points, aribtrarily
 196     j = i+1;
 197     for (k = 0; k < NEIGHBORS; k++) {
 198       if (j >= length) {
 199                                 j = 0;
 200                                 wrap = 1;
 201       }
 202       ineigh[i][k] = j++;
 203     }
 204
 205     // find the max distance
 206     FindMaxDistance(src, ineigh[i], i, &maxdis, &idx);
 207
 208     end = wrap ? i : length;
 209     j = wrap ? i-1 : j;
 210     // search upwards
 211     while (((src[j].x - src[i].x) < maxdis) && (j < end) && (j >= 0)) {
 212       deltax = src[j].x - src[i].x;
 213       deltay = src[j].y - src[i].y;
 214       distance = sqrt(deltax*deltax + deltay*deltay);
 215
 216       if (distance < maxdis) {
 217                                 ineigh[i][idx] = j;
 218                                 FindMaxDistance(src, ineigh[i], i, &maxdis, &idx);
 219       }
 220
 221       if (wrap)
 222                                 j--;
 223       else
 224                                 j++;
 225     }
 226
 227     // only need to search backwards if when picking we didn't wrap around
 228     // if we wrapped around searching upwards will have covered all points
 229     if (!wrap) {
 230       j = i-1;
 231       while (((src[i].x - src[j].x) < maxdis) && (j >= 0)) {
 232                                 deltax = src[i].x - src[j].x;
 233                                 deltay = src[i].y - src[j].y;
 234                                 distance = sqrt(deltax*deltax + deltay*deltay);
 235
 236                                 if (distance < maxdis) {
 237                                   ineigh[i][idx] = j;
 238                                   FindMaxDistance(src, ineigh[i], i, &maxdis, &idx);
 239                                 }
 240                                 j--;
 241       }
 242     }
 243   }
 244 }
 245
 246 static void FindMaxDistance(FvsMinutia_t *src, FvsInt_t *neigh, FvsInt_t ref, FvsFloat_t *d,
 247                             FvsInt_t *idx)
 248 {
 249   FvsInt_t i, maxi = 0, p = neigh[0];
 250   FvsFloat_t deltax = src[p].x - src[ref].x;
 251   FvsFloat_t deltay = src[p].y - src[ref].y;
 252   FvsFloat_t maxd = sqrt(deltax*deltax + deltay*deltay), tempd;
 253
 254   //  printf("origin: %f %f\n", src[ref].x, src[ref].y);
 255   //printf("0: %f %f %f\n", src[p].x, src[p].y, tempd);
 256
 257   for (i = 1; i < NEIGHBORS; i++) {
 258     p = neigh[i];
 259     deltax = src[p].x - src[ref].x;
 260     deltay = src[p].y - src[ref].y;
 261     tempd = sqrt(deltax*deltax + deltay*deltay);
 262     //printf("%d: %f %f %f\n", i, src[p].x, src[p].y, tempd);
 263
 264     if (tempd > maxd) {
 265       maxd = tempd;
 266       maxi = i;
 267     }
 268   }
 269
 270   //printf("max = %f, idx = %d\n\n\n", maxd, maxi);
 271   (*d) = maxd;
 272   (*idx) = maxi;
 273 }
 274
 275 FvsError_t MinutiaSetSort(FvsMinutia_t *src, FvsInt_t length)
 276 {
 277   FvsFloat_t x, y, angle;
 278   FvsMinutiaType_t type;
 279
 280   FvsInt_t i, j, looking;
 281
 282   if (src == NULL)
 283     return FvsMemory;
 284
 285   for (i = 1; i < length; i++) {
 286
 287     x = src[i].x;
 288     y = src[i].y;
 289     angle = src[i].angle;
 290     type = src[i].type;
 291
 292     j = i-1;
 293     looking = 1;
 294
 295     while ((j >= 0) && (looking == 1)) {
 296       // sort with respec to x coordinate
 297       if (x < src[j].x) {
 298                                 src[j + 1].angle = src[j].angle;
 299                                 src[j + 1].x = src[j].x;
 300                                 src[j + 1].y = src[j].y;
 301                                 src[j + 1].type = src[j].type;
 302
 303                                 src[j].angle = angle;
 304                                 src[j].x = x;
 305                                 src[j].y = y;
 306                                 src[j].type = type;
 307
 308                                 --j;
 309                          } else  {
 310                                 looking = 0;
 311                                 src[j + 1].angle = angle;
 312                                 src[j + 1].x = x;
 313                                 src[j + 1].y = y;
 314                                 src[j + 1].type = type;
 315       }
 316     }
 317   }
 318
 319   return FvsOK;
 320 }