function [ ] = analysetransferanne(folder,trch,segmch) 
%ANALYSETRANSFERANNE Summary of this function goes here
%   Detailed explanation goes here

th = 8000;
for i = 1 : numel(folder)
% for i = 3
    runexp(folder{i},th,3,trch,segmch);
end;

th = 15000;
for i = 1 : numel(folder)
% for i = 3
    runexp(folder{i},th,1,trch,segmch);
end;

%--------------------------------------

function [] = runexp(folder,mincellvol,planestart,trch,segmch)

ntrim = numel(trch);
vis = 0;

[a basesave c] = fileparts(folder);
pathsave = [folder '/' 'summary-' basesave  '-' 'planestart' '-' int2str(planestart) '.txt'];


% mincellvol = 20000;
% mincellvol = 15000;
start = 1;
stop = 100;
field = {'control','noncontrol'};
for j = 1 : numel(field)
    data.(field{j}).imagenum = [];        
    data.(field{j}).cellnum = [];                
    data.(field{j}).totimagenum = [];        
    data.(field{j}).totcellnum = [];             
    data.(field{j}).donorvol = [];
    data.(field{j}).volring = [];        
    data.(field{j}).volringcell = [];        
    for i = 1 : ntrim    
        data.(field{j}).transf(i).vol = [];
        data.(field{j}).transf(i).int = [];
        data.(field{j}).transf(i).volincell = [];
        data.(field{j}).transf(i).intincell = [];    
        data.(field{j}).transf(i).meanim = [];
        data.(field{j}).tottransf(i).vol = [];
        data.(field{j}).tottransf(i).int = [];        
        data.(field{j}).donor(i).int = [];
        
    end;
end;

for i = start : stop
    
    name = [folder '/' 'stack' int2str(i) '.mat'];
    try
        D = load(name);
        msg = ['Loading ' name];
        disp(msg);
    catch
        msg = ['Could not load ' name];
        disp(msg);
        continue;
    end;
    
    D.im = D.im(:,:,planestart:end,:);
    imsegm = D.im(:,:,:,segmch);
    h = D.h;
    voxelvol = prod(h);
    dim = size(imsegm);
    if dim(3) < 8
        msg = ['Too few planes, continuing'];
        disp(msg);
        continue;
    end;
    
    se = strel('disk',25);
    imsegm = imtophat(imsegm,se);

    clear prmin;
    prmin.h = h;
    imsegm = cellsegm.smoothim(imsegm,'gaussian','prm',prmin);
    cellbw = gt(imsegm,20);    
    for j = 1 : dim(3)
        cellbw(:,:,j) = bwareaopen(cellbw(:,:,j),100);
    end;

    tr1 = D.im(:,:,:,trch(1));
    if ntrim == 2
        tr2 = D.im(:,:,:,trch(2));
    end;

    %
    % Find donor cells
    %
    donorbw = tr1 > 100;        
%     donorbw = bwareaopen(donorbw,100);
    load ball2;se = getball(ball,2,1);
    donorbw = imclose(donorbw,se);

    for j = 1 : dim(3)
        donorbw(:,:,j) = imfill(donorbw(:,:,j),'holes');
    end;    
%     donorbw = fillholes(donorbw,100,100000,6);
    donorbw = bwareaopen(donorbw,mincellvol);
        
    % remove from border
    a = donorbw;
    donorbw = imclearborderth(donorbw,0.01);
    rem = a - donorbw;    
    
    % the donor cells
    [faserdon,ndoncells] = bwlabeln(donorbw);
    
