imagemanip.c

用VISUAL C++编程实现指纹图像的特征提取以及对指纹图像的识别

/*########################################################################

  The contents of this file are subject to the Mozilla Public License
  Version 1.0(the "License");   You  may  NOT  use this file except in
  compliance with the License. You may obtain a copy of the License at
                http:// www.mozilla.org/MPL/
  Software distributed under the License is distributed on an "AS IS"
  basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See
  the License for the specific language governing rights and limitations
  under the License.

  The Initial Developer of the Original Code is Shivang Patel.

  Copyright(C) 2002. All Rights Reserved.

  Authors: Shivang Patel
           Jaap de Haan(jdh)
  
  This file contains functions that manipulate the image

########################################################################*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

#include "imagemanip.h"

#ifndef min
#define min(a,b) (((a)<(b))?(a):(b))
#endif

/* {{{ Local stretching */

/* helper macro */
#define PIJKL p[i+k + (j+l)*nSizeX]

/* local stretch */
FvsError_t ImageLocalStretch(FvsImage_t image, const FvsInt_t size, const FvsInt_t tolerance)
{
    /* define a bunch of variables */
    int nSizeX = ImageGetWidth(image)  - size + 1;
    int nSizeY = ImageGetHeight(image) - size + 1;
    FvsInt_t i, j, t, l;
    FvsInt_t sum, denom;
    FvsByte_t a = 0;
    FvsInt_t k = 0;
    FvsByte_t b = 255;
    int hist[256];
    FvsByte_t* p = ImageGetBuffer(image);
    if (p==NULL)
       return FvsMemory;
    for (j=0; j < nSizeY; j+=size)
    {
        for (i=0; i < nSizeX; i+=size)
        {
            /* compute local histogram */
            memset(hist, 0, 256*sizeof(int));
            for (l = 0; l<size; l++)
                for (k = 0; k<size; k++)
                    hist[PIJKL]++;
            
            /* stretch locally */
            for (k=0,   sum=0;   k <256; k++)
            {
                sum+=hist[k];
                a = (FvsByte_t)k;
                if (sum>tolerance) break;
            }
            
            for (k=255, sum=0; k >= 0; k--)
            {
                sum+=hist[k];
                b = (FvsByte_t)k;
                if (sum>tolerance) break;
            }

            denom = (FvsInt_t)(b-a);
            if (denom!=0)
            {
                for (l = 0; l<size; l++)
                {
                    for (k = 0; k<size; k++)
                    {
                        if (PIJKL<a) PIJKL = a;
                        if (PIJKL>b) PIJKL = b;
                        t = (FvsInt_t)((((PIJKL)-a)*255)/denom);
                        PIJKL = (FvsByte_t)(t);
                    }
                }
            }
        }
    }

    return FvsOK;
}
/* }}} */
/* {{{ Fingerprint orientation field */

#define P(x,y)      ((int32_t)p[(x)+(y)*pitch])

/*
** In this step, we estimate the ridge orientation field.
** Given a normalized image G, the main steps of the algorithm are as
** follows:
**
** 1 - Divide G into blocks of w x w - (15 x 15)
**
** 2 - Compute the gradients dx(i,j) and dy(i,j) at each pixel (i,j),
**     depending on the computational requirement, the gradient operator
**     may vary from the single Sobel operator to the more complex Marr-
**     Hildreth operator.
**
** 3 - Estimate the local orientation of each block centered at pixel
**     (i,j), using the following operations:
**
**               i+w/2 j+w/2
**               ---   --- 
**               \     \
**     Nx(i,j) =  --    -- 2 dx(u,v) dy(u,v)
**               /     /
**               ---   ---
**            u=i-w/2 v=j-w/2
**
**               i+w/2 j+w/2
**               ---   --- 
**               \     \
**     Ny(i,j) =  --    -- dx²(u,v) - dy²(u,v)
**               /     /
**               ---   ---
**            u=i-w/2 v=j-w/2
**
**                  1    -1  / Nx(i,j) \
**     Theta(i,j) = - tan   |  -------  |
**                  2        \ Ny(i,j) /
**
**     where Theta(i,j) is the least square estimate of the local ridge
**     orientation at the block centered at pixel (i,j). Mathematically,
**     it represents the direction that is orthogonal to the dominant
**     direction of the Fourier spectrum of the w x w window.
**
** 4 - Due to the presence of noise, corrupted ridge and furrow structures,
**     minutiae, etc. in the input image,the estimated local ridge
**     orientation may not always be a correct estimate. Since local ridge
**     orientation varies slowly in a local neighbourhood, where no
**     singular point appears, a low pass filter can be used to modify the
**     incorrect local ridge orientation. In order to perform the low-pass
**     filtering, the orientation image needs to be converted into a
**     continuous vector field, which is defined as follows:
**       Phi_x(i,j) = cos( 2 x theta(i,j) )
**       Phi_y(i,j) = sin( 2 x theta(i,j) )
**     With the resulting vector field, the low-pass filtering can then
**     be performed with a convolution as follows:
**       Phi2_x(i,j) = (W @ Phi_x) (i,j)
**       Phi2_y(i,j) = (W @ Phi_y) (i,j)
**     where W is a 2D low-pass filter.
**
** 5 - Compute the local ridge orientation at (i,j) with
**
**              1    -1  / Phi2_y(i,j) \
**     O(i,j) = - tan   |  -----------  |
**              2        \ Phi2_x(i,j) /
**
** With this algorithm, a fairly smooth orientatin field estimate can be
** obtained.
**
*/

