function einzelreader(handles, h)

global xpix ypix wbox psf_w02;

version = '6/23/08'; %record version (date) of gui being used

rbox=str2double(get(handles.rbox_edit,'String'));

wbox=2*rbox+1;

[xpix,ypix] = meshgrid(-rbox:rbox,-rbox:rbox);

q = str2double(get(handles.cam_pix_edit,'String'));

wvlnth = str2double(get(handles.wvlnth,'String'))/1000; %convert wavelength from nm to um

NA = str2double(get(handles.NA,'String'));

psf_scale = str2double(get(handles.psf_scale_edit,'String'));

psf_w0 = psf_scale*0.55*wvlnth/NA/1.17; % 1/e2 radius of PSF in um, use 1/e2 = FWHM/1.17 from Pawley

% with scale factor (20% for "real objective" due to

% measured PSF from Hess and Webb

% 2002)

psf_std=psf_w0/2; %standard deviation of psf

psf_w02=(psf_w0/q)*(psf_w0/q); %square of 1/e^2 radius in pixels

yfit_psf=exp(-2*((xpix).*(xpix)+(ypix).*(ypix))/psf_w02);

npix=sum(sum(yfit_psf)); % area of molecule in square pixels

total_molecules=0;

xcm_all=zeros(1,1);

ycm_all=zeros(1,1);

framenum_all=zeros(1,1);

n_fail_a0=0;

n_fail_outbox=0;

imagefile = get(handles.im_file_edit, 'String');

outpath = get(handles.out_dir_edit, 'String');

if outpath(length(outpath))~='/'

outpath = [outpath,'/'];

end

preview = get(handles.show_preview,'Value');

fs0 = str2double(get(handles.fs,'String')); %Setup for preview frame skipping

fs1 = fs0;

if get(handles.rolling_ball_cb,'Value') %if using rolling ball subtraction

%Obtain info on image type, etc

file_info = get_file_info2(imagefile, handles);

n_end = file_info.stop;

n_start = file_info.start;

if get(handles.custom_roi,'Value')

x_offset = str2double(get(handles.x_off_edit,'String'));

y_offset = str2double(get(handles.y_off_edit,'String'));

x_size = str2double(get(handles.x_size_edit,'String'));

y_size = str2double(get(handles.y_size_edit,'String'));

else

x_offset = 1;

y_offset = 1;

x_size = file_info.width;

y_size = file_info.height;

end

rball=str2double(get(handles.rb_radius_edit,'String')); %radius of rolling ball

se = strel('ball',rball,rball,0); %structural element, i.e. rolling ball

FWHM=str2double(get(handles.FWHM_edit,'String'));; %FWHM of gaussian smoothing in pixels

rk=(FWHM)/sqrt(2*log(2)); %1/e^2 smoothing radius in pixels

kw=20; %kernal width of smoothing function

[X,Y]=meshgrid(-kw/2:kw/2,-kw/2:kw/2);

kd=sqrt(X.*X+Y.*Y);

gs=exp(-2*kd.*kd/(rk*rk));

gs=gs/sum(sum(gs)); %smoothing function normalized to have area = 1

else %else use integrated widefield sum

%Load wf sum & mean background mat file

meanbkg_mat_file = getappdata(0,'last_mbkg_save'); %strcat(dir,'analysis\',date,'/',base_name,'sum_wf_meanbkg.mat');

load(meanbkg_mat_file,'meanbkg_all','image_sum','x_size','y_size','x_offset','y_offset','field0','n_field','file_info');%,'zero_level'); %mat file from sum_wf_meanbkg.m

n_end = file_info.stop;

n_start = file_info.start;

%give warning if starting/ending frames are not compatible with wf sum & mean background mat file

if(n_start<field0)

warndlg(sprintf('Warning: Starting frame is less than corresponding frame used in Compute WF & Mean Background, need to change starting frame or re-run Compute WF & Mean Background'));

return

end

if(n_end>n_field)

warndlg(sprintf('Warning: Ending frame is greater than corresponding frame used in Compute WF & Mean Background, need to change ending frame or re-run Compute WF & Mean Background'));

return

end

bkg_pc = str2double(get(handles.bkg_percent,'String'))/100;

zero_level = str2double(get(handles.zero_lvl,'String'));

%normalize so that ave pixel value = 1

image_sum_mean=mean(mean(image_sum(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1)));

