function [bw,imsegm,imnucleus,labelim,cmout,nuclei,rim,controlim] = getimageobj(plate,prm)

% prm.lowres.minnumneigh = 15;
% prm.lowres.level = 100;

msg = ['This is ' upper(mfilename) ' finding highres image objects'];
disp(msg);


% take out the image
tic
msg = ['Extracting the images'];
disp(msg);
% [labelim,box] = boundreg(plate.labelim,1,0);
% imsegm = plate.lowresim{prm.segmch}(box(1,1):box(1,2),box(2,1):box(2,2));
% imnucleus = plate.lowresim{prm.nucleusch}(box(1,1):box(1,2),box(2,1):box(2,2));
clear plate:
% labelim = plate.labelim(box(1,1):box(1,2),box(2,1):box(2,2));
% [imnucleus,box] = boundreg(plate.lowresim{prm.nucleusch},1,0);
% [labelim,box] = boundreg(plate.labelim,1,0);
dim = size(imsegm);
toc

% make coordinate image
a = 1:prm.lowres.dim(1):dim(1);
b = 1:prm.lowres.dim(2):dim(2);
na = numel(a);
nb = numel(b);
d = prm.lowres.dim/2-1;

% make centered grid
[cx,cy] = ndgrid(-d(1)-0.5:d(1)+0.5,-d(2)-0.5:d(2)+0.5);
cx = imrotate(cx,-90);
cy = imrotate(cy,-90);
cx = single(cx);
cy = single(cy);
cimx = repmat(cx,na,nb);
cimy = repmat(cy,na,nb);


% find donor cells
tic
msg = ['Finding donor cells'];
disp(msg);
bw = imsegm > prm.detect.leveldid;
se = strel('disk',2);
bw = imclose(bw,se);
% keep only the medium large objects
th = round(prm.detect.smallestcellarea/prm.lowres.voxelvol);
bw = bwareaopen(bw,th);
th = round(prm.detect.largestcellarea/prm.lowres.voxelvol);
bw =logical(bw - bwareaopen(bw,th));
toc
% showall(imsegm,bw)

% find rim
tic
msg = ['Finding rim'];
disp(msg);
rim.bw = imnucleus > prm.detect.levelrim;
th = round(10000/prm.lowres.voxelvol);
rim.bw = bwareaopen(rim.bw,th);
% Keep as many rims as there are wells per line!!
% rim.bw = bwkeep(rim.bw,prm.lowres.nwellsperline,8);
% se = strel('disk',2);
% rim.bw = imclose(rim.bw,se);
[rim.c(:,1),rim.c(:,2)] = find(rim.bw);
rim.nc = size(rim.c,1);
% true voxel size
rim.c = rim.c.* repmat(prm.lowres.h,rim.nc,1);
toc
% showall(imnucleus,rimim)


% find nuclei
tic
msg = ['Finding nuclei'];
disp(msg);
nuclei.bw = imnucleus > prm.detect.levelnuclei;
th = round(prm.detect.smallestnucleiarea/prm.lowres.voxelvol);
nuclei.bw = bwareaopen(nuclei.bw,th);
th = round(prm.detect.largestnucleiarea/prm.lowres.voxelvol);
nuclei.bw = logical(nuclei.bw - bwareaopen(nuclei.bw,th));
toc;
% showall(imnucleus,nuclei.bw)



% find the center of mass of these objects
tic
[faser,L] = bwlabeln(bw);
msg = ['Initially finding ' int2str(L) ' donorcells'];
disp(msg);
cm = NaN(L,2);
for i = 1 : L
    clear c;
    [c(:,1),c(:,2)] = find(faser == i);
    cm(i,:) = mean(c,1);    
end;
toc

