Ale_script.m

%% create MEA object
fileNameMEA='\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X42_20150603\MEA\slab_B_8.mcd';
MC = MCRackRecording(fileNameMEA);

% % Get digital data values and timestamps % %
    % for older recordings: only one digital line used
    [Dig_data,Digital_timestamps] = getDigitalData(MC);
    Dig_data = squeeze(Dig_data);
    Frame_Clock = Dig_data(1,:);
    % for more recent recording: multiple digital lines to communicate with
    % Andor
%     [Dig_data,Digital_timestamps] = getDigitalData(MC,'maxNumberOfDigitalChannels',16);
%     Dig_data = squeeze(Dig_data(1,1,:));
%     Dig_data = squeeze(Dig_data([1,3:5],:,:));
%     Frame_Clock = Dig_data(4,:);
%     Ext_shutter = Dig_data(2,:);

% % Get triggers
    [Trig_ms] = getTrigger(MC);

% % Analog data
    [V_analog, Analog_timestamps] = getAnalogData(MC,1,1,MC.recordingDuration_ms);
    V_analog = squeeze(V_analog);

% % timestamps of frame aquisition (ms) from digital channel
    frameStarts=Digital_timestamps(   find(Frame_Clock(1:end-1)==0 & Frame_Clock(2:end)==1)   );
    frameEnd=Digital_timestamps(   find(Frame_Clock(1:end-1)==1 & Frame_Clock(2:end)==0)   );
    frameCent=(frameEnd+frameStarts)/2; % timestamps for Andor frames
    % % timestamps of frame aquisition (ms) from Analog channel
    % frameStarts=Analog_timestamps(   find(V_analog(1:end-1)<3e6 & V_analog(2:end)>3e6)   );
    % frameEnd=Analog_timestamps(   find(V_analog(1:end-1)>3e6 & V_analog(2:end)<3e6)   );
    % frameCent=(frameEnd+frameStarts)/2; % timestamps for Andor frames

% % start and end of Andor recording
    tStart=Digital_timestamps(find(Frame_Clock>0,1,'first'));
    tEnd=Digital_timestamps(find(Frame_Clock>0,1,'last'));

% % check activity MEA over physical space 1 event
    load layout_100_16x16;
    clear pSegment vsd_t_ms eventStartMs eventWindow
    frame_start = 1; frame_end = 8000;
    pSegment = find( frameCent>=frameCent(frame_start) & frameCent<=frameCent(frame_end)); % time window for analysis
    vsd_t_ms = frameCent(pSegment);
    eventStartMs = frameCent(pSegment(1));
    eventWindow = frameCent(pSegment(end))-frameCent(pSegment(1));
    % all channels
%     [data,timestamps] = MC.getData(En(isfinite(En)),eventStartMs,eventWindow);
    % single channel
    [data,timestamps] = MC.getData(248,eventStartMs,eventWindow);


% % downsample data
    digdata = Dig_data(find(Digital_timestamps>= eventStartMs+timestamps(1),1,'first') : ...
        find(Digital_timestamps<= eventStartMs+timestamps(end),1,'last'));
    % data_downsampled = data(:,1,find(digdata>0));
    % timestamps_downsampled = timestamps(find(digdata>0));
    % figure,
    % [hCbar]=IntensityPhysicalSpacePlot(En(isfinite(En)),mean(squeeze(abs(data)),2),En,'plotElectrodeNumbers',0,'plotGridLines',0,'plotColorBar',1,'markerSize',50);
    % figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,En(isfinite(En)),squeeze(data_downsampled),...
    %     En,'DrawElectrodeNumbers',0);
    % title('MEA traces, aquisition times of frames 4300:4500','FontSize',20,'FontWeight','bold')


%% open Andor file
    fileNameAndor='\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X42_20150603\ANDOR\slab_B_9.sif';
    [~,metatemp] = sifread_ale(fileNameAndor,1);
    aquisition_rate = 1/metatemp.IntegrationCycleTime;
%     [val frame_peak]=min(abs(frameCent - 56000));
%     framestart = round(frame_peak- 8*aquisition_rate);
%     frameend = floor(frame_peak+ 8*aquisition_rate);
    framestart = find(frameCent>6.*10^4,1,'first');
    frameend = find(frameCent > 8.5.*10^4,1,'first');
    [VSData VSDmeta] = sifread_ale(fileNameAndor,framestart,'window',frameend-framestart+1); 
    VSDmeta = VSDmeta(1);
% Add extra fields to metadata
    exp_id = 5;
    % add extra info experiment
    VSDmeta.objective = 4;
    VSDmeta.excitation_wavelength = 550;
    VSDmeta.led_power = 100;
    VSDmeta.dye_uM = 200;
    VSDmeta.exc_filter = '545/30';
    VSDmeta.em_filter = '620/60';
    VSDmeta.frame_start = framestart; 
    VSDmeta.frame_end = frameend;

% Create cell array to save analysis data
%     file_dir = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X42_20150603\ANDOR\MATLAB\SampleB\ROI_MEA_channel_248';
%     filename = 'data_slab_B_wavelength.mat';
%     load([file_dir,filesep,filename]) % load previous analysis data
    % data_analysis{exp_id}.VSD_data = VSData;
    data_analysis{exp_id}.VSD_metadata = VSDmeta;
    % data_analysis{exp_id}.rec_name = 'X42_20150605_SampleB_ROI_MEAch_248';

% apply rotation on data saved with spooling
% rotation of spooling (+90 deg)
rotation_MEA = 0; % = MEA is flipped vertically and horizontaly
VSDmeta.MEA_orientation = rotation_MEA;
angle_rotation = rotation_MEA + 90;
% apply rotation on data saved with spooling
datatemp = imrotate(VSData(:,:,1),angle_rotation);
data = zeros([size(datatemp,1),size(datatemp,2),size(VSData,3)],'uint16');
for i = 1:size(VSData,3)
    data(:,:,i) = imrotate(VSData(:,:,i),angle_rotation);
end
VSData = data; clear data
[pixY,pixX,nFrames]=size(VSData);

