/*############################################################################# * 文件名：imagemanip.c * 功能：  實現了主要的圖像處理操作 * modified by  PRTsinghua@hotmail.com#############################################################################*/#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;