ANDOR_VSD_rise_delay.m

function files = ANDOR_VSD_rise_delay(files,exp_id,rep_id,p_r,p_c,index)

%% visualize delays in VSD signals 
plot_option = 0;

for ii = exp_id
    for jj = rep_id
        
        if isempty(index)
           % valid for ANDOR recordings
           index{2} = 1:files(ii,jj).metadata.height;
           index{1} = 1:files(ii,jj).metadata.width;
           index{3} = 1:files(ii,jj).metadata.NumberImages;
           mode_image = 1; %full image, all pixels included
           single_pixel_analysis = 0;
        elseif numel(index) == 1
            single_pixel_analysis = 1;
        elseif isempty(index{3})
            index{3} = 1:files(ii,jj).metadata.NumberImages;
            single_pixel_analysis = 0;
        end
        
        if single_pixel_analysis
            
            roi_id = index{1};
            p_r = 1; p_c = 1; step_r = 1; step_c = 1;
            size_image = [files(ii,jj).metadata.width,files(ii,jj).metadata.height];
            pxl_id = reshape(1:prod(size_image),size_image(1),size_image(2)); 
        
        else
            
            if rem(max(index{1}),p_r) == 0
                step_r = min(index{1}):p_r:max(index{1});
            else
                step_r = min(index{1}):p_r:max(index{1})-p_r;
            end
            if rem(max(index{2}),p_c) == 0
                step_c = min(index{2}):p_c:max(index{2});
            else
                step_c = min(index{2}):p_c:max(index{2})-p_c;
            end
            roi_id = zeros(p_r * p_c * numel(step_r)*numel(step_c) ,1);
            
            if mode_image % for full image
                pxl_id = reshape(1:length(index{1})*length(index{2}),length(index{1}),length(index{2}));   %reshape(1:256*256,256,256);
            else % for subimage / cropped fov
                width_image = files(ii,jj).metadata.width;
                height_image = files(ii,jj).metadata.height;
                pxl_id = reshape(1:height_image*width_image,height_image,width_image);
            end
            
            counter = 0;
            for c = step_c
                for r = step_r
                    counter = counter + 1;
                    tempid = pxl_id(r: r + p_r-1, c: c + p_c-1);
                    roi_id((1:p_r*p_c) + p_r*p_c*(counter-1) ) = tempid(:);
                end
            end
            
        end
        % options for subtracting dark frame ( TO BE IMPLEMENTED! )