% % Locate epileptiform events in average fluorescence trace over time
% x_pixels = 1:pixX; y_pixels = 1:pixY;
% % avg_I = squeeze(mean(mean(VSData(x_pixels,y_pixels,:),1),2));
% 
% avg_I = reshape(VSData,pixX*pixY,size(VSData,3)); 
% avg_I = squeeze(mean(avg_I,1));
% nFrames = length(avg_I);
% window = round(1/VSDmeta.KineticCycleTime * 15); % 1500 at 100 Hz % median filter (fast implementation compiled in C)
% minSBWnd=300;
% minSBInterval=300;
% avgSignal= avg_I'; %squeeze(mean(mean(VSData,2),3));
% xic = avgSignal(1:floor(window/2));
% xfc = avgSignal(nFrames - floor(window/2)+1:nFrames);
% avgF0 = fastmedfilt1d(avgSignal,window,xic,xfc);
% if ismember(VSDmeta.excitation_wavelength, [405,490,500,525,550])
%    sign_factor = +1;
% elseif ismember(VSDmeta.excitation_wavelength, [435,460,470])
%    sign_factor = +1;
% end
% 
% avgdFoF= sign_factor * (avgSignal-avgF0)./avgF0;
% frameCent = 1:nFrames;
% [BS,BP,BE,BI]=eventLocator(avgdFoF,frameCent,[1;1;1;nFrames],minSBWnd,minSBInterval,'stdThresh',20,'medianWindow',1500);
% 
% % Accept events 
% clear event_peakframes
% select_events = find(BI>0.00002);
% event_starts = BS(select_events);
% event_peaktime = BP(select_events);
% for i = 1:length(event_peaktime)
%     frames_temp = find(frameCent <= event_peaktime(i),1,'last') -48 : find(frameCent >= event_peaktime(i),1,'first')+48  ;
%     [~, ind_frame] = max(avgdFoF(frames_temp));
%     event_peakframes(i) = frames_temp(ind_frame);
% end
% 
% % comparison of average dFoF for different wavelengths
% % or get event peak frames from MEA timestamps
% event_peakframes = peak_frames;
% event_id = 4;
% f = 1/VSDmeta.KineticCycleTime; %f = 1/data_rec(i).metadata.KineticCycleTime;
% frames_window = round(5*f);
% if rem(frames_window,2) == 0
%     frames_window = frames_window+1;
% end
% % n_tests = 6; 
% % avgdFoF_window = zeros(frames_window,n_tests);
% avgdFoF_window(:,exp_id) = squeeze(avgdFoF(event_peakframes(event_id)-floor(frames_window/2):event_peakframes(event_id)+floor(frames_window/2)));

%% Mask for VSD recording
% obtain pixels slab
avg_image = mean(VSData,3); 
figure, subplot(221), imagesc(avg_image), colorbar, title('image of slab')
% automatic threshold based on intensity values
[n,bin] =  hist(avg_image(:));  
subplot(222), bar(bin,n), title('distribution intensities all pixels'), hold on
[n_sorted,sort_order] = sort(n,'descend'); modes = sort(sort_order(1:2)); 
threshold_intensity = bin(sort_order(find(sort_order > modes(1) & sort_order < modes(2),1,'last')));
th_mask = avg_image; th_mask(find(avg_image < threshold_intensity)) = 0;
bw_mask = logical(th_mask); 
lim_axis = ylim; plot([threshold_intensity threshold_intensity],lim_axis,'r','LineWidth',4);
% label object and keep object corresponding to slab
L = bwlabel(bw_mask);
objects = regionprops(bw_mask,'Area'); [~,obj_ind] = max([objects(:).Area]);
slab_mask = bw_mask;  
slab_mask(find(L ~= obj_ind)) = 0; 
pixels = find(logical(slab_mask));
subplot(223), imagesc(slab_mask), colorbar, title('slab mask')
subplot(224), hist(avg_image(pixels)), title('distribution intensity slab pixels')

% find pixels along diagonal of frame
pixel_id = zeros(pixY,pixX);
pixel_id = 1:pixY*pixX;
pixel_id = reshape(pixel_id,pixY,pixX);

% reference image recording
ref_image_filename = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X42_20150603\ANDOR\slab_B_refimage_1.sif';
[ref_image] = sifread_ale(ref_image_filename);
refimage = mean(ref_image,3);
figure, imagesc(refimage), colorbar
ref_rotated = imrotate(refimage,90);
figure, imagesc(ref_rotated), colorbar, title('rotated')
figure, imshow(ref_rotated,[min(ref_rotated(:)) max(ref_rotated(:))]), colorbar, title('rotated')


% Rotation of images to original orientation in focal plane
rotated_image = imrotate(avg_image,90);
figure, imagesc(rotated_image), colorbar, title('10x view of sample 2 rec0005 15April, rotated 90deg')

% reference image masked based on SNR map
image = imrotate(data_rec.downsampled_image.*snr_bw,90);
figure, imagesc(image), colormap('hot')

%% MEA and VSD comparison: single channel, full recording 
channel = 187;
nFrames = VSDmeta.NumberImages; %data_analysis{exp_id}.metadata.NumberImages %VSDmeta.NumberImages
frame_Start = framestart;
frame_End = frameend;
MEA_tStart = frameCent(frame_Start); 
MEA_tEnd = frameCent(frame_End);
time_interval = MEA_tEnd-MEA_tStart;    %frameEnd(nFrames)-tStart; %tEnd-tStart;
[data_1D, timestamps_1D] = MC.getData(channel,MEA_tStart,time_interval);    %MC.getData(channel,tStart,time_interval)
data_1D = squeeze(data_1D); timestamps_1D = squeeze(timestamps_1D);
vsd_ROI = VSData(215:224,100:109,:); 
n_rows = size(vsd_ROI,1); n_cols = size(vsd_ROI,2);
mean_trace = squeeze(mean(mean(vsd_ROI)));
vsd_ROI = reshape(vsd_ROI,size(vsd_ROI,1)*size(vsd_ROI,2),size(vsd_ROI,3));
% vsd_ROI = vsd_ROI(diag(pixel_id),:);
dF_ROI = []; dFoF_ROI = []; F0 = [];
% preallocate :
dF_ROI = zeros(size(vsd_ROI)); dFoF_ROI = zeros(size(vsd_ROI)); F0 = zeros(size(vsd_ROI));
f = 1/VSDmeta.KineticCycleTime; % Hz, frequency data aquisition
window = round(1/VSDmeta.KineticCycleTime * 10); %duration epileptiform event, about 1 sec 
tot_length = size(dF_ROI,2);
if  ismember(VSDmeta.excitation_wavelength, [525,550,435,460,470,490])
    sign_peaks = -1;
elseif ismember(VSDmeta.excitation_wavelength, [405])
    sign_peaks = +1;
end
for p = 1:size(vsd_ROI,1)   %% note:  1 loop takes about 0.48 sec
    x = double(vsd_ROI(p,:));
    xic = x(1:floor(window/2));
    xfc = x((tot_length - floor(window/2)+1):tot_length);
    F0temp = fastmedfilt1d(x, window, xic, xfc);
    F0(p,:) = F0temp;
    temp1 = x(:)-F0temp(:);
    temp2 = temp1./F0temp;
    dF_ROI(p,:) = temp1 *sign_peaks; dFoF_ROI(p,:) = temp2 *sign_peaks;  
