The cycles per degree

主要是解释wiki中的数据 (http://en.wikipedia.org/wiki/Eye)：

Visual acuity, or resolving power, is "theability to distinguish fine detail" and is the property of cone cells.^[34] It is often measured in cyclesper degree (CPD),which measures an angular resolution, or how much aneye can differentiate one object from another in terms of visual angles.Resolution in CPD can be measured by bar charts of different numbers ofwhite/black stripe cycles. For example, if each pattern is 1.75 cm wide andis placed at 1 m distance from the eye, it will subtend an angle of 1 degree, so the numberof white/black bar pairs on the pattern will be a measure of the cycles perdegree of that pattern. The highest such number that the eye can resolve asstripes, or distinguish from a grey block, is then the measurement of visualacuity of the eye.

For a human eye with excellent acuity,the maximumtheoretical resolution is 50 CPD^[35] (1.2 arcminute perline pair, or a 0.35 mm line pair, at 1 m). A rat can resolve only about 1 to 2CPD.^[36] A horsehas higher acuity through most of the visual field of its eyes than a humanhas, but does not match the high acuity of the human eye's central fovearegion.^{[citation needed]}

第一处，可以看做tan(1) = 0.017455，也可以看成2*tan(0.5) = 0.017453。

第二处的50CPD为(1.75 cm)/(0.35 mm)

 

计算的话参照[4]：

% 2) Example  parameters
% Viewing parameters:
% a) image resolution in pixels per degree

  rowpixperdeg = 30 ; colpixperdeg = 30 ;

%  Image width or length is 16/30. =  0.5333 deg
%  Grating frequency is 5.625  cycles per deg
%  ((3 cycles per image/16 pixels per image)) * 30 pixels per deg )

% The pixel area in degrees squared is

  pixelArea = 1/(rowpixperdeg*colpixperdeg) ;

% b) contrast sensitivity filter
% contrast sensitivity function with unity peak gain
%  The Difference of Gaussians filter with unity peak gain has
%    three parameters:
%  i) the frequency cutoff in cycles per degree

  freqcutcpd = 15.4 ;

%  ii) the ratio of surround spread to center spread

  spreadratio = 15.4/1.35 ;

%  iii) the ratio of surround DC amplitude to center DC amplitude

  ampratio = 0.76 ;

%   The above parameter values are a DOG fit to Barten's contrast sensitivity
%    function for a luminance of 50 cd/m^2

%  Using the image size and resolution parameters,
%  the frequency cut in cycles per degree is converted to
%  row and column frequency cut offs in cycles per image.

  rowfreqcut = rows * freqcutcpd / rowpixperdeg ;
  colfreqcut = cols * freqcutcpd / colpixperdeg ;

%  The unscaled contrast sensitivity filter  is
 
  csf = filtdog(rows, cols, rowfreqcut, colfreqcut, ampratio, spreadratio) ;

自己做的：

% ChenVarshney.m
% -------------------------------------------------------------------
% Reference: <A human perception inspired quality metric for image 
%               fusion based on regional information>
%            <Objective assessment of multiresolution image fusion algorithms 
%               for context enhancement in Night vision: A comparative study>
%             http://www.mathworks.com/matlabcentral/fileexchange/3689-wpsnr/content/wpsnr.m
%             http://vision.arc.nasa.gov/personnel/al/code/matlab/filtmod2.htm
% Authors: Sun Li
% Date:    03/12/2014
% Last modified: 03/12/2014
% -------------------------------------------------------------------

% R : the human eye resolution (pixels/deg)
function qcv = ChenVarshney(img1, img2, imgF, R)
    
    img1 = double(img1);
    img2 = double(img2);
    imgF = double(imgF);
    
    alpha = 2;
    wS = 15;
    
    [gx1, gy1] = GradientSobel(img1);
    [gx2, gy2] = GradientSobel(img2);
    
    g1 = sqrt(gx1.^2+gy1.^2);
    g2 = sqrt(gx2.^2+gy2.^2);
    
    fun = @(block_struct) sum(block_struct.data(:).^alpha);
    lambda1 = blockproc(g1, [wS, wS], fun);
    lambda2 = blockproc(g2, [wS, wS], fun);
    
    f1 = img1-imgF;
    f2 = img2-imgF;
    
    funFilter = @(block_struct) MeanSquareCSFFilter(block_struct.data, R);
    
    D1F = blockproc(f1, [wS, wS], funFilter);
    D2F = blockproc(f2, [wS, wS], funFilter);
    
    
    qcv = sum(lambda1(:).*D1F(:) + lambda2(:).*D2F(:))/sum(lambda1(:)+lambda2(:)+eps);
    
end


%% --------------------------------------------
function fc = MannosSkarision(r)     
    fc = 2.6*(0.0192+0.144*r).*exp(-(0.144*r).^1.1);
end

%% -------------------------------------------
function f = CSFFilter(img, R)

    ff = fft2(img);
    ffc = fftshift(ff);
    [m, n] = size(img);
    [xx, yy] = meshgrid(1:n, 1:m);
    xx = xx - (n+1)/2;
    yy = yy - (m+1)/2;
    
    xx = xx/n*R;
    yy = yy/m*R;
    fw = MannosSkarision(sqrt(xx.^2+yy.^2));
    
    f = abs(ifft2(ifftshift(ffc.*fw)));
end

function f = MeanSquareCSFFilter(img, R)
    
    fk = CSFFilter(img, R);
    f = sum(fk(:).^2)/numel(img);% afriad there complex number
end

function [dx, dy] = GradientSobel(img)
    
        Sv = [-1 -2 -1;...
               0  0  0;...
               1  2  1];
        Sh = [-1  0  1;
              -2  0  2;
              -1  0  1];
        dx=conv2(img, Sh,'same'); 
        dy=conv2(img, Sv,'same');

end

 

参考：

【1】http://www.dsprelated.com/showmessage/24539/1.php

【2】http://stackoverflow.com/questions/23564333/define-different-spatial-frequency-thresholds-in-cpd-to-generate-several-image

【3】http://www.normankoren.com/Tutorials/MTF.html

【4】http://vision.arc.nasa.gov/personnel/al/code/matlab/filtmod2.htm