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

vis = 0;
prm.planetrack = 4;

trch = 2;

nfolder = length(folder);
for i = 1 : nfolder
    foldereval1 = [folder{i} '/' 'eval' '/' 'single'];
    [a b c] = mkdir(foldereval1);
    foldereval2 = [folder{i} '/' 'eval' '/' 'multiple'];
    [a b c] = mkdir(foldereval2);
    
    fid = fopen(['stat' '-' folder{i} '.txt'],'wt');
    header{1,1} = folder{i};

    printcell(fid,header,46);
    
    
    nim = 100;
    B0 = [];B1 = [];
    data0 = zeros(0,9);data1 = zeros(0,9);
    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 = 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;
        
        
        % 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 = mean(imtransf(didcell == 1));


        
        % 
        % Find DiD signal etc in segmented cells
        %
        
        clear cells;
        cells.meandid = NaN(ncells,1);
        cells.voldid = NaN(ncells,1);
        cells.ndid = NaN(ncells,1);
        cells.vol = NaN(ncells,1);
        for k = 1 : ncells
            % this cell
            cellhere = cellwat == cellval(k);
            
            % DiD transfer in this cell
            cells.meandid(k,1) = mean(imtransf(cellhere(:)));
            
                                    
            % volume of DiD signals in cell
            a = didbw .* cellhere;
            cells.voldid(k,1) = sum(a(:));
            
            % number of DiD signals
            [faserhere,Lhere] = bwlabeln(a);
            cells.ndid(k,1) = Lhere;
            
            % volume
            cells.vol(k,1) = sum(cellhere(:));            
            
            
        end;%k
        
        % find DID signal in the DiD cell to print to textfile
        reg = didcell == 1;
        did.meandid = mean(imtransf(reg));        
        did.voldid = sum(reg(:));
        
        % find row relations, now there is only one DID cell here
        [cells.rownum,temp] = cell2cellcontact(didcell,cellwat);
        

        %
        % 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);
        a = imtransf.*(didbw | didcell);
        testim = scale(a)*(n+1);       
        
        
        for k = 1 : ncells
            cellhere = cellwat == cellval(k);
            testim(cellhere) = cells.rownum(k);
        end;
                
        
        
        % 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;
        
        
        [cells.rownum,indsort] = sort(cells.rownum);        
        cells.ndid = cells.ndid(indsort);
        cells.voldid = cells.voldid(indsort);
        cells.meandid = cells.meandid(indsort);
        cells.vol = cells.vol(indsort);
        
        
        %
        % For printing
        %
        
        data = NaN(ncells,9);
        header = {['Cell ' int2str(j)],'DiD-cell-mean','DiD-cell-vol','DiD-max-dist','1st-row-vol','1st-row-mean','1st-row-n','2nd-row-vol','2nd-row-mean','2nd-row-n'};
        
        data(1,1) = did.meandid;
        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)*3) = cells.voldid(ind);
            data(1:n,5+(k-1)*3) = cells.meandid(ind);
            data(1:n,6+(k-1)*3) = cells.ndid(ind);
            N = max(n,N);
        end;
        % cut away NaN's
        data = data(1:N,:);
        
        meanval = nanmean(data,1);
        meanval = [{'Mean'} mat2celldirect(meanval)];
        sumval = nansum(data,1);
        sumval = [{'Sum'} mat2celldirect(sumval)];        
        datanum = data;
        data = mat2celldirect(data);
        
        col1 = NaN(N,1);
        col1 = mat2celldirect(col1);
        
        emptyline = NaN(1,10);
        emptyline = mat2celldirect(emptyline);
        
        
        A = [header ; [col1 data]; meanval; sumval; emptyline];
        
        
        if ndid == 0
            B0 = [B0 ; A];
            data0 = [data0 ; datanum];
        elseif ndid >= 1
            B1 = [B1 ; A];
            data1 = [data1 ; datanum];
        end;
        
        
    end;
       
    header = {NaN,'DiD-cell-mean','DiD-cell-vol','DiD-max-dist','1st-row-vol','1st-row-mean','1st-row-n','2nd-row-vol','2nd-row-mean','2nd-row-n'};
    
    % print to file the non-control data
    printcell(fid,B1,20);
    
    % compute mean and sum values for all
    meanval = nanmean(data1,1);
    meanval = [{'Total-mean'} mat2celldirect(meanval)];
    sumval = nansum(data1,1);
    sumval = [{'Total-sum'} mat2celldirect(sumval)];        
    A = [header ; meanval; sumval ; emptyline];
    % print to file
    printcell(fid,A,20);

   
    % print to file the control data
    printcell(fid,B0,20);
           
    
    % compute mean and sum values for all
    meanval = nanmean(data0,1);
    meanval = [{'Total-mean'} mat2celldirect(meanval)];
    sumval = nansum(data0,1);
    sumval = [{'Total-sum'} mat2celldirect(sumval)];  
    save test
    A = [header ; meanval; sumval; emptyline];
    % print to file
    printcell(fid,A,20);
    
    
end;
fclose(fid);

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

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

load ball2;se = getball(ball,2,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;
          