function handles=plot_genes_on_slices(filt_neurons,rolonies,gene_list,plot_individual)
% plot genes from gene_list on slices. Returns figure handles to allow
% saving outside of this function.
% filt_neurons and rolonies: standard output from bARseq
% gene_list: a cell array containing a list of gene names.
% plot individual: plot one gene per image orall genes on the same image. default plot all genes together.

if ~exist('plot_individual','var') 
    plot_individual=0;
end


% if a char array is provided as genes, convert to cells
if ischar(gene_list)
    gene_list={gene_list};
end
%gene_list should be horizontal. If not flip it.
if size(gene_list,1)>1
    gene_list=gene_list';
end

[~,gene_idx]=ismember(gene_list,rolonies.genes(:,1));

uniq_slices=unique(rolonies.slice);
h=round(sqrt(numel(uniq_slices)));
w=ceil(numel(uniq_slices)/h);

handles={};

if plot_individual
    cmap=colormap('lines');
    cmap=cmap(1:7,:);
    for i=1:numel(gene_idx)
        handles{i}=figure('Position',[50 50 1600 800]);
        for n=1:numel(uniq_slices)
            subplot(h,w,n);
            neurons_in_slice=filt_neurons.slice==uniq_slices(n);
            hold on;
            % plot neuron background
            scatter(filt_neurons.pos(neurons_in_slice,1), ...
                filt_neurons.pos(neurons_in_slice,2), ...
                1, ...
                [0.9 0.9 0.9], ...
                'filled' ...
                );
            % plot genes of interest
            rols_in_slice=rolonies.slice==uniq_slices(n);
            scatter(rolonies.pos(rols_in_slice&rolonies.id==gene_idx(i),1), ...
                rolonies.pos(rols_in_slice&rolonies.id==gene_idx(i),2), ...
                3, ...
                cmap(2,:), ...
                'filled', ...
                'MarkerFaceAlpha',0.3,...
                'MarkerEdgeAlpha',0.3...
                );
            title([gene_list{i},', slice',num2str(uniq_slices(n))]);
            set(gca, ...
                'ydir','reverse' ...
                );
            set(gcf, ...
                'color','w' ...
                );
            xl=xlim;yl=ylim;pbaspect([range(xl),range(yl),1]);
            axis off
        end
    end
else
    cmap=hclrainbow(numel(gene_idx)+1);
    handles{1}=figure('Position',[50 50 1200 800]);
    for n=1:numel(uniq_slices)
        subplot(h,w,n);
        neurons_in_slice=filt_neurons.slice==uniq_slices(n);
        hold on;
        % plot neuron background
        scatter(filt_neurons.pos(neurons_in_slice,1), ...
            filt_neurons.pos(neurons_in_slice,2), ...
            1, ...
            [0.9 0.9 0.9], ...
            'filled' ...
            );
        % plot genes of interest
        rols_in_slice=rolonies.slice==uniq_slices(n);

        for i=1:numel(gene_idx)
            scatter(rolonies.pos(rols_in_slice&rolonies.id==gene_idx(i),1), ...
                rolonies.pos(rols_in_slice&rolonies.id==gene_idx(i),2), ...
                2, ...
                cmap(i,:), ...
                'filled', ...
                'MarkerFaceAlpha',0.2,...
                'MarkerEdgeAlpha',0.2...
                );
        end
        num_col=ceil(numel(gene_list)/15);
        legend([{'somas'},gene_list], ...
            'NumColumns',num_col, ...
            'Location','eastoutside' ...
            );
        title([gene_list{i},', slice',num2str(uniq_slices(n))]);
        set(gca, ...
            'ydir','reverse' ...
            );
        set(gcf, ...
            'color','w' ...
            );
        xl=xlim;yl=ylim;pbaspect([range(xl),range(yl),1]);
        axis off

    end

end

end