image_sum_norm = (image_sum(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1))/image_sum_mean;

end

%Load variables from the GUI edit boxes

%--------------------------------------

if get(handles.pixel_photon,'Value')

pix_to_pho = str2double(get(handles.pix_to_pho,'String'));

else

pix_to_pho = 1;

end

iprod_thresh = str2double(get(handles.iprod_thresh_edit,'String'))*pix_to_pho; %initial threshold for a "bright object"

threshold = str2double(get(handles.threshold_edit,'String'))*pix_to_pho; %this value to be included

upper_threshold = str2double(get(handles.upperthresh_edit,'String'))*pix_to_pho; %this value to not be excluded

n_bright_pixel_threshold = str2double(get(handles.min_above_edit,'String')); %min # of pixels above

box_overlap_factor = str2double(get(handles.box_overlap_edit,'String')); %if center to center closer, don't include either

max_pixels_above_threshold = str2double(get(handles.max_above_edit,'String')); %max # of pixels above

%Initialize temporary runtime file for passing variables

%-------------

for x=1:99

temp_runtime = [outpath,'temp_runtime_',num2str(x,'%2.2d'),'.mat'];

if ~exist(temp_runtime,'file')

break

end

end

setappdata(0,'temp_runtime',temp_runtime);

%-------------

%Initialize waitbar and pause component

w = waitbarxmod(0,'Executing "einzelreader.m" ...','CreateCancelBtn','delete(gcf)');

set(w,'Name','Progress Bar');

uicontrol('Style','pushbutton','Parent',w,'String','Pause','Position',[210,10,60,23], ...

'UserData',1,'Callback',@pause_gui);

drawnow; %Draw the extra button immediately

pause(.1); %Pause to ensure window completes drawing

%keepontop('Progress Bar');

wb_norm = n_end-n_start;

if wb_norm == 0 %Avoid divide by zero errors

wb_norm = 1;

end

if h~=0 %If preview window is active, store the handle of the progress bar

setappdata(h,'PrBar_Handle',w); %in the appdata of the figure

end

for fileloop=n_start:n_end

if strcmp(file_info.type,'singlepage')

index = num2str(fileloop,file_info.prec);

infile = [file_info.path,'/',file_info.part_name,index,file_info.ext];

i1=double(imread(infile));

elseif strcmp(file_info.type,'multipage')

i1=double(imread(imagefile,fileloop));

end

compl = (fileloop-n_start)/wb_norm;

drawnow;

waitbarxmod(compl,w,'Executing "einzelreader.m" ...'); %Rename the waitbar

if get(handles.rolling_ball_cb,'Value') %if using rolling ball subtraction

i1_gs = uint16(conv2(i1,gs,'same')); %smoothed original

rb_im = imopen(i1_gs,se);

bkg = zeros(y_size,x_size);

bkg = double(rb_im(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1));

else

bkg = zeros(y_size,x_size);

bkg = image_sum_norm*bkg_pc*meanbkg_all(fileloop)+zero_level;

end

iprod = zeros(y_size,x_size);

iprod = i1(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1)-bkg;

iprod=iprod.*(iprod>0); % set any negative pixel values to zero

high_pixel_mask = zeros(y_size,x_size);

high_pixel_xy = zeros(10000,3); % x,y,active

n_high_pixels = 0;

pix_val=iprod_thresh;

while(pix_val >= iprod_thresh) %continue until reach minimum threshold

pix_val=0;

for i=rbox+1:y_size-rbox-1

for j=rbox+1:x_size-rbox-1

if high_pixel_mask(i,j)==0 && iprod(i,j)>pix_val

pix_val=iprod(i,j); %find brightest pixel

high_pixel_y=i;

high_pixel_x=j;

end

end

end

if(pix_val < iprod_thresh)

break

end

x0box=high_pixel_x-rbox;

y0box=high_pixel_y-rbox;

x1box=high_pixel_x+rbox;

y1box=high_pixel_y+rbox;

high_pixel_mask(y0box:y1box,x0box:x1box)=1;

n_high_pixels=n_high_pixels+1;

high_pixel_xy(n_high_pixels,1)=high_pixel_x;

high_pixel_xy(n_high_pixels,2)=high_pixel_y;

high_pixel_xy(n_high_pixels,3)=1; % active

end

drawnow; %**These added to keep waitbar responsive...

% Without these 'drawnow's, the queue gets large and