end
% % plot MEA and VSD raw trace % %
    x1 = (MEA_tStart+timestamps_1D)./1000;
    y1 = data_1D;
    x2 = frameCent(frame_Start:frame_End)./1000;
    y2 = mean(dFoF_ROI,1);
    % add correction for alignment traces
    y2_mean = mean(y2);
    y2_shift = -60;
    fig = figure,[ax,p1,p2] = plotyy(x1,y1,x2,y2); %mean(DFoF,1) %smooth(DFoF(76,:),25,'sgolay',3)
    for line= 1:length(p2), set(p2(line),'Color','r','LineWidth',2.4), end
    set(ax(2),'YColor','r'), 
    xlabel(ax(2), 'time, s')
    ylabel(ax(1), 'Voltage, uV')
    ylabel(ax(2), 'Raw Fluorescence','Color','r')
    legend('MEA, channel 248','VSD, 10x10 pixel average')
% set axis with predefined scale and limits
    set(ax(1),'YTick',[-20 0 20])
    set(ax(2),'YTick',[5.62 5.64 ].*10^4)
% format image for presentation
    poster_fig_booster('font_weight','bold','legend_position','South',...
        'figure_coordinates',[-1400 800 1200 600])
% computation of mean trace and smoothing:
    f_mea = 20000;
    x = double(mean_trace);
    window = 3;
    tot_lenght = length(x);
    xic = x(1:floor(window/2));
    xfc = x((tot_length - floor(window/2)+1):tot_length);
    trace_smoothed = fastmedfilt1d(x, window, xic, xfc);
    trace_smoothed = smooth(x, 10, 'sgolay',3);
    figure, p = plot(x2,mean_trace), hold on
    set(p,'Color','r','LineWidth',2.4);
    plot(x2,trace_smoothed,'k','LineWidth',2.2),
    xq = [x2(1)*1000:f/f_mea:x2(end)*1000];
    trace_interp = interp1(x2.*1000,trace_smoothed,xq,'linear');
    plot(xq./1000,trace_interp,'g')
% plot uV trace and finite difference of VSD signal
%     fig = figure,[ax,p1,p2] = plotyy(x1,y1,x2(2:end),diff(trace_smoothed));

% % plot MEA and VSD dF/F % %
    x1 = (MEA_tStart+timestamps_1D)./1000;
    y1 = data_1D;
    x2 = frameCent(frame_Start:frame_End)./1000;
    y2 = mean(dFoF_ROI,1);
    fig = figure,[ax,p1,p2] = plotyy(x1,y1,x2,y2); %mean(DFoF,1) %smooth(DFoF(76,:),25,'sgolay',3)
    for line= 1:length(p2), set(p2(line),'Color','r','LineWidth',2.4), end
    set(ax(2),'YColor','r'), 
    xlabel(ax(2), 'time, s')
    ylabel(ax(1), 'Voltage, uV')
    ylabel(ax(2), 'mean dFoF','Color','r')
    legend('MEA, channel 187','VSD, 10x10 pixel average')
% set axis with predefined scale and limits
    centrepoint = 16.930; %find(x1>=15000,1,'first');
    half_interval = 1.500; % second duration * acquisition rate
    set(ax(1:2),'XLim',[67 83]); % xlim = [centrepoint-half_interval centrepoint+half_interval];
    set(ax(1),'YLim',[-600 600],'YTick',[-500 0 500])
    set(ax(2),'YLim',[-0.011 0.011],'YTick',[-0.01 -0.005 0 0.005 0.01]), %'YTick',[-0.01 0.01]
% format image for presentation
    poster_fig_booster('font_weight','bold','legend_position','South',...
        'figure_coordinates',[-1400 1200 1000 600])

% manual thresholding and spike finder
window_analysis = find(timestamps_1D>=1.43e5 & timestamps_1D<=1.62e5);
time_vector_spk = timestamps_1D(window_analysis);
spike_times = time_vector_spk(find(data_1D(window_analysis)<=-20));
spike_trace = data_1D(window_analysis);
spike_dist = pdist(spike_times(:));
spike_link = linkage(spike_dist);
n_clusters = 12;
spike_cluster = cluster(spike_link,'maxclust',n_clusters);
spike_event = zeros(1,n_clusters); spike_frames = zeros(1,n_clusters);
for i = 1:n_clusters
    i
    spike_event(i) = mean(spike_times(find(spike_cluster == i)));
    spike_frames(i) = find(frameCent>=spike_event(i),1,'first');
end
spike_event = sort(spike_event);   spike_frames = sort(spike_frames); 


%% spatial averaging on VSD movie matrix

% using block processing
datatemp = VSData;
block_size = [2 2];
dataout = zeros(size(datatemp,1)/block_size(1),size(datatemp,2)/block_size(2),...
    size(datatemp,3),'uint16');
for i = 1:size(dataout,3)
%     avg_fun = @(block_struct) uint16(mean2(block_struct.data));
%     I = datatemp(:,:,i);
%     I2 = blockproc(I,block_size,avg_fun);   
    I = double(datatemp(:,:,i));
    I2 = uint16(BlockMean(I,2,2)); 
    dataout(:,:,i) = I2;
end

%% timestamps from MEA

% locate events on MEA trace
rec_length = length(tStart+timestamps_1D);
timevector = tStart+timestamps_1D;
[BS,BP,BE,BI]=eventLocator(data_1D,tStart+timestamps_1D,[1;1;1;rec_length],500,200); %'sigma',1,'res',2
% select events
peak_times = BP(find(BI>3e5))
% find frames of peak events
peak_frames = zeros(size(peak_times));
for i = 1:length(peak_times)
    peak_frames(i) = find(frameCent>=peak_times(i),1,'first');
end
data_analysis{exp_id}.peak_frames = peak_frames;

% plot vsd traces with overlaid timestamps from MEA
figure, plot(frameCent,mean_trace,'r','LineWidth',2), hold on
ylimits = ylim;
for i = 1:length(peak_times)
    if i == 1
        plot([peak_times(i) peak_times(i)],[ylimits(1) ylimits(2)],'k','LineWidth',1.9)
        t = [1.5e4 peak_times(i)]; 
        dt(i) = frameCent(find(frameCent<t(2),1,'last')) - frameCent(find(frameCent>t(1),1,'first')) ;
        f = find(frameCent>=t(1) & frameCent<t(2));
        timepoints = frameCent(f);
        mean_F(i) = mean(mean_trace(f));
        dF(i) = mean_trace(f(1)) - mean_trace(f(end));
        plot([timepoints(1) timepoints(end)],[mean_F(i) mean_F(i)],'LineWidth',2.2)
    else
        plot([peak_times(i) peak_times(i)],[ylimits(1) ylimits(2)],'k','LineWidth',1.9)
        t = [peak_times(i-1) peak_times(i)];
        dt(i) = frameCent(find(frameCent<t(2),1,'last')) - frameCent(find(frameCent>t(1),1,'first')) ;
        f = find(frameCent>t(1) & frameCent<t(2));
        dF(i) = mean_trace(f(1)) - mean_trace(f(end));
        timepoints = frameCent(f);
        mean_F(i) = mean(mean_trace(f));
        plot([timepoints(1) timepoints(end)],[mean_F(i) mean_F(i)],'LineWidth',2.2);
    end
