图像增强源代码（Matlab实现）

%-------------------------------------------------------------------------- %fft_enhance_cubs %enhances the fingerprint image %syntax: %[oimg,fimg,bwimg,eimg,enhimg] = fft_enhance_cubs(img) %oimg - [OUT] block orientation image(can be viewed using % view_orientation_image.m) %fimg - [OUT] block frequency image(indicates ridge spacing) %bwimg - [OUT] shows angular bandwidth image(filter bandwidth adapts near the % singular points) %eimg - [OUT] energy image. Indicates the 'ridgeness' of a block (can be % used for fingerprint segmentation) %enhimg- [OUT] enhanced image %img - [IN] input fingerprint image (HAS to be of DOUBLE type) %Contact: % ssc5@cubs.buffalo.edu sharat@mit.edu % http://www.sharat.org %Reference: %1. S. Chikkerur, C.Wu and V. Govindaraju, "Systematic approach for feature % extraction in Fingerprint Images", ICBA 2004 %2. S. Chikkerur and V. Govindaraju, "Fingerprint Image Enhancement using % STFT Analysis", International Workshop on Pattern Recognition for Crime % Prevention, Security and Surveillance, ICAPR 2005 %3. S. Chikkeur, "Online Fingerprint Verification", M. S. Thesis, % University at Buffalo, 2005 %4. T. Jea and V. Govindaraju, "A Minutia-Based Partial Fingerprint Recognition System", % to appear in Pattern Recognition 2005 %5. S. Chikkerur, "K-plet and CBFS: A Graph based Fingerprint % Representation and Matching Algorithm", submitted, ICB 2006 % See also: cubs_visualize_template %-------------------------------------------------------------------------- function [oimg,fimg,bwimg,eimg,enhimg] = fft_enhance_cubs(img) global NFFT; NFFT = 32; %size of FFT BLKSZ = 12; %size of the block OVRLP = 6; %size of overlap ALPHA = 0.4; %root filtering RMIN = 5; %min allowable ridge spacing RMAX = 20; %maximum allowable ridge spacing do_prefiltering = 1; [nHt,nWt] = size(img); img = im2double(img); %convert to DOUBLE nBlkHt = floor((nHt-2*OVRLP)/BLKSZ); nBlkWt = floor((nWt-2*OVRLP)/BLKSZ); fftSrc = zeros(nBlkHt*nBlkWt,NFFT*NFFT); %stores FFT nWndSz = BLKSZ+2*OVRLP; %size of analysis window. warning off MATLAB:divideByZero %------------------------- %allocate outputs %------------------------- oimg = zeros(nBlkHt,nBlkWt); fimg = zeros(nBlkHt,nBlkWt); bwimg = zeros(nBlkHt,nBlkWt); eimg = zeros(nBlkHt,nBlkWt); enhimg = zeros(nHt,nWt); %------------------------- %precomputations %------------------------- [x,y] = meshgrid(0:nWndSz-1,0:nWndSz-1); dMult = (-1).^(x+y); %used to center the FFT [x,y] = meshgrid(-NFFT/2:NFFT/2-1,-NFFT/2:NFFT/2-1); r = sqrt(x.^2+y.^2)+eps; th = atan2(y,x); th(th<0) = th(th<0)+pi; w = raised_cosine_window(BLKSZ,OVRLP); %spectral window %------------------------- %FFT Analysis %------------------------- for i = 0:nBlkHt-1 nRow = i*BLKSZ+OVRLP+1; for j = 0:nBlkWt-1 nCol = j*BLKSZ+OVRLP+1; %extract local block blk = img(nRow-OVRLP:nRow+BLKSZ+OVRLP-1,nCol-OVRLP:nCol+BLKSZ+OVRLP-1); %remove dc dAvg = sum(sum(blk))/(nWndSz*nWndSz); blk = blk-dAvg; %remove DC content blk = blk.*w; %multiply by spectral window %-------------------------- %do pre filtering %-------------------------- blkfft = fft2(blk.*dMult,NFFT,NFFT); if(do_prefiltering) dEnergy = abs(blkfft).^2; blkfft = blkfft.*sqrt(dEnergy); %root filtering(for diffusion) end; fftSrc(nBlkWt*i+j+1,:) = transpose(blkfft(:)); dEnergy = abs(blkfft).^2; %----REDUCE THIS COMPUTATION---- %-------------------------- %compute statistics %-------------------------- dTotal = sum(sum(dEnergy));%/(NFFT*NFFT); fimg(i+1,j+1) = NFFT/(compute_mean_frequency(dEnergy,r)+eps); %ridge separation oimg(i+1,j+1) = compute_mean_angle(dEnergy,th); %ridge angle eimg(i+1,j+1) = log(dTotal+eps); %used for segmentation end;%for j end;%for i %------------------------- %precomputations %------------------------- [x,y] = meshgrid(-NFFT/2:NFFT/2-1,-NFFT/2:NFFT/2-1); dMult = (-1).^(x+y); %used to center the FFT %------------------------- %process the resulting maps %------------------------- for i = 1:3 oimg = smoothen_orientation_image(oimg); %smoothen orientation image end; fimg = smoothen_frequency_image(fimg,RMIN,RMAX,5); %diffuse frequency image cimg = compute_coherence(oimg); %coherence image for bandwidth bwimg = get_angular_bw_image(cimg); %QUANTIZED bandwidth image %------------------------- %FFT reconstruction %------------------------- for i = 0:nBlkHt-1 for j = 0:nBlkWt-1 nRow = i*BLKSZ+OVRLP+1; nCol = j*BLKSZ+OVRLP+1; %-------------------------- %apply the filters %-------------------------- blkfft = reshape(transpose(fftSrc(nBlkWt*i+j+1,:)),NFFT,NFFT); %-------------------------- %reconstruction %-------------------------- af = get_angular_filter(oimg(i+1,j+1),bwimg(i+1,j+1)); rf = get_radial_filter(fimg(i+1,j+1)); blkfft = blkfft.*(af).*(rf); blk = real(ifft2(blkfft).*dMult); enhimg(nRow:nRow+BLKSZ-1,nCol:nCol+BLKSZ-1)=blk(OVRLP+1:OVRLP+BLKSZ,OVRLP+1:OVRLP+BLKSZ); end;%for j end;%for i %end block processing %-------------------------- %contrast enhancement %-------------------------- enhimg =sqrt(abs(enhimg)).*sign(enhimg); enhimg =imscale(enhimg); enhimg =im2uint8(enhimg); %-------------------------- %clean up the image %-------------------------- emsk = segment_print(enhimg,0); enhimg(emsk==0) = 128; %end function fft_enhance_cubs %----------------------------------- %raised_cosine %returns 1D spectral window %syntax: %y = raised_cosine(nBlkSz,nOvrlp) %y - [OUT] 1D raised cosine function %nBlkSz - [IN] the window is constant here %nOvrlp - [IN] the window has transition here %----------------------------------- function y = raised_cosine(nBlkSz,nOvrlp) nWndSz = (nBlkSz+2*nOvrlp); x = abs(-nWndSz/2:nWndSz/2-1); y = 0.5*(cos(pi*(x-nBlkSz/2)/nOvrlp)+1); y(abs(x)<nBlkSz/2)=1; %end function raised_cosine %----------------------------------- %raised_cosine_window %returns 2D spectral window %syntax: %w = raised_cosine_window(blksz,ovrlp) %w - [OUT] 1D raised cosine function %nBlkSz - [IN] the window is constant here %nOvrlp - [IN] the window has transition here %----------------------------------- function w = raised_cosine_window(blksz,ovrlp) y = raised_cosine(blksz,ovrlp); w = y(:)*y(:)'; %end function raised_cosine_window %--------------------------------------------------------------------- %get_angular_filter %generates an angular filter centered around 'th' and with bandwidth 'bw' %the filters in angf_xx are precomputed using angular_filter_bank.m %syntax: %r = get_angular_filter(t0,bw) %r - [OUT