% ignores button presses during execution...

%(this function is fast when there is no

% actual drawing to do)

waitbarxmod(compl,w); % <-check waitbar (fast,see waitbarxmod.m)

for i=1:n_high_pixels

for j=1:i

if (i~=j & high_pixel_xy(i,3)==1 & high_pixel_xy(j,3)==1)

dx=high_pixel_xy(i,1)-high_pixel_xy(j,1);

dy=high_pixel_xy(i,2)-high_pixel_xy(j,2);

dmin_nearest_box=sqrt(dx*dx+dy*dy);

if (dmin_nearest_box<box_overlap_factor*rbox) % these boxes overlap

if j>i

high_pixel_xy(j,3)=-3; % make one of the boxes inactive

else

high_pixel_xy(i,3)=-3; % make one of the boxes inactive

end

end

end

end

end

drawnow;

waitbarxmod(compl,w); %check waitbar

n_boxes=0;

boxes_xy=zeros(10000,3);

for i=1:n_high_pixels

if high_pixel_xy(i,3)==1

n_boxes=n_boxes+1;

boxes_xy(n_boxes,1)=high_pixel_xy(i,1);

boxes_xy(n_boxes,2)=high_pixel_xy(i,2);

boxes_xy(n_boxes,3)=1;

end

end

drawnow;

waitbarxmod(compl,w); %check waitbar

failed=zeros(n_boxes,1);

for i=1:n_boxes

x0_box=boxes_xy(i,1)-rbox;

y0_box=boxes_xy(i,2)-rbox;

x1_box=boxes_xy(i,1)+rbox;

y1_box=boxes_xy(i,2)+rbox;

grab=iprod(y0_box:y1_box,x0_box:x1_box);

minbkg=min(min(grab));

grab=grab-minbkg;

xm_sum=0;

ym_sum=0;

m_sum=0;

for x=x0_box:x1_box

for y=y0_box:y1_box

xind=floor(x);

yind=floor(y);

intens=iprod(yind,xind);

xm_sum=xm_sum+xind*intens;

ym_sum=ym_sum+yind*intens;

m_sum=m_sum+intens;

end

end

x_cm(i)=xm_sum/m_sum;

y_cm(i)=ym_sum/m_sum;

%re-center box around center of mass

if round(x_cm(i)) > 1+rbox && round(x_cm(i)) < x_size-rbox

boxes_xy(i,1)=round(x_cm(i));

end

if round(y_cm(i)) > 1+rbox && round(y_cm(i)) < y_size-rbox

boxes_xy(i,2)=round(y_cm(i));

end

x0_box=boxes_xy(i,1)-rbox;

y0_box=boxes_xy(i,2)-rbox;

x1_box=boxes_xy(i,1)+rbox;

y1_box=boxes_xy(i,2)+rbox;

grab=zeros(wbox,wbox);

grab=iprod(y0_box:y1_box,x0_box:x1_box);

minbkg=min(min(grab));

grab=grab-minbkg;

grab_sum(i)=sum(sum(grab));

image_bright_in_box(:,:,i)=grab;

xc_box=(x0_box+x1_box)*0.5;

yc_box=(y0_box+y1_box)*0.5;

xguess=x_cm(i)-xc_box;

yguess=y_cm(i)-yc_box;

for n=1:wbox

for p=1:wbox

k=(n-1)*wbox+p;

xymerge(k)=0;

zmerge(k)=grab(n,p);

end

end

beta=[xguess,yguess,50];

[betafit,resid,J,COVB,mse] = nlinfit(xymerge,zmerge,@gaussian_merge,beta);

ci = nlparci(betafit,resid,'covar',COVB); %calculate error estimates on parameters

ci_err=(ci(:,2)-ci(:,1))/2;

xf_err(i)=ci_err(1);

yf_err(i)=ci_err(2);

a0_err(i)=ci_err(3);

zfitmerge=gaussian_merge(betafit,xymerge);

for n=1:wbox

for p=1:wbox

k=(n-1)*wbox+p;

zfit(n,p)=zfitmerge(k);

end

end

yf_box(i)=betafit(2)+yc_box;

xf_box(i)=betafit(1)+xc_box;

a0_box(i)=betafit(3);

if(a0_box(i) < 0)

n_fail_a0=n_fail_a0+1;

failed(i)=1;

end