end

%% Create data_rec structure array
i = 1; jj= 1;
data_rec(i,jj).data = VSData;
data_rec(i,jj).metadata = VSDmeta;
clear VSData VSDmeta

%% Check raw data

% select three groups of pixels belonging to the slab (high,medium and low
% intensity pixels)
pixel_analysis = data_analysis{1}.pixels;
npixels = 100;
avg_factor = 10; % multiple integer number of pixels averaged together
baseline = [1000:1800,2000:2800,8400:9400]; % baseline period
event_peakframes = peak_frames;
data_analysis{exp_id} = raw_data_plots( data_rec, npixels, avg_factor, event_peakframes, baseline, pixel_analysis );
data_analysis{exp_id}.metadata = data_rec.metadata;


%% Visualize Spike-triggered projection fields (VSD measurements)

data_baseline = zeros(size(VSData),'uint16');
for i = 1:size(VSData,3)
    data_baseline(:,:,i) = imrotate(VSData(:,:,i),90);
end
metadata_baseline = VSDmeta(1);
clear VSData VSDmeta

data = zeros(size(VSData),'uint16'); 
for i = 1:size(VSData,3)
    data(:,:,i) = imrotate(VSData(:,:,i),90);
end
metadata_STA = VSDmeta(1);
clear VSData VSDmeta

% plotting mean and st
figure, subplot(221), imagesc(mean(data_baseline,3)), colorbar, title('mean fluorescence baseline'), axis square
subplot(222), imagesc(std(double(data_baseline),[],3)), colorbar, title('std fluorescence baseline'), axis square
subplot(223), imagesc(mean(data,3)), colorbar, title('STA raw fluorescence'), axis square
subplot(224), imagesc(std(double(data),[],3)), colorbar, title('std STA fluorescence'), axis square

std_STA = std(double(data),[],3)./std(double(data_baseline),[],3);
std_thresh_STA = std_STA;
std_thresh_STA(find(std_thresh_STA<3)) = 0;
background = uint16(mean(data_baseline,3));
diff_data = zeros(size(data),'uint16'); zscore_data = zeros(size(data),'uint16'); dfof_data = zeros(size(data));
std_baseline = std(double(data_baseline),[],3);
for i = 1:size(data,3)
   diff_data(:,:,i) = data(:,:,i) - background; 
   %dfof_data(:,:,i) = (double(data(:,:,i)) - double(background))./ double(background); 
   zscore_data(:,:,i) = double(data(:,:,i) - background)./std_baseline; 
end

figure, 
subplot(221), imagesc(std_STA), colorbar, title('std STA fluorescence / std baseline'), axis square
subplot(222), %imagesc(std_thresh_STA), colorbar, title('std STA fluorescence / std baseline'), axis square
imagesc(mean(data,3) - mean(data_baseline,3)), colorbar, title('mean STA F - mean baseline F'), axis square
subplot(223), imagesc(mean(diff_data,3)), colorbar, title('mean STA dF'), axis square
subplot(224), imagesc(mean(zscore_data,3)), colorbar, title('mean STA zscored F'), axis square

