function [] = definedonorcell(folder,trch)

msg = ['This is ' upper(mfilename) ' defining donor cells'];
disp(msg);

conn = 18;
nfolder = length(folder);
start = 1;
for i = start : nfolder
        
    nim = 100;
    startim = 1;
    for j = startim : nim
        
        % path to load
        pathload = [folder{i} '/' 'stack' int2str(j) '-segm.mat'];            
        try            
            D = load(pathload);
            msg = ['Loaded ' pathload];
            disp(msg);
        catch         
            msg = ['Could not load ' pathload ', continuing'];
            disp(msg);
            continue;
        end;
        
        im = D.im;
        voxelvol = prod(D.info.prm.h);
        dim = size(im);
        dim3 = dim(1:3);
        cellbw = D.cellbw;
        wat = double(D.wat);

        % segmentation image is first channel
        ntrch = numel(trch);
        imsegm = cell(ntrch,1);
        for k = 1 : ntrch
            imsegm{k} = im(:,:,:,trch(k));
        end;

        
        % segmentation image always thresholded at 100
        bwsegm = cell(ntrch,1);
        for k = 1 : ntrch
            bwsegm{k} = gt(imsegm{k},100);
        end;
        
        % remove small objects
        donorcell = cell(ntrch,1);
        for k = 1 : ntrch
            donorcell{k} = bwareaopen(bwsegm{k},5000,conn);        
            load ball2;se = getball(ball,2,1);
            donorcell{k} = imclose(donorcell{k},se);
            donorcell{k} = logical(donorcell{k} == 1);
            for l = 1 : dim(3)
                donorcell{k}(:,:,l) = imfill(donorcell{k}(:,:,l),'holes');
            end;
        end;
        
        % union of donor cells in the segmentation channels
        donorcellunion = zeros(dim3);
        for k = 1 : ntrch
            donorcellunion = donorcellunion | donorcell{k};
        end;

        % find the largest object to keep that one
        for k = 1 : ntrch
            donorcell{k} = bwkeep(donorcell{k},1,conn);
        end;
        donorcellunion = bwkeep(donorcellunion,1,conn);
        
        % watershed cells where backround is zero
        cellwat = cellbw .* wat;         
        
%         showall(imsegm{1},imsegm{2},donorcellunion)

        %
        % Add the corresponding watershed region as donor cell if its not
        % too big
        %                
        val = cellwat(donorcellunion == 1);
        val(val == 0) = [];
        val = mostfrequent(val);
        % if there is a watershed region here
        donorcellwat = donorcellunion;
        if ~isempty(val)
            donorcellwat = cellwat == val;
            % should not be too big, then its something wrong
            if nnz(donorcellwat) < 2*nnz(donorcellunion)
                donorcellunion = donorcellunion | donorcellwat;
                for k = 1 : ntrch
                    donorcell{k} = donorcell{k} | donorcellwat;
                end;
            end;
        end;
        
        % If we have a control there will be no donor cell, and we will use
        % a watershed region as donor cell (the middle one)
        L = zeros(ntrch,1);        
        [~,ndonorcell] = bwlabeln(donorcellunion,conn);
        
        % number of watershed cells
        cellval = unique(cellwat(cellwat > 0));
        ncells = numel(cellval);
        
        % there is no donor cell!! (control) We take the center watershed
        % cell as donor cell as control.
        if ndonorcell == 0 && ncells > 0
            
            % we add the center cell as the DID cell!
            msg = ['No donor cell, use center watershed cell as donor 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);
            [~,ind] = min(dist);
            reghere = cellwat == cellval(ind);
            donorcellunion = reghere;
            donorcellwat = reghere;
            
        end;


%         % remove donorcell from the watershed definition
%         cellbw(donorcellunion) = 0;
%         wat(donorcellunion) = 0;

        % 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);
        donorcellunion = imdilate(donorcellunion,se);
        for k = 1 : ntrch
            donorcell{k} = imdilate(donorcell{k},se);
        end;
        
            
        % make logical
        donorcellunion = logical(donorcellunion);
        donorcellwat = logical(donorcellwat);
        for k = 1 : ntrch
            donorcell{k} = logical(donorcell{k});
        end;

%         showall(cellbw,imsegm{1},donorcellunion,donorcellwat)
        
        % save the donor cell definition
        pathsave = [folder{i} '/' 'stack' int2str(j) '-donorcell.mat'];
        msg = ['Saving ' pathsave];
        disp(msg);
        save(pathsave,'donorcellunion','donorcell','donorcellwat');
        
    end;
    
end;

          