| seg3d_threshold(xreg,yreg,zsize_int, stk_region, XYm, XZm, YZm, title, gray) |
function [fh_seg3d_threshold, Ldot, Rdot] = seg3d_threshold(xreg,yreg,zsize_int, stk_region, XYm, XZm, YZm, title, gray)
prj_x = YZm;
prj_y = XZm;
prj_z = XYm;
hh.prj_x=nan;
hh.prj_y=nan;
hh.prj_z=nan;
int_lim=65535;
bw_proj_x = [];
bw_proj_y = [];
bw_proj_z = [];
stk_region_int = stk_region;
stkM_region = XYm;
max_int_region=max(stk_region_int(:));
min_int_region=min(stk_region_int(:));
max_edge=max(max(xreg,yreg),zsize_int);
xtmp=round(300*xreg/max_edge);
ytmp=round(300*yreg/max_edge);
ztmp=round(300*zsize_int/max_edge);
fh_seg3d_threshold = figure('position',[100,50,ytmp+ztmp+100,xtmp+ztmp+100]);
set(fh_seg3d_threshold,'Name',title ,...
'NumberTitle','off','MenuBar','none');
handles.xz_prj_axes = axes(...
'Parent',fh_seg3d_threshold,...
'Units','pixels',...
'Position',[50,50,ztmp,xtmp]);
handles.xy_prj_axes = axes(...
'Parent',fh_seg3d_threshold,...
'Units','pixels',...
'Position',[ztmp+60,50,ytmp,xtmp]);
handles.yz_prj_axes = axes(...
'Parent',fh_seg3d_threshold,...
'Units','pixels',...
'Position',[ztmp+60,xtmp+60,ytmp,ztmp]);
handles.slope_axes = axes(...
'Parent',fh_seg3d_threshold,...
'Units','pixels',...
'Position',[50,xtmp+60,ztmp,ztmp]);
handles.thres_next = uicontrol(...
'Parent',fh_seg3d_threshold,...
'Style','pushbutton',...
'String','Next',...
'Position',[ytmp+ztmp-50 20 100 20],...
'Value',0,...
'callback', @thresh_next_Callback);
thresh2d_1=RidlerCalvard(prj_z)*int_lim;
thresh2d_2 = thresh2d_1;
% plot the z-profile curve
prj_y_sort=sort(prj_y,1,'descend');
ave_int_z=mean(prj_y_sort(1:round(xreg/10),:),1);
axes(handles.slope_axes);
plot(ave_int_z);
plot_xlim=[1,zsize_int];
plot_ylim=[min(ave_int_z)-100,max(ave_int_z)+10];
xlim(plot_xlim);
ylim(plot_ylim);
set(gca,'Xtick',[]);
hold on;
% calculate the intensity profile along the z-axis for pixels above
% .old value
thresh2d_ave=(thresh2d_1+thresh2d_2)/2;
profile_z=sum(prj_y>thresh2d_ave,1);
% find the maximum of the intensity profile along the z-axis, and use
% this point as the mid point for the z-range
[zmid_val,zmid]=max(profile_z);
% automatically determing the lower limit of the z-range by looking for
% the 1st frame where the intensity profile value is larger than 5% of
% the maximum value
zflag=0;
zmin=1;
zval=1;
while (zval < zmid && zflag == 0)
if (profile_z(zval) > zmid_val*0.05)
zmin=zval;
zflag=1;
end
zval=zval+1;
end
% move the zmin by minus one to allow some flexibility
zmin=max(zmin-1,1);
% automatically determing the upper limit of the z-range by looking for
% the last frame where the intensity profile value is larger than 5% of
% the maximum value
zflag=0;
zmax=zsize_int;
zval=zsize_int;
while (zval > zmid && zflag == 0)
if (profile_z(zval) > zmid_val*0.05)
zmax=zval;
zflag=1;
end
zval=zval-1;
end
% move the zmax by one to allow some flexibility
zmax=min(zmax+1,zsize_int);
Ldot=[zmin thresh2d_1];
Rdot=[zmax thresh2d_2];
h(1) = line([zmin zmin], plot_ylim, ...
'color' , 'red', ...
'linewidth', 3, ...
'ButtonDownFcn', @startDragFcn);
h(2) = line([zmax zmax], plot_ylim, ...
'color' , 'blue', ...
'linewidth', 3, ...
'ButtonDownFcn', @startDragFcn);
h(3) = line([zmin zmax], [thresh2d_1 thresh2d_2], ...
'color' , 'black', ...
'linewidth', 3, ...
'ButtonDownFcn', @startDragFcn);
h(4) = scatter(zmin,thresh2d_1,100,'yellow','filled','ButtonDownFcn', @startDragFcn);
h(5) = scatter(zmax,thresh2d_2,100,'yellow','filled','ButtonDownFcn', @startDragFcn);
set(fh_seg3d_threshold,'WindowButtonUpFcn', @stopDragFcn);
guidata(fh_seg3d_threshold,handles);
% display the projection image in grayscale
handles.current_stk_region=stkM_region;
guidata(fh_seg3d_threshold,handles);
showfig_region(handles);
% -------------------------------------------------------------------
% Callback functions
function startDragFcn(varargin)
set(fh_seg3d_threshold, 'WindowButtonMotionFcn', @draggingFcn)
end
function draggingFcn(varargin)
handles=guidata(fh_seg3d_threshold);
slope=(Rdot(2)-Ldot(2))/(Rdot(1)-Ldot(1));
pt = get(handles.slope_axes, 'CurrentPoint');
if (gco == h(1))
% drag the left bar
xnew=pt(1,1);
% keep the left bar on the left
if (xnew >= Rdot(1))
xnew=Rdot(1)-1;
end
ynew=Ldot(2)+(xnew-Ldot(1))*slope;
Ldot=[xnew ynew];
elseif (gco == h(2))
% drag the right bar
xnew=pt(1,1);
% keep the right bar on the right
if (xnew <= Ldot(1))
xnew=Ldot(1)+1;
end
ynew=Rdot(2)+(xnew-Rdot(1))*slope;
Rdot=[xnew ynew];
elseif (gco == h(3))
% drag the threshold bar
if (pt(1,1) < Ldot(1))
pt(1,1)=Ldot(1);
elseif (pt(1,1) > Rdot(1))
pt(1,1)=Rdot(1);
end
Ldot(1,2)=pt(1,2)+(Ldot(1)-pt(1,1))*slope;
Rdot(1,2)=pt(1,2)+(Rdot(1)-pt(1,1))*slope;
elseif (gco == h(4))
% drag the left dot
Ldot=pt(1,1:2);
if (Ldot(1) >= Rdot(1))
Ldot(1)=Rdot(1)-1;
end
elseif (gco == h(5))
% drag the right dot
Rdot=pt(1,1:2);
if (Rdot(1) <= Ldot(1))
Rdot(1)=Ldot(1)+1;
end
end
Ldot(1,1)=max(Ldot(1,1),plot_xlim(1)+1);
Ldot(1,1)=min(Ldot(1,1),plot_xlim(2)-1);
Ldot(1,2)=max(Ldot(1,2),plot_ylim(1)+1);
Ldot(1,2)=min(Ldot(1,2),plot_ylim(2)-1);
Ldot=round(Ldot);
Rdot(1,1)=max(Rdot(1,1),plot_xlim(1)+1);
Rdot(1,1)=min(Rdot(1,1),plot_xlim(2)-1);
Rdot(1,2)=max(Rdot(1,2),plot_ylim(1)+1);
Rdot(1,2)=min(Rdot(1,2),plot_ylim(2)-1);
Rdot=round(Rdot);
set(h(1), 'XData', Ldot(1)*[1 1]);
set(h(2), 'XData', Rdot(1)*[1 1]);
set(h(3), 'XData', [Ldot(1) Rdot(1)], 'YData', [Ldot(2) Rdot(2)]);
set(h(4), 'XData', Ldot(1), 'YData', Ldot(2));
set(h(5), 'XData', Rdot(1), 'YData', Rdot(2));
axes(handles.xz_prj_axes);
guidata(fh_seg3d_threshold,handles);
end
function stopDragFcn(varargin)
handles=guidata(fh_seg3d_threshold);
set(fh_seg3d_threshold, 'WindowButtonMotionFcn', '');
% update the outlines
showfig_region(handles);
guidata(fh_seg3d_threshold,handles);
end
function thresh_next_Callback(hObject, eventdata, handles)
handles=guidata(fh_seg3d_threshold);
guidata(fh_seg3d_threshold,handles);
close(gcf);
end
function showfig_region(handles) % displays the image stk
handles=guidata(fh_seg3d_threshold);
axes(handles.xz_prj_axes);
im1=prj_y;
im2=prj_z;
im3=prj_x';
% draw the outlines
slope=(Rdot(2)-Ldot(2))/(Rdot(1)-Ldot(1));
imtmp=uint16(zeros(size(im1)));
for i=Ldot(1):Rdot(1)
imtmp(:,i)=im1(:,i)-(i-Ldot(1))*slope;
end
bw=(im2bw(imtmp,Ldot(2)/int_lim));
bw=bwperim(im2uint16(bw));
bw_proj_y = imfill(bw,'holes');
im1(bw==1)=int_lim;
imtmp=uint16(zeros(ceil(yreg),ceil(xreg),Rdot(1)-Ldot(1)+1));
for i=1:Rdot(1)-Ldot(1)+1
imtmp(:,:,i)=stk_region_int(:,:,i+Ldot(1)-1)-slope*(i-1);
end
imtmp=max(imtmp,[],3);
bw=im2bw(imtmp,Ldot(2)/int_lim);
bw=bwperim(im2uint16(bw));
bw_proj_z = imfill(bw,'holes');
im2(bw==1)=int_lim;
imtmp=uint16(zeros(size(im3)));
for i=Ldot(1):Rdot(1)
imtmp(i,:)=im3(i,:)-(i-Ldot(1))*slope;
end
bw=im2bw(imtmp,Ldot(2)/int_lim);
bw=bwperim(im2uint16(bw));
bw_proj_x = imfill(bw,'holes');
im3(bw==1)=int_lim;
if isempty(get(gca,'Children')) % check whether the figure has been created if not do it
if gray == 1
% colormap gray; % switch to the linear grayscale colormap
end
% obtain the minimum and maximum intensity from the max projection,
% and use it to scale the colormap of the stack
max_int_region=max(stk_region_int(:));
min_int_region=min(stk_region_int(:));
clim=[min_int_region max_int_region];
himage=imagesc(im1,clim);
% % display pixel info and range when mouse over
% hpixelinfopanel = impixelinfo(himage);
set(gca,'Xtick',[]); % remove the X and Y tick marks - makes it faster.
set(gca,'Ytick',[]);
imobject = get(gca, 'Children');
hh.prj_y=findobj(imobject,'type','image');
set(gca,'color','black');
set(imobject, 'Erasemode', 'none');
axes(handles.xy_prj_axes);
himage2=imagesc(im2,clim);
% % display pixel info and range when mouse over
% hpixelinfopanel2 = impixelinfo(himage2);
set(gca,'Xtick',[]); % remove the X and Y tick marks - makes it faster.
set(gca,'Ytick',[]);
imobject = get(gca, 'Children');
hh.prj_z=findobj(imobject,'type','image');
set(gca,'color','black');
set(imobject, 'Erasemode', 'none');
axes(handles.yz_prj_axes);
himage3=imagesc(im3,clim);
% % display pixel info and range when mouse over
% hpixelinfopanel3 = impixelinfo(himage3);
set(gca,'Xtick',[]); % remove the X and Y tick marks - makes it faster.
set(gca,'Ytick',[]);
imobject = get(gca, 'Children');
hh.prj_x=findobj(imobject,'type','image');
set(gca,'color','black');
set(imobject, 'Erasemode', 'none');
else
set(hh.prj_y,'CData',im1);
set(hh.prj_z,'CData',im2);
set(hh.prj_x,'CData',im3);
end
% set(handles.bw_proj_x,'CData', bw_proj_x);
% set(handles.bw_proj_y,'CData', bw_proj_y);
% set(handles.bw_proj_z,'CData', bw_proj_z);
guidata(fh_seg3d_threshold, handles);
end
%make a loop as a hack to keep the pointer from reaching the end such
%that I am able to update the variables
A = 1;
while A == 1
pause(.25);
if findobj('name',title) ~= 0
else
A = 0;
end
end
end
function level = RidlerCalvard(im)
% convert 2D image to 1D vector
im=im2double(im(:));
if max(im) == min(im)
level = im(1);
elseif isempty(im)
%%% im will be empty if the entire image is cropped away by the
%%% CropMask. I am not sure whether it is better to then set the level
%%% to 0 or 1. Setting the level to empty causes problems downstream.
%%% Presumably setting the level to 1 will not cause major problems
%%% because the other blocks will average it out as we get closer to
%%% real objects?
level = 1;
else
%%% We want to limit the dynamic range of the image to 256. Otherwise,
%%% an image with almost all values near zero can give a bad result.
MinVal = max(im)/256;
im(im<MinVal) = MinVal;
im = log(im);
MinVal = min(im);
MaxVal = max(im);
im = (im - MinVal)/(MaxVal - MinVal);
PreThresh = 0;
%%% This method needs an initial value to start iterating. Using
%%% graythresh (Otsu's method) is probably not the best, because the
%%% Ridler Calvard threshold ends up being too close to this one and in
%%% most cases has the same exact value.
NewThresh = graythresh(im);
delta = 0.00001;
while abs(PreThresh - NewThresh)>delta
PreThresh = NewThresh;
Mean1 = mean(im(im<PreThresh));
Mean2 = mean(im(im>=PreThresh));
NewThresh = mean([Mean1,Mean2]);
end
level = exp(MinVal + (MaxVal-MinVal)*NewThresh);
end
end
|
|
