import numpy as np import cv2 from scipy import interpolate def SimpleLowpass(signal, dt, RC): Nsamples = signal.size output = np.zeros(Nsamples) s = dt / (RC + dt) output[0] = s*signal[0] for i in range(1,Nsamples): output[i] = output[i-1] + s * (signal[i] - output[i-1]) return output def SimpleHighpass(signal, dt, RC): Nsamples = signal.size output = np.zeros(Nsamples) s = dt / (RC + dt) output[0] = s*signal[0] for i in range(1,Nsamples): output[i] = output[i-1] + s * (signal[i] - output[i-1]) return signal - output def LuminanceGainControl(signal, dt, RC): Nsamples = signal.size output = np.zeros(Nsamples) # Initialization of output helper = np.zeros(Nsamples) # Internal helper variable, i.e. LP filtered output s = dt / (RC + dt) output[0] = 0 helper[0] = 1 for i in range(1,Nsamples): helper[i] = helper[i-1] + s * (output[i-1] - helper[i-1]) # this LP filters the last output to get new helper output[i] = signal[i] / helper[i] return output def rec(x): return x*(x>0) def nlplus(x, c1): return (1 + rec(x)/c1)**2 def nlc(x, q0, q1, c2): return q0 + q1*x/(c2 + x) def nltan(x, k): return (2/np.pi)*np.arctan(k*x) def ContrastGainControlLoop(signal, dt, RCplus, RCA, c1): Nsamples = signal.size output = np.zeros(Nsamples) cA = np.zeros(Nsamples) d = np.ones(Nsamples) output[0] = 0 d[0] = 1 for i in range(1,Nsamples): # First, calculate the HP filtered input with TC RC s = dt / (RCA + dt) cA[i] = cA[i-1] + s * (signal[i] - cA[i-1]) # In Contrast gain, divisive normalization for output s = dt / (RCplus + dt) d[i] = d[i-1] + s * (output[i-1] - d[i-1]) output[i] = (signal[i] - cA[i]) / nlplus(d[i], c1) return output def delay_signal(signal, dt, delay): """ This delays a signal by "delay" Internally, the signal is (1) resampled at rate RESAMPLED: int (default 100) with cubic spline (2) shifted by the time DELAY (3) resampled at original sampling """ RESAMPLED = 100 times = np.arange(0, len(signal)*dt, dt) interpolationfunction = interpolate.interp1d(times, signal, kind='cubic') times_fine = np.arange(0, np.max(times), dt/RESAMPLED) signal_fine = interpolationfunction(times_fine) new_idx = int( delay / (dt/RESAMPLED) ) delayed_signal = np.roll(signal_fine, new_idx) delayed_signal[0:new_idx] = 0 interpolationfunction2 = interpolate.interp1d(times_fine, delayed_signal, kind='cubic', fill_value="extrapolate") new_signal = interpolationfunction2(times) return new_signal def receptiveField(x, y, xoff, yoff, rcenter): center = np.exp( -((x-xoff)**2 + (y-yoff)**2) / (2*rcenter*rcenter)) /(2*np.pi*rcenter*rcenter) return center def get_RF(xoff, yoff, center): X = np.linspace(-128, 128, 256) Y = np.linspace(-128, 128, 256) x, y = np.meshgrid(X,Y) RF = receptiveField(x, y, xoff, yoff, center) return RF def filter_video_PC(filename, maxframe, fitparameters): xoff = fitparameters['xoff'] yoff = fitparameters['yoff'] rcenterM = fitparameters['rcenterM'] rcenterL = fitparameters['rcenterL'] RFM = get_RF(xoff, yoff, rcenterM) RFL = get_RF(xoff, yoff, rcenterL) vidcap = cv2.VideoCapture(filename) counter = 0 success = True ratesM = [] ratesL = [] while success and (counter < 9750): success, image = vidcap.read() scone, mcone, lcone = image_to_cone(image) rateM = np.sum(mcone*RFM) rateL = np.sum(lcone*RFL) ratesM.append( rateM ) ratesL.append( rateL ) counter += 1 return np.array(ratesM), np.array(ratesL) def filter_video_MC(filename, maxframe, fitparameters): xoff = fitparameters['xoff'] yoff = fitparameters['yoff'] rcenter = fitparameters['rcenter'] rsurr = fitparameters['rsurr'] RFcenter = get_RF(xoff, yoff, rcenter) RFsurround = get_RF(xoff, yoff, rsurr) vidcap = cv2.VideoCapture(filename) counter = 0 success = True rates_center = [] rates_surround = [] while success and (counter < maxframe): success, image = vidcap.read() scone, mcone, lcone = image_to_cone(image) luminance_frame = mcone + lcone rate_center = np.sum(luminance_frame*RFcenter) rates_center.append(rate_center) rate_surround = np.sum(luminance_frame*RFsurround) rates_surround.append(rate_surround) counter += 1 return np.array(rates_center), np.array(rates_surround) def filter_video_BO(filename, maxframe, fitparameters): xoff = fitparameters['xoff'] yoff = fitparameters['yoff'] rcenterY = fitparameters['rcenterY'] rcenterB = fitparameters['rcenterB'] RFY = get_RF(xoff, yoff, rcenterY) RFB = get_RF(xoff, yoff, rcenterB) vidcap = cv2.VideoCapture(filename) counter = 0 success = True ratesL = [] ratesM = [] ratesB = [] while success and (counter < maxframe): success, image = vidcap.read() scone, mcone, lcone = image_to_cone(image) rateL = np.sum(lcone*RFY) rateM = np.sum(mcone*RFY) rateB = np.sum(scone*RFB) ratesL.append( rateL ) ratesM.append( rateM ) ratesB.append( rateB ) counter += 1 return np.array(ratesL), np.array(ratesM), np.array(ratesB) def image_to_cone(image): """ Spectral fits and parameters for this display, measured by Barry and Hans Image has to be BGR! """ #Raw RGB values im_rgb = cv2.cvtColor(image, cv2.COLOR_BGR2RGB) red = np.double(im_rgb[:,:,0]) green = np.double(im_rgb[:,:,1]) blue = np.double(im_rgb[:,:,2]) #Hans' gamma correction rgamma = 0.01451 + 0.9855*((1.0*red/256)**2.3122) ggamma = 0.005123 + 0.9949*((1.0*green/256)**2.2752) bgamma = 0.02612 + 0.9739*((1.0*blue/256)**2.2818) lcone = 1.0*(2.74*rgamma + 3.4*ggamma + 1.34*bgamma) mcone = 1.21*(1.06*rgamma + 3.58*ggamma + 2.07*bgamma) scone = 0.212*rgamma + 8.28*ggamma + 285*bgamma return scone, mcone, lcone def get_MC_response(rates_center, rates_surround, fitparameters, type): rC = SimpleLowpass(fitparameters['lum_factor']*rates_center, fitparameters['dt'], fitparameters['tau1']) rC = SimpleLowpass(rC, fitparameters['dt'], fitparameters['tau1']) rC = SimpleLowpass(rC, fitparameters['dt'], fitparameters['tau1']) rC = LuminanceGainControl(rC, fitparameters['dt'], fitparameters['tau2']) rS = SimpleLowpass(fitparameters['lum_factor']*rates_surround, fitparameters['dt'], fitparameters['tau1']) rS = SimpleLowpass(rS, fitparameters['dt'], fitparameters['tau1']) rS = SimpleLowpass(rS, fitparameters['dt'], fitparameters['tau1']) rS = LuminanceGainControl(rS, fitparameters['dt'], fitparameters['tau2']) #Delay surround delayed_surround = delay_signal(rS, fitparameters['dt'], fitparameters['surround_delay']) #Calculate input signal from center - surround lum_signal = rec( rC - fitparameters['ksurr'] * delayed_surround ) if type == "on-cell": b = nltan(lum_signal, fitparameters['k1']) elif type == "off-cell": b = -1*nltan(lum_signal, fitparameters['k1']) else: error(type + " ...not defined.") filtered_signal = ContrastGainControlLoop(b, fitparameters['dt'], fitparameters['tauplus'], fitparameters['tauA'], fitparameters['c1']) Rpc = rec( nltan(filtered_signal, fitparameters['k2']) ) return Rpc def get_PC_response(lum_signalM, lum_signalL, fitparameters, type): signalM = SimpleLowpass(fitparameters['lum_factor']*lum_signalM, fitparameters['dt'], fitparameters['tau1']) signalM = SimpleLowpass(signalM, fitparameters['dt'], fitparameters['tau1']) signalM = SimpleLowpass(signalM, fitparameters['dt'], fitparameters['tau1']) signalM = LuminanceGainControl(signalM, fitparameters['dt'], fitparameters['tau2']) signalL = SimpleLowpass(fitparameters['lum_factor']*lum_signalL, fitparameters['dt'], fitparameters['tau1']) signalL = SimpleLowpass(signalL, fitparameters['dt'], fitparameters['tau1']) signalL = SimpleLowpass(signalL, fitparameters['dt'], fitparameters['tau1']) signalL = LuminanceGainControl(signalL, fitparameters['dt'], fitparameters['tau2']) bM = nltan(signalM, fitparameters['k1']) bL = nltan(signalL, fitparameters['k1']) if type == "Ron": out = rec(bL - fitparameters['alpha']*bM) elif type == "Gon": out = rec(bM - fitparameters['alpha']*bL) else: error(type + " ...not defined.") out = SimpleHighpass(out, fitparameters['dt'], fitparameters['tau3']) Rpc = rec( fitparameters['k2']*out + fitparameters['o1']) return Rpc def get_BO_response(lum_signalB, lum_signalM, lum_signalL, fitparameters): signalM = SimpleLowpass(lum_signalM, fitparameters['dt'], fitparameters['tau1']) signalM = SimpleLowpass(signalM, fitparameters['dt'], fitparameters['tau1']) signalM = SimpleLowpass(signalM, fitparameters['dt'], fitparameters['tau1']) signalM = LuminanceGainControl(signalM, fitparameters['dt'], fitparameters['tau2']) bM = nltan(signalM, fitparameters['kM']) signalL = SimpleLowpass(lum_signalL, fitparameters['dt'], fitparameters['tau1']) signalL = SimpleLowpass(signalL, fitparameters['dt'], fitparameters['tau1']) signalL = SimpleLowpass(signalL, fitparameters['dt'], fitparameters['tau1']) signalL = LuminanceGainControl(signalL, fitparameters['dt'], fitparameters['tau2']) bL = nltan(signalL, fitparameters['kL']) bY = fitparameters['alpha']*bM + fitparameters['beta']*bL bY = delay_signal(bY, fitparameters['dt'], fitparameters['delay']) signalB = SimpleLowpass(lum_signalB, fitparameters['dt'], fitparameters['tau1']) signalB = SimpleLowpass(signalB, fitparameters['dt'], fitparameters['tau1']) signalB = SimpleLowpass(signalB, fitparameters['dt'], fitparameters['tau1']) signalB = LuminanceGainControl(signalB, fitparameters['dt'], fitparameters['tau2']) bB = nltan(signalB, fitparameters['kB']) out = rec(bB - bY) Rpc = SimpleHighpass(out, fitparameters['dt'], fitparameters['tau3']) Rpc = rec(Rpc - fitparameters['offset']) return Rpc def compare_to_file(filename, Rpc, dt, column, shift): cell_activity = np.loadtxt(filename) time_cell = np.arange(0, len(cell_activity))*dt + 5000 - shift*6.66 cell = cell_activity[(time_cell>5500) * (time_cell<64980), column] time = np.arange(0, len(Rpc)*dt, dt) model = Rpc[(time>=5500) * (time<64980)] r = np.corrcoef(cell, model) return r[0,1]