static FvsError_t FingerprintDirectionLowPass(FvsFloat_t* theta, FvsFloat_t* out,
            FvsInt_t nFilterSize, FvsInt_t w, FvsInt_t h)
{
    FvsError_t nRet = FvsOK;
    FvsFloat_t* filter = NULL;
    FvsFloat_t* phix   = NULL;
    FvsFloat_t* phiy   = NULL;
    FvsFloat_t* phi2x  = NULL;
    FvsFloat_t* phi2y  = NULL;
    FvsInt_t fsize  = nFilterSize*2+1;
    size_t nbytes = (size_t)(w*h*sizeof(FvsFloat_t));
    FvsFloat_t nx, ny;
    FvsInt_t val;
    FvsInt_t i, j, x, y;

    filter= (FvsFloat_t*)malloc((size_t)fsize*fsize*sizeof(FvsFloat_t));
    phix  = (FvsFloat_t*)malloc(nbytes);
    phiy  = (FvsFloat_t*)malloc(nbytes);
    phi2x = (FvsFloat_t*)malloc(nbytes);
    phi2y = (FvsFloat_t*)malloc(nbytes);

    if (filter==NULL || phi2x==NULL || phi2y==NULL || phix==NULL || phiy==NULL)
        nRet = FvsMemory;
    else
    {
        /* reset all fields to 0 */
        memset(filter, 0, (size_t)fsize*fsize*sizeof(FvsFloat_t));
        memset(phix,   0, nbytes);
        memset(phiy,   0, nbytes);
        memset(phi2x,  0, nbytes);
        memset(phi2y,  0, nbytes);

        /* 4 - Compute a continuous field from theta */
        for (y = 0; y < h; y++)
        for (x = 0; x < w; x++)
        {
            val = x+y*w;
            phix[val] = cos(theta[val]);
            phiy[val] = sin(theta[val]);
        }
        /* build the low-pass filter */
        nx = 0.0;
        for (j = 0; j < fsize; j++)
        for (i = 0; i < fsize; i++)
        {
            filter[j*fsize+i] = 1.0;
/*
            filter[j*fsize+i] = (FvsFloat_t)(fsize - (abs(nFilterSize-i)+abs(nFilterSize-j)));
*/
            nx += filter[j*fsize+i]; /* sum of coefficients */
        }
        if (nx>1.0)
        {
            for (j = 0; j < fsize; j++)
            for (i = 0; i < fsize; i++)
                /* normalize the result */
                filter[j*fsize+i] /= nx;
        }
        /* low-pass on the result arrays getting phi2 */
        for (y = 0; y < h-fsize; y++)
        for (x = 0; x < w-fsize; x++)
        {
            nx = 0.0;
            ny = 0.0;
            for (j = 0; j < fsize; j++)
            for (i = 0; i < fsize; i++)
            {
                val = (x+i)+(j+y)*w;
                nx += filter[j*fsize+i]*phix[val];
                ny += filter[j*fsize+i]*phiy[val];
            }
            val = x+y*w;
            phi2x[val] = nx;
            phi2y[val] = ny;
        }
        /* we do not need phix, phiy anymore, delete them */
        if (phix!=NULL) { free(phix); phix=NULL; }
        if (phiy!=NULL) { free(phiy); phiy=NULL; }

        /* 5 - local ridge orientation -> theta */
        for (y = 0; y < h-fsize; y++)
        for (x = 0; x < w-fsize; x++)
        {
            val = x+y*w;
            out[val] = atan2(phi2y[val], phi2x[val])*0.5;
        }
    }

    if (phix!=NULL)  free(phix);
    if (phiy!=NULL)  free(phiy);
    if (phi2x!=NULL) free(phi2x);
    if (phi2y!=NULL) free(phi2y);
    if (filter!=NULL)free(filter);

    return nRet;
}


