function [ ] = analysetransfer( folder )
%ANALYSETRANSFER Summary of this function goes here
%   Detailed explanation goes here


prm.planetrack = 4;

trch = 2;

nfolder = length(folder);
% nfolder = 1;
for i = 1 : nfolder
    foldereval1 = [folder{i} '/' 'eval' '/' 'single'];
    [a b c] = mkdir(foldereval1);
    foldereval2 = [folder{i} '/' 'eval' '/' 'multiple'];
    [a b c] = mkdir(foldereval2);
    
    
    nim = 100;
    B0 = [];B1 = [];
    ndatacol = 11;
    data0 = zeros(0,ndatacol);
    data1 = zeros(0,ndatacol);
    for j = 1 : nim
        dataload = [folder{i} '/' 'im' int2str(j) '-segmaut.mat'];
        try
            
            D = load(dataload);
            msg = ['Loaded ' dataload];
            disp(msg);
        catch            
            continue;
        end;
        cellbw = D.cellbw;
        im = D.im;
        imsegm = im(:,:,:,1);
        wat = double(D.wat);
        dim = size(wat);
        
        % remove from border
        cellbw = imclearborderth(cellbw,0.05);
        
        
        % remove small cells        
        conn = 18;
        cellbw = bwareaopenrange(cellbw,500,conn);
        
        
        % find the DiD cells!        
        imtransf = im(:,:,:,trch);               

                
        % did segmetnation
        didbw = gt(imtransf,100);
        
        % remove small objects to make the did cell image
        didcell = bwareaopen(didbw,5000,conn);
        
        
        % remove the cells from the binary image to do detection of signal
        didbw(didcell == 1) = 0;
        
%         showall(imtransf,didbw,didcell,cellbw)
        
                        
        % find the largest DiD cell to keep that one as the center DiD cell
        a = didcell;
        didcell = bwkeep(didcell,1,conn);
        didrem = a;
        didrem(didcell == 1) = 0;
        [faserdidrem,ndidrem] = bwlabeln(didrem,conn);
               
        