if xf_box(i) > x1_box || xf_box(i) < x0_box || yf_box(i) > y1_box || yf_box(i) < y0_box

n_fail_outbox=n_fail_outbox+1;

failed(i)=1;

end

drawnow %Flush event queue to keep cancel/pause button responsive on every pass.

%This is the slowest function by far and drawnow barely changes its exec time

waitbarxmod(compl,w); %check waitbar

end

n_pixels_above_threshold=zeros(n_boxes,1);

n_pixels_above_upper_threshold=zeros(n_boxes,1);

for i=1:n_boxes

if (boxes_xy(i,3)==1)

n_pixels_above_threshold(i)=0;

for j=1:2*rbox+1

for k=1:2*rbox+1

intens=image_bright_in_box(j,k,i);

if intens>=threshold

n_pixels_above_threshold(i)=n_pixels_above_threshold(i)+1;

end

if intens>=upper_threshold

n_pixels_above_upper_threshold(i)=n_pixels_above_upper_threshold(i)+1;

end

end

end

end

end

drawnow;

waitbarxmod(compl,w); %check waitbar

n_molecules_found=0;

for i=1:n_boxes

if boxes_xy(i,3)==1

if (n_pixels_above_threshold(i)>=n_bright_pixel_threshold && n_pixels_above_upper_threshold(i)<=max_pixels_above_threshold)

n_molecules_found=n_molecules_found+1;

else

if (n_pixels_above_threshold(i)<n_bright_pixel_threshold)

boxes_xy(i,3)=-2; %draw red box

end

if (n_pixels_above_threshold(i)>max_pixels_above_threshold)

boxes_xy(i,3)=-1; %draw yellow box

end

end

end

end

drawnow;

waitbarxmod(compl,w); %check waitbar

if preview %If preview active, update images on screen

set(0,'CurrentFigure',h);

if ~exist('plot2','var')

position2 = [.01,.37,.25,.25];

plot2 = axes('Position', position2);

setappdata(plot2,'pos',position2);

else

set(h,'CurrentAxes',plot2);

end

i1display = i1(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1);

imagesc(i1display,'HitTest','off');

axis image;

colormap gray;

title(['Last Frame Processed: ', num2str(fileloop)])

set(plot2, 'ButtonDownFcn',{@focus_swap,guidata(gcbo)});

clear i1display;

if ~exist('plot1','var')

position1 = [.01,.05,.25,.25];

plot1 = axes('Position', position1,'Color','k','YTickLabel','','XTickLabel','');

axis image;

set(plot1,'Position',position1);

setappdata(plot1,'pos', position1);

end

if get(handles.rolling_ball_cb,'Value')

if ~exist('plot3','var')

position3 = [.01,.70,.25,.25];

plot3 = axes('Position', position3);

setappdata(plot3,'pos',position3);

else

set(h,'CurrentAxes',plot3);

end

imagesc(rb_im,'HitTest','off')

colormap gray, axis image;

title(['Rolling Ball Profile: ',num2str(fileloop)])

set(plot3,'ButtonDownFcn',{@focus_swap,guidata(gcbo)});

end

if ~exist('plot4','var')

if x_size < y_size

positionf = [.18,.065,.87,.87];

else

positionf = [.32,.3,.6,.65];

end

plot4 = axes('Position', positionf);

setappdata(0,'focused',plot4);

else

set(h,'CurrentAxes',plot4);

end

if get(handles.mnl_scale_check,'Value')

clim = str2double(get(handles.mnl_scale_edit,'String'));

imagesc(iprod,'HitTest','off',[0 clim]);

else

imagesc(iprod,'HitTest','off');

end

axis image;

colormap gray;

title(['Einzel Reader: Frame: ',num2str(fileloop)]);

set(plot4, 'ButtonDownFcn',{@focus_swap,guidata(gcbo)});

if (get(handles.show_colorbar,'Value') == 1)

focused = getappdata(0, 'focused');

set(h,'CurrentAxes', focused);

pos = get(focused, 'Position');

colorbar('Position',[.947, pos(2), .015, pos(4)],'DrawMode','fast')

end

set(h,'CurrentAxes', plot4);

end

if n_boxes>0

if preview %Only draw if preview is active

hold on

draw_boxes(n_boxes,boxes_xy,rbox);

end

for i=1:n_boxes

if boxes_xy(i,3)==1 && failed(i) ~= 1

