minutia.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
  27 /* For doing time calculations */
  28 #include <time.h>
  29 #include <sys/time.h>
  30 #include <sys/resource.h>
  31
  32 #include <string.h>
  33 #include <stdlib.h>
  34 #include <stdio.h>
  35 #include <math.h>
  36
  37 #include "minutia.h"
  38 #include "matching.h"
  39 #include "fvs.h"
  40
  41
  42 typedef struct iFvsMinutiaSet_t
  43 {
  44     FvsInt_t     nbminutia; /* number of minutia set in the table */
  45     FvsInt_t     tablesize; /* max number of minutia that can be stored in the table */
  46     FvsMinutia_t ptable[1]; /* the minutia table starts here */
  47 } iFvsMinutiaSet_t;
  48
  49
  50
  51
  52 FvsMinutiaSet_t MinutiaSetCreate(const FvsInt_t size)
  53 {
  54     iFvsMinutiaSet_t* p = NULL;
  55     p = (iFvsMinutiaSet_t*)malloc(sizeof(iFvsMinutiaSet_t)+size*sizeof(FvsMinutia_t));
  56     if (p!=NULL)
  57     {   /* no minutia in table yet */
  58         p->nbminutia = 0;
  59         p->tablesize = size;
  60     }
  61     return (FvsMinutiaSet_t)p;
  62 }
  63
  64
  65 void MinutiaSetDestroy(FvsMinutiaSet_t minutia)
  66 {
  67     iFvsMinutiaSet_t* p = NULL;
  68     if (minutia==NULL)
  69         return;
  70     p = minutia;
  71     free(p);
  72 }
  73
  74
  75 /* returns the maximum number of minutia the table can contain */
  76 FvsInt_t MinutiaSetGetSize(const FvsMinutiaSet_t min)
  77 {
  78     const iFvsMinutiaSet_t* minutia = (iFvsMinutiaSet_t*)min;
  79     FvsInt_t nret = 0;
  80
  81     if (minutia!=NULL)
  82         nret = minutia->tablesize;
  83
  84     return nret;
  85 }
  86
  87
  88 /* returns the number of minutia in the set */
  89 FvsInt_t MinutiaSetGetCount(const FvsMinutiaSet_t min)
  90 {
  91     const iFvsMinutiaSet_t* minutia = (iFvsMinutiaSet_t*)min;
  92     FvsInt_t nret = 0;
  93
  94     if (minutia!=NULL)
  95         nret = minutia->nbminutia;
  96
  97     return nret;
  98 }
  99
 100
 101 /* returns a pointer to the table of minutia */
 102 FvsMinutia_t* MinutiaSetGetBuffer(FvsMinutiaSet_t min)
 103 {
 104     iFvsMinutiaSet_t* minutia = (iFvsMinutiaSet_t*)min;
 105     FvsMinutia_t* pret = NULL;
 106
 107     if (minutia!=NULL)
 108         pret = minutia->ptable;
 109
 110     return pret;
 111 }
 112
 113
 114 /* empty the minutia set */
 115 FvsError_t MinutiaSetEmpty(FvsMinutiaSet_t min)
 116 {
 117     iFvsMinutiaSet_t* minutia = (iFvsMinutiaSet_t*)min;
 118     FvsError_t nRet = FvsOK;
 119
 120     if (minutia!=NULL)
 121         minutia->nbminutia = 0;
 122     else
 123         nRet = FvsMemory;
 124
 125     return nRet;
 126 }
 127
 128
 129 FvsError_t MinutiaSetAdd(FvsMinutiaSet_t min,
 130        const FvsFloat_t x, const FvsFloat_t y,
 131        const FvsMinutiaType_t type, const FvsFloat_t angle)
 132 {
 133     iFvsMinutiaSet_t* minutia = (iFvsMinutiaSet_t*)min;
 134     FvsError_t nRet = FvsOK;
 135
 136     if (minutia->nbminutia < minutia->tablesize)
 137     {
 138         minutia->ptable[minutia->nbminutia].x       = x;
 139         minutia->ptable[minutia->nbminutia].y       = y;
 140         minutia->ptable[minutia->nbminutia].type    = type;
 141         minutia->ptable[minutia->nbminutia].angle   = angle;
 142         minutia->nbminutia++;
 143     }
 144     else
 145         /* no place left in table */
 146         nRet = FvsMemory;
 147
 148     return nRet;
 149 }
 150
 151
 152 //taken from matching.c by petertu and adapted for angles from -M_PI/2 to M_PI/2
 153 static inline FvsFloat_t anglediff(FvsFloat_t a, FvsFloat_t b)
 154 {
 155   FvsFloat_t dif;
 156
 157   dif = a-b;
 158
 159   if (dif > -M_PI/2 && dif <= M_PI/2)
 160     return dif;
 161   else if (dif <= -M_PI/2)
 162     return (M_PI + dif);
 163   else
 164     return (M_PI - dif);
 165 }
 166
 167 static FvsError_t MinutiaSetCheckClean(FvsMinutiaSet_t min)
 168 {
 169     iFvsMinutiaSet_t* minutia = (iFvsMinutiaSet_t*)min;
 170     FvsError_t    nRet = FvsOK;
 171     FvsFloat_t    tx, ty; /* tolerances */
 172     FvsInt_t      i, j;
 173     FvsMinutia_t* mi, *mj;
 174
 175     tx = 8.0;
 176     ty = 8.0;
 177
 178     if (minutia!=NULL && minutia->nbminutia > 1)
 179     {
 180                 // test if a minutia with these coordinates, type and
 181         //       angle is already in the table
 182         for (j = 0;   j < minutia->nbminutia; j++)
 183             for (i = j+1; i < minutia->nbminutia; i++)
 184             {
 185                 mi = minutia->ptable + i;
 186                 mj = minutia->ptable + j;
 187
 188                 // Minutia that are near each other and opposite directions -> Remove both
 189
 190                 if ( (fabs(mi->x     - mj->x    ) < tx     ) &&
 191                      (fabs(mi->y     - mj->y    ) < ty     ) &&
 192                      mi->type == FvsMinutiaTypeEnding &&
 193                      mj->type == FvsMinutiaTypeEnding &&
 194                      anglediff(mi->angle,mj->angle) < M_PI/8
 195                     )
 196                 {
 197                     // found two line endings that are close to each other -> remove one at a time
 198
 199                     if(i < minutia->nbminutia-1)
 200                     {
 201                         // swap mi with last non-dupilcate minutia
 202                         mi->x = minutia->ptable[minutia->nbminutia-1].x;
 203                         mi->y = minutia->ptable[minutia->nbminutia-1].y;
 204                         mi->angle = minutia->ptable[minutia->nbminutia-1].angle;
 205                         mi->type = minutia->ptable[minutia->nbminutia-1].type;
 206                     }
 207                     minutia->nbminutia--;
 208
 209                     // do same for nj
 210                     if(j < minutia->nbminutia-1)
 211                     {
 212                         mj->x = minutia->ptable[minutia->nbminutia-1].x;
 213                         mj->y = minutia->ptable[minutia->nbminutia-1].y;
 214                         mj->angle = minutia->ptable[minutia->nbminutia-1].angle;
 215                         mj->type = minutia->ptable[minutia->nbminutia-1].type;
 216                     }
 217                     minutia->nbminutia--;
 218                 }
 219             }
 220     }
 221     else
 222         /* no place left in table */
 223         nRet = FvsMemory;
 224
 225     return nRet;
 226 }
 227
 228
 229 #define P(x,y)      p[(x)+(y)*pitch]
 230
 231
 232 /* draw the mintuia set into the image */
 233 FvsError_t MinutiaSetDraw(const FvsMinutiaSet_t min, FvsImage_t image)
 234 {
 235     FvsInt_t w       = ImageGetWidth(image);
 236     FvsInt_t h       = ImageGetHeight(image);
 237     FvsInt_t pitch   = ImageGetPitch(image);
 238     FvsByte_t* p     = ImageGetBuffer(image);
 239
 240     FvsInt_t n, x, y;
 241     FvsFloat_t fx, fy;
 242     FvsMinutia_t* minutia = MinutiaSetGetBuffer(min);
 243     FvsInt_t nbminutia    = MinutiaSetGetCount(min);
 244
 245     if (minutia==NULL || p==NULL)
 246         return FvsMemory;
 247
 248     /* draw each minutia */
 249     for (n = 0; n < nbminutia; n++)
 250     {
 251         x = (FvsInt_t)minutia[n].x;
 252         y = (FvsInt_t)minutia[n].y;
 253         if (x<w-5 && x>4)
 254         {
 255             if (y<h-5 && y>4)
 256             {
 257                 switch (minutia[n].type)
 258                 {
 259                 case FvsMinutiaTypeEnding:
 260                     P(x,y)    = 0xFF;
 261                     P(x-1, y) = 0xA0;
 262                     P(x+1, y) = 0xA0;
 263                     P(x, y-1) = 0xA0;
 264                     P(x, y+1) = 0xA0;
 265                     break;
 266                 case FvsMinutiaTypeBranching:
 267                     P(x,y)    = 0xFF;
 268                     P(x-1, y-1) = 0xA0;
 269                     P(x+1, y-1) = 0xA0;
 270                     P(x-1, y+1) = 0xA0;
 271                     P(x+1, y+1) = 0xA0;
 272                     break;
 273                 default:
 274                     continue;
 275                 }
 276                 fx = sin(minutia[n].angle);
 277                 fy = -cos(minutia[n].angle);
 278                 P(x+(int32_t)(fx)    ,y+(int32_t)(fy)    ) = 0xFF;
 279                 P(x+(int32_t)(fx*2.0),y+(int32_t)(fy*2.0)) = 0xFF;
 280                 P(x+(int32_t)(fx*3.0),y+(int32_t)(fy*3.0)) = 0xFF;
 281                 P(x+(int32_t)(fx*4.0),y+(int32_t)(fy*4.0)) = 0xFF;
 282                 P(x+(int32_t)(fx*5.0),y+(int32_t)(fy*5.0)) = 0xFF;
 283             }
 284         }
 285     }
 286     return FvsOK;
 287 }
 288
 289
 290 /* defines to facilitate reading */
 291 // note by petertu: caution: these defines differ from the ones in imagemanip.c
 292 #define P1  P(x  ,y-1)
 293 #define P2  P(x+1,y-1)
 294 #define P3  P(x+1,y  )
 295 #define P4  P(x+1,y+1)
 296 #define P5  P(x  ,y+1)
 297 #define P6  P(x-1,y+1)
 298 #define P7  P(x-1,y  )
 299 #define P8  P(x-1,y-1)
 300
 301
 302 /*
 303 ** This function locates the minutia and compute their properties
 304 */
 305 FvsError_t MinutiaSetExtract
 306     (
 307     FvsMinutiaSet_t       minutia,
 308     const FvsImage_t      image,
 309     const FvsFloatField_t direction,
 310     const FvsImage_t      mask
 311     )
 312 {
 313     FvsInt_t w      = ImageGetWidth(image);
 314     FvsInt_t h      = ImageGetHeight(image);
 315     FvsInt_t pitch  = ImageGetPitch(image);
 316     FvsInt_t pitchm = ImageGetPitch(mask);
 317     FvsByte_t* p    = ImageGetBuffer(image);
 318     FvsByte_t* m    = ImageGetBuffer(mask);
 319     FvsInt_t   x, y;
 320     FvsFloat_t angle = 0.0;
 321     FvsInt_t   whitecount;
 322
 323     if (m==NULL || p==NULL)
 324         return FvsMemory;
 325
 326     (void)MinutiaSetEmpty(minutia);
 327
 328     /* loop through the image and find the minutas */
 329     for (y = 1; y < h-1; y++)
 330     for (x = 1; x < w-1; x++)
 331     {
 332         if (m[x+y*pitchm]==0)
 333              continue;
 334         if (P(x,y)==0xFF)
 335         {
 336             whitecount = 0;
 337             if (P1!=0) whitecount++;
 338             if (P2!=0) whitecount++;
 339             if (P3!=0) whitecount++;
 340             if (P4!=0) whitecount++;
 341             if (P5!=0) whitecount++;
 342             if (P6!=0) whitecount++;
 343             if (P7!=0) whitecount++;
 344             if (P8!=0) whitecount++;
 345
 346             switch(whitecount)
 347             {
 348             case 0:
 349                 /* isolated point, ignore it */
 350                 break;
 351             case 1:
 352                 /* detect angle */
 353                     angle = FloatFieldGetValue(direction, x, y);
 354                 (void)MinutiaSetAdd(minutia, (FvsFloat_t)x, (FvsFloat_t)y,
 355                               FvsMinutiaTypeEnding, (FvsFloat_t)angle);
 356                 break;
 357             case 2:
 358                 break;
 359             default:
 360                 /* normal line or branching */
 361                 /* it may be a branching but we must make sure */
 362 /*
 363                 if ((P(x-1, y+1) && P(x-1, y-1)) ||
 364                     (P(x-1, y-1) && P(x+1, y-1)) ||
 365                     (P(x+1, y-1) && P(x+1, y+1)) ||
 366                     (P(x+1, y+1) && P(x-1, y+1)))
 367 */
 368                                                         if ((P1 && P2) || (P2 && P3) || (P3 && P4) || (P4 && P5) || (P5 && P6) || (P6 && P7) || (P7 && P8) || (P8 && P1))
 369                                                                 {
 370                                                                         //do nothing
 371                                                                 }
 372                                                                 else
 373                 {
 374                     angle = FloatFieldGetValue(direction, x, y);
 375                     (void)MinutiaSetAdd(minutia, (FvsFloat_t)x, (FvsFloat_t)y,
 376                                   FvsMinutiaTypeBranching, (FvsFloat_t)angle);
 377                 }
 378                 /* is branch double check */
 379                 break;
 380             }
 381         }
 382     }
 383     (void)MinutiaSetCheckClean(minutia);
 384
 385     return FvsOK;
 386 }
 387
 388 static float subtractTimeSpec(struct timespec value1, struct timespec value2)
 389 {
 390     long time1 = value1.tv_sec * 1000000000 + value1.tv_nsec;
 391     long time2 = value2.tv_sec * 1000000000 + value2.tv_nsec;
 392
 393     return (float)(time2 - time1) / 1000000000;
 394 }
 395
 396 static float sumFloats(float *array, unsigned int size)
 397 {
 398     unsigned int i;
 399     float sum = 0;
 400     for (i = 0; i < size; i++)
 401     {
 402         sum += array[i];
 403     }
 404     return sum;
 405 }
 406
 407 FvsMinutiaSet_t* loadMinutiaSetsFromFiles(int numberOfFiles, const char **listOfFiles)
 408 {
 409     FvsMinutiaSet_t *minutiaArray = (FvsMinutiaSet_t*)malloc(numberOfFiles * sizeof(FvsMinutiaSet_t));
 410     struct timespec time1, time2;
 411     float elaspedTime, averageThreadTime, averageWallTime, *perThreadTime;
 412     perThreadTime = (float*)malloc(numberOfFiles * sizeof(float));
 413     clock_gettime(CLOCK_MONOTONIC_RAW, &time1);
 414
 415 #pragma omp parallel for
 416     for (int i = 0; i < numberOfFiles; ++i)
 417     {
 418         minutiaArray[i] = calculateMinutiaSet(listOfFiles[i], &perThreadTime[i]);
 419     }
 420
 421     /* sort the sets, this code was removed from MatchingCompareMinutiaSets */
 422 #pragma omp parallel for
 423     for (int i = 0; i < numberOfFiles; ++i)
 424         MinutiaSetSort(MinutiaSetGetBuffer(minutiaArray[i]), MinutiaSetGetCount(minutiaArray[i]));
 425
 426     clock_gettime(CLOCK_MONOTONIC_RAW, &time2);
 427     elaspedTime = subtractTimeSpec(time1, time2);
 428     averageWallTime = elaspedTime / numberOfFiles;
 429     averageThreadTime = sumFloats(perThreadTime, numberOfFiles) / numberOfFiles;
 430     printf("Load %d files: %02.3fsec = %02.3fsec/file, cpu time: %02.3fsec, speedup: %02.2fx\n",
 431            numberOfFiles, elaspedTime, averageWallTime, averageThreadTime, averageThreadTime / averageWallTime);
 432     free(perThreadTime);
 433     return minutiaArray;
 434 }
 435
 436 FvsFloat_t* compareMinutiaSets(int numberOfSets, const FvsMinutiaSet_t *minutiaArray)
 437 {
 438     FvsFloat_t *goodness = (FvsFloat_t*)malloc(numberOfSets * sizeof(FvsFloat_t));
 439     goodness[0] = 1;            /* base set is always equal to itself */
 440     /* start loop at 1! */
 441 #pragma omp parallel for
 442     for (int i = 1; i < numberOfSets; ++i)
 443     {
 444         MatchingCompareMinutiaSets(minutiaArray[0], minutiaArray[i], goodness + i);
 445     }
 446     return goodness;
 447 }
 448
 449 void cleanupMinutiaSetArray(int numberOfSets, FvsMinutiaSet_t *minutiaArray)
 450 {
 451     for (int i = 0; i < numberOfSets; ++i)
 452     {
 453         MinutiaSetDestroy(minutiaArray[i]);
 454     }
 455     free(minutiaArray);
 456 }
 457
 458 static float subtractTimeval(struct timeval value2, struct timeval value1)
 459 {
 460     long time1 = value1.tv_sec * 1000000 + value1.tv_usec;
 461     long time2 = value2.tv_sec * 1000000 + value2.tv_usec;
 462
 463     return (float)(time2 - time1) / 1000000;
 464 }
 465
 466 FvsMinutiaSet_t calculateMinutiaSet(const char *file, float *thread_time)
 467 {
 468     FvsImage_t image;
 469     FvsImage_t mask;
 470     FvsFloatField_t direction;
 471     FvsMinutiaSet_t minutia;
 472     FvsFloatField_t frequency;
 473     struct rusage usage1, usage2;
 474     getrusage(RUSAGE_THREAD, &usage1);
 475
 476     mask      = ImageCreate();
 477     image     = ImageCreate();
 478     direction = FloatFieldCreate();
 479     frequency = FloatFieldCreate();
 480     minutia   = MinutiaSetCreate(100);
 481
 482     if (mask!=NULL && image!=NULL && direction!=NULL && minutia!=NULL && frequency!=NULL)
 483     {
 484         FvsImageImport(image, file);
 485         ImageSoftenMean(image, 3);
 486         ImagePixelNormalize(image, 9);
 487         // creates a fixed-size mask, which is needed to ignore the mostly
 488         // unusable border areas in an image
 489         FingerprintGetFixedMask(image, mask, 6);
 490         FingerprintGetDirection(image, direction, 5, 4);
 491         // apply mask to image
 492         ImageLogical(image,mask, FvsLogicalAnd);
 493         FingerprintGetFrequency(image, direction, frequency);
 494         ImageEnhanceGabor(image, direction, frequency, mask, 4.0);
 495         // Binarize the image
 496         ImageBinarize(image, (FvsByte_t)0x80);
 497         // Thinning: Fist the standard fvs algorithm (ImageThinHitMiss) is
 498         // used, then the ImageEnhanceThinned algorithm written by petertu is
 499         // applied
 500         ImageThinHitMiss(image);
 501         ImageEnhanceThinned(image, direction);
 502         // detecting the minutia in the fingerprint and writing them to minutia
 503         // list
 504         MinutiaSetExtract(minutia, image, direction, mask);
 505         // this function draws the minutias as small vectors. the contents of
 506         // "image" are erased first
 507         ImageClear(image);
 508         MinutiaSetDraw(minutia, image);
 509     }
 510     ImageDestroy(image);
 511     ImageDestroy(mask);
 512     FloatFieldDestroy(direction);
 513
 514     getrusage(RUSAGE_THREAD, &usage2);
 515     *thread_time = subtractTimeval(usage2.ru_utime, usage1.ru_utime);
 516     return minutia;
 517 }
 518
 519 FvsFloat_t* calculateGoodnesses(int numberOfFiles, const char **listOfFiles)
 520 {
 521     FvsMinutiaSet_t *minutiaArray = loadMinutiaSetsFromFiles(numberOfFiles, listOfFiles);
 522     FvsFloat_t *goodness = compareMinutiaSets(numberOfFiles, minutiaArray);
 523     cleanupMinutiaSetArray(numberOfFiles, minutiaArray);
 524     return goodness;
 525 }