dft.c 16.2 KB
Newer Older
1 2
/*******************************************************************************

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
License:
This software and/or related materials was developed at the National Institute
of Standards and Technology (NIST) by employees of the Federal Government
in the course of their official duties. Pursuant to title 17 Section 105
of the United States Code, this software is not subject to copyright
protection and is in the public domain.

This software and/or related materials have been determined to be not subject
to the EAR (see Part 734.3 of the EAR for exact details) because it is
a publicly available technology and software, and is freely distributed
to any interested party with no licensing requirements.  Therefore, it is
permissible to distribute this software as a free download from the internet.

Disclaimer:
This software and/or related materials was developed to promote biometric
standards and biometric technology testing for the Federal Government
in accordance with the USA PATRIOT Act and the Enhanced Border Security
and Visa Entry Reform Act. Specific hardware and software products identified
in this software were used in order to perform the software development.
In no case does such identification imply recommendation or endorsement
by the National Institute of Standards and Technology, nor does it imply that
the products and equipment identified are necessarily the best available
for the purpose.

This software and/or related materials are provided "AS-IS" without warranty
of any kind including NO WARRANTY OF PERFORMANCE, MERCHANTABILITY,
NO WARRANTY OF NON-INFRINGEMENT OF ANY 3RD PARTY INTELLECTUAL PROPERTY
or FITNESS FOR A PARTICULAR PURPOSE or for any purpose whatsoever, for the
licensed product, however used. In no event shall NIST be liable for any
damages and/or costs, including but not limited to incidental or consequential
damages of any kind, including economic damage or injury to property and lost
profits, regardless of whether NIST shall be advised, have reason to know,
or in fact shall know of the possibility.

By using this software, you agree to bear all risk relating to quality,
use and performance of the software and/or related materials.  You agree
to hold the Government harmless from any claim arising from your use
of the software.
41 42 43

*******************************************************************************/

44

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
/***********************************************************************
      LIBRARY: LFS - NIST Latent Fingerprint System

      FILE:    DFT.C
      AUTHOR:  Michael D. Garris
      DATE:    03/16/1999
      UPDATED: 03/16/2005 by MDG

      Contains routines responsible for conducting Discrete Fourier
      Transforms (DFT) analysis on a block of image data as part of
      the NIST Latent Fingerprint System (LFS).

***********************************************************************
               ROUTINES:
                        dft_dir_powers()
                        sum_rot_block_rows()
                        dft_power()
                        dft_power_stats()
                        get_max_norm()
                        sort_dft_waves()
***********************************************************************/

#include <stdio.h>
#include <lfs.h>

70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
/*************************************************************************
**************************************************************************
#cat: dft_dir_powers - Conducts the DFT analysis on a block of image data.
#cat:         The image block is sampled across a range of orientations
#cat:         (directions) and multiple wave forms of varying frequency are
#cat:         applied at each orientation.  At each orentation, pixels are
#cat:         accumulated along each rotated pixel row, creating a vector
#cat:         of pixel row sums.  Each DFT wave form is then applied
#cat:         individually to this vector of pixel row sums.  A DFT power
#cat:         value is computed for each wave form (frequency0 at each
#cat:         orientaion within the image block.  Therefore, the resulting DFT
#cat:         power vectors are of dimension (N Waves X M Directions).
#cat:         The power signatures derived form this process are used to
#cat:         determine dominant direction flow within the image block.

   Input:
      pdata     - the padded input image.  It is important that the image
                  be properly padded, or else the sampling at various block
                  orientations may result in accessing unkown memory.
      blkoffset - the pixel offset form the origin of the padded image to
                  the origin of the current block in the image
      pw        - the width (in pixels) of the padded input image
      ph        - the height (in pixels) of the padded input image
      dftwaves  - structure containing the DFT wave forms
      dftgrids  - structure containing the rotated pixel grid offsets
   Output:
      powers    - DFT power computed from each wave form frequencies at each
                  orientation (direction) in the current image block
   Return Code:
      Zero     - successful completion
      Negative - system error
**************************************************************************/
int dft_dir_powers(double **powers, unsigned char *pdata,
               const int blkoffset, const int pw, const int ph,
               const DFTWAVES *dftwaves, const ROTGRIDS *dftgrids)
{
   int w, dir;
   int *rowsums;
   unsigned char *blkptr;

   /* Allocate line sum vector, and initialize to zeros */
   /* This routine requires square block (grid), so ERROR otherwise. */
   if(dftgrids->grid_w != dftgrids->grid_h){
      fprintf(stderr, "ERROR : dft_dir_powers : DFT grids must be square\n");
      return(-90);
   }
   rowsums = (int *)malloc(dftgrids->grid_w * sizeof(int));
   if(rowsums == (int *)NULL){
      fprintf(stderr, "ERROR : dft_dir_powers : malloc : rowsums\n");
      return(-91);
   }
121
   memset(rowsums, 0, dftgrids->grid_w * sizeof(int));
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

   /* Foreach direction ... */
   for(dir = 0; dir < dftgrids->ngrids; dir++){
      /* Compute vector of line sums from rotated grid */
      blkptr = pdata + blkoffset;
      sum_rot_block_rows(rowsums, blkptr,
                         dftgrids->grids[dir], dftgrids->grid_w);

      /* Foreach DFT wave ... */
      for(w = 0; w < dftwaves->nwaves; w++){
         dft_power(&(powers[w][dir]), rowsums,
                   dftwaves->waves[w], dftwaves->wavelen);
      }
   }

   /* Deallocate working memory. */
   free(rowsums);

   return(0);
}

