function varargout = CellMeasure(varargin)
% CELLMEASURE M-file for CellMeasure.fig
% CELLMEASURE, by itself, creates a new CELLMEASURE or raises the existing
% singleton*.
%
% H = CELLMEASURE returns the handle to a new CELLMEASURE or the handle to
% the existing singleton*.
%
% CELLMEASURE('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in CELLMEASURE.M with the given input arguments.
%
% CELLMEASURE('Property','Value',...) creates a new CELLMEASURE or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before view_stack_OpeningFunction gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to CellMeasure_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help CellMeasure
% Last Modified by GUIDE v2.5 08-Aug-2009 20:50:05
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @CellMeasure_OpeningFcn, ...
'gui_OutputFcn', @CellMeasure_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin & isstr(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before CellMeasure is made visible.
function CellMeasure_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to CellMeasure (see VARARGIN)
% Choose default command line output for CellMeasure
% current axis position
% Update handles structure
handles.path_exist = 0;
handles.loaded = 0; % flag for determining whether stack file is loaded
guidata(hObject, handles);
% UIWAIT makes CellMeasure wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = CellMeasure_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
%varargout{1} = handles.output;
% --- Executes on button press in loadch2.
function loadch1_Callback(hObject, eventdata, handles)
% acctually loadch2 the file
CurrentDir=pwd;
B = 0;
A = B;
data=tiffread;
name = data.filename;
handles.name = name;
stk = data;
set(handles.stackname,'String','name');
%[stk,stacklength]=FNtiffread([data filename]);
%stk = data.imageData;
%find the physical size
prompt={'Enter the pixel-size for xy pixels in um','Enter the pixel-size for z pixels in um'};
name='Input for Peaks function';
numlines=1;
defaultanswer={'0.33','0.20'};
pixelSize=inputdlg(prompt,name,numlines,defaultanswer);
pixelSize = str2double(pixelSize);
sz = size(data);
stacklength = sz(2);
% set maximum plane for kymograph
% reset length of slider
smin = get(handles.plane_slider,'Min');
smax = stacklength;
set(handles.plane_slider,'Max',stacklength,'Sliderstep', [1/(smax-smin),10/(smax-smin)]);
% extract projection and other parameters from stack
X = data(1).width;
Y = data(1).height;
stkM = zeros(Y,X);
stk = struct('data',{});
for i = 1:stacklength
stk(i).data = reshape(data(i).data,X,Y);
t_s(i) = 1;
stkM = max(cat(3,stkM,stk(i).data),[],3);
end
% display projection
showfig(stkM, -1,1);
% save important variables
handles.loaded = 1; % file is loaded
handles.pointed = 0; % points have been identified
handles.tracked = 0; % points have been tracked
handles.stacklength = stacklength; % number of planes in the stack
handles.X = X; % width of image in pixels
handles.Y = Y; % height of image in pixels
handles.stkM = stkM; % maximum projection of stack
handles.stk = stk; % entire stack structure written by metamorph
handles.stk2 = []; % entire stack structure written by metamorph
handles.current_stk = handles.stkM; % current_stk
handles.current_t = -1; % current time, -1 for maximum projection
handles.t = t_s; % time point
handles.shiftx = zeros(length(stk),1);
handles.shifty = zeros(length(stk),1);
handles.regions = 0;
handles.pixelSize = pixelSize;
'channel one loaded: '
handles.name
guidata(hObject, handles);
return
function loadch2_Callback(hObject, eventdata, handles)
data=tiffread;
stk = data;
name = data.filename;
handles.name2 = name;
%[stk,stacklength]=FNtiffread([data filename]);
%stk = data.imageData;
sz = size(data);
stacklength = sz(2);
% set maximum plane for kymograph
% reset length of slider
smin = get(handles.plane_slider,'Min');
smax = stacklength;
% extract projection and other parameters from stack
X = data(1).width;
Y = data(1).height;
stkM = zeros(Y,X);
stk = struct('data',{});
for i = 1:stacklength
stk(i).data = reshape(data(i).data,X,Y);
stkM = max(cat(3,stkM,stk(i).data),[],3);
end
handles.stk2 = stk;
handles.stkM2 = stkM;
guidata(hObject, handles);
'channel two loaded: '
handles.name2
return
function view_projection_Callback(hObject, eventdata, handles)
handles.current_stk = handles.stkM; handles.current_t = -1;
showfig(handles.current_stk, handles.current_t,1);
guidata(hObject,handles);
% --- Executes during object creation, after setting all properties.
function plane_slider_CreateFcn(hObject, eventdata, handles)
% hObject handle to plane_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
% 'usewhitebg' to 0 to use default. See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%set(1, 'DefaultTextColor', 'yellow');
set(hObject,'Value',1);
smin = 1;
smax = 100;
sstep = [1/(smax-smin) 10/(smax-smin)];
set(hObject,'Min',1,'Max',100,'SliderStep',sstep);
% --- Executes on slider movement.
function plane_slider_Callback(hObject, eventdata, handles)
plane = round(get(hObject,'Value'));
set(handles.plane_current,'String', num2str(plane));
handles.current_stk = handles.stk(plane).data; handles.current_t = plane;
showfig(handles.current_stk, handles.current_t,1);
return
% --- Executes on button press in view_single.
function view_single_Callback(hObject, eventdata, handles)
plane = str2num(get(handles.plane_current, 'String'));
handles.current_stk = handles.stk(plane).data; handles.current_t = plane;
showfig(handles.current_stk, handles.current_t,1);
guidata(hObject, handles);
return
function plane_current_Callback(hObject, eventdata, handles)
% hObject handle to plane_current (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of plane_current as text
% str2double(get(hObject,'String')) returns contents of plane_current as a double
% obtain number of stacks, etc...
smin = 1;
smax = handles.stacklength; %round(get(handles.plane_slider,'Max'));
plane = round(str2num(get(hObject, 'String')));
% make sure limits are not breached
if (plane > smax)
plane = smax;
elseif (plane < smin)
plane = smin;
end
% update the slider
set(hObject,'String', plane);
set(handles.plane_slider, 'Value', plane);
% display the figure
handles.current_stk = handles.stk(plane).data; handles.current_t = plane;
showfig(handles.current_stk, handles.current_t,1);
guidata(hObject, handles);
return
%%%%% MY OWN FUNCTIONS %%%%%
function showfig(im,t,fignum) % displays the image stk
figure(fignum);
if isempty(get(gcf,'Children')) % check whether the figure has been created if not do it
colormap gray; % switch to the linear grayscale colormap
% imagesc(im);
imagesc(im);
set(gca,'Xtick',[]); % remove the X and Y tick marks - makes it faster.
set(gca,'Ytick',[]);
imobject = get(gca, 'Children');
set(gca,'color','black');
set(imobject, 'Erasemode', 'none');
else
imobjects = get(gca, 'Children');
imobject = findobj(imobjects,'type','image');
set(imobject, 'CData', im);
end
if (t ~= -1)
title(['zsection = ' num2str(t) '']);
else
title('maximum projection')
end
return
function newroi_Callback(hObject, eventdata, handles)
loc = str2num(get(handles.plane_current, 'String'));
N = length(handles.stk);
figure(1);
[x1,y1,BW,RECT] = imcrop;
regions = handles.regions + 1;
handles.processed(regions) = 0;
handles.roi(regions).RECT = RECT;
handles.regions = regions;
guidata(hObject, handles);
showrectangle(handles.roi);
return
function showregions_Callback(hObject, eventdata, handles)
loc = str2num(get(handles.plane_current, 'String'));
flag = get(hObject,'Value');
if (flag==1)
figure(1);
showrectangle(handles.roi);
else
close(1); figure(1);
showfig(handles.current_stk, handles.current_t,1);
end
return
function showrectangle(rois)
N = length(rois);
figure(1); hold on;
for i = 1:N
rectangle('Position',rois(i).RECT, 'EdgeColor','r');
end
return
function clearroi_Callback(hObject, eventdata, handles)
handles.regions = 0;
handles.roi = [];
guidata(hObject,handles);
% --- Executes on button press in segment.
function segment_Callback(hObject, eventdata, handles)
for roi = 1:handles.regions
%xz plot
clearvars X2;
clearvars X3;
for i = 1:handles.stacklength
X2(i,:,:) = imcrop(handles.stk(i).data(:,:), handles.roi(roi).RECT);
X3(i,:,:) = imcrop(handles.stk2(i).data(:,:), handles.roi(roi).RECT);
end
[zmax,xmax,ymax]=size(X2);
%Run the thresholding program
z_scale = handles.pixelSize(2)/handles.pixelSize(1);
stk_region = permute(X2, [2 3 1]);
[yreg,xreg,zsize]=size(stk_region);
[X1, XZm, YZm] = threeWayProject(handles.stkM, handles, roi, X2);
[X12, XZm2, YZm2] = threeWayProject(handles.stkM2, handles, roi, X3);
[fh_seg3d_threshold, Ldot, Rdot] = seg3d_threshold(xreg,yreg,handles.stacklength, stk_region, X1, XZm, YZm, 'Adjust the threshold to find the cell outline', 1);
close();
stk_region2 = permute(X3, [2 3 1]);
[fh_seg3d_threshold2, Ldot2, Rdot2] = seg3d_threshold(xreg,yreg,handles.stacklength, stk_region2, X12, XZm2, YZm2, 'Adjust the threshold find the outline of the organelle', 0);
close();
zi=Ldot(1);
zf=Rdot(1);
zi2=Ldot2(1);
zf2=Rdot2(1);
thresh2d_val=[Ldot(2) Rdot(2)];
%create the waitbar
h = waitbar(0,['Processing Region number ' num2str(roi) ' of ' num2str(handles.regions)]);
%do the segmentation of the cell outline
ellipsevar = [];
int_lim=65535;
slope=(Rdot(2)-Ldot(2))/(Rdot(1)-Ldot(1));
imtmp=uint16(zeros(ceil(yreg),ceil(xreg),Rdot(1)-Ldot(1)+1));
for i=zi:zf
imtmp=stk_region(:,:,i);
bw=im2uint16(im2bw(imtmp,(Ldot(2) + slope*(i-zi))/int_lim));
bw=bwperim(imfill(imfill(bwperim(bw), 'holes'), 'holes'));
ellipsevar(:,:,i) = bw;
waitbar((i)/(3*(zf-zi)),h);
end
%do the segmentation of the organelles
ellipsevar2 = [];
clearvars imtmp;
slope=(Rdot2(2)-Ldot2(2))/(Rdot2(1)-Ldot2(1));
imtmp=uint16(zeros(ceil(yreg),ceil(xreg),Rdot2(1)-Ldot2(1)+1));
for i=zi:zf
imtmp=stk_region2(:,:,i);
%get the broad outline
bw=im2uint16(im2bw(imtmp,(Ldot2(2)+slope*(i-zi))/int_lim));
bw=bwperim(bw);
%
ellipsevar2(:,:,i) = bw;
waitbar((zf-zi + i)/(3*(zf-zi)),h);
end
%translate the black and white (1 and 0) image maps into coordinates
pointsK = [];
pointsSelect = [];
pointsK2 = [];
pointsSelect2 = [];
% figure;
for k = zi:zf
pointsK = [];
for j = 1:ymax
for i = 1:xmax
if ellipsevar(i,j,k) == 1
pointsK = [pointsK; i j k];
end
if ellipsevar2(i,j,k) == 1
pointsK2 = [pointsK2; i j k];
end
end
end
sizeator = size(pointsK);
pointsSelect = [pointsSelect; pointsK];
pointsSelect2 = [pointsSelect2; pointsK2];
waitbar((2*(zf-zi) + k)/(3*(zf-zi)),h);
end
%shift the coordinates to corrispond to the actual size of the cell
pointsSelect(:,:,:) = [pointsSelect(:,1).*handles.pixelSize(1) pointsSelect(:,2).*handles.pixelSize(1) pointsSelect(:,3).*handles.pixelSize(2)];
pointsSelect2(:,:,:) = [pointsSelect2(:,1).*handles.pixelSize(1) pointsSelect2(:,2).*handles.pixelSize(1) pointsSelect2(:,3).*handles.pixelSize(2)];
%take a random 10% of the surface points
%pSSize = size(pointsSelect);
%pSindex = ceil(rand(ceil(.1*pSSize(1)), 1).*pSSize(1));
%pointsSelect = [pointsSelect(pSindex, 1) pointsSelect(pSindex, 2) pointsSelect(pSindex, 3)];
handles.roi(roi).pointsSelectCell = pointsSelect;
handles.roi(roi).pointsSelectOrganelle = pointsSelect2;
% figure(3 + roi);
% scatter3(pointsSelect(:,1), pointsSelect(:,2), pointsSelect(:,3), 'blue');
% hold on;
% scatter3(pointsSelect2(:,1), pointsSelect2(:,2), pointsSelect2(:,3), 'red');
%figure(3 + roi);
%scatter(pointsSelect(:,1), pointsSelect(:,3), 'blue');
%hold on;
%scatter(pointsSelect2(:,1), pointsSelect2(:,3), 'red');
%figure(3 + roi);
%K = convhulln(pointsSelect);
%trisurf(K,pointsSelect(:,1),pointsSelect(:,2),pointsSelect(:,3));
close(h);
end
[a b] = fileparts(handles.name);
filename = strcat(b, '.mat');
ch1_filename = handles.name;
ch2_filename = handles.name2;
out_points = handles.roi;
save(filename, 'ch1_filename', 'ch2_filename', 'Ldot', 'Rdot', 'out_points');
return
function [h, POINTS] = plot_ellipsoid(center, radii, angles)
theta = 0:.25:pi;
phi = 0:.25:2*pi;
k = 1;
for i = 1:length(phi)
for j = 1:length(theta)
y(k) = center(1) + radii(1)*sin(phi(i)-angles(1))*cos(theta(j)-angles(2));
x(k) = center(2) + radii(2)*sin(phi(i)-angles(1))*sin(theta(j)-angles(2));
z(k) = center(3) + radii(3)*cos(phi(i)-angles(1));
k = k + 1;
end
end
scatter3(x,y,z, 'red')
return
function [image] = seg_img(I, thresh)
BWs = edge(I, 'sobel', thresh);
se90 = strel('line', 3, 90);
se0 = strel('line', 3, 0);
BWsdil = imdilate(BWs, [se90 se0]);
BWdfill = imfill(BWsdil, 'holes');
BWnobord = imclearborder(BWdfill, 4);
seD = strel('diamond',1);
BWfinal = imerode(BWnobord,seD);
BWfinal = imerode(BWfinal,seD);
image = bwperim(BWfinal);
return
function [X1, XZm, YZm] = threeWayProject(stkM, handles, roi, X2)
%xy plot
X1 = imcrop(stkM, handles.roi(roi).RECT);
XZm = zeros(handles.stacklength,round(handles.roi(roi).RECT(4)));
for i = 1:handles.roi(roi).RECT(3)
XZm = max(cat(3,XZm,X2(:,:,i)),[],3);
end
XZm = transpose(XZm);
for i = 1:round(handles.roi(roi).RECT(3))
for j = 1:round(handles.roi(roi).RECT(4))
Y2(:,i,j) = X2(:,j,i);
end
end
%YZ
YZm = zeros(handles.stacklength,round(handles.roi(roi).RECT(3)));
for i = 1:handles.roi(roi).RECT(4)
YZm = max(cat(3,YZm,Y2(:,:,i)),[],3);
end
YZm = transpose(YZm);
return
function stackname_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function stackname2_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
% function to calculate 3D guassian filter
function h=gauss3D_filter(fsize,alpha)
r=fsize/2;
for i=1:fsize
for j=1:fsize
for k=1:fsize
a(i,j,k)=(i-r)^2+(j-r)^2+(k-r)^2;
b(i,j,k)=exp(a(i,j,k)*(-1/(2*alpha^2)));
end
end
end
h=b/sum(b(:));
return