%     showall(tr1,donorbw)
    
    if ndoncells == 0
        celltype = 'control';
    else
        celltype = 'noncontrol';
    end;

    % make fake control cells to compare 
    if isequal(celltype,'control')
        c = [dim(1)*0.33 dim(2)*0.33];
        c = round(c);
        p = [];
        for j = 1 : 2
            for k = 1 : 2
                p = [p ; j*c(j) k*c(k)];
            end;
        end;
        a = zeros(dim);
        for j = 1 : size(p,1)
            a(p(j,1),p(j,2)) = 1;
        end;
        b = bwdist(a);
        donorbw = b < 40;
        

        % the donor cells
        [faserdon,ndoncells] = bwlabeln(donorbw);

    end;
    
    load ball2;se = getball(ball,2,1);
    donorbwdil = imdilate(donorbw,se);
    
    % remove donor cell in cellbw
    cellbw(donorbwdil == 1) = 0;        
      
    % transfer rates
    tr1bw = tr1 > 130;    
    if ntrim == 2
        tr2bw = tr2 > 130;        
    end;
    
    % remove transfer in donor cells
    tr1bw(donorbwdil == 1 | rem == 1) = 0;    
    if ntrim == 2
        tr2bw(donorbwdil == 1 | rem == 1) = 0;
    end;
    

    for j = 1 : ndoncells
    
        % this donor cell
        reg = faserdon == j;
                        
        load ball2;se = getball(ball,2,1);
        regdil = imdilate(reg,se);    
    
        % find distance from the donor cell
        d = bwdist(regdil);
        dbw = d < 50;

        % find volume of ring inside cells
        bwringcell = dbw .* cellbw;
        volringcell = voxelvol*nnz(bwringcell);                
        data.(celltype).volringcell = [data.(celltype).volringcell ; volringcell];
        
        % find transfer in channel 1, inside cells  
        vol = tr1bw .* dbw .* cellbw;
        int = vol .* tr1;        
        vol = voxelvol*sum(vol(:));
        int = sum(int(:));
        data.(celltype).transf(1).volincell = [data.(celltype).transf(1).volincell ; vol];
        data.(celltype).transf(1).intincell = [data.(celltype).transf(1).intincell ; int];
        
        if ntrim == 2
            % find transfer 2, inside cells    
            vol = tr2bw .* dbw .* cellbw;        
            int = vol .* tr2;
            vol = voxelvol*sum(vol(:));        
            int = sum(int(:));
            data.(celltype).transf(2).volincell = [data.(celltype).transf(2).volincell ; vol];
            data.(celltype).transf(2).intincell = [data.(celltype).transf(2).intincell ; int];
        end;
        
        if vis
            [j ndoncells]
            int1
            vol1
            if ntrim > 1
                int2
                vol2
                showall(imsegm,cellbw,reg,intim1,intim2)
            else
                showall(imsegm,cellbw,reg,intim1)
            end;
            
        end;
        
        % find volume of total ring
        bwring = dbw;
        volring = voxelvol*nnz(bwring);
        data.(celltype).volring = [data.(celltype).volring ; volring];
        
        % find transfer in channel 1       
        vol = tr1bw .* dbw;        
        int = vol .* tr1;
        vol = voxelvol*sum(vol(:));
        int = sum(int(:));
        data.(celltype).transf(1).vol = [data.(celltype).transf(1).vol ; vol];
        data.(celltype).transf(1).int = [data.(celltype).transf(1).int ; int];
        
        if ntrim == 2
            % find transfer 2    
            vol = tr2bw .* dbw;        
            int = vol .* tr2;
            vol = voxelvol*sum(vol(:));        
            int = sum(int(:));            
            data.(celltype).transf(2).vol = [data.(celltype).transf(2).vol ; vol];
            data.(celltype).transf(2).int = [data.(celltype).transf(2).int ; int];
        end;
        
        % labeling
        data.(celltype).imagenum = [data.(celltype).imagenum ; i];
        data.(celltype).cellnum = [data.(celltype).cellnum ; j];
                
        % donor cell volume
        val = nnz(reg)*voxelvol;
        data.(celltype).donorvol = [data.(celltype).donorvol ; val];
        
        val = tr1(reg == 1);
        val = sum(val(:));
        data.(celltype).donor(1).int = [data.(celltype).donor(1).int ; val];
    
        if ntrim == 2
            val = tr2(reg == 1);
            val = sum(val(:));
            data.(celltype).donor(2).int = [data.(celltype).donor(2).int ; val];
        end;

    end;
    

    % total transfer 1
    vol = tr1bw;
    int = tr1bw .* tr1;
    vol = round(voxelvol*sum(vol(:)));    
    int = sum(int(:));
    data.(celltype).tottransf(1).vol = [data.(celltype).tottransf(1).vol ; vol];
    data.(celltype).tottransf(1).int = [data.(celltype).tottransf(1).int ; int];
        
   
    if ntrim == 2
        % total transfer 2
        vol = tr2bw;
        int = tr2bw .* tr2;
        vol = round(voxelvol*sum(vol(:)));    
        int = sum(int(:));
        data.(celltype).tottransf(2).vol = [data.(celltype).tottransf(2).vol ; vol];
        data.(celltype).tottransf(2).int = [data.(celltype).tottransf(2).int ; int];
    end;

    % mean image value
    data.(celltype).transf(1).meanim = [data.(celltype).transf(1).meanim ; mean(tr1(:))];
    if ntrim == 2
        data.(celltype).transf(2).meanim = [data.(celltype).transf(2).meanim ; mean(tr2(:))];
    end;
    
    % image number
    data.(celltype).totimagenum = [data.(celltype).totimagenum ; i];    
    

    filesave = fullfile(folder,['stack' int2str(i) '-imsegm.tiff']);
    imwritemulttif(filesave,imsegm);

    filesave = fullfile(folder,['stack' int2str(i) '-tr1bw.tiff']);
    imwritemulttif(filesave,imsegm);

    filesave = fullfile(folder,['stack' int2str(i) '-tr2bw.tiff']);
    imwritemulttif(filesave,imsegm);

