#
# visible contrast energy with masking by a homogeneous background
# program begins with subroutines (last one is dBVM)
# fltgaus1 generates a 1D Gaussian low pass filter
# n: spatial image length (must be even)
# f: 1/e frequency cutoff in cycles per image
fltgaus1<-function( n, f){
n2 = round(n/2);
if (n2 == n/2) {
freq = c(0:n2,(1:(n2-1))-n2);
} else {
freq = c(0:n2,(0:(n2-1))-n2);
}
exp(-(freq*freq)/(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------------------------- #
#
# 2D Gaussian filter (outer product of two 1D Gaussians) function
filtgaus<-function(r, c, rf, cf){
fltgaus1(r,rf)%o%fltgaus1(c,cf) }
#---------end filtgaus-------------------------
# filtdog makes a Difference of Gaussians filter
# r,c: # of rows and columns
# rf,cf: row and column high frequency cutoffs
# sr (spread ratio): wide spatial spread / narrow one
# ar (amplitude ratio): surround DC level / center DC level
filtdog <- function(r, c, rf, cf, ar, sr){
filtgaus(r,c,rf,cf)-ar*filtgaus(r,c,rf/sr,cf/sr) }
#---------end filtdog-------------------------
# dBB, contrast energy in dB re 10-6 deg^2 sec
dBB <- function(conImage,pixelArea=1/120^2,duration=0.25){
60 + 10*log10(pixelArea*duration*sum(conImage^2)) }
#-----end of dBB--------------
# dBV, the visible contrast energy subroutine
dBV <- function(conImage,csf,pixelArea=1/120^2,duration=0.25){
power = sum(abs(fft(conImage)*csf)^2) ;
60 + 10*log10(pixelArea*duration*power/prod(dim(conImage))) }
#-----end of dBV--------------
# dBMask, the masking factor subroutine function
dBMask <- function(conImage,maskConThresh= 0.05){
10*log10(1+mean(conImage^2)/(maskConThresh*maskConThresh)) }
#-----end of dBVMask--------------
# dBVMask, the visible masking factor subroutine function
dBVMask <- function(conImage,csf,maskConThresh= 0.05){
mscon = mean(abs(fft(conImage)*csf)^2) ;
10*log10(1+mscon/(maskConThresh*maskConThresh*prod(dim(conImage)))) }
#-----end of dBVMask--------------
# dBVM, the masked visible contrast energy subroutine
dBVM <- function(conImage1,conImage2,csf,pixelArea=1/(120*120),duration=0.25,maskConThresh=0.05){
dBV(conImage1-conImage2,csf,pixelArea,duration)-dBVMask(conImage1,csf,maskConThresh)}
#-----end of dBVM--------------
# Observer contrast sensitivity
# Here we set the observer's best possible contrast energy sensitivity to 0 dBB.
# This is unusually good, but not exceeding rare, sensitivity (Watson, 2000 Optics Express).
# For a duration of 0.25 sec and an area of 1 deg^2, it corresponds to a contrast sensitivity
# of 500 (a contrast threshold of 0.002).
# Relative to this standard, the ModelFest observers had approximately
# a mean sensitivity of -6 dB with a standard deviation of 3 dB.
dBObs = 0 ;
# images are assumed to be space-time separable with an equivalent energy duration of
duration = 0.25 ;
# Outline
# 1) Example input images
# 2) Example parameters # 3) Example computation of visible masked contrast energy
# 1) Example images # Example: grating increment detection
# Background or pedestal (image1) with a contrast of 0.1,
# backgound + increment (image2) with a contrast of 0.2,
# 16 pixels in each dimension and a # grating frequency of 3 cycles per image.
# Background level (irrelevant) set to 128
# image size
rows = 16; cols = 16;
image1 = array(128*(1.0 + 0.1*sin(2*pi*3*(1:cols)/cols)),c(rows,cols)) ;
image2 = array(128*(1.0 + 0.2*sin(2*pi*3*(1:cols)/cols)),c(rows,cols)) ;
# 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 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 = 8 ;
# iii) the ratio of surround DC amplitude to center DC amplitude
ampratio = 0.685 ;
# 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) ;
# The unity peak gain contrast sensitivity filter is
csf = csf/max(max(csf)) ;
# c) contrast masking threshold
maskconthresh=0.05 ;
# 5# rms contrast raises threshold by 3 dB
# 3) Example computation of visible masked contrast energy
# convert images to contrast images
meanlum = mean(image1) ;
con1 = image1/meanlum - 1 ;
con2 = image2/meanlum - 1 ;
# contrast energy in difference signal is 25.5 dB
dBB12 = dBB(con1-con2,pixelArea,duration)
# visible contrast energy in difference signal is 25.2 dB
dBV12 = dBObs + dBV(con1-con2,csf,pixelArea,duration)
# contrast energy in difference signal masked by background is 20.6 dB
dBVM12 = dBObs + dBVM(con1,con2,csf,pixelArea,duration,maskconthresh)
# Example dprime value = 10.7563
dprime = 10^(dBVM12/20) ;
# Example conversion from d' to threshold increment.
# image1 had a contrast of 0.1, image2 had a contrast of 0.2,
# so the model predicts the contrast difference threshold
# for a base contrast of 0.1 is the actual difference/dprime = 0.0093 ;
diffthresh = (0.2 - 0.1)/ dprime ;
# show results
c(dBVM12,dBV12,dBB12,dprime,diffthresh)