imagemanip.c

很好的一套指纹识别系统C源码

/*#############################################################################
 * 文件名：imagemanip.c
 * 功能：  实现了主要的图像处理操作
 * modified by  [e-mail]
#############################################################################*/

#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


/* 宏定义 */
#define PIJKL p[i+k + (j+l)*nSizeX]


/******************************************************************************
  * 功能：图像缩放操作
  * 参数：image       指纹图像
  *       size        缩放的图像块大小
  *       tolerance   消去直方图的边界
  * 返回：错误编号
******************************************************************************/
FvsError_t ImageLocalStretch(FvsImage_t image, const FvsInt_t size, 
						const FvsInt_t tolerance)
{
    /* 定义一些变量 */
    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)
        {
            /* 计算直方图 */
            memset(hist, 0, 256*sizeof(int));
            for (l = 0; l<size; l++)
                for (k = 0; k<size; k++)
                    hist[PIJKL]++;
            
            /* 伸缩 */
            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;
}


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

/******************************************************************************
** 估算脊线的方向
** 给定一个归一化的指纹图像，算法的主要步骤如下：
**
** 1 - 将G分成大小为 w x w - (15 x 15) 的块；
**
** 2 - 计算每个象素 (i,j)的梯度 dx(i,j) 和 dy(i,j) ,
**     根据计算的需求，梯度算子可以从简单的Sobel算子到复杂的Marr-Hildreth 算子。
**
** 3 - 估算优势方向(i,j), 使用如下的操作：
**
**               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) /
**
**     这里，Theta(i,j)是局部脊线方向的最小方差估计，以像素 (i,j) 为中心。
**     从数学的角度看，它代表傅立叶频谱中直角占有时的方向。
**
** 4 - 由于有噪声，脊线的中断，细节点等等的存在，在输入图像中，对局部脊线
**     方向的估计并不总是正确的。由于局部脊线方向变化缓慢，所以可以用低通
**     滤波器来修正不正确的脊线方向。为了运用低通滤波器，方向图必须转换成
**     连续的矢量域，定义如下：
**       Phi_x(i,j) = cos( 2 x theta(i,j) )
**       Phi_y(i,j) = sin( 2 x theta(i,j) )
**     在矢量域，可以用如下的卷积低通滤波：
**       Phi2_x(i,j) = (W @ Phi_x) (i,j)
**       Phi2_y(i,j) = (W @ Phi_y) (i,j)
**     W是一个二维的低通滤波器。
**
** 5 - 用如下公式计算 (i,j) 处的方向：
**
**              1    -1  / Phi2_y(i,j) \
**     O(i,j) = - tan   |  -----------  |
**              2        \ Phi2_x(i,j) /
**
** 用这个算法可以得到相当平滑的方向图
**
*/

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
    {
        /* 置 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 */
        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]);
        }
        /* 构造低通滤波器 */
        nx = 0.0;
        for (j = 0; j < fsize; j++)
        for (i = 0; i < fsize; i++)
        {
            filter[j*fsize+i] = 1.0;

            nx += filter[j*fsize+i]; /* 系数和 */
        }
        if (nx>1.0)
        {
            for (j = 0; j < fsize; j++)
            for (i = 0; i < fsize; i++)
                /* 归一化结果 */
                filter[j*fsize+i] /= nx;
        }
        /* 低通滤波 */
        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;
        }
        /* 销毁 phix, phiy */
        if (phix!=NULL) 
        { 
        	free(phix); 
        	phix=NULL; 
        }
        if (phiy!=NULL) 
        { 
        	free(phiy); 
        	phiy=NULL; 
        }

        /* 步骤5 */
        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;
}


