gtAnalysePencils_B1.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assemble B1 pencil scan %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gtbb=analyse_B1(prefix)

dir=pwd;

% spot bounding box [xmin ymin xsize ysize]

% bounding box for each blob
gtbb(1,:)=[890 655 620 505];
gtbb(2,:)=[483 695 690 405];
gtbb(3,:)=[803 701 560 485];
gtbb(4,:)=[516 675 786 460];
gtbb(5,:)=[940 1151 530 484];
gtbb(6,:)=[454 1200 699 450];

% pmo position for each blob
pmo(1,:)=[-89.2 -87.7];
pmo(2,:)=[90.8  92.3];
pmo(3,:)=[-27.1 -25.6];
pmo(4,:)=[152.9 154.4];
pmo(5,:)=[38.4  39.9];
pmo(6,:)=[218.4 219.9];

omegaimages=15;
gridzimages=53;
nblobs=length(gtbb);
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;

    %tmp=zeros(rows,cols);
    stack=zeros(rows,cols,gridzimages*omegaimages);
    vol=zeros(rows,cols,gridzimages);

    % collect the omega stepping images
    %%%%%%%%%%
    % step 1 %
    %%%%%%%%%%
    for slice=1:gridzimages
        disp(slice)
        for i=1:omegaimages
            imagendx=(omegaimages*(slice-1))+i+omegaimages*gridzimages*(spot-1);
            %imagendx=(15*(slice-1))+i;
            stackndx=(omegaimages*(slice-1))+i;
            stack(:,:,stackndx)=edf_read(sprintf('B1pencilscan_%04d.edf', 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:end, omegaimages:omegaimages: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);
    for i=1:gridzimages
        rndx=(1:omegaimages)+((i-1)*omegaimages);
        minivol=filtstack(:,:,rndx);
        minivol(minivol<0)=0;
        minivolmask=imdilate(minivol>5, ones(5,5,5));
        vol(:,:,i)=sum(minivol.*minivolmask, 3);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 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*gridzimages*nblobs)+(spot-1)*53+((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;
    edf_write(volf,sprintf('spot%d_vol2.edf',spot))

    clear vol vol2 volf filttmlvol minivol filtstack stack

end % for j


end % end function