%         % remove small objects
%         didcell = bwareaopenrange(didcell,500,conn);
           
        % cells on WAT        
        cellwat = cellbw .* wat;

        % label DID cells to konw the total number of DiD cells (number of
        % objects large enough to be possible cells)
        [didfaser,ndid] = bwlabeln(didcell,conn);        
                

        % remove the watershed regions around these false DiD cells            
        for k = 1 : ndidrem
            % this DiD cell
            reghere = faserdidrem == k;
            
            % the watershed region of this DiD cell we remove since it
            % shall not be analysed
            val = cellwat(reghere);
            val(val == 0) = [];
            val = mostfrequent(val);
            % it may not exist in watershed image
            if ~isempty(val)
                a = wat == val;
                wat(a) = 0;
                cellwat(a) = 0;
            end;
            
            % find the neighbours of the DiD cell
            [rownum,~] = cell2cellcontact(reghere,cellwat);
            
            % remove these cells from wat and cellwat
            ind = find(rownum == 1);
            for l = 1 : numel(ind)
                r = wat == ind(l);
                wat(r) = 0;
                cellwat(r) = 0;
                cellbw(r) = 0;
            end;
            
        end;% k
        
                                                       
        % update cell indicator
        cellval = unique(cellwat(:));
        cellval(cellval == 0) = [];
        ncells = numel(cellval);


        % we have not detected any cells, just continue
        if ncells == 0
            continue;
        end
        
        
        if ndid == 0
            
            % we add the center CT cell as the DID cell!
            msg = ['No DiD cells, use center cell as DID cell'];
            disp(msg);
                             
            cm = NaN(ncells,3);
            for k = 1 : ncells
                reghere = cellwat == cellval(k);
                clear c;
                [c(:,1) c(:,2) c(:,3)] = ind2sub(dim,find(reghere));
                cm(k,:) = mean(c,1); 
            end;

            cmall = mean(cm,1);
            dist = distrep(cm,cmall);
            [minval,ind] = min(dist);
            reghere = cellwat == cellval(ind);
            didcell = reghere;
            cellwat(reghere) = 0;
            wat(reghere) = 0;
            cellbw(reghere) = 0;
                   
        end;
        
        %
        % Now we have exactly one DiD cell anyway. Remove the corresonding
        % watershed region if there is any
        %
        
        val = cellwat(didcell == 1);
        val(val == 0) = [];
        val = mostfrequent(val); 
        if ~isempty(val)
            a = cellwat == val;
            % remove
            cellwat(a) = 0;
            wat(a) = 0;
        end;
        
        % make sure there is no signal from large did cells in the 
        % segmented cells, we only want the small transferred dots!        
        se = ones(3,3,1);
        reg = imdilate(didcell,se);
        cellwat(reg == 1) = 0;
        cellbw(reg == 1) = 0;
        wat(reg == 1) = 0;
                
        
        % update cell indicator
        cellval = unique(cellwat(:));
        cellval(cellval == 0) = [];
        ncells = numel(cellval);
        
       
        % we may have removed the only cell we had when we turned it into a
        % DiD cell!
        if ncells == 0
            continue;
        end


        % only do if there is a did cell at all        
        did.did = sum(imtransf(didcell == 1));


        
        % 
        % Find DiD signal etc in segmented cells
        %
        
        clear cells;
        cells.voldid = NaN(ncells,1);
        cells.pos = NaN(ncells,1);
        cells.intdid = NaN(ncells,1);
        
        % threshold for a signal to be measured for a cell to be positive
        thpos = 5;
        for k = 1 : ncells
            % this cell
            cellhere = cellwat == cellval(k);
                                       
            % NB if you add a term here also change the sorting below!!            
            
            % volume of DiD signals inside cell
            a = logical(didbw .* cellhere);
            cells.voldid(k,1) = sum(a(:));
                         
            % sum of DiD signal inside cell   
            val = imtransf(a == 1);
            cells.intdid(k,1) = sum(val);

            
            % label cell as positive if it has at least 5 spots
            [faserhere,L] = bwlabeln(a);
            cells.pos(k,1) = ge(L,thpos);
                                    
        end;%k
        
        
        % find DID signal in the DiD cell to print to textfile
        reg = didcell == 1;
        did.intdid = sum(imtransf(reg == 1));        
        did.voldid = sum(reg(:));
        
        % find total transferred DiD signal outside centered DiD cell
        tottr.intdid = sum(imtransf(didbw == 1));
        tottr.voldid = sum(didbw(:));
        
        
        % find row relations, now there is only one DID cell here
        [cells.rownum,temp] = cell2cellcontact(didcell,cellwat);
        
%         [cells.voldid        cells.pos   cells.rownum cellval]
%         showall(cellwat,imtransf,didbw,imsegm)
        
        

        %
        % Find the longest distance in the didcell
        %
        
        maxprojdidcell = maxprojimage(didcell);
        ind = find(maxprojdidcell);
        clear c;
        [c(:,1) c(:,2)] = ind2sub(dim,ind);
        n = numel(ind);
        dist = 0;
        for k = 1 : n
            disthere = bsxfun(@minus,c(k:end,:),c(k,:));
            disthere = sqrt(sum(disthere.^2,2));
            disthere = max(disthere);
            dist = max(dist,disthere);
        end            
        did.maxdist = dist;
        
        
        
        % label the test image with the row number of the cells
        % for evaluation
%         n = max(cells.rownum);
        testim = zeros(dim);
                
        for k = 1 : ncells
            cellhere = cellwat == cellval(k);
            testim(cellhere) = cells.rownum(k);
        end;
        ind = didcell == 1;
        a = scale(imtransf)*2;
        testim(ind) = a(ind);
        testim(didbw) = 0;