/******************************************************************************
  * 功能：计算指纹图像脊线的方向。
          该算法在许多论文中都有描述，如果图像做了归一化，并且对比度较高，
          则最后的处理效果也较好。
          方向的值在-PI/2和PI/2之间，弧度和脊并不相同。
          选取的块越大，分析的效果也越好，但所需的处理计算时间也越长。
          由于指纹图像中脊线方向的变化比较缓慢，所以低通滤波器可以较好的
          过虑掉方向中的噪声和错误。
  * 参数：image          指向图像对象的指针
  *       field          指向浮点域对象的指针，保存结果
  *       nBlockSize     块大小
  *       nFilterSize    滤波器大小
  * 返回：错误编号
******************************************************************************/
FvsError_t FingerprintGetDirection(const FvsImage_t image, 
						FvsFloatField_t field, const FvsInt_t nBlockSize, 
						const FvsInt_t nFilterSize)
{
    /* 输入图像的宽度和高度 */
    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;
    FvsFloat_t dx[(nBlockSize*2+1)][(nBlockSize*2+1)];
    FvsFloat_t dy[(nBlockSize*2+1)][(nBlockSize*2+1)];
    FvsFloat_t nx, ny;
    FvsFloat_t* out;
    FvsFloat_t* theta  = NULL;
    FvsError_t nRet = FvsOK;

    /* 输出图像 */
    nRet = FloatFieldSetSize(field, w, h);
    if (nRet!=FvsOK) return nRet;
    nRet = FloatFieldClear(field);
    if (nRet!=FvsOK) return nRet;
    out = FloatFieldGetBuffer(field);

    /* 为方向数组申请内存 */
    if (nFilterSize>0)
    {
        theta = (FvsFloat_t*)malloc(w * h * sizeof(FvsFloat_t));
        if (theta!=NULL)
            memset(theta, 0, (w * h * sizeof(FvsFloat_t)));
    }

    /* 内存错误，返回 */
    if (out==NULL || (nFilterSize>0 && theta==NULL))
        nRet = FvsMemory;
    else
    {
        /* 1 - 图像分块 */
        for (y = nBlockSize+1; y < h-nBlockSize-1; y++)
        for (x = nBlockSize+1; x < w-nBlockSize-1; x++)
        {
            /* 2 - 计算梯度 */
            for (j = 0; j < (nBlockSize*2+1); j++)
            for (i = 0; i < (nBlockSize*2+1); i++)
            {
                dx[i][j] = (FvsFloat_t)
                           (P(x+i-nBlockSize,   y+j-nBlockSize) -
                            P(x+i-nBlockSize-1, y+j-nBlockSize));
                dy[i][j] = (FvsFloat_t)
                           (P(x+i-nBlockSize,   y+j-nBlockSize) -
                            P(x+i-nBlockSize,   y+j-nBlockSize-1));
            }

            /* 3 - 计算方向 */
            nx = 0.0;
            ny = 0.0;
            for (v = 0; v < (nBlockSize*2+1); v++)
            for (u = 0; u < (nBlockSize*2+1); u++)
            {
                nx += 2 * dx[u][v] * dy[u][v];
                ny += dx[u][v]*dx[u][v] - dy[u][v]*dy[u][v];
            }
            /* 计算角度 (-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);
    return nRet;
}


/* 指纹频率域 */

/******************************************************************************
** 这个步骤里，我们估计指纹脊线的频率。在局部邻域里，没有凸现的细节点或者孤点，
** 沿着脊线和谷底，可以用一个正弦曲线波形作为模型，因此，局部脊线频率是指纹图
** 像的另一个本质的特征。对指纹图像G进行归一化，O是其方向图，估算局部脊线频率
** 的步骤如下：
**
** 1 - 图像分块 w x w - (16 x 16)
**
** 2 - 对每块，计算大小为l x w (32 x 16)的方向图窗口
**
** 3 - 对中心在 (i,j) 的每块, 计算脊线和谷底的 x-signature
**     X[0], X[1], ... X[l-1] 采用如下公式：
**
**              --- 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)
**
**     如果方向图窗口中没有细节点和孤立的点，则x-signature形成了一个离散
**     的正弦曲线波，与方向图中脊线和谷底的频率一样。因此，脊线和谷底的
**     频率可以由x-signature来估计。设T(i,j)是两个峰顶的平均距离，则频率
**     OHM(i,j)可以这样计算：OHM(i,j) = 1 / T(i,j)。
**
**     如果没有两个连续的峰顶，则频率置为-1，说明其无效。
**
** 4 - 对于一个指纹图像而言，脊线频率的值在一个范围之内变动，比如说对于500
**     dpi的图像，变动范围为[1/3, 1/25]，因此，如果估计出的频率不在这个范
**     围内，说明频率估计无效，同意置为-1。
**
** 5 - 如果某块有断点或者细节点，则不会有正弦曲线，其频率可以由邻块的频率
**     插值估计（比如说高斯函数，均值为0，方差为9，宽度为7）。
**
** 6 - 脊线内部距离变化缓慢，可以用低通滤波器
**
*/

/* 宽度 */
#define BLOCK_W     16
#define BLOCK_W2     8

/* 长度 */
#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)
{
    /* 输入图像的宽度和高度 */
    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;

    /* 输出图像的内存申请 */
    nRet = FloatFieldSetSize(frequency, w, h);
    if (nRet!=FvsOK) return nRet;
    (void)FloatFieldClear(frequency);
    freq = FloatFieldGetBuffer(frequency);
    if (freq==NULL)
        return FvsMemory;

    /* 输出的内存申请 */
    size = w*h*sizeof(FvsFloat_t);
    out  = (FvsFloat_t*)malloc(size);
    if (out!=NULL)
    {
        FvsFloat_t dir = 0.0;
        FvsFloat_t cosdir = 0.0;
        FvsFloat_t sindir = 0.0;