143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
/*************************************************************************
**************************************************************************
#cat: sum_rot_block_rows - Computes a vector or pixel row sums by sampling
#cat:               the current image block at a given orientation.  The
#cat:               sampling is conducted using a precomputed set of rotated
#cat:               pixel offsets (called a grid) relative to the orgin of
#cat:               the image block.

   Input:
      blkptr    - the pixel address of the origin of the current image block
      grid_offsets - the rotated pixel offsets for a block-sized grid
                  rotated according to a specific orientation
      blocksize - the width and height of the image block and thus the size
                  of the rotated grid
   Output:
      rowsums   - the resulting vector of pixel row sums
**************************************************************************/
160
void sum_rot_block_rows(int *rowsums, const unsigned char *blkptr,
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
                        const int *grid_offsets, const int blocksize)
{
   int ix, iy, gi;

   /* Initialize rotation offset index. */
   gi = 0;

   /* For each row in block ... */
   for(iy = 0; iy < blocksize; iy++){
      /* The sums are accumlated along the rotated rows of the grid, */
      /* so initialize row sum to 0.                                 */
      rowsums[iy] = 0;
      /* Foreach column in block ... */
      for(ix = 0; ix < blocksize; ix++){
         /* Accumulate pixel value at rotated grid position in image */
         rowsums[iy] += *(blkptr + grid_offsets[gi]);
         gi++;
      }
   }
}

/*************************************************************************
**************************************************************************
#cat: dft_power - Computes the DFT power by applying a specific wave form
#cat:             frequency to a vector of pixel row sums computed from a
#cat:             specific orientation of the block image

   Input:
      rowsums - accumulated rows of pixels from within a rotated grid
                overlaying an input image block
      wave    - the wave form (cosine and sine components) at a specific
                frequency
      wavelen - the length of the wave form (must match the height of the
                image block which is the length of the rowsum vector)
   Output:
      power   - the computed DFT power for the given wave form at the
                given orientation within the image block
**************************************************************************/
199
void dft_power(double *power, const int *rowsums,
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
               const DFTWAVE *wave, const int wavelen)
{
   int i;
   double cospart, sinpart;

   /* Initialize accumulators */
   cospart = 0.0;
   sinpart = 0.0;

   /* Accumulate cos and sin components of DFT. */
   for(i = 0; i < wavelen; i++){
      /* Multiply each rotated row sum by its        */
      /* corresponding cos or sin point in DFT wave. */
      cospart += (rowsums[i] * wave->cos[i]);
      sinpart += (rowsums[i] * wave->sin[i]);
   }

   /* Power is the sum of the squared cos and sin components */
   *power = (cospart * cospart) + (sinpart * sinpart);
}

/*************************************************************************
**************************************************************************
223 224 225 226 227 228 229 230 231
#cat: dft_power_stats - Derives statistics from a set of DFT power vectors.
#cat:           Statistics are computed for all but the lowest frequency
#cat:           wave form, including the Maximum power for each wave form,
#cat:           the direction at which the maximum power occured, and a
#cat:           normalized value for the maximum power.  In addition, the
#cat:           statistics are ranked in descending order based on normalized
#cat:           squared maximum power.  These statistics are fundamental
#cat:           to selecting a dominant direction flow for the current
#cat:           input image block.
232 233

   Input:
234 235 236 237 238 239 240 241 242
      powers   - DFT power vectors (N Waves X M Directions) computed for
                 the current image block from which the values in the
                 statistics arrays are derived
      fw       - the beginning of the range of wave form indices from which
                 the statistcs are to derived
      tw       - the ending of the range of wave form indices from which
                 the statistcs are to derived (last index is tw-1)
      ndirs    - number of orientations (directions) at which the DFT
                 analysis was conducted
243
   Output:
244 245 246 247 248 249 250 251 252 253
      wis      - list of ranked wave form indicies of the corresponding
                 statistics based on normalized squared maximum power. These
                 indices will be used as indirect addresses when processing
                 the power statistics in descending order of "dominance"
      powmaxs  - array holding the maximum DFT power for each wave form
                 (other than the lowest frequecy)
      powmax_dirs - array to holding the direction corresponding to
                  each maximum power value in powmaxs
      pownorms - array to holding the normalized maximum powers corresponding
                 to each value in powmaxs
254 255 256 257
   Return Code:
      Zero     - successful completion
      Negative - system error
**************************************************************************/
258 259 260
int dft_power_stats(int *wis, double *powmaxs, int *powmax_dirs,
                     double *pownorms, double **powers,
                     const int fw, const int tw, const int ndirs)