%     save test
%     pause
%     if vis
%         showall(imsegm,tr1,tr2,donorbw,tr1bw,tr2bw)   
%     end; 
end;

[a basesave c] = fileparts(folder);
pathsave = [folder '/' 'summary-' basesave  '-planestart' int2str(planestart) '.txt'];

msg = ['Opening file ' pathsave];
disp(msg);
fid = fopen(pathsave,'wt');
p = 20;

fieldv = {'control','noncontrol'};
fields = {'CONTROL-DATA','NON-CONTROL-DATA'};
for i = 1 : 2
    field = fieldv{i};
    clear A;
    A{1} = fields{i};
    printcell(fid,A,p);
    clear A;
    if ntrim == 2
        A = {'Image','cell','donorvol','volringcell','volring'...
            'voltransf1','voltransfincell1','inttransf1','inttransfincell1','donorint1',...
            'voltransf2','voltransfincell2','inttransf2','inttransfincell2','donorint2'};
    else
        A = {'Image','cell','donorvol','volringcell','volring'...
            'voltransf1','voltransfincell1','inttransf1','inttransfincell1','donorint1'};
    end;
    printcell(fid,A,p);

    clear A;
    if ntrim == 2
        val = [data.(field).imagenum data.(field).cellnum data.(field).donorvol data.(field).volringcell data.(field).volring ...
           data.(field).transf(1).vol data.(field).transf(1).volincell data.(field).transf(1).int data.(field).transf(1).intincell data.(field).donor(1).int...
           data.(field).transf(2).vol data.(field).transf(2).volincell data.(field).transf(2).int data.(field).transf(2).intincell data.(field).donor(2).int];
    elseif ntrim == 1
        val = [data.(field).imagenum data.(field).cellnum data.(field).donorvol data.(field).volringcell data.(field).volring ...
           data.(field).transf(1).vol data.(field).transf(1).volincell data.(field).transf(1).int data.(field).transf(1).intincell data.(field).donor(1).int];        
    end
    A = mat2celldirect(val);
    printcell(fid,A,p);
    fprintf(fid,'\n');

    if ntrim == 2
        A = {'Image',...
            'tottransfvol1','tottransfint1','meanim1',...
            'tottransfvol2','tottransfint2','meanim2'};
        printcell(fid,A,p);
        clear A;
        val = [data.(field).totimagenum ...
               data.(field).tottransf(1).vol data.(field).tottransf(1).int  data.(field).transf(1).meanim ...
               data.(field).tottransf(2).vol data.(field).tottransf(2).int  data.(field).transf(2).meanim];
    else
        A = {'Image','tottransfvol1','tottransfint1','meanim1'};
        printcell(fid,A,p);
        clear A;
        val = [data.(field).totimagenum ...
               data.(field).tottransf(1).vol data.(field).tottransf(1).int  data.(field).transf(1).meanim];    
    end;
    A = mat2celldirect(val);
    printcell(fid,A,p);
    fprintf(fid,'\n');

end;

fclose(fid);


