raw_data_plots.m

function [ data_out ] = raw_data_plots( data_rec, npixels, avg_factor, event_peakframes, baseline, pixels_select )
%Raw_data_plots: Plots intensity of pixels of slab, visualize traces for
% different pixels intensities, plot SNR of single pixels with different
% intensities
% Function takes in: 
                    % data_rec: raw data and metadat
                    % npixels (multiple of 2) = number of representative pixels for each
                    % category of light intensity
% Function gives out:
                    % data_out.pixels = selected pixels
                    % data_out.DF = F(t) - F0 traces for pixels
                    % data_out.DFoF = delta F over F
                    % data_out.SNR = signal amplitude divided by RMS
                    % baseline

% npixels = 10;
i = 1; % case of only one recording

%for i = 1:size(data_rec,2)    
    disp(['experiment filename >>  ',data_rec(i).metadata.FileName]);
    data = data_rec(i).data;
    % obtain pixels slab
    avg_image = mean(data_rec(i).data,3); 
    figure, subplot(221), imagesc(avg_image), colorbar, title('image of slab')
    % automatic threshold
%     [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);
%         %   figure, imagesc(bw_mask), colorbar 
%     % 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')
    % high, medium, low intensity pixels
    if isempty(pixels_select)
        [pixels_intensities, pixel_ind] = sort(avg_image(pixels));
        tot_pixels = numel(pixels_intensities);
        pixels_low = pixels(pixel_ind(randi([1 round(tot_pixels/100)*10],npixels,1)));
        pixels_middle = pixels(pixel_ind(randi([round(tot_pixels/100)*45 round(tot_pixels/100)*55],npixels,1)));
        pixels_high = pixels(pixel_ind(randi([round(tot_pixels/100)*90 tot_pixels],npixels,1)));
        pixels_select = [pixels_low, pixels_middle, pixels_high];
    else
        pixels_low = pixels_select(:,1);
        pixels_middle = pixels_select(:,2);
        pixels_high = pixels_select(:,3);
    end
    % compute dF/F
    f = 1/data_rec(i).metadata.KineticCycleTime; % Hz, frequency data aquisition
    window = round(1/data_rec(i).metadata.KineticCycleTime * 15); % 1500 at 100 Hz % median filter (fast implementation compiled in C)
    DF = zeros([size(data,3),size(pixels_select)]);
    DFoF = zeros([size(data,3),size(pixels_select)]);
    tot_length = size(data,3);
    data1 = reshape(data,[size(data,1)*size(data,2),size(data,3)]);
    % sign of dF/F
    if ismember(data_rec(i).metadata.excitation_wavelength, [525,550,460,470,490])
        sign_peaks = -1;
    elseif ismember(data_rec(i).metadata.excitation_wavelength, [500,405,435])
        sign_peaks = +1;
    end
    for ii = 1:size(DF,3)
        for jj = 1:size(DF,2)
            x = double(data1(pixels_select(jj,ii),:));
            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;
            DF(:,jj,ii) = temp1 *sign_peaks; DFoF(:,jj,ii) = temp2 *sign_peaks; 
        end
    end
    clear data1
   % plots
    n_events = numel(event_peakframes);
    for t = 1:n_events
        timewindow(:,t) = round(event_peakframes(t)-5*f):round(event_peakframes(t)+5*f); % frames in time window peak-2 sec:peak+2sec
    end
    if timewindow(1)<0
        timewindow(1) = 1
    elseif timewindow(end) > size(DF,1)
        timewindow(find(timewindow> size(DF,1))) = size(DF,1);
    end
    
    if isempty(baseline)
       baseline = setdiff(floor(window/2):size(traces,1)-floor(window/2),timewindow(:));
    end
    traces = DF; 
    avg_bsl_DF = squeeze(mean(traces(baseline,:,:),1));
    sigma = squeeze(sqrt(mean((traces(baseline,:,:) -  ...
        permute(repmat(avg_bsl_DF,[1,1,length(baseline)]),[3,1,2])).^2,1)));
    amplitudes = squeeze(prctile(traces(timewindow(:),:,:),95,1));
%     noise = squeeze(sqrt(mean(traces(baseline,:,:).^2,1)));
    SNR = amplitudes ./ noise;