261
{
262 263
   int w, i;
   int ret; /* return code */
Daniel Drake's avatar
Daniel Drake committed
264

265 266 267
   for(w = fw, i = 0; w < tw; w++, i++){
      get_max_norm(&(powmaxs[i]), &(powmax_dirs[i]),
                    &(pownorms[i]), powers[w], ndirs);
Daniel Drake's avatar
Daniel Drake committed
268 269
   }

270 271 272
   /* Get sorted order of applied DFT waves based on normalized power */
   if((ret = sort_dft_waves(wis, powmaxs, pownorms, tw-fw)))
      return(ret);
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297

   return(0);
}

/*************************************************************************
**************************************************************************
#cat: get_max_norm - Analyses a DFT power vector for a specific wave form
#cat:                applied at different orientations (directions) to the
#cat:                current image block.  The routine retuns the maximum
#cat:                power value in the vector, the direction at which the
#cat:                maximum occurs, and a normalized power value.  The
#cat:                normalized power is computed as the maximum power divided
#cat:                by the average power across all the directions.  These
#cat:                simple statistics are fundamental to the selection of
#cat:                a dominant direction flow for the image block.

   Input:
      power_vector - the DFT power values derived form a specific wave form
                     applied at different directions
      ndirs      - the number of directions to which the wave form was applied
   Output:
      powmax     - the maximum power value in the DFT power vector
      powmax_dir - the direciton at which the maximum power value occured
      pownorm    - the normalized power corresponding to the maximum power
**************************************************************************/
298
void get_max_norm(double *powmax, int *powmax_dir,
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
               double *pownorm, const double *power_vector, const int ndirs)
{
   int dir;
   double max_v, powsum;
   int max_i;
   double powmean;

   /* Find max power value and store corresponding direction */
   max_v = power_vector[0];
   max_i = 0;

   /* Sum the total power in a block at a given direction */
   powsum = power_vector[0];

   /* For each direction ... */
   for(dir = 1; dir < ndirs; dir++){
      powsum += power_vector[dir];
      if(power_vector[dir] > max_v){
         max_v = power_vector[dir];
         max_i = dir;
      }
   }

   *powmax = max_v;
   *powmax_dir = max_i;

   /* Powmean is used as denominator for pownorm, so setting  */
   /* a non-zero minimum avoids possible division by zero.    */
   powmean = max(powsum, MIN_POWER_SUM)/(double)ndirs;

   *pownorm = *powmax / powmean;
}

/*************************************************************************
**************************************************************************
#cat: sort_dft_waves - Creates a ranked list of DFT wave form statistics
#cat:                  by sorting on the normalized squared maximum power.

   Input:
      powmaxs  - maximum DFT power for each wave form used to derive
                 statistics
      pownorms - normalized maximum power corresponding to values in powmaxs
      nstats   - number of wave forms used to derive statistics (N Wave - 1)
   Output:
      wis      - sorted list of indices corresponding to the ranked set of
                 wave form statistics.  These indices will be used as
                 indirect addresses when processing the power statistics
                 in descending order of "dominance"
   Return Code:
      Zero     - successful completion
      Negative - system error
**************************************************************************/
351
int sort_dft_waves(int *wis, const double *powmaxs, const double *pownorms,
352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379
                   const int nstats)
{
   int i;
   double *pownorms2;

   /* Allocate normalized power^2 array */
   pownorms2 = (double *)malloc(nstats * sizeof(double));
   if(pownorms2 == (double *)NULL){
      fprintf(stderr, "ERROR : sort_dft_waves : malloc : pownorms2\n");
      return(-100);
   }

   for(i = 0; i < nstats; i++){
      /* Wis will hold the sorted statistic indices when all is done. */
      wis[i] = i;
      /* This is normalized squared max power. */
      pownorms2[i] = powmaxs[i] * pownorms[i];
   }

   /* Sort the statistic indices on the normalized squared power. */
   bubble_sort_double_dec_2(pownorms2, wis, nstats);

   /* Deallocate the working memory. */
   free(pownorms2);

   return(0);
}