%
function dprime = symboldiff1()
% Version of single filter model
% with masking by a non-homogeneous background.
% Adapted to assess discriminability of symbol pairs
% by including translation invariance
% and a range of size invariance.
% In this simple version the possibly large symbol must be in image 2.
% Local contrast is computed using a local luminance image.
% computed as a blurred mixture of the image and a background.
% Masking is computed using a local contrast energy image.
% Outline
% 1) Example input symbols
% 2) Example model parameters
% 3) Example computation of model output
% 1) Example inputs
% image 1 is a small o image 2 is a bigger u.
im1 = [ 0 0 0 0 0 0 0 0 0
0 1 1 1 1 1 1 0 0
0 1 1 1 1 1 1 0 0
0 1 1 0 0 1 1 0 0
0 1 1 0 0 1 1 0 0
0 1 1 1 1 1 1 0 0
0 1 1 1 1 1 1 0 0
0 0 0 0 0 0 0 0 0 ] ;
im2 = [ 0 1 1 1 0 0 0 1 1 1
0 1 1 1 0 0 0 1 1 1
0 1 1 1 0 0 0 1 1 1
0 1 1 1 0 0 0 1 1 1
0 1 1 1 0 0 0 1 1 1
0 1 1 1 1 1 1 1 1 1
0 1 1 1 1 1 1 1 1 1
0 1 1 1 1 1 1 1 1 1] ;
% Image size
rows = 16 ; cols = rows ;
image1 = surround(im1,0,rows,cols) ;
image2 = surround(im2,0,rows,cols) ;
% Image conditioning parameters
repfactor = 4; % allows subpixel accuracy on translation
% and increased accuracy on size ratio estimation
borderwidth = 0;
sizeratio = 0.5 ; % size minification range limit
% sizeratio = 1 ; % size minification range limit set to not vary size
image1 = double(replicate(image1,repfactor,repfactor));
image2 = double(replicate(image2,repfactor,repfactor));
% 2) Example model parameters
% Viewing parameters:
% a) image resolution in pixels per degree
rowpixperdeg = 30 ; colpixperdeg = 30 ;
rowpixperdeg = repfactor*rowpixperdeg;
colpixperdeg = repfactor*colpixperdeg ;
% b) contrast sensitivity filter parameters
% The Difference of Gaussians filter with unity peak gain has
% three parameters:
% i) the center spread in minutes of arc
centerspread = 2.0 ;
% or the frequency cutoff in cycles per degree
freqcutcpd = 60 /(sqrt(pi) * centerspread) ;
% ii) the ratio of surround spread to center spread
spreadratio = 0.42/(2/60) ;
% iii) the ratio of surround DC amplitude to center DC amplitude
ampratio = 0.8 ;
% The above parameter values are from Ahumada (2005 ECVP).
% Convert these parameters to parameters of the non-homogeneous model
% An 'optical' blur filter cutoff
opticalcut = freqcutcpd ;
% A luminance spread filter cutoff
lumspreadcut = opticalcut/spreadratio ;
% A contrast energy masking spread filter cutoff
maskspreadcut = lumspreadcut ;
% The background level or image preceding the images
imBackLum1 = image1(1,1) ;
imBackLum2 = image2(1,1) ;
% The proportion of image vs. background determining local luminance
pIm = ampratio ;
% Contrast sensitivity, masking, and spatial integration parameters
consensmax=500 ;% lowest threshold near 0 dBB for 0.25 sec duration
% see (Watson, 2000 Optics Express)
maskconthresh=0.05 ; % 5% rms contrast raises threshold by 3 dB
beta=2 ; %contrast energy summation
% 3) Example model output
% just noticeable differences between images
dprime = symboldiff(image1,image2, sizeratio, repfactor,...
rowpixperdeg, colpixperdeg, ...
imBackLum1, imBackLum2, pIm, ...
opticalcut, lumspreadcut, maskspreadcut, ...
consensmax, maskconthresh, beta) ;
% Example o vs u
% dprime best_size_ratio yoffset xoffset
% 4.5301 0.6875 0 1.5000
% Example dprime without size invariance
% 11.2822 1.0000 0.5000 1.5000
% --------------- begin symboldiff, the main subroutine
function dprime = symboldiff(image1,image2, sizeratio, repfactor,...
rowpixperdeg, colpixperdeg, ...
imBackLum1, imBackLum2, pIm, ...
opticalcut, lumspreadcut, maskspreadcut, ...
consensmax, maskconthresh, beta) ;
cols = size(image1,2) ;
rows = size(image1,1) ;
% convert image 1 to a visible contrast image
con1 = localconmask(image1,rowpixperdeg, colpixperdeg, ...
imBackLum1, pIm, ...
opticalcut, lumspreadcut, maskspreadcut, ...
maskconthresh) ;
dat = []; % will collect the dprimes as function of minification
for irows = rows:-2:floor(sizeratio*rows)
imr = shrink2(image2,irows,irows);
imr = surround(imr, imBackLum2, rows, cols);
con2 = localconmask(imr,rowpixperdeg, colpixperdeg, ...
imBackLum2, pIm, ...
opticalcut, lumspreadcut, maskspreadcut, ...
maskconthresh) ;
[dist, offset] = mindist(con1,con2) ;
dat = [dat ; [dist offset]] ;
end
% dat
[dist mindex] = min(dat(:,1)) ;
dprime = consensmax*dist/(rowpixperdeg*colpixperdeg)^(1/beta) ;
bestsizeratio = 1-2*(mindex-1)/rows ;
ypixoffset = dat(mindex,2)/repfactor ;
xpixoffset = dat(mindex,3)/repfactor ;
dprime = [dprime bestsizeratio ypixoffset xpixoffset] ;
%---- end of symboldiff, the main subroutine --------------
% localconmask, the visible contrast image subroutine
function con = localconmask(imLum, ...
rowpixperdeg, colpixperdeg, ...
imBackLum, pIm, ...
opticalcut, lumspreadcut, maskspreadcut, ...
maskconthresh)
%list of local variables:
% rows, cols, rowfrequcut, colfreqcut, maskfactor
% blurfilter, im0f, im1f, im0, im1
rows = size(imLum,1) ; cols = size(imLum,2);
% Blur with 'optical' blur to get 'visible' luminance image
% 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 * opticalcut / rowpixperdeg ;
colfreqcut = cols * opticalcut / colpixperdeg ;
% The optical blur filter is
blurfilter = filtgaus(rows, cols, rowfreqcut, colfreqcut) ;
imf = fft2(imLum) ;
if prod(size(imBackLum))> 1
imBackf = fft2(imBackLum) ;
end
im = real(ifft2(imf.*blurfilter)) ;
% Blur with the luminance spread function to get local luminance image
rowfreqcut = rows * lumspreadcut / rowpixperdeg; % cycles per image
colfreqcut = cols * lumspreadcut / colpixperdeg;
% The local luminance blur filter is
blurfilter = filtgaus(rows, cols, rowfreqcut, colfreqcut) ;
% Local luminance image
pBack = 1 - pIm ;
if prod(size(imBackLum))> 1
imf = real(ifft2((pBack*imBackf+pIm*im0).*blurfilter)) ;
else
imf = pBack*imBackLum+pIm*real(ifft2(imf.*blurfilter)) ;
end
% Contrast image
im = im./imf - 1 ;
% Masking contrast energy image
imf = real(ifft2(fft2(im.*im).*blurfilter)) ;
% Masked contrast images
maskfactor = 1/maskconthresh^2;
con = im./sqrt(1 + imf*maskfactor) ;
%-----end of localconmask--------------
function filter = filtgaus(r, c, rf, cf)
% 2D Gaussian filter (outer product of two 1D Gaussians)
filter = fltgaus1(r,rf)'*fltgaus1(c,cf);
%---------end filtgaus-------------------------
function filter = fltgaus1( n, f)
% fltgaus1 generates a 1D Gaussian low pass filter
% n: spatial image length (may be odd or even)
% f: 1/e frequency cutoff in cycles per image
n2 = round(n/2);
if n2 == n/2
filter = [[0:n2] [1:n2-1]-n2];
else
filter = [[0:n2] [0:n2-1]-n2];
end
filter = exp(-(filter.*filter)/(f*f)) ;
% For a sequence of length N (N even) Matlab orders
% the Discrete Fourier Transform coefficients
% by their frequencies in cycles per image
% [0 (DC) 1 2 3 ... N/2 -N/2+1 ... -2 -1]
%---------end fltgaus1-------------------------
%-----begin surround ---------------
function as = surround(a, s, nys, nxs)
% Outputs an nys by nxs array, surrounding array a with the value s.
% FUNCTION CALL (Example)
% surround(repmat(1,[2,2]), 0, 4, 4)
ny = size(a,1);
nx = size(a,2);
as = repmat(s,[nys,nxs]) ;
ny0 = floor((nys-ny)/2) ;
nx0 = floor((nxs-nx)/2) ;
as(ny0+1:ny0+ny,nx0+1:nx0+nx) = a ;
%---------------end surround ------------------
%-----begin replicate ---------------
function ar = replicate(a,nyf,nxf)
% Replicates the values of the array by integer factors nyf and nxf
% EXAMPLE FUNCTION CALL
% replicate([1 2; 2 3],2,2)
ny = size(a,1);
nx = size(a,2);
nyr = ny*nyf;
nxr = nx*nxf;
indexy = 1+floor(([1:nyr]-1)/nyf) ;
indexx = 1+floor(([1:nxr]-1)/nxf) ;
ar=a(indexy,indexx) ;
%---------------end replicate ------------------
%--------------begin shrink2 --------------
function imr = shrink2(im,nyr,nxr)
% frequency domain image shrink to given dimensions
% after Watson(1986)
% shrinks from even to even sizes only (e.g., to nyr-1 if nyr is odd).
% EXAMPLE CALL: shrink2(ones(8,8),6,6)
ny = size(im,1);
nx = size(im,2);
nyr2 = floor(nyr/2);
nxr2 = floor(nxr/2);
nyr21 = nyr2+1;
nxr21 = nxr2+1;
nyr22 = ny-nyr2+2;
nxr22 = nx-nxr2+2;
ny1 = ny-nyr2+1;
nx1 = nx-nxr2+1;
imf = fft2(im);
imrf = imf([1:nyr21 nyr22:ny],[1:nxr21 nxr22:nx]);
imrf(nyr21,:)=0.5*(imf(nyr21,[1:nxr21 nxr22:nx]) ...
+ imf( ny1,[1:nxr21 nxr22:nx]));
imrf(:,nxr21)=0.5*(imf([1:nyr21 nyr22:ny],nxr21) ...
+ imf([1:nyr21 nyr22:ny],nx1));
imrf(nyr21,nxr21)=0.25*(imf(nyr21,nxr21) ...
+imf(ny1,nxr21) +imf(nyr21,nx1)+imf(ny1,nx1));
imr = ((4*nyr2*nxr2)/(ny*nx))*real(ifft2(imrf));
%-------------end shrink2 -----------------
%-------------begin mindist -----------------
function [dmin, offset]= mindist(a,b)
% minimum distance from a to b under translations
% example call: mindist([1 0; 0 0],[0 1.1; 0 0])
afft = fft2(a);
bfft = fft2(b);
c = real(ifft2(conj(afft).*bfft)) ;
[cmax, indexyx] = max(c(:)) ;
ssa = sum(a(:).*a(:)) ;
ssb = sum(b(:).*b(:)) ;
dmin = sqrt(ssa + ssb - (cmax + cmax));
rows = size(a,1);
ind1 = indexyx -1;
offset = [ mod(ind1,rows) fix(ind1/rows) ] ;
%-------------end mindist -----------------
%