%         if darkframe
%             temp = files{ii,jj}.RawData.traces{1,2};
%         else
%             temp = files{ii,jj}.RawData.traces{1,1};
%             DF = mean(temp(:,:,1:100),3);
%             temp = temp - repmat(int16(DF),1,1,size(temp,3));
%         end

        % get the time series for all pixels in rois
        temp = files(ii,jj).data;
        v_record = reshape(temp,size(temp,1)*size(temp,2),size(temp,3));
        roi_traces = v_record(roi_id,:);
        

        % determine average pixel fluorescence of each roi over time
        n_rois = size(roi_traces,1)/(p_r*p_c);
        if single_pixel_analysis
            avg_roi_traces = roi_traces'; %double(roi_traces')
        else
            avg_roi_traces = zeros(size(roi_traces,2),n_rois,'uint16');
            for i = 1:size(roi_traces,1)/(p_r*p_c)
               indexes = (1:p_r*p_c) + p_r*p_c*(i-1);
               tracetemp = mean(roi_traces(indexes,:),1);
               avg_roi_traces(:,i) =  uint16(tracetemp(:));
            end
        end
        if single_pixel_analysis
            ref_avg_image = zeros(size_image,'uint16');
            ref_avg_image(index{1}) = uint16(mean(avg_roi_traces,1));
        else
            ref_avg_image = uint16(reshape(mean(avg_roi_traces',2),numel(step_r),numel(step_c)));
        end
        
        % median filter (fast implementation compiled in C)
        F0 = zeros(size(avg_roi_traces),'uint16');
        W = 1/files(ii,jj).metadata.KineticCycleTime * 15; % 1500 at 100 Hz 
        if W >= size(temp,3)/2
            W = size(temp,3)/2 - 1;
        end
        tot_length = size(avg_roi_traces,1);
        tic;
        disp('start of F0 computation');
        for i = 1:size(avg_roi_traces,2)
            x = avg_roi_traces(:,i);
            xic = x(1:floor(W/2));
            xfc = x(tot_length - floor(W/2)+1:tot_length);
            F0(:,i) = uint16(fastmedfilt1d(x, W, xic, xfc));
        end
        telapsed = toc
        
        dFoF = (avg_roi_traces - F0) ./ F0;
        
        if ismember(files(ii,jj).metadata.excitation_wavelength, [525,550])
            sign_peaks = -1;
        elseif ismember(files(ii,jj).metadata.excitation_wavelength, [490,500,405,435,460])
            sign_peaks = +1;
        end
            dFoF = sign_peaks .* dFoF;
        
        % fit a polynomial to smooth traces (alternatives!)
%         interval = 866:size(dFoF,1);
        s_dFoF = zeros(size(dFoF),'uint16');
        for i = 1:size(dFoF,2)
            trace = double(dFoF(:,i));
            s1 = smooth(trace,25,'sgolay',3); 
            s_dFoF(:,i) = uint16(s1);
        end

        % up-sample and interpolate the smoothed traces on finer time scale (ms)
        f = 1/files(ii,jj).metadata.KineticCycleTime; % Hz, frequency data aquisition
        f1 = 1000; % 1 frame every 1 ms

        x = [1:size(s_dFoF,1)];
        xq = [x(1):f/f1:x(end)];
       
        dFoF_interp = uint16(interp1(x,single( s_dFoF(:,i)),xq,'spline'));
        dFoF_interp = uint16(interp1(x,single(s_dFoF),xq,'spline'));
    
        % determine time-frame of events from average dFoF 
        trace = mean(s_dFoF,2);
        

        
% % algorithm 2 to find epileptiform events (Mark's function)
        ic=[1;1;1;numel(trace)];
        minSBWnd = 1000; % scale of duration events (ms)
        minSBInterval = 500;% estimanted 
        [BS,BP,BE,BI]=eventLocator(trace,1:length(trace),ic,minSBWnd,minSBInterval);
        
        % keep events with integrated activity above 2 and of duration
        % above 2 sec
        event_id = find(BI > 1 & BE-BS > 800);
        peakfinal = trace(BP(event_id));
        locfinal = BP(event_id);


%         crossing1 = zeros(length(step_r),length(step_c),length(peakfinal));
%         crossing2 = zeros(length(step_r),length(step_c),length(peakfinal));
%         location = NaN(size(dFoF_interp,2),length(peakfinal));  
%         peaks = NaN(size(dFoF_interp,2),length(peakfinal));         %     peaks = cell(size(dFoF_interp,2),2);
%         half_height = NaN(size(dFoF_interp,2),length(peakfinal));  %     half_height = cell(size(dFoF_interp,2),2);
%         H_height_delay = NaN(size(dFoF_interp,2),length(peakfinal));%     peaks{1,1} = 'peak values';
                                                                    %     peaks{1,2} = 'locations peaks';
                                                                    %     half_height{1,1} = 'half height';
                                                                    %     half_height{1,2} = 'delays rising phase';
%         for i = 1:size(dFoF_interp,2)
%             trace = dFoF_interp(:,i);
        if single_pixel_analysis
            crossing3 = zeros(size_image(1),size_image(2),length(peakfinal));
            delay_peak = zeros(size_image(1),size_image(2),length(peakfinal));
        else
            crossing3 = zeros(length(step_r),length(step_c),length(peakfinal));
            delay_peak = zeros(length(step_r),length(step_c),length(peakfinal));
        end
            for e = 1:length(peakfinal)
                
                if locfinal(e) - 629 > 0
                    window = locfinal(e)-629 : locfinal(e)+150;           
        %             window = flipud(window(:));
                else
                    window = 1 : locfinal(j)+150;
        %             window = flipud(window(:));
                end
                trace = dFoF_interp(window,:);
           %% crossing during rising phase of VSD signal based on descen algorithm that localizes half-the-maximum-response
           %  delays to peak are 
                [trace_peak peakpos] = sort(trace(end-249:end,:),'descend');
                peaks = median(trace_peak(1:5,:));
                temp5 = abs( trace-repmat(peaks./2,size(trace,1),1));
                [~,min_crossing3] = min(temp5(300:650,:)); %min(abs(temp5 - 0.5));
                min_crossing3 = min_crossing3+299;
                peak_CoM = sum(trace_peak(1:5,:).* peakpos(1:5,:),1) ./ ...
                    sum(trace_peak(1:5,:),1);
                peak_delays = round(peak_CoM + length(window)-250);    
                if single_pixel_analysis
                    crossing3temp = zeros(size_image);
                    crossing3temp(ismember(pxl_id,roi_id)) = min_crossing3;
                    crossing3(:,:,e) = crossing3temp;
                    peakdelaytemp = zeros(size_image);
                    peakdelaytemp(ismember(pxl_id,roi_id)) = peak_delays;
                    delay_peak(:,:,e) = peakdelaytemp;
                else
                    crossing3(:,:,e) =  reshape(min_crossing3,length(step_r),length(step_c));
                    delay_peak(:,:,e) = reshape(peak_delays,length(step_r),length(step_c));
                end
            end  
%                 trace11 = mean(trace,2);
%                 com = sum([1:400].*trace11(1:400)')/sum([1:400]);
%                 trace12 = trace - repmat(com,size(trace,1),size(trace,2)); 
%                 trace = com - trace12;
                %                 temp2 = trace(window);
%                 b_temp2 = (trace - repmat(min(trace),size(trace,1),1));
%                 nn_temp2 = b_temp2 ./ repmat(max(b_temp2),size(b_temp2,1),1);    
%                 cutoff = mean(nn_temp2(1:350,:)) + 3*std(nn_temp2(1:350,:));
% %                 rois_null = find(cutoff > 0.5);
%                 temp3 = nn_temp2 - repmat(cutoff,size(nn_temp2,1),1);
%                 [~,min_crossing] = min(abs(temp3));
%                 
%                 crossing1(:,:,e) =  reshape(min_crossing,length(step_r),length(step_c));
%                 
%                 temp4 = nn_temp2 - 0.5; %repmat(cutoff,size(nn_temp2,1),1);
%                 [~,min_crossing2] = min(abs(temp4));
%                 
%                 crossing2(:,:,e) =  reshape(min_crossing2,length(step_r),length(step_c));

          
                
%                 rng = range(temp2);
%                 crossing = max(temp2) - rng/2;
%                 [~, min_crossing] = min(abs(temp2 - crossing));
%                 H_height_delay(i,jj) = window(min_crossing);
%                 H_height(i,jj) = temp2(min_crossing);
%                 [ptemp ltemp] = max(temp2);
%                 peaks(i,jj) = ptemp;
%                 location(i,jj) = window(ltemp);
 
%         end

        files(ii,jj).downsampled_image = ref_avg_image;
        files(ii,jj).dFoF = dFoF;
        files(ii,jj).dFoF_interpolated = dFoF_interp;
        files(ii,jj).event_HalfRiseTimes = crossing3;
        files(ii,jj).event_timestamps = [BS(event_id);locfinal;BE(event_id)];
        files(ii,jj).event_PeakTimes = delay_peak;

%     delays = H_height_delay - repmat(min(H_height_delay),size(H_height_delay,1),1);
    
    end
end

%% diagnostics plots
if plot_option
    x = [25:27;57:59;28:30]'; y = [21:23;25:27;5:7]';
    w = max(x(:,1)) - min(x(:,1)); h = max(y(:,1)) - min(y(:,1));
    linear_ind = reshape(1:size(step_r,2)*size(step_c,2),size(step_r,2),size(step_c,2));
    pairs = [];
    for r = 1:size(x,2)
        [p,q] = meshgrid(y(:,r), x(:,r));
        if r == 1
            pairs = [p(:) q(:)];
        else
            pairs = cat(1,pairs,[p(:) q(:)]);
        end
    end
    avg_downsampled_image = reshape(mean(avg_roi_traces,1),size(step_r,2),size(step_c,2));
    idx = sub2ind(size(avg_downsampled_image), pairs(:,1),pairs(:,2));
    bw_roi = zeros(size(avg_downsampled_image));
    bw_roi(idx) = 1;
    figure, 
    subplot(121), imagesc(avg_downsampled_image), colorbar, axis square, hold on
    for r = 1:size(x,2)
        rectangle('Position',[x(1,r) y(1,r) w h], 'LineWidth',2, 'EdgeColor','k');
    end
    title('spatial averaging 3x5 pixels') 
    subplot(122), imagesc(avg_downsampled_image.*bw_roi), axis square, colorbar, 


    figure, hold on
    cols = ['b','g','r'];
    for r = 1:size(x,2)
       indtemp = linear_ind(y(:,r),x(:,r));
       indtemp = indtemp(:);
       plot(trace(:,indtemp),cols(r)),
       m(r) = plot(mean(trace(:,indtemp),2),cols(r),'LineWidth',4);
    end
        intensity_threshold = 2999;
        mask_plotting = reshape(mean(avg_roi_traces,1),45,83);
        mask_plotting(find(mask_plotting<= intensity_threshold)) = 0; 
        mask_plotting = logical(mask_plotting);
        figure,
        for f = 1:5
            image = crossing3(:,:,f).*mask_plotting;
    %         min_scale = prctile( image( find(image(:)>0) ),[5]);
    %         max_scale = max( image( find(image(:)>0) ),[95]);
            subplot(2,3,f), imagesc(image, [400 630]), colorbar, axis tight
            title('half-max delays, rise VSD signal')
        end
        figure,
        for f = 1:5
            image = delay_peak(:,:,f).*mask_plotting;
    %         min_scale = prctile( image( find(image(:)>0) ),[5]);
    %         max_scale = max( image( find(image(:)>0) ),[95]);
            subplot(2,3,f), imagesc(image, [630 780]), colorbar, axis tight
            title('delays peak VSD signal')
        end
end
end