%         testim = scale(testim)*(n+1); 
%         [cells.pos cells.rownum]
%         showall(testim)
                
        
        
        % print images to disc for easy evaluation
        filesave = [foldereval2 '/' 'im' int2str(j) '-eval' '.tif'];
        mat2tifdirect(testim,filesave);
        filesave = [foldereval1 '/' 'im' int2str(j) '-eval-oneplane' '.tif'];
        imwrite(prepimwrite(testim(:,:,prm.planetrack)),filesave,'tiff');

        filesave = [foldereval2 '/' 'im' int2str(j) '-imsegm' '.tif'];
        mat2tifdirect(imsegm,filesave);
        filesave = [foldereval1 '/' 'im' int2str(j) '-imsegm-oneplane' '.tif'];
        imwrite(prepimwrite(imsegm(:,:,prm.planetrack)),filesave,'tiff');
        
        
        vis = 0;
        if vis == 1
            
            ndid
            figure(1);
            subplot(2,2,1);imagesc(imsegm(:,:,prm.planetrack));colormap(gray);axis image;axis off;
            subplot(2,2,2);imagesc(imtransf(:,:,prm.planetrack));colormap(gray);axis image;axis off;
            subplot(2,2,3);imagesc(cellwat(:,:,prm.planetrack));colormap(gray);axis image;axis off;
            subplot(2,2,4);imagesc(testim(:,:,prm.planetrack));colormap(gray);axis image;axis off;
            
            pause
        end;
        

        % NB sort the array according row number!!
        [cells.rownum,indsort] = sort(cells.rownum);        
        cells.voldid = cells.voldid(indsort);
        cells.intdid = cells.intdid(indsort);
        cells.pos = cells.pos(indsort);
        

        % store information in struct to print to textfile
        stats.folder(i).name = folder{i};
        stats.folder(i).im(j).num = j;
        stats.folder(i).im(j).did = did;
        stats.folder(i).im(j).cells = cells;
        n = numel(cells.rownum);
        stats.folder(i).im(j).n = n;
        stats.folder(i).im(j).ndid = ndid;
        stats.folder(i).im(j).tottr = tottr;
        stats.folder(i).im(j).ind1 = cells.rownum == 1;
        stats.folder(i).im(j).ind2 = cells.rownum == 2;
        stats.folder(i).im(j).ind = (1:n)';
        
        %
        % For printing
        %
        
        data = NaN(ncells,ndatacol);
        header1 = {'DiD-cell-int','DiD-cell-vol','DiD-max-dist','1st-row-vol','1st-row-int',['1st-row-pos(>' int2str(thpos) ')'],'2nd-row-vol','2nd-row-int',['2nd-row-pos(>' int2str(thpos) ')'],'tot-transf-vol','tot-transf-int'};
        header = [['Cell ' int2str(j)] header1];
        

        data(1,1) = did.intdid;
        data(1,2) = did.voldid;
        data(1,3) = did.maxdist;
        nrow = 2;
        N = 0;
        for k = 1 : nrow

            ind = find(cells.rownum == k);
            n = numel(ind);
            data(1:n,4+(k-1)*2) = cells.voldid(ind);
            data(1:n,5+(k-1)*2) = cells.intdid(ind);
            data(1:n,6+(k-1)*2) = cells.pos(ind);
            N = max(n,N);
        end;
        data(1,6+(k-1)*2 + 1) = tottr.voldid;
        data(1,6+(k-1)*2 + 2) = tottr.intdid;
        
        % cut away NaN's
        data = data(1:N,:);
        
        sumval = nansum(data,1);
        sumval = [{'Sum'} mat2celldirect(sumval)];        
        datanum = data;
        data = mat2celldirect(data);
        
        col1 = NaN(N,1);
        col1 = mat2celldirect(col1);
        
        emptyline = NaN(1,ndatacol+1);
        emptyline = mat2celldirect(emptyline);
        
        % combine the data to print to file 
        A = [header ; [col1 data]; sumval; emptyline];
        
        
        if ndid == 0
            B0 = [B0 ; A];
            data0 = [data0 ; datanum];
%             summary0 = [summary0 ; nansum(datanum,1) sum(cells.rownum == 1) sum(cells.rownum == 2) ...
%                 nansum(datanum(:,4))+nansum(datanum(:,6))  ...
%                 nansum(datanum(:,5))+nansum(datanum(:,7)) nansum(datanum(:,8))+nansum(datanum(:,9))];             
             

        elseif ndid >= 1
            B1 = [B1 ; A];
            data1 = [data1 ; datanum];