FvsError_t FingerprintGetDirection(const FvsImage_t image, FvsFloatField_t field,
            const FvsInt_t nBlockSize, const FvsInt_t nFilterSize)
{
    /* width & height of the input image */
    FvsInt_t w       = ImageGetWidth (image);
    FvsInt_t h       = ImageGetHeight(image);
    FvsInt_t pitch   = ImageGetPitch (image);
    FvsByte_t* p     = ImageGetBuffer(image);
    FvsInt_t i, j, u, v, x, y;
    FvsInt_t s = nBlockSize*2+1;
    FvsFloat_t *dx = malloc(sizeof(FvsFloat_t) * s*s);
    FvsFloat_t *dy = malloc(sizeof(FvsFloat_t) * s*s);
    FvsFloat_t nx, ny;
    FvsFloat_t* out;
    FvsFloat_t* theta  = NULL;
    FvsError_t nRet = FvsOK;

    /* (re-)allocate the output image */
    nRet = FloatFieldSetSize(field, w, h);
	if (nRet!=FvsOK)
	{
		free(dx); free(dy);
		return nRet;
	}
    nRet = FloatFieldClear(field);
    if (nRet!=FvsOK)
	{
		free(dx); free(dy);
		return nRet;
	}
    out = FloatFieldGetBuffer(field);

    /* allocate memory for the orientation values */
    if (nFilterSize>0)
    {
        theta = (FvsFloat_t*)malloc(w * h * sizeof(FvsFloat_t));
        if (theta!=NULL)
            memset(theta, 0, (w * h * sizeof(FvsFloat_t)));
    }

    /* detect any allocation error */
    if (out==NULL || (nFilterSize>0 && theta==NULL))
        nRet = FvsMemory;
    else
    {
        /* 1 - divide the image in blocks */
        for (y = nBlockSize+1; y < h-nBlockSize-1; y++)
        for (x = nBlockSize+1; x < w-nBlockSize-1; x++)
        {
            /* 2 - for the block centered at x,y compute the gradient */
            for (j = 0; j < s; j++)
            for (i = 0; i < s; i++)
            {
                dx[i*s+j] = (FvsFloat_t)
                           (P(x+i-nBlockSize,   y+j-nBlockSize) -
                            P(x+i-nBlockSize-1, y+j-nBlockSize));
                dy[i*s+j] = (FvsFloat_t)
                           (P(x+i-nBlockSize,   y+j-nBlockSize) -
                            P(x+i-nBlockSize,   y+j-nBlockSize-1));
            }

            /* 3 - compute orientation */
            nx = 0.0;
            ny = 0.0;
            for (v = 0; v < s; v++)
            for (u = 0; u < s; u++)
            {
                nx += 2 * dx[u*s+v] * dy[u*s+v];
                ny += dx[u*s+v]*dx[u*s+v] - dy[u*s+v]*dy[u*s+v];
            }
            /* compute angle (-pi/2 .. pi/2) */
            if (nFilterSize>0)
                theta[x+y*w] = atan2(nx, ny);
            else
                out[x+y*w] = atan2(nx, ny)*0.5;
        }
        if (nFilterSize>0)
            nRet = FingerprintDirectionLowPass(theta, out, nFilterSize, w, h);
    }

    if (theta!=NULL) free(theta);

	free(dx); free(dy);
	
	return nRet;
}

/* }}} */
/* {{{ Fingerprint frequency field */

