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

vis = 0;


trch = 2;

nfolder = length(folder);
for i = 1 : nfolder
    foldereval = [folder{i} '/' 'eval'];
    mkdir(foldereval);
    fid = fopen([foldereval '/' 'stat' '-' folder{i} '.txt'],'wt');
    header{1,1} = folder{1};
    header{2,2} = 'Cell number';
    header{2,3} = 'DiD intensity';
    header{2,4} = 'Row number';
    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;
        
        % 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(:));
        didbw = gt(imtransf,10*meanval);
        load ball3;se = getball(ball,3,3);
        didbw = imopen(didbw,se);
        didbw = imclose(didbw,se);
        didbw = bwareaopenrange(didbw,500,6);
        a = didbw;
        didbw = bwkeep(didbw,1,6);                
        d = a - didbw;
        
        
        % we have at least another cell and must remove the neighbours in
        % the first row
        if ~isempty(find(d,1))
            b = wat;
            b(d == 1) = -1;
            
            neigh = neighreg(b,-1,10);            
            for k = 1 : numel(neigh)
                cellbw(wat == neigh(k)) = 0;
            end;
            
            clear b a;
        end

        % only do if there is a did cell at all
        isdid = ~isempty(find(didbw,1));
        if isdid
            did.did = mean(imtransf(didbw == 1));
            ind = find(didbw);
            clear c;
            [c(:,1) c(:,2) c(:,3)] = ind2sub(size(didbw),ind);
            meancmdid = mean(c,1);
        end;
        
        % cellwat                
        cellwat = cellbw .* wat;
        
        % no did cell in CT cells!        
        cellbw(didbw == 1) = 0;
        
        val = unique(cellwat(:));
        val(val == 0) = [];
        ncells = numel(val);
        
        clear cells;
        
        % find DiD in CT 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(:)));
                        
        end;

        
        
%         show(imtransf,4)
%         show(didbw,3)
%         pause
%         continue;

        % for evaluation
        planetrack = 4;
        testim = double(cellbw) + scale(imtransf);
        testim(1,1,planetrack) = 2;

        cells.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 = wat == 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;
                                                
            dim = size(cellbw);
            cellval = [];
            r = meancmct;
            if isdid
                while 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) = planetrack;
                    ind = sub2ind(size(didbw),c(1),c(2),c(3));
                    cellval = [cellval ; cellwat(ind)];
    
    
                    testim(ind) = 2;                            
                    
                    atdid = didbw(ind);
                    if atdid
                        break;
                    end;        
                end;
            
            
                cellval(cellval == 0) = [];
                cellval(cellval == valhere) = [];
                cellval = unique(cellval);
                cells.rownum(k,1) = 1 + numel(cellval);
            end;
            
            filesave = [foldereval '/' 'im' int2str(j) '-eval' '.tif'];
            imwrite(prepimwrite(testim(:,:,planetrack)),filesave,'tiff');

%             if j >= 4
%             end;
        end;
        vis = 0;
        if vis == 1
            
            show(imsegm(:,:,planetrack),1)
            show(testim(:,:,planetrack),2)
            show(cellwat(:,:,planetrack),3)
            show(imtransf(:,:,planetrack),4)
            pause
        end;
            
        [cells.rownum,indsort] = sort(cells.rownum);        
        
        cells.did = cells.did(indsort);
        
             
        clear B;
        B{1,1} = dataload;
        B{2,1} = ['DiD cell intensity: ' num2str(did.did)];
        printcell(fid,B,46);
        A = [cell(ncells,1) mat2celldirect((1:ncells)') mat2celldirect(cells.did) mat2celldirect(cells.rownum)];        
        printcell(fid,A,46);
    end;
    
end;
fclose(fid);