if preview %Only plot if preview is active

%plot(x_cm(i),y_cm(i),'.m','HitTest','off'); %centroid

%plot(xf_box(i),yf_box(i),'.b','HitTest','off'); %gaussian fit

end

total_molecules=total_molecules+1;

xcm_all(total_molecules)=x_cm(i);

ycm_all(total_molecules)=y_cm(i);

xf_all(total_molecules)=xf_box(i);

yf_all(total_molecules)=yf_box(i);

a0_all(total_molecules)=a0_box(i);

grab_sum_all(total_molecules)=grab_sum(i);

framenum_all(total_molecules)=fileloop;

xf_err_all(total_molecules)=xf_err(i);

yf_err_all(total_molecules)=yf_err(i);

a0_err_all(total_molecules)=a0_err(i);

end

end

nmol_all(fileloop)=n_boxes;%%

if preview %Remove the hold

hold off

end

drawnow; %Draw any boxes and update einzelreader image (flush queue and check waitbar if preview not active)

waitbarxmod(compl,w); %Update the waitbar

end

if get(handles.frame_delay_check,'Value')

delay = str2double(get(handles.frame_delay_edit,'String'))/1000;

if delay ~= inf

pause(delay);

end

end

if preview %If preview active, compute wf & merge on first, last, and nth frames (always)

if fileloop == n_end

setappdata(0,'Save',1); %Trigger image save on last frame

save(temp_runtime); %Save variables to a temporary file for function passing

elseif fileloop == n_start+fs1-1

fs1=fs1+fs0;

save(temp_runtime);

elseif fileloop == n_start

save(temp_runtime);

end

end

end

if ~preview %If preview is off, continue to wf & merge at end and save the result

plot1 = 0;

save(temp_runtime);

setappdata(0,'Save',1);

fpalm_render_einzelreader(handles,h);

else %Otherwise delete progress bar and handle stored in appdata

delete(w);

rmappdata(h,'PrBar_Handle');

end

if get(handles.oc_wf,'Value') %save widefield image if selected

if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators

wf_file = [outpath,'/',file_info.part_name,'wf.tif'];

else

wf_file = [outpath,'/',file_info.part_name,'_wf.tif'];

end

imwrite(uint8((image_sum/max(image_sum(:)))*255),wf_file,'Compression','none');

end

if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators

out_file = [outpath,file_info.part_name,num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];

else

out_file = [outpath,file_info.part_name,'_',num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];

end

if ~exist(outpath,'dir') %Make sure the directory exists

answer = questdlg('Output directory could not be accessed. Choose a new save location?','Warning!','Yes','Discard','Yes');

if strcmp(answer,'Yes')

outpath = uigetdir();

if ~outpath

return

end

if outpath(length(outpath))~='/'

outpath = [outpath,'/'];

end

if (file_info.part_name(length(file_info.part_name)) == '_')

out_file = [outpath,file_info.part_name,num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];

else

out_file = [outpath,file_info.part_name,'_',num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'.mat'];

end

else

return

end

end

prec_edit = get(handles.prec_edit,'String'); %<- These to be saved with the

bkgn = str2double(get(handles.bkgn_noise,'String')); %<- mat file for loading purposes

ppp = str2double(get(handles.pix_to_pho,'String')); %<-

jshift = str2double(get(handles.jshift,'String')); %<-

ishift = str2double(get(handles.ishift,'String')); %<-

exf = str2double(get(handles.exf_edit,'String')); %<- expansion factor

a0_phot=a0_all/ppp; % peak amplitude of each molecule in photons

N=npix*a0_phot; % number of photons for each molecule

%localization precision in um

lp2=((psf_std^2)+(q^2)/12)*1./N+8*pi*(psf_std^4)*(bkgn^2)/(q^2)*1./(N.*N);

lp=sqrt(lp2);

list = who; %Use regular expressions to ensure handles are not saved

match = regexp(list,'^plot|handles|\<h\>|\<w\>');

for i=length(match):-1:1

if match{i}

list(i)=[];

end

end

save(out_file,list{:}); %Save variables used in einzelreader.m

setappdata(0,'last_einzel_save',out_file); %Save location of einzelreader mat file for other functions to use

%Clean up temp files and appdata

if exist(temp_runtime,'file')

delete(temp_runtime)

end

if isappdata(0,'temp_runtime')

rmappdata(0,'temp_runtime')

end