tic
% only keep objects that are at a certain distance to all other objects
remove = zeros(L,1);
remhere = zeros(L,1);
for i = 1 : L
    for j = i+1 : L
        d = cm(j,:) - cm(i,:);
        d = d.*prm.lowres.h;
        absd = norm(d);
        if absd < prm.detect.smallestdistcells
            remove(i) = 1;            
            remove(j) = 1;
            remhere(i) = 1;            
            remhere(j) = 1;            
        end;
    end;
end;
msg = ['Cell-to-cell distance removed ' int2str(sum(remhere)) ' objects'];
disp(msg);
toc


tic
% only keep cells that are in a dense collection of cells
d = prm.detect.averagecelldiameter;
remhere = zeros(L,1);
for i = 1 : L
    cmhere = cm(i,:);    
    cmhere = round(cmhere);
    m1 = max(1,cmhere(1)-d);m2 = min(dim(1),cmhere(1)+d);
    n1 = max(1,cmhere(2)-d);n2 = min(dim(2),cmhere(2)+d);
    imhere = nuclei.bw(m1:m2,n1:n2);
    [a,n] = bwlabeln(imhere);
    if n < prm.detect.minnumneigh
        remove(i) = 1;
        remhere(i) = 1;
    end;
%     showall(imheredid,imhere)
end;
msg = ['Density requirement removed ' int2str(sum(remhere)) ' objects'];
disp(msg);
toc

tic
% only keep the objects that have a certain distance to the rim
remhere = zeros(L,1);
for i = 1 : L
    cmhere = cm(i,:);
    d = rim.c - repmat(cmhere.*prm.lowres.h,rim.nc,1);    
    d = d.^2;
    d = sum(d,2);
    d = sqrt(d);   
    if d < prm.detect.smallestdistrim
        remove(i) = 1;            
        remhere(i) = 1;
    end;
end;
msg = ['Cell-to-rim distance removed ' int2str(sum(remhere)) ' objects'];
disp(msg);
toc


% showall(imsegm,faser)
% remove
% showall(faser,bw)
ind = find(remove);
if ~isempty(ind)
    msg = ['Removing objects '];
    disp(msg);
    
    % remove from CM
    cm(ind,:) = [];    
    for i = 1 : numel(ind)            
        % remove from BW image
        bw(faser == ind(i)) = 0;
    end;
end;

% H = figure('Position',[100 100 1000 800]);
% subplot(2,3,1);imagesc(imsegm);colormap(gray);axis image;axis off;
% subplot(2,3,2);imagesc(imnucleus);colormap(gray);axis image;axis off;
% subplot(2,3,3);imagesc(nuclei.bw);colormap(gray);axis image;axis off;
% subplot(2,3,4);imagesc(rim.bw);colormap(gray);axis image;axis off;
% subplot(2,3,5);imagesc(bw);colormap(gray);axis image;axis off;

% Find the centered coordinates of the hits
cm = round(cm);
n = size(cm,1);
cmout.did.label = cm;
cm2 = zeros(n,2);
% showall(imsegm,bw)
for i = 1 : size(cm,1)
    cm2(i,1) = cimx(cm(i,1),cm(i,2));
    cm2(i,2) = cimy(cm(i,1),cm(i,2));
end;
cmout.did.global = cm2;
cmout.did.global(:,1) = cmout.did.global(:,1) + prm.lowres.offset;

%
% Get control images
%
[cmout.control.label,controlim] = controlroi(imsegm,labelim,nuclei,rim,prm);
cm = cmout.control.label;
n = size(cmout.control.label,1);
cm2 = zeros(n,2);
for i = 1 : n
    cm2(i,1) = cimx(cm(i,1),cm(i,2));
    cm2(i,2) = cimy(cm(i,1),cm(i,2));
end;
cmout.control.global = cm2;
cmout.control.global(:,1) = cmout.control.global(:,1) + prm.lowres.offset;

% fix the offset
% for i = 1 : size(cm,1)
%     val = plate.labelim(cm(i,1),cm(i,2));
%     cm(i,1) = cm(i,1) + val;
%     cm(i,2) = cm(i,2) + val;
% end;