% variability STA fluorescence
avg_STA_F = mean(data,3) - mean(data_baseline,3);
lowcut = -100; highcut = 300;
lower = find(avg_STA_F <= lowcut & logical(slab_mask)); 
higher = find(avg_STA_F >= highcut & logical(slab_mask)); 
datatemp = double(reshape(data,size(data,1)*size(data,2),size(data,3)));
baselinetemp = double(reshape(data_baseline,size(data_baseline,1)*size(data_baseline,2),size(data_baseline,3)));
minimum = double(min([min(min(datatemp(pixels,:))) min(min(baselinetemp(pixels,:)))]));
maximum = double(max([max(max(datatemp(pixels,:))) max(max(baselinetemp(pixels,:)))]));
x = linspace(minimum,maximum,10);
[l_bsl,~] = hist(reshape(baselinetemp(lower,:),1,numel(baselinetemp(lower,:))), x);
[l_sta,~] = hist(reshape(datatemp(lower,:),1,numel(datatemp(lower,:))), x );
[h_bsl,~] = hist(reshape(baselinetemp(higher,:),1,numel(baselinetemp(higher,:))), x );
[h_sta,~] = hist(reshape(datatemp(higher,:),1,numel(datatemp(higher,:))), x );
figure, subplot(121), bar(x,l_bsl./sum(l_bsl),'BarWidth',0.9), hold on, bar(x,l_sta./sum(l_sta),'c','BarWidth',0.7)
subplot(122), bar(x,h_bsl./sum(h_bsl),'BarWidth',0.9,'FaceColor',[0.3 0 0]), hold on, bar(x,h_sta./sum(h_sta),'r','BarWidth',0.7)
figure,
plot([squeeze(mean(baselinetemp(lower,:),2))./squeeze(mean(baselinetemp(lower,:),2)) squeeze(mean(datatemp(lower,:),2))./squeeze(mean(baselinetemp(lower,:),2))]','o-','Color',[0 0 0.7])
hold on
plot([squeeze(mean(baselinetemp(higher,:),2))./squeeze(mean(baselinetemp(higher,:),2)) squeeze(mean(datatemp(higher,:),2))./squeeze(mean(baselinetemp(higher,:),2))]','o-','Color',[0.7 0 0])

% countour plot
avg_STA_norm = (mean(data,3) ./ mean(data_baseline,3)) - 1;
avg_STA_norm(find(~logical(slab_mask))) = NaN;
figure, contour(avg_STA_norm,[-5.5:1.5:10.5].*1e-3)

% MEA grid
x_grid = zeros(1,16); y_grid = zeros(1,16);
x_grid(1) = 5; y_grid(1) = 14;
step = 15.5;
for i = 2:16
   x_grid(i) = round(x_grid(i-1) +  step-(0.005*i));
   y_grid(i) = round(y_grid(i-1) +  step-(0.005*i));
end
[X,Y] = meshgrid(x_grid,y_grid);
mea_contact = avg_image;
mea_contact(Y,X) = NaN; 
clf, imagesc(mea_contact)


%% Segment MEA contact sites


n_hood = ones(3); n_hood([1,3,7,9]) = 0;
n_hood = ones(5); n_hood([1,2,4,5,6,10,16,20,21,22,24,25]) = 0;
% apply a standard deviation based spatial filter to enhance MEA contacts
n_hood = ones(7); n_hood(sort([1,2,3,5,6,7,8,9,13,14,15,21,50-[1,2,3,5,6,7,8,9,13,14,15,21]])) = 0;
filt_image = stdfilt(avg_image,n_hood);
% blur the MEA contacts to find boundaries better
dx = 1; dy = 1; width = 2;
kernel_con = GenGauss(dx, dy, width);
blur_image = conv2(filt_image,kernel_con,'same');
%  Segment image into foreground and background using active contour
se = strel('disk',2);
% slab_mask1 = avg_image; slab_mask1(find(slab_mask1 < 9000)) = 0; slab_mask1 = logical(slab_mask1); 
mask = imclose(slab_mask,se);
bw = activecontour(blur_image,mask,300);
% now erode out the MEA sites
se = strel('disk',1);   se = strel('line', 6,90)     
erodedBW = imerode(bw,se);
IM2 = imclearborder(erodedBW,4);
se = strel('disk',3); IM3 = imdilate(IM2,se);
% label connected elements
MEA = regionprops(IM3, 'Area', 'Centroid');
site_selected = find([MEA.Area]<150);
centroidtemp = reshape([MEA(site_selected(:)).Centroid],2,length(site_selected))';
x_centroids = centroidtemp(:,1); y_centroids = centroidtemp(:,2);
r_dist = pdist(y_centroids); c_dist = pdist(x_centroids);
r_tree = linkage(r_dist); c_tree = linkage(c_dist);
n_electrode_dim = 16; 
r_clusters = cluster(r_tree,'maxclust',n_electrode_dim); c_clusters = cluster(c_tree,'maxclust',n_electrode_dim);

% plot centroid positions over image 
% set the range of the axes: the image will be stretched to this
min_x = 1;
max_x = size(avg_image,2);
min_y = 1;
max_y = size(avg_image,1);
figure, imagesc([min_x max_x], [min_y max_y], avg_image);
hold on;
plot(x_centroids, y_centroids, 'ko','MarkerSize',5,'MarkerFaceColor','k');


%% dFoF MEA ROIs

x_MEA_centres = round(x_centroids);
y_MEA_centres = round(y_centroids);
dF_MEA_centres = [];
dFoF_MEA_centres = [];
for x = 1:length(x_MEA_centres)
    for y = 1:length(y_MEA_centres)
        
    end
end

%% dFoF for selected pixels
VSData=permute(VSData,[3,1,2]);
f = 1/VSDmeta.KineticCycleTime; % Hz, frequency data aquisition
f1 = 1000; % 1 frame every 1 ms
rows = [40,90,123]
cols = [106,106,106]
DF = zeros(size(VSData,1),length(rows));
DFoF = zeros(size(VSData,1),length(rows));
sDFoF = zeros(size(VSData,1),length(rows));
int_dFoF = zeros(length([1:f/f1:size(VSData,1)]),length(rows));
tot_length = size(VSData,1);
for i = 1:length(rows)
    x = double(VSData(:,rows(i),cols(i)));
    xic = x(1:floor(window/2));
    xfc = x((tot_length - floor(window/2)+1):tot_length);
    F0temp = fastmedfilt1d(x, window, xic, xfc);
    temp1 = x-F0temp;
    temp2 = temp1./F0temp;
    temp3 = smooth(temp2./F0temp,25,'sgolay',3); 
    xpoints = [1:size(temp3,1)];
    xq = [xpoints(1):f/f1:xpoints(end)];
    temp4 = interp1(xpoints,temp3,xq,'spline');
    DF(:,i) = temp1; DFoF(:,i) = temp2; sDFoF(:,i)=temp3; int_dFoF(:,i) = temp4';
end

%% Compute dFoF and obtain delays of events for all selected pixels
p_r = 1; p_c = 1; darkframe = 0; % spatial binning of pixels (steps along rows and columns) 
data_rec(1).data = VSData;
data_rec(1).metadata = VSDmeta;
% % specify information about imaging
data_rec(1).metadata.objective = 4;
data_rec(1).metadata.excitation_wavelength = 405;
data_rec(1).metadata.led_power = 35;
data_rec(1).metadata.dye_uM = 150;
% INDEX is a 3x1 cell array where you specify: {rows,cols,frames} to use for finding delays in DF/F traces
% 1) if index{1} = pixels; you enter the single pixels of the image you want to restrict your analysis to
% 2) if index is empty {rows,cols,frames} are taken from the metadata
% 3) if isempty(index{3}), {rows,cols,frames} with frames = total number of frames
index = {pixels};
data_rec = ANDOR_VSD_rise_delay(data_rec,1,1,p_r,p_c,index);

%% Compute pixel-wise SNR
pixels = 1:size(data_rec.dFoF,2);
plotflag = 1;
size_image = size(data_rec.downsampled_image);
[ data_rec ] = VSD_SNR( data_rec, pixels, plotflag, size_image);
SNR_map = reshape(data_rec.SNR_rois,size_image);
figure, imagesc(imrotate(SNR_map,90)), colorbar
snr_threshold = 4.5;
snr_mask = SNR_map;
snr_mask(find(SNR_map < snr_threshold)) = 0;
snr_bw = logical(snr_mask); figure, imagesc(imrotate(snr_bw.*SNR_map,90)), colorbar 
snr_pixels = find(logical(snr_mask));


%% MEA data out for a defined event epoch
tStart=Digital_timestamps(find(Dig_data>0,1,'first'));
tEnd=Digital_timestamps(find(Dig_data>0,1,'last'));
%
clear pSegment vsd_t_ms eventStartMs eventWindow
frame_start = 650; frame_end = 5650;
pSegment=find( frameCent>frameCent(frame_start) & frameCent<frameCent(frame_end)); % time window for analysis
vsd_t_ms = frameCent(pSegment);
eventStartMs=frameCent(pSegment(1));
eventWindow=frameCent(pSegment(end))-frameCent(pSegment(1));
[M1 event_ms1]=MC.getData(MC.channelNumbers,eventStartMs,eventWindow);
%
pSegment=find( frameCent>3.2e4 & frameCent<3.5e4 ); % time window for analysis
vsd_t_ms = frameCent(pSegment);
eventStartMs=frameCent(pSegment(1));
eventWindow=frameCent(pSegment(end))-frameCent(pSegment(1));
[M2 event_ms2]=MC.getData(MC.channelNumbers,eventStartMs,eventWindow);

% comparison of average VSD signal (all pixels) to average MEA signal (all
% channels)
vsd_trace = avgdFoF;
figure,[ax,p1,p2] = plotyy(frameCent(frame_start) +event_ms1',squeeze(mean(M1,1)),frameCent(frame_start:frame_end),vsd_trace(frame_start:frame_end))
for line= 1:length(p2), set(p2(line),'Color','r','LineWidth',2.4), end
set(ax(2),'YColor','r')
xlabel(ax(2), 'time, ms')
ylabel(ax(1), 'Voltage, uV')
ylabel(ax(2), 'dFoF','Color','r')
legend('MEA, avg all channels','VSD, average all slab pixels')


% plot MEA raw traces over physical space of array
load layout_100_16x16; 
channels = En(isfinite(En));
figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,channels(:),squeeze(M1(channels(:),:,:)),En,'DrawElectrodeNumbers',0);
title('MEA raw traces, event 1','FontSize',20,'FontWeight','bold')
figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,channels(:),squeeze(M2(channels(:),:,:)),En,'DrawElectrodeNumbers',0);
title('MEA raw traces, event 4','FontSize',20,'FontWeight','bold')
%
En1 = En; %fliplr(flipud(En));
figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,1:252,squeeze(M1),En1,'DrawElectrodeNumbers',0);
title('MEA raw traces, event 4','FontSize',20,'FontWeight','bold')
figure, [hCbar]=IntensityPhysicalSpacePlot(1:252,(max(squeeze(M1),[],2)-min(squeeze(M1),[],2))./std(squeeze(M1),[],2),En1,'plotElectrodeNumbers',0,'plotGridLines',0);
figure, [hCbar]=IntensityPhysicalSpacePlot(1:252,std(squeeze(M1),[],2),En1,'plotElectrodeNumbers',0,'plotGridLines',0);