%             summary1 = [summary1 ; nansum(datanum,1) sum(cells.rownum == 1) sum(cells.rownum == 2) ...
%                 nansum(datanum(:,4))+nansum(datanum(:,6))  ...
%                 nansum(datanum(:,5))+nansum(datanum(:,7)) nansum(datanum(:,8))+nansum(datanum(:,9))];             
            
        end;
        
        
    end;
     
%     % slight change in header
%     header = ['Summary' header1];
%     
    %
    % print to file the non-control data
    %
    

    fid = fopen(['stat' '-' folder{i} '.txt'],'wt');
    header{1,1} = folder{i};
    printcell(fid,header,46);
    fprintf(fid,'\n');

    printcell(fid,B1,20);

    % summary
    sumval = nansum(data1,1);
    sumval = [{'Total-sum'} mat2celldirect(sumval)];        
    A = [header ; sumval ; emptyline];
    
    % print to file
    printcell(fid,A,20);

   
    %
    % print to file the control data
    %    
    
    printcell(fid,B0,20);
    
    % summary
    sumval = nansum(data0,1);
    sumval = [{'Total-sum'} mat2celldirect(sumval)];        
    A = [header ; sumval; emptyline];
    
    % print to file
    printcell(fid,A,20);
    
end;

fclose(fid);

printsummary(stats);
save('summary.mat','stats');

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

function [] = printsummary(stats)

fid = fopen('summary.txt','wt');
p = 30;

A{1}  =['Summary from ' mfilename];
printcell(fid,A,p);
fprintf(fid,'\n');
description{1,1} = 'DiD-cell-int: The total DiD signal inside the DiD cell';
description{2,1} = 'DiD-cell-vol: Volume of DiD cell';
description{3,1} = 'DiD-max-dist: Maximum distance within the DiD cell';
description{4,1} = '1st-row-vol: Number of pixels in the first row cells above the threshold';
description{5,1} = '1st-row-int: Total intensity in the first row cells above the threshold';
description{6,1} = '1st-row-pos: Number of positive cells (with number of spots above a threshold)';
description{7,1} = 'tot-transf-vol: Total volume transferred DiD from donor cell';
description{8,1} = 'tot-transf-int: Total intensities transferred DiD from donor cell';

printcell(fid,description,p);
fprintf(fid,'\n');
varname = {'DiD-cell-int','DiD-cell-vol','DiD-max-dist','1st-row-vol','1st-row-int','1st-row-pos(>5)','1st-row-n-cells','2nd-row-vol','2nd-row-int','2nd-row-pos(>5)','2nd-row-n-cell','1and2-row-vol','1and2-row-pos(>5)','1and2-row-n-cells','tot-transf-vol','tot-transf-int'};
% nvar = numel(varname);

header = [{'Stack'},varname];
nfolder = numel(stats.folder);
for i = 1 : nfolder

    meanvalnc{i} = [];
    meanvalc{i} = [];
    
    folder = stats.folder(i).name;
    clear A;
    A{1} = folder;
    printcell(fid,A,p);

    printcell(fid,header,p);
    nim = numel(stats.folder(i).im);
    B = [];C = [];
    for j = 1 : nim
        imhere = stats.folder(i).im(j);
        if isempty(imhere.n)
            continue;
        end;
        
        ind1 = imhere.ind1;
        ind2 = imhere.ind2;
        ind = imhere.ind;

%         ind1
%         ind2
%         imhere.cells.intdid
%         imhere.cells.voldid
%         pause

        addval = [imhere.num imhere.did.intdid imhere.did.voldid imhere.did.maxdist ...
            sum(imhere.cells.voldid(ind1)) sum(imhere.cells.intdid(ind1)) sum(imhere.cells.pos(ind1)) sum(ind1)...
            sum(imhere.cells.voldid(ind2)) sum(imhere.cells.intdid(ind2)) sum(imhere.cells.pos(ind2)) sum(ind2)...
            sum(imhere.cells.voldid(ind)) sum(imhere.cells.pos(ind)) sum(ind) imhere.tottr.voldid imhere.tottr.intdid];
        
        
       
        
        if imhere.ndid >= 1
            B = [B ; addval];
                
        elseif imhere.ndid == 0
            C = [C ; addval];
            
        end;
        
    end;

    % non-control data
    A = {'Non-control data'};
    printcell(fid,A,p);
    
    A = mat2celldirect(B);
    printcell(fid,A,p);    
    A = mean(B,1);
    meanvalnc{i} = [meanvalnc{i} ; A(:,2:end)];
    
    A = mat2celldirect(A);  
    A{1} = 'Mean';
    printcell(fid,A,p);

    % control data
    A = {'Control data'};
    printcell(fid,A,p);

    A = mat2celldirect(C);
    printcell(fid,A,p);    
    A = mean(C,1);
    meanvalc{i} = [meanvalc{i} ; A(:,2:end)];
    
    A = mat2celldirect(A);  
    A{1} = 'Mean';
    printcell(fid,A,p);

    fprintf(fid,'\n');
    
