diff --git a/QuickCSF/QuickCSF.py b/QuickCSF/QuickCSF.py index 474e424..445e746 100644 --- a/QuickCSF/QuickCSF.py +++ b/QuickCSF/QuickCSF.py @@ -4,10 +4,28 @@ See Lesmes, L. A., Lu, Z. L., Baek, J., & Albright, T. D. (2010). Bayesian adaptive estimation of the contrast sensitivity function: The quick CSF method. Journal of vision, 10(3), 17-17. Assumes a 4-parameter model of CSF: - Peak sensitivity: highest contrast level across all spatial frequencies - Peak frequency: the spatial frequency at which peak sensitivity occurs - Bandwidth: full-width at half-maximum - Delta: difference between peak sensitivity and truncation at low-frequencies + a = Peak sensitivity: highest contrast level across all spatial frequencies + b = Peak frequency: the spatial frequency at which peak sensitivity occurs + c = Bandwidth: full-width at half-maximum + d = Delta: difference between peak sensitivity and truncation at low-frequencies + + log10(CSF(f)) = log10(a) - 4 * log10(2) * [log10(f/b) / log10(2 * c)]^2 if f >= b else log10(a) - d + + + The peripheral version is based on + + Robert Rosén; Linda Lundström; Abinaya Priya Venkataraman; Simon Winter; Peter Unsbo (2014). + Quick contrast sensitivity measurements in the periphery. + Journal of Vision July 2014, Vol.14, 3. doi:https://doi.org/10.1167/14.8.3 + + The peripheral version assumes a different 4-parameter model of CSF: + a = Peak sensitivity: highest contrast level across all spatial frequencies. start=1, stop=256, num=32 (log) + b = Peak frequency: the spatial frequency at which peak sensitivity occurs. start=2, stop=12, num=8 (log) + c = Bandwidth: full-width at half-maximum start=1, stop=9, num=27 (log) + d = The high frequency truncation. start=2, stop=50, num=41 (log) + + log10(CSF(f)) = log10(a) - 4 * log10(2) * [log10(f/b) / log10(2 * c)]^2 if f <= d else 0 + ''' import logging @@ -62,101 +80,12 @@ def makeFrequencySpace(min=.2, max=36, count=20): return frequencySpace -def csf_unmapped(parameters, frequency): - '''The truncated log-parabola model for human contrast sensitivity - - Expects UNMAPPED parameters - Param order = peak sensitivity, peak frequency, bandwidth, log delta - ''' - # Get everything into log-units - [peakSensitivity, peakFrequency, logBandwidth, delta] = mapCSFParams(parameters) - - return csf(peakSensitivity, peakFrequency, logBandwidth, delta, frequency) - -def csf(peakSensitivity, peakFrequency, logBandwidth, delta, frequency): - frequency = numpy.log10(frequency) - - if not isinstance(peakSensitivity, Iterable): - delta = numpy.log10(peakSensitivity) - numpy.log10(peakSensitivity-delta) - delta = numpy.array([delta]) - - peakSensitivity = numpy.array([numpy.log10(peakSensitivity)]) - peakFrequency = numpy.array([numpy.log10(peakFrequency)]) - logBandwidth = numpy.array([numpy.log10(logBandwidth)]) - - frequency = numpy.array([[frequency]]) - - n = len(peakSensitivity) - m = len(frequency[0]) - - frequency = frequency.repeat(n, 0) - - peakFrequency = peakFrequency[:,numpy.newaxis].repeat(m,1) - peakSensitivity = peakSensitivity[:,numpy.newaxis].repeat(m,1) - delta = delta[:,numpy.newaxis].repeat(m,1) - - divisor = numpy.log10(2)+logBandwidth - divisor = divisor[:,numpy.newaxis].repeat(m,1) - truncation = (4 * numpy.log10(2) * numpy.power(numpy.divide(frequency-peakFrequency, divisor), 2)) - - logSensitivity = numpy.maximum(0, peakSensitivity - truncation) - Scutoff = numpy.maximum(logSensitivity, peakSensitivity-delta) - logSensitivity[frequency 0: - gonePositive = True - - if gonePositive and height <= 0: - done = True - else: - area += height * bucketWidth - - return area - -def mapCSFParams(params, exponify=False): - ''' - Maps parameter indices to log values - - Exponify will de-log them, leaving the following units: - Peak Sensitivity: 1/contrast - Peak Frequency: cycles per degree - Bandwidth: octaves - Delta: 1/contrast (Difference between Peak Sensitivity and the truncation) - ''' - peakSensitivity = 0.1*params[:,0] + 0.3 - peakFrequency = -0.7 + 0.1*params[:,1] - bandwidth = 0.05 * params[:,2] - logDelta = -1.7 + 0.1 * params[:,3] - delta = numpy.power(10, logDelta) - - if exponify: - deltaDiff = numpy.power(10, peakSensitivity-delta) - - peakSensitivity = numpy.power(10, peakSensitivity) - peakFrequency = numpy.power(10, peakFrequency) - bandwidth = numpy.power(10, bandwidth) - delta = peakSensitivity - deltaDiff - - return numpy.stack((peakSensitivity, peakFrequency, bandwidth, delta)) - def entropy(p): return numpy.multiply(-p, numpy.log(p)) - numpy.multiply(1-p, numpy.log(1-p)) + class QuickCSFEstimator(): - def __init__(self, stimulusSpace=None): + def __init__(self, stimulusSpace=None, periphery=False): '''Create a new QuickCSF estimator with the specified input/output spaces Args: @@ -169,12 +98,20 @@ def __init__(self, stimulusSpace=None): makeFrequencySpace() ] - parameterSpace = [ - numpy.arange(0, 28), # Peak sensitivity - numpy.arange(0, 21), # Peak frequency - numpy.arange(0, 21), # Log bandwidth - numpy.arange(0, 21) # Low frequency truncation (log delta) - ] + if not periphery: + parameterSpace = [ + numpy.arange(0, 28), # Peak sensitivity + numpy.arange(0, 21), # Peak frequency + numpy.arange(0, 21), # Log bandwidth + numpy.arange(0, 21) # Low frequency truncation (log delta) + ] + else: + parameterSpace = [ + numpy.arange(0, 32), # Peak sensitivity (32 values from 2..256) + numpy.arange(0, 28) , # Peak frequency + numpy.arange(0, 27), # Log bandwidth + numpy.arange(0, 41) # High frequency truncation limit + ] logger.info('Initializing QuickCSFEStimator') logger.debug('Initializing QuickCSFEstimator stimSpace='+str(stimulusSpace).replace('\n','')+', paramSpace='+str(parameterSpace).replace('\n','')) @@ -182,6 +119,8 @@ def __init__(self, stimulusSpace=None): self.stimulusSpace = stimulusSpace self.parameterSpace = parameterSpace + self.periphery = periphery + self.stimulusRanges = [len(sSpace) for sSpace in self.stimulusSpace] self.stimComboCount = numpy.prod(self.stimulusRanges) @@ -191,13 +130,137 @@ def __init__(self, stimulusSpace=None): self.d = 0.5 self.sig = 0.25 - # Probabilities (initialize all of them to equal values that sum to 1) + # Prior probabilities (initialize all of them to equal values that sum to 1) self.probabilities = numpy.ones((self.paramComboCount,1))/self.paramComboCount self.currentStimulusIndex = None self.currentStimParamIndices = None self.responseHistory = [] + def mapCSFParams(self, params, exponify=False): + ''' + Maps parameter indices to log values + + Exponify will de-log them, leaving the following units: + Peak Sensitivity: 1/contrast + Peak Frequency: cycles per degree + Bandwidth: octaves + Delta: if not periphery 1/contrast (Difference between Peak Sensitivity and the truncation) + if periphery cycles/degree where truncation occurs + ''' + if not self.periphery: + peakSensitivity = 0.1*params[:,0] + 0.3 + peakFrequency = -0.7 + 0.1*params[:,1] + bandwidth = 0.05 * params[:,2] + logDelta = -1.7 + 0.1 * params[:,3] + delta = numpy.power(10, logDelta) + else: + peakSensitivity = 0.30103 + 0.06797452 * params[:,0] # log10(2) + (log10(256) - log10(2)) / (32 -1) * i + peakFrequency = 0.0 + 0.04 * params[:,1] # log10(1) + (log10(12) - log10(1)) / (28 - 1) * i + bandwidth = 0.0 + 0.03670163 * params[:,2] # log10(1) + (log10(9) - log10(1)) / (27 - 1) * i + delta = 0.30103 + 0.0349485 * params[:,3] # log10(2) + (log10(50) - log10(2)) / (41 - 1) * i + + if exponify: + if not self.periphery: + deltaDiff = numpy.power(10, peakSensitivity-delta) + + peakSensitivity = numpy.power(10, peakSensitivity) + peakFrequency = numpy.power(10, peakFrequency) + bandwidth = numpy.power(10, bandwidth) + + if not self.periphery: + delta = peakSensitivity - deltaDiff + else: + delta = numpy.power(10, delta) # I am not sure about this one + + return numpy.stack((peakSensitivity, peakFrequency, bandwidth, delta)) + + def csf_unmapped(self, parameters, frequency): + '''The truncated log-parabola model for human contrast sensitivity + + Expects UNMAPPED parameters + Param order = peak sensitivity, peak frequency, bandwidth, log delta + ''' + # Get everything into log-units + [peakSensitivity, peakFrequency, logBandwidth, delta] = self.mapCSFParams(parameters) + + return self.csf(peakSensitivity, peakFrequency, logBandwidth, delta, frequency) + + def csf(self, peakSensitivity, peakFrequency, logBandwidth, delta, frequency): + '''Return values of log10(csf(frequency)). + + If periphery is true, then the peripheral version of the CSF function is used. + + Parameters + ---------- + peakSensitivity : float + peakFrequency : float + logBandwidth : float + delta : float + If periphery is true, this is the left-truncation point of the CSF function (d) in cycles/degree. + If periphery is false, this is the truncation level at low spatial frequencies log(sensitivity). + + Returns + ------- + logSensitivity: [float] + ''' + frequency = numpy.log10(frequency) + + if not isinstance(peakSensitivity, Iterable): + if not self.periphery: + delta = numpy.log10(peakSensitivity) - numpy.log10(peakSensitivity-delta) + delta = numpy.array([delta]) + + peakSensitivity = numpy.array([numpy.log10(peakSensitivity)]) + peakFrequency = numpy.array([numpy.log10(peakFrequency)]) + logBandwidth = numpy.array([numpy.log10(logBandwidth)]) + + frequency = numpy.array([[frequency]]) + + n = len(peakSensitivity) + m = len(frequency[0]) + + frequency = frequency.repeat(n, 0) + + peakFrequency = peakFrequency[:,numpy.newaxis].repeat(m,1) + peakSensitivity = peakSensitivity[:,numpy.newaxis].repeat(m,1) + delta = delta[:,numpy.newaxis].repeat(m,1) + + divisor = numpy.log10(2)+logBandwidth + divisor = divisor[:,numpy.newaxis].repeat(m,1) + truncation = (4 * numpy.log10(2) * numpy.power(numpy.divide(frequency-peakFrequency, divisor), 2)) + + logSensitivity = numpy.maximum(0, peakSensitivity - truncation) + + if not self.periphery: + Scutoff = numpy.maximum(logSensitivity, peakSensitivity-delta) + logSensitivity[frequency delta] = 0 + + return logSensitivity + + def aulcsf(self, peakSensitivity, peakFrequency, logBandwidth, delta, bucketWidth=.1): + def myCSF(frequency): + return self.csf(peakSensitivity, peakFrequency, logBandwidth, delta, frequency)[0][0] + + gonePositive = False + done = False + frequency = 0 + area = 0 + while not done: + frequency += bucketWidth + height = myCSF(frequency) + if height > 0: + gonePositive = True + + if gonePositive and height <= 0: + done = True + else: + area += height * bucketWidth + + return area + def next(self): '''Determine the next stimulus to be tested''' @@ -272,7 +335,7 @@ def _pmeas(self, parameterIndex, stimulusIndex=None): stimulusIndices = self.inflateStimulusIndex(stimulusIndex) frequencies = self.stimulusSpace[1][stimulusIndices[:,1]].reshape(1,-1) - csfValues = csf_unmapped(parameters, frequencies) + csfValues = self.csf_unmapped(parameters, frequencies) # Make vector of sensitivities contrast = self.stimulusSpace[0][stimulusIndices[:,0]] @@ -353,7 +416,7 @@ def getResults(self, leaveAsIndices=False): results = estimatedParamMeans.reshape(1,len(self.parameterRanges)) if not leaveAsIndices: - results = mapCSFParams(results, True).T + results = self.mapCSFParams(results, True).T results = results.reshape(4).tolist() @@ -362,5 +425,5 @@ def getResults(self, leaveAsIndices=False): 'peakFrequency': results[1], 'bandwidth': results[2], 'delta': results[3], - 'aulcsf': aulcsf(*results) + 'aulcsf': self.aulcsf(*results) } diff --git a/QuickCSF/plot.py b/QuickCSF/plot.py index 2c033cc..f650ae6 100644 --- a/QuickCSF/plot.py +++ b/QuickCSF/plot.py @@ -21,7 +21,7 @@ def plot(qCSFEstimator, graph=None, unmappedTrueParams=None, showNumbers=True, s plt.show() if unmappedTrueParams is not None: - truthData = QuickCSF.csf_unmapped(unmappedTrueParams.reshape(1, -1), frequencyDomain) + truthData = qCSFEstimator.csf_unmapped(unmappedTrueParams.reshape(1, -1), frequencyDomain) truthData = numpy.power(10, truthData) truthLine = graph.fill_between( frequencyDomain.reshape(-1), @@ -38,7 +38,7 @@ def plot(qCSFEstimator, graph=None, unmappedTrueParams=None, showNumbers=True, s estimatedParamMeans['bandwidth'], estimatedParamMeans['delta'], ]]) - estimatedData = QuickCSF.csf_unmapped(estimatedParamMeans.reshape(1, -1), frequencyDomain) + estimatedData = qCSFEstimator.csf_unmapped(estimatedParamMeans.reshape(1, -1), frequencyDomain) estimatedData = numpy.power(10, estimatedData) estimatedLine = graph.fill_between( @@ -74,13 +74,13 @@ def plot(qCSFEstimator, graph=None, unmappedTrueParams=None, showNumbers=True, s graph.grid() if showNumbers: - estimatedParamMeans = QuickCSF.mapCSFParams(estimatedParamMeans, exponify=True) + estimatedParamMeans = qCSFEstimator.mapCSFParams(estimatedParamMeans, exponify=True) estimatedParamMeans = estimatedParamMeans.reshape(1,-1).tolist()[0] paramEstimates = '%03.2f, %.4f, %.4f, %.4f' % tuple(estimatedParamMeans) estimatedLine.set_label(f'Estim: {paramEstimates}') if truthData is not None: - trueParams = QuickCSF.mapCSFParams(unmappedTrueParams, True).T.tolist()[0] + trueParams = qCSFEstimator.mapCSFParams(unmappedTrueParams, True).T.tolist()[0] trueParams = '%03.2f, %.4f, %.4f, %.4f' % tuple(trueParams) truthLine.set_label(f'Truth: {trueParams}') diff --git a/QuickCSF/simulate.py b/QuickCSF/simulate.py index 01197e5..13bbd3d 100644 --- a/QuickCSF/simulate.py +++ b/QuickCSF/simulate.py @@ -10,8 +10,6 @@ import matplotlib.pyplot as plt import pathlib -import argparseqt.gui -import argparseqt.groupingTools from . import QuickCSF from .plot import plot @@ -30,6 +28,7 @@ def runSimulation( 'truePeakSensitivity':18, 'truePeakFrequency':11, 'trueBandwidth':12, 'trueDelta':11, }, + periphery=False ): logger.info('Starting simulation') @@ -49,9 +48,9 @@ def runSimulation( parameters['trueBandwidth'], parameters['trueDelta'], ]]) - qcsf = QuickCSF.QuickCSFEstimator(stimulusSpace) + qcsf = QuickCSF.QuickCSFEstimator(stimulusSpace, periphery=periphery) - graph = plot(qcsf, unmappedTrueParams=unmappedTrueParams) + graph = plot(qcsf, unmappedTrueParams=unmappedTrueParams, show = imagePath is None) # Trial loop for i in range(trials): @@ -63,7 +62,7 @@ def runSimulation( if usePerfectResponses: logger.debug('Simulating perfect response') frequency = newStimValues[:,1] - trueSens = numpy.power(10, QuickCSF.csf_unmapped(unmappedTrueParams, numpy.array([frequency]))) + trueSens = numpy.power(10, qcsf.csf_unmapped(unmappedTrueParams, numpy.array([frequency]))) testContrast = newStimValues[:,0] testSens = 1 / testContrast @@ -78,7 +77,7 @@ def runSimulation( # Update the plot graph.clear() graph.set_title(f'Estimated Contrast Sensitivity Function ({i+1})') - plot(qcsf, graph, unmappedTrueParams) + plot(qcsf, graph, unmappedTrueParams, show = imagePath is None) if imagePath is not None: plt.savefig(pathlib.Path(imagePath+'/%f.png' % time.time()).resolve()) @@ -93,13 +92,14 @@ def runSimulation( paramEstimates = qcsf.getResults() logger.info('Results: ' + str(paramEstimates)) - trueParams = QuickCSF.mapCSFParams(unmappedTrueParams, True).T + trueParams = qcsf.mapCSFParams(unmappedTrueParams, True).T print('******* Results *******') print(f'\tEstimates = {paramEstimates}') print(f'\tActuals = {trueParams}') print('***********************') plt.ioff() - plt.show() + if imagePath is None: + plt.show() def entropyPlot(qcsf): params = numpy.arange(qcsf.paramComboCount).reshape(-1, 1) @@ -128,8 +128,10 @@ def entropyPlot(qcsf): parser = argparse.ArgumentParser() + parser.add_argument('--noQT', default=False, action='store_true', help='Do not use QT ui for input') + parser.add_argument('-n', '--trials', type=int, help='Number of trials to simulate') - parser.add_argument('--imagePath', default=None, help='If specified, path to save images') + parser.add_argument('--imagePath', default=None, help='If specified, path to save images (does not plot to screen).') parser.add_argument('-perfect', '--usePerfectResponses', default=False, action='store_true', help='Whether to simulate perfect responses, rather than probablistic ones') stimuliSettings = parser.add_argument_group('stimuli') @@ -147,12 +149,36 @@ def entropyPlot(qcsf): parameterSettings.add_argument('-b', '--trueBandwidth', type=int, default=12, help='True bandwidth (index)') parameterSettings.add_argument('-d', '--trueDelta', type=int, default=11, help='True delta truncation (index)') - settings = argparseqt.groupingTools.parseIntoGroups(parser) - if settings['trials'] is None: - from qtpy import QtWidgets - from . import ui - app = QtWidgets.QApplication() - settings = ui.getSettings(parser, settings, ['trials']) + parser.add_argument('--periphery', default=False, action='store_true', help='If present, use peripheral quickCSF. Rosen suggests -minf 1 -maxf 50 -fr 50 -cr 64.') + + args = parser.parse_args() + if args.noQT: + settings = {} + if args.trials: + settings['trials'] = args.trials + if args.imagePath: + settings['imagePath'] = args.imagePath + if args.usePerfectResponses: + settings['usePerfectResponses'] = args.usePerfectResponses + settings['stimuli'] = { + 'minContrast': args.minContrast, 'maxContrast': args.maxContrast, 'contrastResolution': args.contrastResolution, + 'minFrequency': args.minFrequency, 'maxFrequency': args.maxFrequency, 'frequencyResolution': args.frequencyResolution, + } + settings['parameters'] = { + 'truePeakSensitivity': args.truePeakSensitivity, 'truePeakFrequency': args.truePeakFrequency, + 'trueBandwidth': args.trueBandwidth, 'trueDelta': args.trueDelta, + } + if args.periphery: + settings['periphery'] = args.periphery + else: + import argparseqt.gui + import argparseqt.groupingTools + settings = argparseqt.groupingTools.parseIntoGroups(parser) + if settings['trials'] is None: + from qtpy import QtWidgets + from . import ui + app = QtWidgets.QApplication() + settings = ui.getSettings(parser, settings, ['trials']) if settings is not None: runSimulation(**settings) diff --git a/README.md b/README.md index 39d81c6..c11c7db 100644 --- a/README.md +++ b/README.md @@ -2,11 +2,19 @@ A fast, adaptive approach to estimating contrast sensitivity function parameters. -This implmentation is based on: +This implementation is based on: > Lesmes, L. A., Lu, Z. L., Baek, J., & Albright, T. D. (2010). Bayesian adaptive estimation of the contrast sensitivity function: The quick CSF method. *Journal of vision*, 10(3), 17-17. -Special thanks to Dr. Tianshi Lu at Wichita State University for providing a Matlab implemenation of the fundamental algorithm and Dr. Rui Ni at Wichita State University for the motivation. +Special thanks to Dr. Tianshi Lu at Wichita State University for providing a Matlab implementation of the fundamental algorithm and Dr. Rui Ni at Wichita State University for the motivation. + +Peripheral CSF based on + +> Rosén, R., Lundström, L., Venkataraman A.P., Winter, S., Unsbo, P. (2014). + Quick contrast sensitivity measurements in the periphery. + Journal of Vision July 2014, Vol.14, 3. doi:https://doi.org/10.1167/14.8.3 + +Thanks to Prof. Andrew Turpin at Lions Eye Institute, Australia, for incorporating the peripheral parameterisation into our code and providing the noQT option.. ## Dependencies ~~~~ @@ -14,10 +22,8 @@ $ pip3 install -e . ~~~~ Requires: * `numpy` -* `qtpy` +* `qtpy` (not necessary for simulation) * Qt bindings (via `PySide2`, `PyQt5`, `PySide`, or `PyQt`) - -Optional (for simulation visuals): * `matplotlib` ## Usage @@ -30,6 +36,7 @@ A settings dialog will appear; session ID and viewing distance are required. Arg ~~~bash $ python -m QuickCSF.app --help ~~~ + ### Simulate and visualize an evaluation Run: ~~~bash @@ -38,4 +45,9 @@ $ python -m QuickCSF.simulate A settings dialog will appear; the number of trials is required. Arguments can also be specified on the command line. use the `--help` flag to see all options: ~~~bash $ python -m QuickCSF.simulate --help +~~~ + +As an example for a command line version that does not use `qtpy` +~~~bash +$ python -m QuickCSF.simulate --noQT --periphery -minf 1 -maxf 50 -fr 50 -cr 64 -s 27 -d 20 -b 10 -perfect -n 200 ~~~ \ No newline at end of file