%% low frequency propagation on MEA
En_ref = En(isfinite(En));
subchannels = En_ref(:);%[En_ref(:)',30,2,5,248,245,63,104,74,172,163];
% subchannels = setdiff(subchannels(:),[133,35,30,178,218,143,173,246]);

figure, [delaysMean_ms]=lowFreqPropagationPattern_Ale(M2(subchannels(:),:,:),...
    MC.samplingFrequency,subchannels(:),En,'minMaxDelay',[1000 2000]); %'minMaxDelay',[-150 +300]
% alternatively plot: 
F=figure('Position',[620   485   519   383]);h=axes;
[hCbar]=IntensityPhysicalSpacePlot(subchannels(:),squeeze(delaysMean_ms),...
    En,'plotElectrodeNumbers',0,'plotGridLines',0,'plotColorBar',1,'markerSize',50,...
    'h',h,'Ilim',[1300 1600]);

% plot multiple events of wave propagation from MEA
delaysMean_allevents = NaN(length(subchannels),length(BS)-1);
for i = 1:length(BS)-1
    pSegment=find( frameCent>BS(i)-500 & frameCent<BE(i)+500 ); % time window for analysis
    vsd_t_ms = frameCent(pSegment);
    eventStartMs=frameCent(pSegment(1));
    eventWindow=frameCent(pSegment(end))-frameCent(pSegment(1));
    [M event_ms]=MC.getData(MC.channelNumbers,eventStartMs,eventWindow);
    [delaysMean_ms]=lowFreqPropagationPattern_Ale(M(subchannels(:),:,:),MC.samplingFrequency,subchannels(:),En1(7:13,6:11));
    delaysMean_allevents(:,i) = squeeze(delaysMean_ms);
end
% plot map of average delay
F=figure('Position',[620   485   519   383]);
avg_delay_MEA = mean(delaysMean_allevents,2) - min(mean(delaysMean_allevents,2));
h=axes; minMaxDelay = [0 max(avg_delay_MEA(:))];
[hCbar]=IntensityPhysicalSpacePlot(subchannels(:),avg_delay_MEA,En1(7:13,6:11),'plotElectrodeNumbers',0,'plotGridLines',0,'plotColorBar',1,'markerSize',50,'h',h,'Ilim',minMaxDelay);


%% delay maps from VSD 
delays_rise = data_rec.event_HalfRiseTimes;
figure, imagesc(delays_rise(:,:,4)), colorbar
% % mask for display
ds_image = data_rec.downsampled_image; figure, imagesc(imrotate(ds_image,90)), colorbar, title('ref image delays,rotated')
threshold_intensity = 9000;
th_mask = ds_image;
th_mask(find(ds_image < threshold_intensity)) = 0;
mask = logical(th_mask); figure, imagesc(mask), colorbar, title('mask delays')
% or %
mask = snr_bw;
% % diplay masked delay map (here either use avg_image mask or SNR mask)
HRDelays = data_rec.event_HalfRiseTimes;
figure,
events_selected = [1,3:12];
c = 0;
for i = 1:length(events_selected)+1 
    if i<= length(events_selected) 
        image = HRDelays(:,:,events_selected(i)).* mask;
        min_delay = min(image(find(mask==1)));  %lowcut_delay = prctile(image(find(mask==1)),0.05);
        image = image - min_delay; %image = image - lowcut_delay; 
        max_delay = max(image(find(mask==1)));
        image(find(~mask)) = NaN; image = imrotate(image,90);  
        subplot(3,4,i), imagesc(image), axis square %highcut_delay = prctile( image( find(isfinite(image)) ) , 99.95), imagesc(image,color_lim) 
        h=colorbar; freezeColors %set(h, 'YTick', [0 max_delay],'YTickMode','manual','YTickLabel',{'0',mat2str(max_delay)})
        c = c+1;
        titlestring = strcat('event # ',num2str(c),' ; max delay: ',' ',num2str(max_delay),' ms' )
        title(titlestring,'FontSize',16,'FontWeight','bold');
    else
        image =  mean(HRDelays(:,:,events_selected(:)),3).*mask;
        mindelay = min(image(find(mask==1)));
        image = image - mindelay; 
        max_delay = max(image(find(mask==1)));
        image(find(~mask)) = NaN; image = imrotate(image,90);
        subplot(3,4,i), imagesc(image), axis square
        %color_lim = prctile((image(find(mask))),[0.05 99.95]); imagesc(image, color_lim), 
        h=colorbar;  freezeColors %set(h, 'YTick', [0 max_delay],'YTickMode','manual'),
        titlestring = strcat('mean all events ',' ; max avg delay: ',' ',num2str(round(max_delay)),' ms' )
        title(titlestring,'FontSize',16,'FontWeight','bold','Color','r');
    end
end

%% Overlap MEA and VSD in physical space

% rotate VSD data to the right orientation of the fov
sz_downsamp_image = size(data_rec.downsampled_image);
nFrames = size(data_rec.dFoF,1);
VSD_dFoF_mat = reshape(data_rec.dFoF,[nFrames,sz_downsamp_image]);
VSD_sIdFoF_mat = reshape(data_rec.dFoF_interpolated,[size(data_rec.dFoF_interpolated,1),sz_downsamp_image]);
VSD_dFoF_rotated = zeros(size(VSD_dFoF_mat));
VSD_sIdFoF_rotated = zeros(size(VSD_sIdFoF_mat));
for t = 1:size(VSD_dFoF_mat,1)
    temp_rotation = imrotate(squeeze(VSD_dFoF_mat(t,:,:)),90);
    VSD_dFoF_rotated(t,:,:) = temp_rotation;
end
for tt = 1:size(VSD_sIdFoF_mat,1)
    temp_rotation = imrotate(squeeze(VSD_sIdFoF_mat(tt,:,:)),90);
    VSD_sIdFoF_rotated(tt,:,:) = temp_rotation;
end

grid_x = [19 38 57 77 96 115]
grid_y = [8 27 46 66 85 104 124]

VDS_dFoF_griddata = zeros(size(VSD_dFoF_rotated,1),length(grid_y),length(grid_x));
VDS_sIdFoF_griddata = zeros(size(VSD_sIdFoF_rotated,1),length(grid_y),length(grid_x));
% for y = 1:length(grid_y)
%     for x = 1:length(grid_x)
%         VDS_dFoF_griddata(:,y,x) = mean(mean( VSD_dFoF_rotated(:,grid_y(y)-1:grid_y(y)+1,grid_x(x)-1:grid_x(x)+1) ,2),3);
%         VDS_sIdFoF_griddata(:,y,x) = mean(mean( VSD_sIdFoF_rotated(:,grid_y(y)-1:grid_y(y)+1,grid_x(x)-1:grid_x(x)+1) ,2),3);
%     end
% end
for y = 1:length(grid_y)
    for x = 1:length(grid_x)
        VDS_dFoF_griddata(:,y,x) =  VSD_dFoF_rotated(:,grid_y(y):grid_y(y),grid_x(x):grid_x(x)) ;
        VDS_sIdFoF_griddata(:,y,x) =  VSD_sIdFoF_rotated(:,grid_y(y):grid_y(y),grid_x(x):grid_x(x)) ;
    end
end

% define vector of event timestamps
tStart=Digital_timestamps(find(Dig_data>0,1,'first'));
tEnd=Digital_timestamps(find(Dig_data>0,1,'last'));
clear pSegment vsd_t_ms eventStartMs eventWindow
for i = 1:length(BS)
    temp = find( frameCent>BS(i)-500 & frameCent<BE(i)+500 ); % time window for analysis
    pSegment{i} = temp;
    pSegment_int{i} = length(1:f/f1:temp(1)):length(1:f/f1:temp(end));
    vsd_t_ms{i} = frameCent(temp);
    eventStartMs{i} = frameCent(temp(1));
    eventWindow{i} = frameCent(temp(end))-frameCent(temp(1));
end

% single VSD trace / single MEA channel comparison
[Mea_trace Mea_timestamps] = MC.getData(148,tStart,tEnd-tStart);
Vsd_single = VDS_dFoF_griddata(:,4,4);
figure,[ax,p1,p2] = plotyy(tStart+Mea_timestamps',squeeze(Mea_trace),frameCent,Vsd_single)
for line= 1:length(p2), set(p2(line),'Color','r','LineWidth',2.4), end
set(ax(2),'YColor','r')
xlabel(ax(2), 'time, ms')
ylabel(ax(1), 'Voltage, uV')
ylabel(ax(2), 'dFoF','Color','r')
legend('MEA (4,4)','VSD, average pixels in ROI (4,4)')


% activity MEA over physical space
event = 5;
En_ref = En1(7:13,6:11);
subchannels = En_ref(:);
[M M_ms] = MC.getData(subchannels,eventStartMs{event},eventWindow{event});
figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,subchannels,squeeze(M),En_ref,'DrawElectrodeNumbers',0);
title('MEA raw traces, event 5','FontSize',20,'FontWeight','bold')
figure, [hCbar]=IntensityPhysicalSpacePlot(subchannels,range(squeeze(M),2),En_ref,'plotElectrodeNumbers',0,'plotGridLines',0);
% activity VSD over physical space
index_points = 1:6*7;
%
vsd_traces = reshape(VDS_dFoF_griddata(pSegment{event},:,:),length(pSegment{event}),length(grid_y)*length(grid_x));
vsd_traces = vsd_traces(:,index_points)';
channel_vector = flipud(En_ref);
figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,channel_vector(:),vsd_traces,En_ref,'DrawElectrodeNumbers',0,...
 'scaling','std','traceColor',[0.9 0.1 0.1],'LineWidth',1.8);
