function [res_idx, stage] = gmesh(conn, emesh, gene, orientation, stagesel, show, img_dir)
% [res_idx, stage] = gmesh(conn, emesh, gene, orientation, stagesel, show, img_dir)
% displays all images for a gene
% returns res_idx = index of files in emesh.files returned, stage = stages
% gene = single gene or cell array of genes
% orientation = e.g. lateral
% stagesel = stage number or 0 for all stages
% show = 0 for no figure, otherwise shows figure
% img_dir = image directory for raw images (displayed when specified)

if nargin < 6,
    fprintf('Missing arguments');
end

% Query restriction for ALL_ANNOTATED
st_qrestr = 'and list_el.list_id = 66782 and m.gene_id = list_el.element_id and m.date <= ''2007-02-31'' ';
% Ordering
st_order = 'order by a.stage';
% Orientation
st_orient = 'and i.dv = ''%s'' ';
st_orientwc = 'and i.dv like ''%s''';


st_info = ['select distinct(i.image_path), a.stage ' ... 
    'from main m, annot a, image i, list_el where ' ...
    'i.annot_id = a.id and a.main_id = m.id and i.image_path not like ''production_images%%'' and ' ...
    'm.%s %s ''%s'' '];

st_infostage = ['select distinct(i.image_path), a.stage ' ... 
    'from main m, annot a, image i, list_el where ' ...
    'i.annot_id = a.id and a.main_id = m.id and i.image_path not like ''production_images%%'' and ' ...
    'm.%s %s ''%s'' and a.stage = %d '];


% Assemble query
st_info = [st_info, st_qrestr];
st_infostage = [st_infostage, st_qrestr];
if strcmp(orientation, '%'),
    st_info = [st_info, st_orientwc];
    st_infostage = [st_infostage, st_orientwc];
else
    st_info = [st_info, st_orient];
    st_infostage = [st_infostage, st_orient];
end
st_info = [st_info, st_order];
st_infostage = [st_infostage, st_order];


if iscell(gene),
    gene_cell = gene;
else
    gene_cell{1} = gene;
end

res_idx = [];

for gene_no = 1:length(gene_cell),
    % figure out if to query for symbol or cgname
    if strncmp(gene_cell{gene_no}, 'CG', 2),
        qi_gene = 'cgname';
    else
        qi_gene = 'flybase_name';
    end
    % figure out wildcard
    if isempty(strfind(gene_cell{gene_no}, '%')),
        gene_match = '=';
    else
        gene_match = 'like';
    end
    
    % query for single stage or all stages
    if stagesel == 0,
        info = pg_query(conn, 0, st_info, qi_gene, gene_match, gene_cell{gene_no}, orientation)
    else
        info = pg_query(conn, 0, st_infostage, qi_gene, gene_match, gene_cell{gene_no}, stagesel, orientation)
    end

    if strcmp(info{1}, 'No Data'),
        fprintf('No genes found for %s', gene_cell{gene_no});
        continue;
    end
    
    stage = [];
    for i = 1:size(info, 1)
        stage(i) = info{i, 2};
    end

    ct_stage = histc(stage, [min(stage):max(stage)]);

    ph = max(ct_stage);
    pv = length(ct_stage);

    if show > 0, figure; end;
    
    p = 1; pos_info = 1;
    
    for i = 1:pv
        for j = 1:ct_stage(i)
            
            % Extract filename
            fullpath = info{pos_info,1};
            fn = regexp(fullpath, '\/', 'split');
            pos_info = pos_info + 1;
            if length(fn) ~= 2, fprintf('Skipping unparseable %s\n', fullpath); continue; end;
            fn = fn{2};
            
            eidx = strmatch(fn, emesh.files);
            res_idx = [res_idx; eidx];
            % fprintf('%d=%s\n', eidx, info{pos_info,1});
            if show > 0 && ~isempty(eidx),
                subplot(pv, ph, p); displaymesh(emesh.stain(:, eidx), emesh.p_scale, emesh.t);
                title(fn);
            end
            p = p + 1;

        end
        p = p + (ph - ct_stage(i));
    end
    if show > 0, suptitle(gene); end;

    if show > 0 && nargin > 6,

        figure;
        p = 1; pos_info = 1;
        for i = 1:pv
            for j = 1:ct_stage(i)

                filename = info{pos_info,1};
                dir = pg_query(conn,0,pg_insitus('path'),filename);
                if strcmp(dir{1}, 'No Data')
                    error('Path not found');
                end

                fullname = strcat(img_dir,'/',dir{1},'/',filename);
                fprintf('Loading %s ... ', filename);
                im = imread(fullname);
                im2 = imresize(im, 0.5);
                fprintf('done\n');

                eidx = strmatch(info{pos_info,1}, emesh.files);
                subplot(pv, ph, p);
                image(im2); axis off;
                title(info{pos_info,1});
                pos_info = pos_info + 1;
                p = p + 1;
            end
            p = p + (ph - ct_stage(i));
        end
        suptitle(gene);
    end
end