/* 
** In this step, we estimate the ridge frequency. In a local neighbour-
** hood where no minutiae and singular points appear, the gray levels
** along ridges and furrows can be modelled as a sinusoidal-shaped wave
** along a direction normal to the local ridge orientation. Therefore,
** local ridge frequency is another intrinsic property of a fingerprint
** image. Let G be a normalized image G, and O be the orientation image
** (computed at step B). Then the steps involved in local ridge 
** frequency estimation are as follows:
**
** 1 - Divide G into blocks of w x w - (16 x 16)
**
** 2 - For each block centered at pixel (i,j), compute an oriented
**     window of size l x w (32 x 16) that is defined in the ridge
**     coordinates system.
**
** 3 - For each block centered at pixel (i,j), compute the x-signature
**     X[0], X[1], ... X[l-1] of the ridges and furrows within the
**     oriented window where:
**
**              --- w-1
**            1 \
**     X[k] = -  --  G (u, v), k = 0, 1, ..., l-1
**            w /
**              --- d=0
**
**     u = i + (d - w/2).cos O(i,j) + (k - l/2).sin O(i,j)
**
**     v = j + (d - w/2).sin O(i,j) - (k - l/2).cos O(i,j)
**
**     If no minutiae and singular points appear in the oriented window,
**     the x-signature forms a discrete sinusoidal-shape wave, which has
**     the same frequency as that of the ridges and furrows in the
**     oriented window. Therefore, the frequency of ridges and furrows
**     can be estimated from the x-signature. Let T(i,j) be the average
**     number of pixels between two consecutive peaks in the x-signature,
**     then the frequency, OHM(i,j) is computed as OHM(i,j) = 1 / T(i,j).
**
**     If no consecutive peaks can be detected from the x-signature, then
**     the frequency is assigned a value of -1, to differenciate it from
**     the valid frequency values.
**
** 4 - For a fingerprint image scanned at a fixed resolution, the value
**     of the frequency of the ridges and furrows in a local neighbour-
**     hood lies in a certain range. For a 500 dpi image, this range is
**     [1/3, 1/25]. Therefore, if the estimated value of the frequency
**     is out of this range, then the frequency is assigned a value of
**     -1 to indicate that an valid frequency cannot be obtained.
**
** 5 - The blocks in which minutiae and/or singular points appear and/or
**     ridges and furrows are corrupted do not form a well defined
**     sinusoidal-shaped wave. The frequency value for these blocks need
**     to be interpolated from the frequency of the neighbouring blocks
**     which have a well-defined frequency. (Ex: Gaussian kernel mean 0,
**     and variance 9 and of size 7).
**
** 6 - Inter-ridges distance change slowly in a local neighbourhood. A
**     low-pass filter can be used to remove the outliers.
**
*/

/* width */
#define BLOCK_W     16
#define BLOCK_W2     8

/* length */
#define BLOCK_L     32
#define BLOCK_L2    16

#define EPSILON     0.0001
#define LPSIZE      3

#define LPFACTOR    (1.0/((LPSIZE*2+1)*(LPSIZE*2+1)))


FvsError_t FingerprintGetFrequency(const FvsImage_t image, const FvsFloatField_t direction,
           FvsFloatField_t frequency)
{
    /* width & height of the input image */
    FvsError_t nRet = FvsOK;
    FvsInt_t w      = ImageGetWidth (image);
    FvsInt_t h      = ImageGetHeight(image);
    FvsInt_t pitchi = ImageGetPitch (image);
    FvsByte_t* p    = ImageGetBuffer(image);
    FvsFloat_t* out;
    FvsFloat_t* freq;
    FvsFloat_t* orientation = FloatFieldGetBuffer(direction);

    FvsInt_t x, y, u, v, d, k;
    size_t size;

    if (p==NULL)
        return FvsMemory;

    /* (re-)allocate the output image */
    nRet = FloatFieldSetSize(frequency, w, h);
    if (nRet!=FvsOK) return nRet;
    (void)FloatFieldClear(frequency);
    freq = FloatFieldGetBuffer(frequency);
    if (freq==NULL)
        return FvsMemory;

    /* allocate memory for the output */
    size = w*h*sizeof(FvsFloat_t);
    out  = (FvsFloat_t*)malloc(size);
    if (out!=NULL)