title('VSD dFoF traces, event 5','FontSize',20,'FontWeight','bold')
%
vsd_sI_traces = reshape(VDS_sIdFoF_griddata(pSegment_int{event},:,:),length(pSegment_int{event}),length(grid_y)*length(grid_x));
vsd_sI_traces = vsd_sI_traces(:,index_points)';
figure, h=axes;[hPlot]=activityTracePhysicalSpacePlot(h,channel_vector(:),vsd_sI_traces,En_ref,'DrawElectrodeNumbers',0,...
 'scaling','std','traceColor',[0.9 0.1 0.1],'LineWidth',1.8);
title('VSD smoothed upsampled dFoF, event 5','FontSize',20,'FontWeight','bold')
figure, [hCbar]=IntensityPhysicalSpacePlot(channel_vector(:),range(vsd_sI_traces,2),En_ref,'plotElectrodeNumbers',0,'plotGridLines',0);

% compare delays from raw traces
test_electrodes = [185,151,183,148,213]
ind = find(ismember(subchannels,test_electrodes));
mea_traces = flip(squeeze(M(ind,:,:)),1);
colmat = colormap('jet');
colors = colmat(1:floor(length(colmat)./size(mea_traces,1)):length(colmat),:);
figure, hold on
for i = 1:size(mea_traces,1)
    ph(i) = plot(eventStartMs{event}+M_ms,mea_traces(i,:),'LineWidth',1.8,'Color',colors(i,:));
end
legend('MEA 1','MEA 2','MEA 3','MEA 4','MEA 5')
xlim([5e4 5.007e4]), ylim([-250 400])
xlabel('time, ms','FontSize',20,'FontWeight','bold'), ylabel('raw MEA trace, uV','FontSize',20,'FontWeight','bold')

rows = [1:5], cols = 4;
vsdtraces = VDS_sIdFoF_griddata(pSegment_int{event},rows,cols);
starttime = vsd_t_ms{event}; starttime = starttime(1);
endtime = vsd_t_ms{event};
timevector = starttime:starttime+length(vsdtraces)-1;
figure, hold on
for i = 1:size(mea_traces,1)
    ph(i) = plot(timevector,vsdtraces(:,i),'LineWidth',1.8,'Color',colors(i,:));
end
legend('VSD pixel 1','VSD pixel 2','VSD pixel 3','VSD pixel 4','VSD pixel 5')
xlim([4.985e4 5.03e4]), ylim([-2.5e-3 14e-3])
xlabel('time, ms','FontSize',20,'FontWeight','bold'), ylabel('dFoF, smoothed & upsampled','FontSize',20,'FontWeight','bold')

% algorithm for MEA delays based on threshold crossing
for i = 1:size(mea_traces,1)
    trace = mea_traces(i,:);
    datapoints = length(trace);
    [val ind_sort] = sort(trace,'descend');
    Rmax = median(val(1:20));
    max_COM = sum(val(1:20).*ind_sort(1:20)) / sum(val(1:20));
    rising_phase = find([1:datapoints]<=max_COM);
    [val1 ind_sort1] = sort(trace(rising_phase));
    R_trough = median(val1(1:20));
    trough_COM = sum(val1(1:20).*ind_sort1(1:20)) / sum(val1(1:20));
    window_crossing = find([1:datapoints]>=trough_COM & [1:datapoints]<=max_COM);
    threshold_crossing_MEA(i) = window_crossing( find(trace(window_crossing)>= R_trough + (Rmax-R_trough)/2 ,1, 'first') );
