Macaque-ganglion-cells/retinatools/library.py

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]