%     SNR = (amplitudes - avg_bsl_DF) ./ sigma;
    figure, 
    plot(traces(baseline,:,1)), hold on, plot(mean(traces(baseline,:,1),2),'k','LineWidth',2)
    plot(traces(baseline,:,2)+300), plot(mean(traces(baseline,:,2)+300,2),'k','LineWidth',2)
    plot(traces(baseline,:,3)+600), plot(mean(traces(baseline,:,3)+600,2),'k','LineWidth',2), 
    xlabel('frames baseline'),ylabel('F(t)-F0'),title('fluorescence during baseline for increasing intensities, bottom to top')
    % first figure : compare traces at different intensities
    figure, subplot(221), plot(squeeze(traces(timewindow(:),:,1)),'Color',[0.2 0.5 0.9]), hold on
        plot(mean(traces(timewindow(:),:,1),2),'b','LineWidth',4), 
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('pixel with low intensity','FontSize',18,'FontWeight','bold')
    subplot(222), plot(squeeze(traces(timewindow(:),:,2)),'Color',[0.2 0.9 0.5]), hold on
        plot(mean(traces(timewindow(:),:,2),2),'Color',[0 0.45 0],'LineWidth',4), 
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('pixel with medium intensity','FontSize',18,'FontWeight','bold')
    subplot(223), plot(squeeze(traces(timewindow(:),:,3)),'Color',[0.9 0.4 0.4]); hold on
        plot(mean(traces(timewindow(:),:,3),2),'Color',[0.5 0 0],'LineWidth',4), 
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('pixel with max intensity','FontSize',18,'FontWeight','bold')
    subplot(224), p=plot(squeeze(mean(traces(timewindow(:),:,:),2)),'LineWidth',2), hold on
    set(p(1),'Color','b');set(p(2),'Color',[0 0.45 0]);set(p(3),'Color',[0.5 0 0]);
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('average trace for different intensities','FontSize',18,'FontWeight','bold'),
    legend('low intensity pixels','medium intensity pixels','high intensity pixels')
    figure_handle = gcf;
    allAxesInFigure = findall(figure_handle,'type','axes');
    h_axLabels = get(allAxesInFigure,{'XLabel' 'YLabel'}); 
    set(allAxesInFigure,'FontSize',14,'FontWeight','bold')
    set([h_axLabels{:}],'FontSize',14,'FontWeight','bold')
    % second figure : compare traces with averaging of increasing number of
    % pixels
%     avg_factor = 10;
    n_traces = round(size(traces,2)/avg_factor);
    traces1 = zeros([size(traces,1),n_traces,size(traces,3)]);
    string_legend = cell(1,n_traces);
    for b = 1:n_traces
        string_legend{b}  = strcat(num2str(avg_factor*b),' pixel average');
        for bb = 1:size(traces,3)
            traces1(:,b,bb) = mean(traces(:,1:b*avg_factor,bb),2);
        end
    end
    figure, subplot(221), p=plot(squeeze(traces1(timewindow(:),:,1))); hold on
    colors = distinguishable_colors(round(size(traces,2)/avg_factor));
    for j = 1:length(p)
        set(p(j),'Color',colors(j,:))
    end
    legend(string_legend);
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('avg pixel with low intensity','FontSize',18,'FontWeight','bold')
    subplot(222), p=plot(squeeze(traces1(timewindow(:),:,2))); hold on
    for j = 1:length(p)
        set(p(j),'Color',colors(j,:))
    end
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('avg pixel with medium intensity','FontSize',18,'FontWeight','bold')
    subplot(223),  p=plot(squeeze(traces1(timewindow(:),:,3))); hold on
    for j = 1:length(p)
        set(p(j),'Color',colors(j,:))
    end
    xlabel('frames','FontSize',18,'FontWeight','bold'), ylabel('F(t)-F0','FontSize',18,'FontWeight','bold')
    title('avg pixel with max intensity','FontSize',18,'FontWeight','bold')
    subplot(224), 
    % plot SNR versus intensities
    scatter(log10(avg_image(pixels_low(:))),SNR(:,1),'b'), hold on
    plot(mean(log10(avg_image(pixels_low(:)))),mean(SNR(:,1)),'sb','MarkerFaceColor','b','MarkerSize',10)
    scatter(log10(avg_image(pixels_middle(:))),SNR(:,2),[],[0 0.45 0]),
    plot(mean(log10(avg_image(pixels_middle(:)))),mean(SNR(:,2)),'s','Color',[0 0.45 0],'MarkerFaceColor',[0 0.45 0],'MarkerSize',10)
    scatter(log10(avg_image(pixels_high(:))),SNR(:,3),[],'r'),
    plot(mean(log10(avg_image(pixels_high(:)))),mean(SNR(:,3)),'rs','MarkerFaceColor','r','MarkerSize',10)
    xlabel('log10(avg pixel intensity)','FontSize',18,'FontWeight','bold'), ylabel('SNR','FontSize',18,'FontWeight','bold')
    title('SNR of single pixels with different intensities','FontSize',18,'FontWeight','bold')
    figure_handle = gcf; xLIM = get(gca,'XLim'); set(gca,'xLim',[xLIM(1)-range(xLIM)/2 xLIM(2)+range(xLIM)/2]);
    allAxesInFigure = findall(figure_handle,'type','axes');
    h_axLabels = get(allAxesInFigure,{'XLabel' 'YLabel'}); 
    set(allAxesInFigure,'FontSize',12,'FontWeight','bold')
    set([h_axLabels{:}],'FontSize',12,'FontWeight','bold')   
% end

data_out.pixels = pixels_select;
data_out.DF = DF;
data_out.DFoF = DFoF;
data_out.SNR = SNR;
end