gtAnalysePencils_B1post.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assemble B1post_ scan %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gtbb=analyse_B1post(prefix,gridyfill)

dir=pwd;

% spot bounding box [xmin ymin xsize ysize]

% bounding box for each blob
% gtbb=[883 1069 697 526]; % B1post_pencils_ scan
% pmo position for each blob
% pmo=[38.4  39.9]; % B1post_pencils_ scan
% omegaimages=1;
% gridzimages=53;

% bounding box for each blob
gtbb(5,:)=[450 850 600 520]; % B1post_linescan_far_ scan # 5

% pmo position for each blob
pmo(5,:)=[37.7 41.1]; % B1post_linescan_far_ scan # 5

omegaimages(1)=25;
omegaimages(2)=25;
omegaimages(3)=25;
omegaimages(4)=25;
omegaimages(5)=34;
omegaimages(6)=34;

gridzimages=55;



nblobs=size(gtbb,1);
maxx=max(gtbb(:,4));
maxy=max(gtbb(:,3));


for j=1:nblobs
    
    % number of rows (v direction)
    rows=gtbb(j,4);
    % number of columns (u direction)
    cols=gtbb(j,3);
    spot=j;
    tot(j)=sum(omegaimages(1:j-1));
    %tmp=zeros(rows,cols);
    stack=zeros(rows,cols,gridzimages*omegaimages(j));
    vol=zeros(rows,cols,gridzimages);

    % collect the omega stepping images
    %%%%%%%%%%
    % step 1 %
    %%%%%%%%%%
    for slice=1:gridzimages
        disp(slice)
        for i=1:omegaimages(j)
            imagendx=(omegaimages(j)*(slice-1))+i+tot(j)*gridzimages
            %imagendx=(15*(slice-1))+i;
            stackndx=(omegaimages(j)*(slice-1))+i
            stack(:,:,stackndx)=edf_read(sprintf('%s/%s%04d.edf', dir, prefix, imagendx), gtbb(j,:));
        end
    end

    % need to de-noise...
    filtstack=zeros(size(stack));
    for i=1:size(stack, 3);
        filtstack(:,:,i)=gtImgMedianFilter2(stack(:,:,i));
    end
    % ...and need a dark
    dark=median(filtstack(:,:,[1:omegaimages(j):end, omegaimages(j):omegaimages(j):end]), 3);
    % subtract from stack
    for i=1:size(stack, 3)
        filtstack(:,:,i)=filtstack(:,:,i)-dark;
    end

    % sum these to make an integrated volume
    %minivol=zeros(rows,cols, omegaimages(j));
    for i=1:gridzimages
        rndx=(1:omegaimages(j))+((i-1)*omegaimages(j));
        minivol=filtstack(:,:,rndx);
        minivol(minivol<0)=0;
        minivolmask=imdilate(minivol>5, ones(5,5,5));
        vol(:,:,i)=sum(minivol.*minivolmask, 3);
    end

    if gridyfill
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % collect the other pencil steps 2-7 %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    vol2=zeros(rows,cols,gridzimages);
    for step=2:7
        disp(step)
        tmpvol=zeros(rows,cols,gridzimages);
        for slice=1:gridzimages
            % which image
            ndx=(omegaimages(j)*gridzimages*nblobs)+(spot-1)*gridzimages+((step-2)*(gridzimages*nblobs))+slice;
            tmpvol(:,:,slice)=edf_read(sprintf('%s/%s%04d.edf', dir, prefix, ndx), gtbb(j,:));
        end
        % need to de-noise...
        filttmpvol=zeros(size(tmpvol));
        for i=1:size(tmpvol, 3);
            filttmpvol(:,:,i)=gtImgMedianFilter2(tmpvol(:,:,i));
        end
        % prepare a dark
        dark=min(filttmpvol(:,:,[1 end]),[], 3);
        b=dark(1:5, :);
        b=sum(b(:));
        % subtract
        for i=1:size(tmpvol, 3);
            a=filttmpvol(1:5,:,i);
            scalefactor=sum(a(:))/b;
            scaleddark=dark*scalefactor;
            filttmpvol(:,:,i)=filttmpvol(:,:,i)-scaleddark;
        end

        vol2=vol2+filttmpvol;
    end


    volf=vol2+vol;

    else
        volf=vol;
    end
    
    edf_write(volf,sprintf('spot%d_vol2.edf',spot))

end % for j


end % end function