end;

save test1

%
% Print all the means from above for each folder
%

headerhere{1} = 'Total-means';
printcell(fid,headerhere,p);
clear folder;
for j = 1 : nfolder
    folder{1,j} = stats.folder(j).name;
end;
A = [];
for i = 1 : nfolder
    A = [A ;folder{i} mat2celldirect(meanvalnc{i})];
end;
headerhere = cell(1,size(A,2));
headerhere{1} = 'Non-control-means';
% A = mat2celldirect(A);  
printcell(fid,headerhere,p);
headerhere = [{'Folder'} varname];
printcell(fid,headerhere,p);
printcell(fid,A,p);
fprintf(fid,'\n');

A = [];
for i = 1 : nfolder
    A = [A; folder{i} mat2celldirect(meanvalnc{i})];
end;
headerhere = cell(1,size(A,2));
headerhere{1} = 'Control-means';
A = mat2celldirect(A);  
printcell(fid,headerhere,p);
headerhere = [{'Folder'} varname];
printcell(fid,headerhere,p);
printcell(fid,A,p);
fprintf(fid,'\n');

%
% Print cell information to excel friendly format
%

varnameexcel = {'voldid','intdid','pos','voldid','intdid','pos','voldid','intdid','pos'};
indexcel = {'ind1','ind1','ind1','ind2','ind2','ind2','ind','ind','ind'};
rowexcel = {'row1','row1','row1','row2','row2','row2','all','all','all'};
n = numel(varnameexcel);

save test2

headervarini = cell(1,nfolder);
headernc = cell(1,nfolder);
headernc{1} = 'Non-control-data';
headerc = cell(1,nfolder);
headerc{1} = 'Control-data';
headerfolder = folder;

% print to file EACH VARIABLE at a time, for all folders as column argument
for i = 1 : n    
    % get values
    [b,c] = getvarcells(stats,varnameexcel{i},indexcel{i});
    headervar = headervarini;
    headervar{1} = [varnameexcel{i} '-' rowexcel{i}];
    b = [b; mean(b,1)];
    c = [c; mean(c,1)];
        
    A = [headervar ; headerfolder ; headernc ; mat2celldirect(b) ; headerc ; mat2celldirect(c)];
    printcell(fid,A,p);
    fprintf(fid,'\n');
end;
fprintf(fid,'\n');

%
% Print DiD information to excel friendly format
%

save test3

varnameexcel = {'voldid','intdid','maxdist'};
n = numel(varnameexcel);
% print to file EACH VARIABLE at a time, for all folders as column argument
for i = 1 : n    
    % get values
    [b,c] = getvardid(stats,varnameexcel{i});
    headervar = headervarini;
    headervar{1} = [varnameexcel{i} '-' rowexcel{i}];
    b = [b; mean(b,1)];
    c = [c; mean(c,1)];
    
    A = [headervar ; headerfolder ; headernc ; mat2celldirect(b) ; headerc ; mat2celldirect(c)];
    printcell(fid,A,p);
    fprintf(fid,'\n');
end;

% close file
fclose(fid);

%-----------------------------------------------
function [valB,valC] = getvarcells(stats,varhere,inds)

