function [ output_args ] = 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{1};
    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;
    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);
        
%         aa = cellbw;

        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);
               
%         showall(didbw,imtransf)
%         ndid
%         ndidrem
        
        % 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) = [];
            ncells = numel(val);
            
        
        % remove the CT neighbours of the other DID cells
        elseif ndid > 1            
            temp = zeros(dim);
            
            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);

%         showall(imsegm,didbw,cellbw,aa)
        % 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);        

        
        
        % find DiD in CT cells
        clear cells;
        cells.did = 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;
%         save test
        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;
            
        [cells.rownum,indsort] = sort(cells.rownum);        
        cells.did = cells.did(indsort);
        cells.vol = cells.vol(indsort);
             
        clear B;
        B{1,1} = dataload;        
        B{1,2} = num2str(did.did);
        printcell(fid,B,46);
        A = [cell(ncells,1) cell(ncells,1) mat2celldirect((1:ncells)') mat2celldirect(cells.rownum) mat2celldirect(cells.did)  mat2celldirect(cells.vol)];        
        printcell(fid,A,46);
    end;
    
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);

            
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;
    while flag == 0         

        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 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;
          