end
% algorithm for VSD delays based on threshold crossing
for i = 1:size(vsdtraces,2)
    trace = vsdtraces(:,i);
    datapoints = length(trace);
    [trace_peak peakpos] = sort(trace,'descend');
    peaks = median(trace_peak(1:5));
    peak_COM = sum(trace_peak(1:5).*peakpos(1:5)) / sum(trace_peak(1:5));
    temp5 = abs( trace-repmat(peaks./2,size(trace,1),1));
    window_crossing = find([1:datapoints]<peak_COM);
    [~,min_crossing] = min(temp5(window_crossing));
    threshold_crossing_VSD(i) = min_crossing;
end

time_delays_MEA = eventStartMs{event} + M_ms(threshold_crossing_MEA);
time_delays_VSD = timevector(threshold_crossing_VSD);

figure, plot(time_delays_MEA-min(time_delays_MEA),'--','LineWidth',2), hold on
plot(time_delays_VSD-min(time_delays_VSD),'--k','LineWidth',2)
for i = 1:size(mea_traces,1)
   plot( i, time_delays_MEA(i)-min(time_delays_MEA),'o','Color',colors(i,:),'MarkerSize',8,'MarkerFaceColor',colors(i,:))
   plot( i, time_delays_VSD(i)-min(time_delays_VSD),'s','Color',colors(i,:),'MarkerSize',8,'MarkerFaceColor',colors(i,:)) 
end


%% Delays VSD on grid points

grid_x = [19 38 57 77 96 115]
grid_y = [8 27 46 66 85 104 124]
VSD_delays_grid = NaN(size(data_rec.event_HalfRiseTimes,3),length(grid_y),length(grid_x));
vsd_delaymap = NaN(size(data_rec.event_HalfRiseTimes));
for i = 1:size(data_rec.event_HalfRiseTimes,3)
    vsd_delaymap(:,:,i) = imrotate(data_rec.event_HalfRiseTimes(:,:,i),90);
end
for y = 1:length(grid_y)
    for x = 1:length(grid_x)
        VSD_delays_grid(:,y,x) =  vsd_delaymap(grid_y(y),grid_x(x),:) ;
    end
end

figure,
events_selected = [1,3:12];
c = 0;
mask_roi = ismember(flipud(En1(7:13,6:11)),subchannels);
for i = 1:length(events_selected)+1 
    if i<= length(events_selected) 
        image = squeeze(VSD_delays_grid(i,:,:)).*mask_roi;
        min_delay = min(image(find(mask_roi==1)));  
        image = image - min_delay; 
        max_delay = max(image(find(mask_roi==1)));
        image(find(~mask_roi)) = NaN; 
        subplot(3,4,i), imagesc(image), axis square 
        h=colorbar; freezeColors 
        c = c+1;
        titlestring = strcat('event # ',num2str(c),' ; max delay: ',' ',num2str(max_delay),' ms' )
        title(titlestring,'FontSize',16,'FontWeight','bold');
    else
        image =  squeeze(mean(VSD_delays_grid,1)).*mask_roi;
        mindelay = min( image(find(mask_roi==1)) );
        image = image - mindelay; 
        max_delay = max(image(find(mask_roi==1)));
        image(find(~mask_roi)) = NaN; 
        subplot(3,4,i), imagesc(image), axis square
        h=colorbar;  freezeColors 
        titlestring = strcat('mean all events ',' ; max avg delay: ',' ',num2str(round(max_delay)),' ms' )
        title(titlestring,'FontSize',16,'FontWeight','bold','Color','r');
    end
end


% plot over physical space 
F=figure('Position',[620   485   519   383]);
avg_delay_VSD = mean(VSD_delays_grid(:,find(ismember(En1(7:13,6:11),subchannels))),1) - min(mean(VSD_delays_grid(:,find(ismember(En1(7:13,6:11),subchannels))),1));
h=axes; minMaxDelay = [0 max(avg_delay_VSD(:))];
[hCbar]=IntensityPhysicalSpacePlot(subchannels(:),avg_delay_VSD,En1(7:13,6:11),'plotElectrodeNumbers',0,'plotGridLines',0,'plotColorBar',1,'markerSize',50,'h',h,'Ilim',minMaxDelay);

%% clear variables

clearvars -except data_analysis

%% save files

% cell array to save data
data_analysis{exp_id}.VSD_ROI_data = {'raw','dF','dFoF';vsd_ROI,dF_ROI,dFoF_ROI};
data_analysis{exp_id}.VSD_metadata = VSDmeta;
data_analysis{exp_id}.rec_name = 'X42_20150605_SampleB_ROI_MEAch_248';

data_saved = data_marcel;
data_saved = rmfield(data_saved,{'data','dFoF_interpolated'});
rec_dir = 'C:\Users\Superuser\Documents';
filename = 'data_marcel.mat';
save(strcat(rec_dir,filesep,filename),'data_saved','-v7.3')

% save some variables of workspace
file_dir = '\\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X42_20150603\ANDOR\MATLAB\SampleB\ROI_MEA_channel_248';
filename = 'data_slab_B_wavelength.mat';
save(strcat(file_dir,filesep,filename),'data_analysis')

%% figure booster
current_position = get(gcf,'Position');
line_flag = 0;
poster_fig_booster(line_flag, 'font_weight','bold','legend_position','South',...
        'figure_coordinates',current_position,'font_legend',13,'weight_legend','bold')
% remove tick from axes
set(gca,'XTickLabelMode','manual','XTickLabel',[])

%% save figures
 
% edit figures
% edit axes in figure
% figure_handle = gcf;
% set(figure_handle, 'Position', [100, 100, 900, 600]);
% allAxesInFigure = findall(figure_handle,'type','axes');
% allTitlesInFigure = findall(figure_handle,'type','text');
% h_axLabels = get(allAxesInFigure,{'XLabel' 'YLabel'}); 
% set(allAxesInFigure,'FontSize',14,'FontWeight','bold')
% set([h_axLabels{:}],'FontSize',14,'FontWeight','bold')
% set(allTitlesInFigure,'FontSize',14,'FontWeight','bold')
% legend('Location','South')


set(gcf,'PaperPositionMode','auto')
% cd \\storage.laur.corp.brain.mpg.de\Data_2\CoattiAlex\Experiments_2015\X42_20150603\ANDOR\MATLAB\SampleB\Mat_figures
cd \\storage.laur.corp.brain.mpg.de\Personal\CoattiAlex\Presentations\Lab_meetings\2015_11_09
print X42_slabB_rec9_smallLFP_VSDvsMEA  -dpng -r300;