nfolder = numel(stats.folder);
for j = 1 : nfolder

    nim = numel(stats.folder(j).im);
    valB{j} = [];
    valC{j} = [];
    for k = 1 : nim
        imhere = stats.folder(j).im(k);        
        ind = imhere.(inds);        
        % non-control cells
        if imhere.ndid >= 1
            valB{j} = [valB{j}; sum(imhere.cells.(varhere)(ind))];
        % control cells
        else
            valC{j} = [valC{j}; sum(imhere.cells.(varhere)(ind))];
        end;                
            
    end;
            
end;

% convert to matrix
valB = cell2mat(valB);
valC = cell2mat(valC);

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

function [valB,valC] = getvardid(stats,varhere)

nfolder = numel(stats.folder);
for j = 1 : nfolder

    nim = numel(stats.folder(j).im);
    valB{j} = [];
    valC{j} = [];    
    for k = 1 : nim
        imhere = stats.folder(j).im(k);                
        if imhere.ndid >= 1
            valB{j} = [valB{j}; sum(imhere.did.(varhere))];
        else
            valC{j} = [valC{j}; sum(imhere.did.(varhere))];
        end;                

            
    end;
            
end;
% convert to matrix
valB = cell2mat(valB);
valC = cell2mat(valC);




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

function [rownum,valcellwat] = cell2cellcontact(didbw,cellwat)

load ball3;se = getball(ball,3,1);
dildidbw = imdilate(didbw,se);


valcellwat = unique(cellwat(cellwat > 0));

% find the neighbours
val = cellwat(dildidbw == 1);
val(val == 0) = [];
val = unique(val);
[inters,inda,indb] = intersect(val,valcellwat);


% 1 is direct neighbour and 2 is distant neighbour
% default is 2, distant neighbour
rownum = 2*ones(numel(valcellwat),1);
rownum(indb) = 1;
% index = valcellwat(indb);

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

function [rownum,index,testim] = cell2cell(didbw,cellwat,testim,prm)

dim = size(cellwat);
ind = find(didbw);
[c(:,1) c(:,2) c(:,3)] = ind2sub(dim,ind);
% meancmdid = mean(c,1);
% amke sure that the point is within the cell!!!
mid = round(size(c,1)/2);
meancmdid = c(mid,:);
            
val = unique(cellwat(cellwat > 0));
ncells = numel(val);
rownum = NaN(ncells,1);        
% find how many cells we are crossing on our way to the DiD cell
for k = 1 : ncells
    valhere = val(k);
    cellhere = cellwat == valhere;
    ind = find(cellhere);
    clear c;
    [c(:,1) c(:,2) c(:,3)] = ind2sub(size(didbw),ind);
    meancmct = mean(c,1);

    % slope
    a = meancmct - meancmdid;
    a = a / norm(a);


    % for evaluation
%             testim(cellhere == 1) = 2;
%             ind = didbw == 1;
%             testim(ind) = testim(ind) + 1;

    
    cellval = [];
    r = meancmct;
    flag = 0;
    niter = 0;
    while flag == 0         
        niter = niter + 1;
        r = r - a;

        c = round(r);
        c(1) = max(c(1),1);c(1) = min(c(1),dim(1));
        c(2) = max(c(2),1);c(2) = min(c(2),dim(2));
        c(3) = max(c(3),1);c(3) = min(c(3),dim(3));
        % we need to stay low in the planes to get the true
        % neighbor relations!
%         c(3) = prm.planetrack;
        ind = sub2ind(dim,c(1),c(2),c(3));
        cellval = [cellval ; cellwat(ind)];


        % for visualization in one plane!
%         c(3) = prm.planetrack;
%         indtestim = sub2ind(dim,c(1),c(2),c(3));
        
%         testim(indtestim) = 2;                            
        flag = didbw(ind);
        
%         if k == 10
%             r
%             a
%             k
%         show(testim(:,:,prm.planetrack),2)
%         show(cellwat,3)
%         show(didbw,4)
%         pause
%         end;
        
        if niter > 100000
%             save test
%             pause
            warning('Could not hit cell');
            break;
        end;
        
        if flag == 1
            break;
        end;        
        
    end;

    cellval(cellval == 0) = [];
    cellval = unique(cellval);
    % remove itself
    cellval(cellval == valhere) = [];
    rownum(k,1) = 1 + numel(cellval);

end;

index = val;
          