% This will return a conductance-as-function-of-time for treating dendrites
% as single compartment
%
% class: 1 = 10x10 box
% 2 = 10x80 bar
% 3 = looming
% 4 = accelarating box
% 5 = loomlike box
% velocity - in deg/s -- for class = 4, specify start and end velocity in 2-element vector
% velocity - in deg/s -- for class = 5, specify the angle
% direction: 1 = AP, -1 = PA, 2=DV, -2=VD
% dt: time step for *simulation*
% synaptic_map_path: which map to use? uniform is current dflt
% n_syns_per_input: how many synapses *per facet*?
%
% The following variables alter timing to account for feedforward inhibition and the observation
% that larger inputs are more synchronized.
%
% dvn: parameters of delay-versus-nfacets function - delay = v(1)./(1+exp((nsyns-v(2))/v(3)))+v(4)
% sdvn: parameters of sd of delay-versus-nsyns function - delay_jitter = v(1)./(1+exp((nsyns-v(2))/v(3)))+v(4)
% note that the actual delay is jittered by drawing from a normal distro with the above sd;
% you read it correctly -- dvn is based on the number of facets, while sdvn is based on the
% number of synapses
% vel_sf: this accounts for the observation that there is speed modulation speed_sf = v(1)./(1+exp((nsyns-v(2))/v(3)))+v(4)
%
function g_syn_of_t = get_synaptic_pattern(class, direction, velocity, syn_tau, syn_gmax, ...
duration, dt, synaptic_map_path, n_syns_per_input, dvn, sdvn, vel_sf)
% --- Definitions
% synaptic map -- to determine how many synapses get triggered
if (exist('synaptic_map_path') == 0 | synaptic_map_path == -1)
synaptic_map_path = 'uniform_synmap.mat'; % The synaptic mapping [cmpt frac az el]
end
dt_image = 5; % screen inter-step time; monitor specific
half_angle = 2.5; % acceptance half-angle in degrees; add
% for translating stimuli:
edge_sf = [1 0]; % apply 1 to the leading edge, 2 to the trailing edge [1 0], for instance, means leading only
% --- Preliminaries
load(synaptic_map_path);
n_syns = length(synmap);
% --- derive the timing pattern
switch (class)
case 1 % 10x10 box
vert_lims = [-5-half_angle 5+half_angle];
hor_lims = [85-half_angle 95+half_angle];
bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
switch(direction)
case 1 % AP
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [50 130], dt_image, bar_thickness, vert_lims, edge_sf);
case -1 % PA
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [130 50], dt_image, bar_thickness, vert_lims, edge_sf);
case 2 % DV
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [40 -40], dt_image, bar_thickness, hor_lims, edge_sf);
case -2 % VD
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [-40 40], dt_image, bar_thickness, hor_lims, edge_sf);
end;
case 1.1 % 10x10 box, but with bounds that are used for the SFA 'review'/model paper
vert_lims = [-5-half_angle 5+half_angle];
hor_lims = [85-half_angle 95+half_angle];
bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
switch(direction)
case 1 % AP
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [60 120], dt_image, bar_thickness, vert_lims, edge_sf);
case -1 % PA
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [120 60], dt_image, bar_thickness, vert_lims, edge_sf);
case 2 % DV
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [30 -30], dt_image, bar_thickness, hor_lims, edge_sf);
case -2 % VD
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [-30 30], dt_image, bar_thickness, hor_lims, edge_sf);
end;
case 2 % 10x80 bar
vert_lims = [-40-half_angle 40+half_angle];
hor_lims = [50-half_angle 130+half_angle];
bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
switch(direction)
case 1 % AP
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [50 130], dt_image, bar_thickness, vert_lims, edge_sf);
case -1 % PA
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [130 50], dt_image, bar_thickness, vert_lims, edge_sf);
case 2 % DV
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [40 -40], dt_image, bar_thickness, hor_lims, edge_sf);
case -2 % VD
[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [-40 40], dt_image, bar_thickness, hor_lims, edge_sf);
end;
case 3 % loom l/v = velocity
% calculate so that you finish with 80 deg, then go for 100 ms more
end_time = velocity/(tand(80/2));
start_time = duration - (end_time + 100);
theta_start = 2*atand(velocity/start_time);
[synpos, timing, inst_vel] = pattern_circular_looming_stimulus(synmap, velocity, [0 90], theta_start, dt_image);
case 4 % 10x10 box, accelarating linearly
vert_lims = [-5-half_angle 5+half_angle];
hor_lims = [85-half_angle 95+half_angle];
bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
switch(direction)
case 1 % AP
[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'v', velocity, [60 120], 5, bar_thickness, vert_lims, edge_sf);
case -1 % PA
[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'v', velocity, [120 60], 5, bar_thickness, vert_lims, edge_sf);
case 2 % DV
[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'h', velocity, [30 -30], 5, bar_thickness, hor_lims, edge_sf);
case -2 % VD
[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'h', velocity, [-30 30], 5, bar_thickness, hor_lims, edge_sf);
end;
case 5 % 10x10 box, loomlike linearly
vert_lims = [-5-half_angle 5+half_angle];
hor_lims = [85-half_angle 95+half_angle];
bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
switch(direction)
case 1 % AP
[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'v', velocity, [60 120], 5, bar_thickness, vert_lims, edge_sf);
case -1 % PA
[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'v', velocity, [120 60], 5, bar_thickness, vert_lims, edge_sf);
case 2 % DV
[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'h', velocity, [30 -30], 5, bar_thickness, hor_lims, edge_sf);
case -2 % VD
[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'h', velocity, [-30 30], 5, bar_thickness, hor_lims, edge_sf);
end;
end;
% --- the preceding logic returns a timing vector, which tells you the timing of ommatidial
% --- activation. however, response at the level of the LGMD is non-linear; specifically,
% --- synaptic resposne increases superlinearly and occurs faster for larger luminance changes;
% --- this is true for both inhibitory and excitatory inputs. We will do this by making the
% --- assumption that intensity change is approximated by the fraction of synapses activated per timestep
% --- for each element of timing, calculate the scaling factor based on velocity ...
% --- this must be done BEFORE any other timing alterations as it depends on the OMMATIDIAL
% --- transition time, which is not the same as the timing vetor after considering delay
vsf = ones(size(timing));
tv = inst_vel.t_vec;
for t=0:dt:duration
idx = find (timing >= t & timing < t + dt_image);
ivi = min(find(tv >= t-1));
if (length(ivi) == 0)
%disp(['warning: ivi not found; t = ' num2str(t)]) ;
continue;
end
vsf(idx) = vel_sf(1)./(1+exp((inst_vel.vel_vec(ivi)-vel_sf(2))/vel_sf(3)))+vel_sf(4);
%if(length(idx) > 0) ; disp(['vsf: ' num2str(vsf(idx(1))) ' vel: ' num2str(inst_vel.vel_vec(ivi)) ' t: ' num2str(t) ' ivi: ' num2str(ivi)]); end
end
% --- first, introduce mean delay via dvn based on number of facets hit in a given timestep
n_timing = timing;
for t=0:dt_image:duration
idx = find(timing >= t & timing < t + dt_image);
nfacets = length(idx);
n_timing(idx) = timing(idx) + dvn(1)./(1+exp((nfacets-dvn(2))/dvn(3)))+dvn(4);
end
timing = n_timing;
% --- Now account for the fact that you do not necessarily have the same number
% --- of synapses as facests
n_timing = zeros(1,floor(length(synpos)*n_syns_per_input));
cf = 1;
for t=0:dt_image:duration
idx = find(timing >= t & timing < t + dt_image);
nfacets = length(idx);
nsyns = nfacets*n_syns_per_input;
% For a non-integer number of synapses, scale based on a normal pdf
if (nsyns ~= round(nsyns))
rem = nsyns - floor(nsyns);
if (rem > rand)
nsyns = floor(nsyns) + 1;
else
nsyns = floor(nsyns);
end
end
if (nsyns ~= round(nsyns)) ; disp (['wtf:' num2str(nsyns) ' rem: ' num2str(rem)]); end
% Now make nsyns in n_timing have the time
n_timing(cf:min(length(n_timing),cf+nsyns)) = t;
if (length(idx) > 0)
n_vsf(cf:min(length(n_timing),cf+nsyns)) = mean (vsf(idx));
else
n_vsf(cf:min(length(n_timing),cf+nsyns)) = 1;
end
cf = cf + nsyns;
if (cf > length(n_timing))
disp('Warning: vector too big');
cf = length(n_timing);
end
end
timing = n_timing(find(n_timing > 0));
vsf = n_vsf;
% --- it is still unknown if there is an actual increase in input strength from stronger stimuli;
% --- what is certain is that jitter declines. do jitter here. *still based on number of facets*
n_timing = timing;
for t=0:dt_image:duration
idx = find(timing >= t & timing < t + dt_image);
nsyns = length(idx);
jitter = normrnd(0, sdvn(1)./(1+exp(nsyns-sdvn(2))/sdvn(3))+sdvn(4), length(idx), 1);
if (length(idx) > 0) ; n_timing(idx) = timing(idx) + jitter'; end
end
timing = n_timing;
% --- compute the conductance-as-function-of-time
g_syn_of_t = zeros(1,round(duration/dt)+1);
for t=0:dt:duration
active_idx = find (timing >= t-syn_tau*20 & timing <= t);
for a = 1:length(active_idx)
T = t-timing(active_idx(a));
S = vsf(active_idx(a));
g_syn_of_t(1+round(t/dt)) = g_syn_of_t(1+round(t/dt)) + S*syn_gmax*(T/syn_tau)*exp(1-T/syn_tau);
end
end
% --- sanity check
if (0 == 1)
plot(0:dt:duration, g_syn_of_t);
end
%
% This should return timing and synapses for loomlike bars of CONSTANT
% ANGULAR SIZE AND CHANGING ANGULAR VELOCITY. Each position is activated
% TWICE - once for the ON and once for the OFF.
%
% RETURNS:
% synpos - [cmpt_idx cmpt_fraction]
% timing - in ms
%
% ASSIGNED:
% synmap - [cmpt cmpt_frac az el]
% direction - 'h' 'v' for horizontal, vertical (DIRECTION OF MOTION, not bar
% orientation!) - default is vertical
% velocity - in degrees / second -- l/v ; negative --> receding
% limits - [start end] where to start; if horizontal motion, give azimuths;
% if vertical, elevations.
% dt - temporal resolution, in ms
% thick - thickness of the bar, in degrees (i.e., dimension along axis PARALLEL
% to motion vector)
% length_limits - edges of the bar, in degrees for dimension along axis PERPENDICULAR to
% motion vector
%
%
function [synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, direction, velocity, ...
limits, dt, thick, length_limits, ...
edge_sf)
% equivalent disc RADIUS at which to start
theta_start = 5;
% Preliminary
timing = [];
% Orientation-dependent component
if (strcmp('h', direction))
positional_array = synmap(:,4);
length_array = synmap(:,3);
else
positional_array = synmap(:,3);
length_array = synmap(:,4);
end
used = zeros(1,length(positional_array));
% Introduce variation into the eye -- vary each ommatidial view by sd of 5 degrees
eye_deformation = normrnd(0,10,length(positional_array),1);
positional_array = positional_array + eye_deformation;
% Determine valid based on length limits
length_valids = find(length_array > min(length_limits) & length_array < max(length_limits));
% determine total time based on the limits and the accelaration
% loom time
t_o = velocity/tand(theta_start);
t_f = velocity/tand(diff(limits) + theta_start);
if (velocity < 0) ; t_o = -1*t_o; t_f = -1*t_f ; end
total_time = 5*(1+floor(abs(t_f-t_o)/5));
% instantaneous velocity vector
inst_vel.vel_vec = zeros(length(0:dt:total_time),1);
ivi = 1;
% loop thru time, and add accordingly
position = limits(1); % LEADING edge
for time=0:dt:total_time
last_position = position;
if (velocity < 0) % 'receding'
theta_t = atand(velocity/(t_f+time));
theta_tp1 = atand(velocity/(t_f+(time+dt)));
else % normal
theta_t = atand(velocity/(t_o-time));
theta_tp1 = atand(velocity/(t_o-(time+dt)));
end
dp = abs(theta_tp1-theta_t);
% change in position/time is velocity!
int_vel.vel_vec = dp;
if (time+dt > t_o) ; dp = 0; end
stimulated_leading = [];
stimulated_trailing = [];
% Increasing position value (e.g., rightward)
if (limits(1) < limits(2))
position = position + dp;
% Leading edge
if (position < limits(2))
stimulated_leading = find(positional_array >= last_position & positional_array < position);
end
% Trailing edge
if (position-thick > limits(1))
stimulated_trailing = find(positional_array >= last_position-thick & positional_array < position-thick);
end
% For decreasing position value
else
position = position - dp;
% Leading edge
if (position < limits(1))
stimulated_leading = find(positional_array < last_position & positional_array >= position);
end
% Trailing edge
if (position+thick > limits(2))
stimulated_trailing = find(positional_array < last_position+thick & positional_array >= position+thick);
end
end
% Constrain by length
stimulated_trailing = intersect(stimulated_trailing,length_valids);
stimulated_leading = intersect(stimulated_leading,length_valids);
if (edge_sf(1) > 0)
stimulated_leading = stimulated_leading(1:round(length(stimulated_leading)*edge_sf(1)));
else
stimulated_leading = [];
end
if (edge_sf(2) > 0)
stimulated_trailing = stimulated_trailing(1:round(length(stimulated_trailing)*edge_sf(2)));
else
stimulated_trailing = [];
end
% And append returned arrays if needbe
if (length(stimulated_leading) > 0)
for s=1:length(stimulated_leading)
synpos(length(timing) + 1, :) = [synmap(stimulated_leading(s),1) synmap(stimulated_leading(s),2)];
timing(length(timing) + 1) = time;
end
end
if (length(stimulated_trailing) > 0)
for s=1:length(stimulated_trailing)
synpos(length(timing) + 1, :) = [synmap(stimulated_trailing(s),1) synmap(stimulated_trailing(s),2)];
timing(length(timing) + 1) = time;
end
end
end
% Scale instantaneous velocity by dt -- in s
inst_vel.t_vec = 0:dt:total_time;
inst_vel.vel_vec = abs([inst_vel.vel_vec/(.001*dt)])';
%
% This should return timing and synapses for accelarating bars of CONSTANT
% ANGULAR SIZE AND CHANGING ANGULAR VELOCITY. Each position is activated
% TWICE - once for the ON and once for the OFF.
%
% RETURNS:
% synpos - [cmpt_idx cmpt_fraction]
% timing - in ms
%
% ASSIGNED:
% synmap - [cmpt cmpt_frac az el]
% direction - 'h' 'v' for horizontal, vertical (DIRECTION OF MOTION, not bar
% orientation!) - default is vertical
% velocity - in degrees / second (1): start (2): end
% limits - [start end] where to start; if horizontal motion, give azimuths;
% if vertical, elevations. Bar will 'edge in' from one side and edge
% out the other (i.e., it does not start wholly formed).
% dt - temporal resolution, in ms
% thick - thickness of the bar, in degrees (i.e., dimension along axis PARALLEL
% to motion vector)
% length_limits - edges of the bar, in degrees for dimension along axis PERPENDICULAR to
% motion vector
%
%
function [synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, direction, velocity, ...
limits, dt, thick, length_limits, ...
edge_sf)
% Preliminary
timing = [];
% Orientation-dependent component
if (strcmp('h', direction))
positional_array = synmap(:,4);
length_array = synmap(:,3);
else
positional_array = synmap(:,3);
length_array = synmap(:,4);
end
used = zeros(1,length(positional_array));
% Introduce variation into the eye -- vary each ommatidial view by sd of 5 degrees
eye_deformation = normrnd(0,10,length(positional_array),1);
positional_array = positional_array + eye_deformation;
% Determine valid based on length limits
length_valids = find(length_array > min(length_limits) & length_array < max(length_limits));
% determine total time based on the limits and the accelaration
% total time = total_distance/mean_velocity [assuming LINEAR accelaration]
mean_velocity = mean(velocity);
total_time = abs(abs(diff(limits)))/(mean_velocity/1000);
% instantaneous velocity vector
inst_vel.vel_vec = zeros(length(0:dt:total_time),1);
ivi = 1;
% loop thru time, and add accordingly
position = limits(1); % LEADING edge
dv_dt = diff(velocity)/1000/total_time; % again, convert to deg/ms2
v_o = velocity(1)/1000; % convert to d/ms
for time=0:dt:total_time
last_position = position;
vel_t = dv_dt*time + v_o;
dp = vel_t*dt;
% change in position/time is velocity!
int_vel.vel_vec = dp;
stimulated_leading = [];
stimulated_trailing = [];
% Increasing position value (e.g., rightward)
if (limits(1) < limits(2))
position = position + dp;
% Leading edge
if (position < limits(2))
stimulated_leading = find(positional_array >= last_position & positional_array < position);
end
% Trailing edge
if (position-thick > limits(1))
stimulated_trailing = find(positional_array >= last_position-thick & positional_array < position-thick);
end
% For decreasing position value
else
position = position - dp;
% Leading edge
if (position < limits(1))
stimulated_leading = find(positional_array < last_position & positional_array >= position);
end
% Trailing edge
if (position+thick > limits(2))
stimulated_trailing = find(positional_array < last_position+thick & positional_array >= position+thick);
end
end
% Constrain by length
stimulated_trailing = intersect(stimulated_trailing,length_valids);
stimulated_leading = intersect(stimulated_leading,length_valids);
if (edge_sf(1) > 0)
stimulated_leading = stimulated_leading(1:round(length(stimulated_leading)*edge_sf(1)));
else
stimulated_leading = [];
end
if (edge_sf(2) > 0)
stimulated_trailing = stimulated_trailing(1:round(length(stimulated_trailing)*edge_sf(2)));
else
stimulated_trailing = [];
end
% And append returned arrays if needbe
if (length(stimulated_leading) > 0)
for s=1:length(stimulated_leading)
synpos(length(timing) + 1, :) = [synmap(stimulated_leading(s),1) synmap(stimulated_leading(s),2)];
timing(length(timing) + 1) = time;
end
end
if (length(stimulated_trailing) > 0)
for s=1:length(stimulated_trailing)
synpos(length(timing) + 1, :) = [synmap(stimulated_trailing(s),1) synmap(stimulated_trailing(s),2)];
timing(length(timing) + 1) = time;
end
end
end
% Scale instantaneous velocity by dt -- in s
inst_vel.t_vec = 0:dt:total_time;
inst_vel.vel_vec = abs([inst_vel.vel_vec/(.001*dt)])';
%
% This should return timing and synapses for translating bars of CONSTANT
% ANGULAR SIZE AND ANGULAR VELOCITY. Each position is activated TWICE - once
% for the ON and once for the OFF.
%
% RETURNS:
% synpos - [cmpt_idx cmpt_fraction]
% timing - in ms
%
% ASSIGNED:
% synmap - [cmpt cmpt_frac az el]
% direction - 'h' 'v' for horizontal, vertical (DIRECTION OF MOTION, not bar
% orientation!) - default is vertical
% velocity - in degrees / second
% limits - [start end] where to start; if horizontal motion, give azimuths;
% if vertical, elevations. Bar will 'edge in' from one side and edge
% out the other (i.e., it does not start wholly formed).
% dt - temporal resolution, in ms
% thick - thickness of the bar, in degrees (i.e., dimension along axis PARALLEL
% to motion vector)
% length_limits - edges of the bar, in degrees for dimension along axis PERPENDICULAR to
% motion vector
% edge_sf - scales the input strength of each edge
%
%
function [synpos, timing, inst_vel] = pattern_translating_bar(synmap, direction, velocity, ...
limits, dt, thick, length_limits, ...
edge_sf )
% Preliminary
timing = [];
% Orientation-dependent component
if (strcmp('h', direction))
positional_array = synmap(:,4);
length_array = synmap(:,3);
else
positional_array = synmap(:,3);
length_array = synmap(:,4);
end
used = zeros(1,length(positional_array));
% Introduce variation into the eye -- vary each ommatidial view by sd of 5 degrees
eye_deformation = normrnd(0,10,length(positional_array),1);
positional_array = positional_array + eye_deformation;
% Determine valid based on length limits
length_valids = find(length_array > min(length_limits) & length_array < max(length_limits));
% determine total time ... limits + thickness scaled by velocity
total_time = abs(abs(diff(limits)) + thick)/(velocity/1000);
% loop thru time, and add accordingly
position = limits(1); % LEADING edge
dp = (velocity/1000)*dt;
for time=0:dt:total_time
last_position = position;
stimulated_leading = [];
stimulated_trailing = [];
% Increasing position value (e.g., rightward)
if (limits(1) < limits(2))
position = position + dp;
% Leading edge
if (position < limits(2))
stimulated_leading = find(positional_array >= last_position & positional_array < position);
end
% Trailing edge
if (position-thick > limits(1))
stimulated_trailing = find(positional_array >= last_position-thick & positional_array < position-thick);
end
% For decreasing position value
else
position = position - dp;
% Leading edge
if (position < limits(1))
stimulated_leading = find(positional_array < last_position & positional_array >= position);
end
% Trailing edge
if (position+thick > limits(2))
stimulated_trailing = find(positional_array < last_position+thick & positional_array >= position+thick);
end
end
% Constrain by length
stimulated_trailing = intersect(stimulated_trailing,length_valids);
stimulated_leading = intersect(stimulated_leading,length_valids);
if (edge_sf(1) > 0)
stimulated_leading = stimulated_leading(1:round(length(stimulated_leading)*edge_sf(1)));
else
stimulated_leading = [];
end
if (edge_sf(2) > 0)
stimulated_trailing = stimulated_trailing(1:round(length(stimulated_trailing)*edge_sf(2)));
else
stimulated_trailing = [];
end
% And append returned arrays if needbe
if (length(stimulated_leading) > 0)
for s=1:length(stimulated_leading)
synpos(length(timing) + 1, :) = [synmap(stimulated_leading(s),1) synmap(stimulated_leading(s),2)];
timing(length(timing) + 1) = time;
end
end
if (length(stimulated_trailing) > 0)
for s=1:length(stimulated_trailing)
synpos(length(timing) + 1, :) = [synmap(stimulated_trailing(s),1) synmap(stimulated_trailing(s),2)];
timing(length(timing) + 1) = time;
end
end
end
% Scale instantaneous velocity by dt -- in s
inst_vel.t_vec = 0:dt:total_time;
inst_vel.vel_vec = abs(velocity*ones(size(inst_vel.t_vec)));
%
% This should return timing and synapses for looming stimulus. Each position
% is activated only once. The stimulus is square.
% synpos - [cmpt_idx cmpt_fraction]
% timing - in ms
% synmap - [cmpt cmpt_frac az el]
% l_over_v - the l/v parameter
% center - [el az] where to start
% theta_start - size at start; zero would take infinite time. in degrees.
% dt - temporal resolution, in ms
%
function [synpos, timing] = pattern_square_looming_stimulus(synmap, l_over_v, center, theta_start, dt)
% Preliminary
timing = [];
theta_final = 80; % Degrees - final size, as in screen
% Start, end time - compute it
start_time = l_over_v/(tand(theta_start/2));
end_time = l_over_v/(tand(theta_final/2));
% GEnerate a center-centric coordinate set
centered_el = synmap(:,4)-center(1);
centered_az = synmap(:,3)-center(2);
% determine total time ...
total_time = start_time - end_time;
% loop thru time, and add accordingly
theta_0 = theta_start; % Theta at START of step
for time=0:dt:total_time
theta_f = 2*atand(l_over_v/(total_time + end_time - time));
% Find any points between theta_0 and theta_f object size
stimulated_1 = find(centered_el < theta_f & centered_el >= (-1*theta_f) & centered_az > theta_0 & centered_az <= theta_f);
stimulated_2 = find(centered_el < theta_f & centered_el >= (-1*theta_f) & centered_az < -1*theta_0 & centered_az >= -1*theta_f);
stimulated_3 = find(centered_az > -1*theta_0 & centered_az < theta_0 & centered_el >= -1*theta_f & centered_el < -1*theta_0);
stimulated_4 = find(centered_az > -1*theta_0 & centered_az < theta_0 & centered_el <= theta_f & centered_el > theta_0);
stimulated = unique([stimulated_1' stimulated_2' stimulated_3' stimulated_4']);
% And append returned arrays if needbe
if (length(stimulated) > 0)
for s=1:length(stimulated)
synpos(length(timing) + 1, :) = [synmap(stimulated(s),1) synmap(stimulated(s),2)];
timing(length(timing) + 1) = time;
end
end
theta_0 = theta_f;
end
%
% This should return timing and synapses for looming stimulus. Each position
% is activated only once. The stimulus is square.
% synpos - [cmpt_idx cmpt_fraction]
% timing - in ms
% inst_vel - the instantaneous velocity of the stimulus
%
% synmap - [cmpt cmpt_frac az el]
% l_over_v - the l/v parameter
% center - [el az] where to start
% theta_start - size at start; zero would take infinite time. in degrees.
% dt - temporal resolution, in ms
%
function [synpos, timing, inst_vel] = pattern_circular_looming_stimulus(synmap, l_over_v, center, theta_start, dt)
% Preliminary
timing = [];
theta_final = 80; % Degrees - final size, as in screen
recede = 0;
if (l_over_v < 0) % recede
recede = 1;
l_over_v = -1*l_over_v;
theta_start = -1*theta_start;
end
% Start, end time - compute it
start_time = l_over_v/(tand(theta_start/2));
end_time = l_over_v/(tand(theta_final/2));
% GEnerate a center-centric coordinate set
centered_el = synmap(:,4)-center(1);
centered_az = synmap(:,3)-center(2);
% determine total time ...
total_time = start_time - end_time;
% For each point, determine the angular deflection from the center
for s=1:length(centered_el)
alpha(s) = acosd(cosd(centered_el(s))*cosd(centered_az(s)));
end
% instantaneous velocity vector
inst_vel.vel_vec = zeros(length(0:dt:total_time),1);
ivi = 1;
% loop thru time, and add accordingly
if (recede == 0)
theta_0 = theta_start; % Theta at START of step
for time=0:dt:total_time
theta_f = 2*atand(l_over_v/(total_time + end_time - time));
% Find any points between theta_0 and theta_f object size
stimulated = find(alpha < theta_f & alpha > theta_0);
% And append returned arrays if needbe
if (length(stimulated) > 0)
for s=1:length(stimulated)
synpos(length(timing) + 1, :) = [synmap(stimulated(s),1) synmap(stimulated(s),2)];
timing(length(timing) + 1) = time;
end
end
inst_vel.vel_vec(ivi) = theta_f - theta_0; ivi = ivi + 1;
theta_0 = theta_f;
end
else
theta_0 = theta_final;
for time=total_time:-1*dt:0
theta_f = 2*atand(l_over_v/(total_time + end_time - time));
% Find any points between theta_0 and theta_f object size
stimulated = find(alpha > theta_f & alpha < theta_0);
% And append returned arrays if needbe
if (length(stimulated) > 0)
for s=1:length(stimulated)
synpos(length(timing) + 1, :) = [synmap(stimulated(s),1) synmap(stimulated(s),2)];
timing(length(timing) + 1) = total_time-time;
end
end
inst_vel.vel_vec(ivi) = theta_f - theta_0; ivi = ivi + 1;
theta_0 = theta_f;
end
end
% Scale instantaneous velocity by dt -- in s
inst_vel.t_vec = 0:dt:total_time;
inst_vel.vel_vec = abs([inst_vel.vel_vec/(.001*dt)])';