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
    foldereval = [folder{i} '/' 'eval'];
    mkdir(foldereval);
    fid = fopen(['stat' '-' folder{i} '.txt'],'wt');
    header{1,1} = folder{i};
%     header{2,1} = 'Name';
%     header{2,3} = 'Cell number';
%     header{2,5} = 'DiD intensity';
%     header{2,6} = 'Volume';
%     header{2,4} = 'Row number';
%     header{2,2} = ['DiD cell intensity'];
    printcell(fid,header,46);
    
    
    nim = 100;
    A = [];    
    firstcol = [];
    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        
        cellbw = bwareaopenrange(cellbw,500,6);
        
        % find the DiD cells!        
        imtransf = im(:,:,:,trch);        
        meanval = mean(imtransf(:));
        
        highsegm = gt(imtransf,50);
        highsegm = bwareaopen(highsegm,100);
        ratio = sum(highsegm(:)) / prod(dim);
        
        if ratio > 0.001    
            didbw = gt(imtransf,5*meanval);

            % remove small objects
            didbw = bwareaopen(didbw,5000);
            load ball1;se = getball(ball,1,1);
            didbw = imdilate(didbw,se);

            load ball1;se = getball(ball,1,1);
            didbw = imclose(didbw,se);
            
            didbw = bwareaopenrange(didbw,500,6);
        else
            % later we add  CT cell
            didbw = zeros(dim);
        end;
             

        % no did cell in CT cells!        
        cellbw(didbw == 1) = 0;

        % count DID cells
        [faser,ndid] = bwlabeln(didbw);
        
        a = didbw;
        didbw = bwkeep(didbw,1,6);                
        % the did cells to remove neighbors around
        didbwrem = a - didbw;
        [faserdidrem,ndidrem] = bwlabeln(didbwrem);
        
        % cellwat                
        cellwat = cellbw .* wat;
        
        % the cell indicator
        val = unique(cellwat(:));
        val(val == 0) = [];
        ncells = numel(val);
               
        if ncells == 0
            continue;
        end
        

        % we add the center CT cell as the DID cell!
        if ndid == 0
            
            cm = NaN(ncells,3);
            for k = 1 : ncells
                reghere = cellwat == val(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 == val(ind);
            didbw = reghere;
            cellwat(reghere) = 0;
        
            % remove from list
            val(ind) = [];            
        
        % remove the CT neighbours of the other DID cells
        elseif ndid >= 1            
            temp = zeros(dim);
            
%             ndidrem
%             pause 
            for k = 1 : ndidrem
                reghere = faserdidrem == k;
                
                % find row relations
                [rownum,index,temp] = cell2cell(reghere,cellwat,temp,prm);
                
                % first row cells
                val = index(rownum == 1);
                for l = 1 : numel(val)
                    cellwat(cellwat == val(l)) = 0;
                end;
            end;    
                        
        end;
            
                
        % new cellbw
        cellbw = gt(cellwat,0);


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

        % repeat since we now have fewer cells potentially
        val = unique(cellwat(:));
        val(val == 0) = [];
        ncells = numel(val);        


        % we may have removed the only cell we had when we turned it into a
        % DiD cell!
        if ncells == 0
            continue;
        end

        
        % find DiD in CT cells
        clear cells;
        cells.did = NaN(ncells,1);
        cells.vol = NaN(ncells,1);
        for k = 1 : ncells
            % this cell
            cellhere = wat == val(k);
            
            % DiD transfer in this cell
            cells.did(k,1) = mean(imtransf(cellhere(:)));
             
            % volume
            cells.vol(k,1) = sum(cellhere(:));
            
        end;

        
        
        % for evaluation
        testim = double(cellbw) + 2*scale(imtransf);
        
        
        % find row relations, now there is only one DID cell here
        [cells.rownum,index,testim] = cell2cell(didbw,cellwat,testim,prm);
        
        % label the test image with the row number of the cells
        nhere = max(cells.rownum);
        indhere = testim == 2;
        for k = 1 : ncells
            cellhere = cellwat == val(k);
            testim(cellhere) = cells.rownum(k);
        end;
        
        testim(indhere) = nhere + 1;
%         showall(testim,imtransf)        
        
        filesave = [foldereval '/' 'im' int2str(j) '-eval' '.tif'];
        mat2tifdirect(testim,filesave);
        filesave = [foldereval '/' 'im' int2str(j) '-eval-oneplane' '.tif'];
        imwrite(prepimwrite(testim(:,:,prm.planetrack)),filesave,'tiff');
        
        
        vis = 0;
        if vis == 1
            
            show(imsegm(:,:,prm.planetrack),1)
            show(testim(:,:,prm.planetrack),2)
            show(cellwat(:,:,prm.planetrack),3)
            show(imtransf(:,:,prm.planetrack),4)
            pause
        end;
            
%         save test cells
%         cells
%         pause

        [cells.rownum,indsort] = sort(cells.rownum);        
        cells.did = cells.did(indsort);
        cells.vol = cells.vol(indsort);
                     
        
        %
        % For printing
        %
        
        add = NaN(ncells+1,13);
        % did cell
        add(1,1) = did.did;
        
        % loop over each cell row
        nmax = 0;
        th = 5;
        for k = 1 : 3
            ind = find(cells.rownum == k);
            didhere = cells.did(ind);
            volhere = cells.vol(ind);
            n = numel(ind);           
            nmax = max(n,nmax);
            add(1:n,2*k) = didhere;
            add(1:n,2*k + 1) = volhere;
%             % the sum
%             add(1,10 + k) = sum(didhere);
%             % the percentage over th
%             gt(didhere,th);
        end;
             
        % the summing up for percentage
        c = 0;        
        
        for k = 2 : 2 : 6
            c = c + 1;
            didhere = add(:,k);
            n = sum(gt(didhere,0));
            indth = find(didhere > th);
            perth = 100*numel(indth) / n;
%             if isnan(perth)
%                 perth = 0;
%             end;
            add(1,7 + c) = perth;
                                    
        end;
        
        % the summing up        
        for k = 2 : 2 : 6
            c = c + 1;
            didhere = add(:,k);
            val = nansum(didhere);        
            % only NaN's, no values to sum up!
            if val == 0
                val = NaN;
            end;
            add(1,7 + c) = val;
                                    
        end;

        
        % cut away the NaN that we dont need        
        add = add(1:nmax+1,:);
        
        n = size(add,1);
        a = cell(n,1);
        a{1,1} = ['Cell' int2str(j)];
        firstcol = [firstcol; a];
        
        % data for printing
        A = [A ; add];
                               
    end;
    
    

    % find averages
    a = size(A,1);
    for k = 1 : size(A,2)
        datahere = A(:,k);
        val = datahere(~isnan(datahere));
        n = numel(val);
        meanval = mean(val);
        stdval = std(val);
        semval = stdval/sqrt(n);
        summary(1,k) = meanval;
        summary(2,k) = stdval;
        summary(3,k) = semval;
        
    end;
    A = [A; summary];

    header = {'','DiD', '1st-row-did', '1st-row-vol', '2nd-row-did', '2nd-row-vol', '3rd-row-did', '3rd-row-vol', ...
        ['1st-row-perc-gt' int2str(th)], ['2nd-row-perc-gt' int2str(th)], ['3rd-row-perc-gt' int2str(th)], '1st-row-sum', '2nd-row-sum', '3rd-row-sum'};

    firstcol = [firstcol ; {'Average'} ; {'Stdev'} ; {'SEM'}];
        

    A = [header ; firstcol mat2celldirect(A)];
%     % remove zeros for nicer printing
%     for k = 1 : size(A,1)
%         for l = 1 : size(A,2)
%             if (A{k,l} == 0) 
%                 A{k,l} = 'NaN';
%             end;
%             if isempty(A{k,l})
%                 A{k,l} = 'NaN';                
%             end;
%         end;
%     end
        
    printcell(fid,A,20);
    
end;
fclose(fid);

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

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;
          