diff --git a/2-Channel Codes/channel_calibration.m b/2-Channel Codes/channel_calibration.m index 6ee4ce6..db6c9d5 100644 --- a/2-Channel Codes/channel_calibration.m +++ b/2-Channel Codes/channel_calibration.m @@ -8,10 +8,21 @@ % load calibration data and store in conglomorate array for i = 1:numel(files) load(files(i).name) +<<<<<<< HEAD + xf = [xf; cdata.red.ncoords(:,1)]; + yf = [yf; cdata.red.ncoords(:,2)]; + N = [N; cdata.red.ncoords(:,3)]; + xf = [xf; cdata.orange.ncoords(:,1)]; + yf = [yf; cdata.orange.ncoords(:,2)]; + N = [N; cdata.orange.ncoords(:,3)]; + fnum = [fnum; cdata.red.framenumber + max([max(fnum),0])]; + fnum = [fnum; cdata.orange.framenumber + max([max(fnum),0])]; +======= xf = [xf; fits(:,1)]; yf = [yf; fits(:,2)]; N = [N; fits(:,3)]; fnum = [fnum; framenumber + max([max(fnum),0])]; +>>>>>>> master end % display calibration parameters of splitting channels via horizontal diff --git a/2-Channel Codes/two_channel_drift_with_ruler.mlx b/2-Channel Codes/two_channel_drift_with_ruler.mlx index 89ee0bf..ff8782f 100644 Binary files a/2-Channel Codes/two_channel_drift_with_ruler.mlx and b/2-Channel Codes/two_channel_drift_with_ruler.mlx differ diff --git a/2-Channel Codes/xy_feature.m b/2-Channel Codes/xy_feature.m index c8fea3a..b3792a2 100644 --- a/2-Channel Codes/xy_feature.m +++ b/2-Channel Codes/xy_feature.m @@ -2,6 +2,10 @@ x = x(:); y=y(:); % vector = [x.^3, x.^2.*y, x.*y.^2, x.^2, y.^2, x.*y,x, y, x*0+1]; +<<<<<<< HEAD +% vector = [x.^2, y.^2, x.*y,x, y, x*0+1]; +======= vector = [x.^2, y.^2, x.*y,x, y, x*0+1]; +>>>>>>> master vector = [x.^3, y.^3, x.^2.*y, x.*y.^2, x.^2, y.^2, x.*y,x, y, x*0+1]; \ No newline at end of file diff --git a/3D-Localization/total_calibration.m b/3D-Localization/total_calibration.m index 2204a31..7b9da2c 100644 --- a/3D-Localization/total_calibration.m +++ b/3D-Localization/total_calibration.m @@ -284,18 +284,18 @@ end % Find plane of least confusion -ds = ssx - ssy; % Subtract sigma-y from sigma-x +ds = ssy - ssx; % Subtract sigma-y from sigma-x tds = uitab(tg3,'Title','DS curve'); ax = axes(tds); plot(ax,abs(ds)) xlabel(ax,'Index') ylabel(ax,'Difference between average sx and sy') tdx = input('Choose an Index to grab minimum of the DS curve'); -ind = find(abs(ds(tdx:end)) == min(abs(ds(tdx:end)))); % Find minimum of absolute value +ind = find(abs(ds(tdx:tdx + 30)) == min(abs(ds(tdx:tdx + 30)))); % Find minimum of absolute value z0s = zus - zus(ind(1)+tdx); % Call the found index the 0 point % Data Representation of Sigma and Z space -t3 = uitab(tg,'Title','Sigmas'); +t3 = uitab(tg,'Title','Sigmas2'); ax = axes(t3); plot(ax,z0(indy)-zus(ind(1)+tdx),sx(indy),'.') % Plot Toleranced Sigma-x data hold on @@ -304,13 +304,17 @@ plot(ax,z0s,ssy) % Plot average sigma-y % Determining Z parameters -ind = abs(z0s) < 1.2; % Limit height over which to fit sigmas +ind = abs(z0s) < 0.5; % Limit height over which to fit sigmas % z_net = train_neural_z(z0s(indy),ims(:,indy)); -z_cal = get_z_params(z0s(ind),ssx(ind),ssy(ind)); % Fit sigma curves to data and return result - +z_cal = get_z_params(z0s(ind),ssx(ind),ssy(ind)); +% Fit sigma curves to data and return result +% Redfeining z_call and % Display results of Z-calibration yx = z_cal_fit(z0s(ind),z_cal(1:5)); % Determine ideal fitted Sig-x Values -yy = z_cal_fit(z0s(ind),z_cal(6:end)); % Determine ideal fitted Sig-y values +yy = z_cal_fit(z0s(ind),z_cal(6:end)); + + +% Determine ideal fitted Sig-y values % Overlay result on scattered / average sig-x plot plot(ax,z0s(ind),yx,'gx') plot(ax,z0s(ind),yy,'gx') @@ -321,7 +325,7 @@ % It's known that in astigmatism there may be a slight slant that % exists over the Z direction, in this section we'll address that -zf_um = getdz(sx,sy,z_cal); % Get Z-values +zf_um = getdz(sx,sy,z_cal,2); % Get Z-values upz = 0.5; dnz = -0.5; indy = indy & zf_um dnz; % Limit view to fitted region listed above @@ -341,17 +345,22 @@ % Overlay stage v fit curves zc = z0; -for i = 1:numel(psfs) - id = psf == i & abs(zf_um) < 1; % ident by psf and position +figure +for i = 4:numel(psfs) + id = psf == i & abs(z0 - 0.88) < 0.22; % ident by psf and position zs = z0(id); % subset of stage positions zfs = zf_um(id); % subset of fitted positions + plot(zs,zfs,'.') a = polyfit(zfs,zs,1); % linear fit to obtain x-intercept of stage zc(id) = zc(id) - a(2); % adjust stage position to 'center' by subtracting x-interecept of all psfs + hold on end + % plot(ax,z0(indy),zf_um(indy),'.') % a = polyfit(z0(indy),zf_um(indy),1); % The slope of this distribution gives us the 'correction' for absolute Z +%purge clumped psfs for orange % grab subsets dist = 0.5; @@ -360,6 +369,11 @@ yfs = yf(indy); zfs = zf_um(indy)/q; pfs = psf(indy); +done_id = psf <4; +xfs(done_id) = []; +yfs(done_id) = []; +zfs(done_id) = []; +pfs(done_id) = []; zs = (min(zfs*q):0.040:max(zfs*q))/q; % zfs is in pixels, are values are in pixels right now for i = 1:numel(zs)-1 % Attempt to grab the magnification curve of stage versus fit id = zf_um >= zs(i)*q & zf_um < zs(i+1)*q; @@ -369,7 +383,15 @@ z0m(i) = mean(zc(id2)); zfm(i) = mean(zf_um(id2)); end -a = polyfit(z0m,zfm,1); +id = zf_um < 0.18; +z0m = zc(id); +% a = polyfit(z0m,zfm,1); +ind = psf >3 & abs(z0 - 0.88) < 0.22; +plot(z0(ind),zf_um(ind),'.') +a_actual = polyfit(z0(ind),zf_um(ind),1); +hold on +plot(z0(ind),z0(ind)*a_actual(1) + a_actual(2),'r') + clear zfm plot(ax,zc(indy),zf_um(indy),'.'); hold on @@ -384,11 +406,14 @@ xt = []; yt = []; zt = []; -for i = 1:max(pfs) - ind = pfs == i; - xt = [xt;xfs(ind)-mean(xfs(ind))]; - yt = [yt;yfs(ind)-mean(yfs(ind))]; - zt = [zt;zfs(ind)]; +for i = 4:max(psf) + ind = psf == i; + try + xt = [xt;xfs(ind)-mean(xfs(ind))]; + yt = [yt;yfs(ind)-mean(yfs(ind))]; + zt = [zt;zfs(ind)]; + catch + end end % indy = indy & xfc.^0.5*q < 0.005 & yfc.^0.5*q < 0.005; % indy = indy & abs(xf) < 0.5 & abs(yf) < 0.5; @@ -398,12 +423,12 @@ zsel = []; for i = 1:numel(zs) - 1 % Loop over a variety of height windows for averaging of points ind1 = zt >= zs(i) & zt <=zs(i+1); % grab all points within height window - ind1 = ind1 & (xt.^2+yt.^2).^0.5 < 0.5; % grab all points who lie within a half pixel radius of center (at this point all scans have a 0 mean in XY) + ind1 = ind1 & ((xt - mean(xt(ind1))).^2+(yt - mean(yt(ind1))).^2).^0.5 < 0.5; % grab all points who lie within a half pixel radius of center (at this point all scans have a 0 mean in XY) % Subsets xts = xt(ind1); yts = yt(ind1); zts = zt(ind1); - rts = (xts.^2+yts.^2).^0.5; % convert subsetted xy to R + rts = ((xts-mean(xts)).^2+(yts-mean(yts)).^2).^0.5; % convert subsetted xy to R mr = mean(rts); % grab mean of R str = std(rts); % grab stdev of R % my = mean(yts); @@ -419,7 +444,7 @@ zfm(i) = mean(zts(ind2)); end -splz = (-1:0.001:1)/q; +splz = (-0.5:0.001:0.5)/q; xtilt = gausssmooth(xfm,5,10); ytilt = gausssmooth(yfm,5,10); splx = spline(zs(1:end-1) + mean(diff(zs))/2,xtilt,splz); @@ -442,6 +467,7 @@ zf = zf_um/q; xc = spline(zs(1:end-1) + mean(diff(zs))/2,xtilt,zf); yc = spline(zs(1:end-1) + mean(diff(zs))/2,ytilt,zf); +% new drift correction attempt xf_c = xf - xc; yf_c = yf - yc; @@ -464,4 +490,4 @@ %% Calibration Handling disp(['Distance between Minima is ', num2str(1000*(cal.z_cal(2)-cal.z_cal(7))),'nm']); -save('z_calib.mat','cal') +save('z_calib_orange.mat','cal') diff --git a/3D-Localization/total_calibration_1.m b/3D-Localization/total_calibration_1.m new file mode 100644 index 0000000..1a57629 --- /dev/null +++ b/3D-Localization/total_calibration_1.m @@ -0,0 +1,468 @@ +% Calibration for 3D +% This script is intended to be a one-stop-shop for astigmatic calibration +% of scan data. Comments are heavy to help guide user through programming +% AJN 10-8-18 Ryan Lab + +% Preclean UI +clearvars; +close all; +clc; +%% User variables +pixw = 6; % ROI 'radius' +q = 0.120; % Pixel size in um +step = 20; % Steps between frames in nm +CCs = 10; % number of frames to x-correlate over +wind = -pixw:pixw; % create window for segmentation + +%% END USER INPUT +% Create display figure as a tab group +% f = figure('units','Normalized','OuterPosition',[0 0 0.5 0.5]); % Initialize figure +f = figure; +tg = uitabgroup(f); % tg is the tabgroup to help reduce clutter of figures + +%% File loading and image segmentation +files = dir('*tif'); +t1 = uitab(tg,'Title','Mean Images'); % Create tab for image representation +tg1 = uitabgroup(t1); % Make tabgroup object +A= {}; % Saving images as a double cell variable +psfs = {}; % Saving scans of PSF as double cell variable +pind = 1; % Counting variable for populating psfs +for i = 1:numel(files) % Loop over all files +% A{i} = roball(readtiff(files(i).name),6,4); % store total images into variable +% A{i} = bandpass(readtiff(files(i).name),1.15,5); +% A{i} = imgaussfilt(readtiff(files(i).name), + A{i} = gpu_rball(readtiff(files(i).name)); +% A{i} = readtiff(files(i).name); + A{i}(:,1:4,:) = 0; + A{i}(end-4:end,:,:) = 0; + A{i}(:,end-4:end,:) = 0; + A{i}(1:4,:,:) = 0; +% for j = 1:numel(A{i}(1,1,:)) +% A{i}(:,:,j) = imgaussfilt(A{i}(:,:,j),0.18/q); +% A{i}(:,:,j) = band_pass(A{i}(:,:,j),2); +% end +% A{i} = bp_subtract(A{i}); + ax = axes(uitab(tg1,'Title',files(i).name(1:end-4))); % get axes for appropriate tab + imagesc(ax,max(A{i},[],3)) % Represent maximal image + psf = denoise_psf(max(A{i},[],3),2); % use wavelet transform to identify molecules + dps = das_peaks(psf,50); % Peak finder + [row,col] = find(dps == 1); % Find peaks in dps +% ind = find_dupes([col,row],1.5*pixw); % remove overlapping molecules + % Remove duplicate entries +% row(ind) = []; +% col(ind) = []; + for j = 1:numel(row) % loop over all locations, but attempt to center the maximum pixel on each sub image + draw_boxes([col,row],pixw); % Show location of found regions on representation + try + for k = 1:numel(A{i}(1,1,:)) + [m,n,o] = size(A{i}); + if row(j)-pixw > 0 && row(j)+ pixw < m && col(j)-pixw > 0 && col(j)+ pixw < n + [r,c] = find(A{i}(row(j) + wind, col(j) + wind,k) == max(max(A{i}(row(j) + wind, col(j) + wind,k)))); % find maxima pixel in subregion + psf_off{pind}(k,:) = [c(1),r(1)]; + psfs{pind}(:,:,k) = A{i}((row(j) - pixw - 1) + r(1) + wind, col(j) - pixw - 1 + c(1) + wind,k); % Store psfs in their own variable + end + end + pind = pind +1; + catch lsterr + end + end +end + +% Memory Management +% Here all interested regions are saved in psfs +clear A dps psf row col files pind + +%% Cross Correlate PSFs +[m,n,o] = size(psfs{1}); % Grab size of PSF regions +% t2 = uitab(tg,'Title','X-Correlation'); % Make tab for X-Corr representation +% ax = axes(t2); % Get axes for X-Corr plot +% +% % Correlate all PSFS to the first frame scan +% for i = 1:numel(psfs) +% [ind(i),mc] = cross3d(psfs{1},psfs{i}(:,:,round(o/2) + (-CCs:CCs))); % Perform Correlation +% plot(ax,(1:numel(mc(:))),mc(:)/max(mc(:))) % Represent Result +% hold on +% end +% hold off +% disp = ind - ind(1); % Make Correction to x-corrs by subtracting the identity index + +ang = []; +%% Fitting Analysis +% Get Eliptical Angle +for i = round(o/2) + (-CCs:-10) % Perform eliptical Angle determination over several frames + try + [a] = get_elip_ang(psfs{1}(:,:,i),2.5,1.5); % Find optimal eliptical angle + ang = [ang;a]; % Save result to an array + catch + end +end +ang = mean(ang); % take mean angle + +% Preallocate fitting variables +xf = []; +yf = []; +N = []; +sx = []; +sy = []; +O = []; +xfc = []; +yfx = []; +Oc =[]; +Nc = []; +sxc = []; +syc = []; +yfc = []; +llv = []; +fnum = []; +psf = []; +ims = []; +% gpu fit of PSFs +% NOTE TO USER: because we assume an axis of elipticity that +% need not be parallel to the axis of the camera we rotate our +% grid to perform the fit. This rotation is corrected for on +% the back end and the X-Y positions returned are understood to +% be on the axis of the camera. However, the sigma values exist +% along the rotated or prime axis. This is representated in the +% following comments once and is to be understood that +% Sigma-x/y corresponds to the prime axis and not the camera +% axis +fms = 1:o; % input framenum +for i = 1:numel(psfs) % Loop over all identified ROIs + if i ~= 1 && i ~= 6 + [fits, crlb, lv,fnout,im0s] = slim_locs(psfs{i},fms,zeros(o,2),ang,50,100); % Perform Fit + fnout = fnout.'; % Save frame number which corresponds to Z-position + xa = fits(:,1)+psf_off{i}(fnout,1); % define 'assignment' variable for x + ya = fits(:,2)+psf_off{i}(fnout,2); % define 'assignment' variable for y this allows centering around distribution + % xa = fits(:,1); + % ya = fits(:,2); + xf = [xf;xa-mean(xa)]; % X-Position + yf = [yf;ya-mean(ya)]; % Y-Position + % xf = [xf;xa]; % X-Position + % yf = [yf;ya]; % Y-Position + N = [N; fits(:,3)]; % Number of Photons + sx = [sx;fits(:,4)]; % Sigma in x' direction + sy = [sy;fits(:,5)]; % Sigma in y' direction + O = [O; fits(:,6)]; % Offset + % Lower bound on Variance of fitted variables + xfc = [xfc;crlb(:,1)]; + yfc = [yfc;crlb(:,2)]; + Nc = [Nc;crlb(:,3)]; + Oc = [Oc;crlb(:,6)]; + sxc = [sxc;crlb(:,4)]; + syc = [syc;crlb(:,5)]; + psf = [psf; fits(:,1)*i./fits(:,1)]; % Log the PSF the fit is associated with + llv = [llv;-abs(lv)]; % Log Likelihood Value + % fnum = [fnum;fnout - disp(i)]; % Correct the Frame number based off correlation result + ims = [ims,im0s]; + fnum = [fnum;fnout]; + end +end +z1 = fnum*step/1000; % Populate Z-positions +% Tolerance Step +%(User can adjust these values to change fidelity of data used to determine calibration results) +% For use with 10-3-18 scans +% indy = llv./N > -0.1 & N >0 & N < 10000 & abs(syc) < 0.01 & abs(sxc) < 0.01 ; +% indy = llv./N > -0.15 & N >0 & N < 100000 & abs(syc).^0.5*q < 0.006 & abs(sxc).^0.5*q < 0.006 ; +indy = sx < 5 & sx >0 & sy < 5 & sy > 0; + +%% Scan Correlation +refp = mode(psf(indy)); % find psf scan w/ the most number of passed fits, this is our reference psf +rid = psf == refp; % build an index for the reference ID +sigr = [fnum(indy&rid),sx(indy &rid),sy(indy & rid)]; % build reference sigma scan +fnumc = fnum; +for i = 1:numel(psfs) % loop over all psfs + if i ~= 1 && i ~= 6 + cid = psf == i; % grab Id for psf scan to correct + sigc = [fnum(indy&cid),sx(indy & cid),sy(indy & cid)]; + dff = xcorrsig(sigr,sigc); + + fnumc(cid) = fnumc(cid) + dff; + end +end +z0 = fnumc*step/1000; + +%% Data Representation +% ind = psf == 1 || psf == 6; +% fnumc(ind) = []; + + +% x corr corrections +t3 = uitab(tg,'Title','Correlation Correction'); +tg3 = uitabgroup(t3); +tunc = uitab(tg3,'Title','Without Height Adjustment'); % start with uncorrected +ax = axes(tunc); +plot(ax,z1,sx,'.'); +hold on +plot(ax,z1,sy,'.'); +ylim(ax,[0,5]) +xlim(ax,[0,3]); +hold off +legend('Sigma X', 'Sigma Y','Location','North'); +xlabel(ax,'Height (um)'); +ylabel(ax,'Width (pix)'); + +tc = uitab(tg3, 'Title','With Height Adjustment'); +ax = axes(tc); +plot(ax,z0,sx,'.'); +hold on +plot(ax,z0,sy,'.'); +ylim(ax,[0,5]) +xlim(ax,[0,3]); +hold off +legend('Sigma X', 'Sigma Y','Location','North'); +xlabel(ax,'Height (um)'); +ylabel(ax,'Width (pix)'); + +% Display X-Y Position +t4 = uitab(tg,'Title','Fitting Outputs'); +tg4 = uitabgroup(t4); +txy = uitab(tg4,'Title','X-Y Position'); +ax = axes(txy); +plot(ax,xf(indy),yf(indy),'.'); +axis equal +xlabel('X-Fit') +ylabel('Y-Fit'); + +% Display X-Y Uncertainties +tlp = uitab(tg4,'Title','X-Y Uncertainty'); +ax = axes(tlp); +histogram(ax,xfc(indy).^0.5*q,'Normalization','Probability'); +hold on +histogram(ax,yfc(indy).^0.5*q,'Normalization','Probability'); +legend('X-unc','Y-Unc'); +xlabel('Uncertainty in um') +ylabel('Probability'); +hold off + +% Display Uncertainties in Sigma Values +tsxy = uitab(tg4,'Title','Sigma Uncertainties'); +ax = axes(tsxy); +histogram(ax,sxc(indy).^0.5*q,'Normalization','Probability'); +hold on +histogram(ax,syc(indy).^0.5*q,'Normalization','Probability'); +legend('X-unc','Y-Unc'); +xlabel('Uncertainty in um') +ylabel('Probability'); + +% Display Histogram of Number of Photons +tN = uitab(tg4,'Title','Number of Photons'); +ax = axes(tN); +histogram(ax,N(indy),'Normalization','Probability'); +xlabel('Photons Detected') +ylabel('Probability'); + +% Display Histogram of Log-Likelihood Values +tllv = uitab(tg4,'Title','LLV'); +ax = axes(tllv); +histogram(ax,llv(indy),'Normalization','Probability'); +xlabel('Log Likelihood Value') +ylabel('Probability'); + +% Display 'goodness of fit' metric LLV/N +tiln = uitab(tg4,'Title','LLV./N'); +ax = axes(tiln); +histogram(ax,llv(indy)./N(indy),'Normalization','Probability'); +xlabel('Log Likelihood Value') +ylabel('Probability'); + +%% Z-Calibration +% select subset of data that passed tolerances for sx/sy fitting +sxs = sx(indy); % Sigma X +sys = sy(indy); % Sigma Y +z0s = z0(indy); % Z-values +zus = unique(z0s); % Create a Unique list of Z-positions for indexing purposes +for i = 1:numel(zus) % Loop over unique Z-Positions + ind = z0s == zus(i); % select the fits corresponding to a particular Z-Position + subsx = sxs(ind); % Create subset of Sigma-x corresponding to this Z-pos + subsy = sys(ind); % Create subset of Sigma-y corresponding to this Z-pos + msx = mean(subsx); % Determine Average Value of Sigma-x at this Z-pos + stx = std(subsx); % Determine standard deviation of Sigma-x at this Z-pos + msy = mean(subsy); % Determine Average Value of Sigma-y at this Z-pos + sty = std(subsy); % Determine standard deviation of Sigma-x at this Z-pos + ssx(i) = mean(subsx(subsx >= msx - 2* stx & subsx <= msx + 2* stx)); % Select only values within 2 standard deviations of the average for X' + ssy(i) = mean(subsy(subsy >= msy - 2* sty & subsy <= msy + 2* sty)); % Select only values within 2 standard deviations of the average for Y' + +end + +% Find plane of least confusion +ds = ssx - ssy; % Subtract sigma-y from sigma-x +tds = uitab(tg3,'Title','DS curve'); +ax = axes(tds); +plot(ax,abs(ds)) +xlabel(ax,'Index') +ylabel(ax,'Difference between average sx and sy') +tdx = input('Choose an Index to grab minimum of the DS curve'); +ind = find(abs(ds(tdx:end)) == min(abs(ds(tdx:end)))); % Find minimum of absolute value +z0s = zus - zus(ind(1)+tdx); % Call the found index the 0 point + +% Data Representation of Sigma and Z space +t3 = uitab(tg,'Title','Sigmas'); +ax = axes(t3); +plot(ax,z0(indy)-zus(ind(1)+tdx),sx(indy),'.') % Plot Toleranced Sigma-x data +hold on +plot(ax,z0(indy)-zus(ind(1)+tdx),sy(indy),'.') % Plot Toleranced Sigma-y data +plot(ax,z0s,ssx) % Plot average sigma-x +plot(ax,z0s,ssy) % Plot average sigma-y + +% Determining Z parameters +ind = abs(z0s) < 1.2; % Limit height over which to fit sigmas +% z_net = train_neural_z(z0s(indy),ims(:,indy)); +z_cal = get_z_params(z0s(ind),ssx(ind),ssy(ind)); % Fit sigma curves to data and return result + +% Display results of Z-calibration +yx = z_cal_fit(z0s(ind),z_cal(1:5)); % Determine ideal fitted Sig-x Values +yy = z_cal_fit(z0s(ind),z_cal(6:end)); % Determine ideal fitted Sig-y values +% Overlay result on scattered / average sig-x plot +plot(ax,z0s(ind),yx,'gx') +plot(ax,z0s(ind),yy,'gx') +hold off +xlim(ax,[-1.2,1.2]); +legend('Sig-X','Sigy-Y','Location','North'); +%% Determine Necessary x-y-z correction +% It's known that in astigmatism there may be a slight slant that +% exists over the Z direction, in this section we'll address that + +zf_um = getdz(sx,sy,z_cal); % Get Z-values +upz = 0.5; +dnz = -0.5; +indy = indy & zf_um dnz; % Limit view to fitted region listed above +d3 = uitab(tg4,'Title','3-D Positions'); +ax = axes(d3); +scatter3(xf(indy),yf(indy),zf_um(indy)/q,[],psf(indy)); +axis equal +xlabel('Lateral-X') +ylabel('Lateral-Y'); +zlabel('Axial-Z'); + +% Display z0 true v zf +tt4 = uitab(tg,'Title','Final Corrections'); +tg5 = uitabgroup(tt4); +tf = uitab(tg5,'Title','Z0 v. Zf'); +ax = axes(tf); + +% Overlay stage v fit curves +zc = z0; +for i = 1:numel(psfs) + id = psf == i & abs(zf_um) < 1; % ident by psf and position + zs = z0(id); % subset of stage positions + zfs = zf_um(id); % subset of fitted positions + a = polyfit(zfs,zs,1); % linear fit to obtain x-intercept of stage + zc(id) = zc(id) - a(2); % adjust stage position to 'center' by subtracting x-interecept of all psfs +end + +% plot(ax,z0(indy),zf_um(indy),'.') +% a = polyfit(z0(indy),zf_um(indy),1); % The slope of this distribution gives us the 'correction' for absolute Z + + +% grab subsets +dist = 0.5; +% indy = indy & abs(zf_um) < dist; +xfs = xf(indy); +yfs = yf(indy); +zfs = zf_um(indy)/q; +pfs = psf(indy); +zs = (min(zfs*q):0.040:max(zfs*q))/q; % zfs is in pixels, are values are in pixels right now +for i = 1:numel(zs)-1 % Attempt to grab the magnification curve of stage versus fit + id = zf_um >= zs(i)*q & zf_um < zs(i+1)*q; + zcm = mean(zc(id)); + zstd = std(zc(id)); + id2 = id & zc < zcm + zstd & zc > zcm - zstd; + z0m(i) = mean(zc(id2)); + zfm(i) = mean(zf_um(id2)); +end +a = polyfit(z0m,zfm,1); +clear zfm +plot(ax,zc(indy),zf_um(indy),'.'); +hold on +plot(ax,zc(indy),a(1)*zc(indy)+a(2),'g') +legend('Data','Fit','Location','North'); +ylim(ax,[dnz,upz]) +hold off +xlabel('Stage Position'); +ylabel('Found Position'); + +% Align the PSFs by their average positions +xt = []; +yt = []; +zt = []; +for i = 1:max(pfs) + ind = pfs == i; + xt = [xt;xfs(ind)-mean(xfs(ind))]; + yt = [yt;yfs(ind)-mean(yfs(ind))]; + zt = [zt;zfs(ind)]; +end +% indy = indy & xfc.^0.5*q < 0.005 & yfc.^0.5*q < 0.005; +% indy = indy & abs(xf) < 0.5 & abs(yf) < 0.5; +% Grab average position in the aligned psfs +xsel = []; +ysel = []; +zsel = []; +for i = 1:numel(zs) - 1 % Loop over a variety of height windows for averaging of points + ind1 = zt >= zs(i) & zt <=zs(i+1); % grab all points within height window + ind1 = ind1 & (xt.^2+yt.^2).^0.5 < 0.5; % grab all points who lie within a half pixel radius of center (at this point all scans have a 0 mean in XY) + % Subsets + xts = xt(ind1); + yts = yt(ind1); + zts = zt(ind1); + rts = (xts.^2+yts.^2).^0.5; % convert subsetted xy to R + mr = mean(rts); % grab mean of R + str = std(rts); % grab stdev of R +% my = mean(yts); +% sty = std(yts); + ind2 = rts <= (mr + 0.325*str); % Inclusion criteria is based on mean and std deviation + % record subset + xsel = [xsel;xts(ind2)]; + ysel = [ysel;yts(ind2)]; + zsel = [zsel;zts(ind2)]; + % make final measurement of averaged points + xfm(i) = mean(xts(ind2)); + yfm(i) = mean(yts(ind2)); + zfm(i) = mean(zts(ind2)); +end + +splz = (-1:0.001:1)/q; +xtilt = gausssmooth(xfm,5,10); +ytilt = gausssmooth(yfm,5,10); +splx = spline(zs(1:end-1) + mean(diff(zs))/2,xtilt,splz); +sply = spline(zs(1:end-1) + mean(diff(zs))/2,ytilt,splz); + +ts = uitab(tg5,'Title','Z-Tilt w/ Spline'); +ax = axes(ts); +plot3(ax,xfm,yfm,zfm,'.','Color',[0,1,0],'MarkerSize',20); +hold on +plot3(ax,xt,yt,zt,'.','Color',[0.1,0.1,0.1],'MarkerSize',2); +plot3(ax,xsel,ysel,zsel,'.','Color',[0,0,1],'MarkerSize',8); +plot3(ax,splx,sply,splz,'Color',[1,0,0],'LineWidth',10); +hold off +xlabel('Lat-X'); +ylabel('Lat-Y'); +zlabel('Axi-Z'); +legend('Avg Pts','Spline Curve'); + +% Make Drift Correction +zf = zf_um/q; +xc = spline(zs(1:end-1) + mean(diff(zs))/2,xtilt,zf); +yc = spline(zs(1:end-1) + mean(diff(zs))/2,ytilt,zf); + +xf_c = xf - xc; +yf_c = yf - yc; +zf_c = zf/(a(1)); +tc = uitab(tg5,'Title','Corrected localizations'); +ax = axes(tc); +scatter3(ax,xf_c(indy)*q,q*yf_c(indy),q*zf_c(indy),[],psf(indy)); +colormap('jet') +xlabel('Lat-X'); +ylabel('Lat-Y'); +zlabel('Axi-Z'); +axis equal +hold off +cal.z_cal = z_cal; +cal.a = a; +cal.tilt.x = xtilt; +cal.tilt.y = ytilt; +cal.ang = ang; +cal.tilt.zs = zs; + +%% Calibration Handling +disp(['Distance between Minima is ', num2str(1000*(cal.z_cal(2)-cal.z_cal(7))),'nm']); +save('z_calib.mat','cal') diff --git a/3D-Localization/xcorrsig.m b/3D-Localization/xcorrsig.m index 08c7f84..e9faeca 100644 --- a/3D-Localization/xcorrsig.m +++ b/3D-Localization/xcorrsig.m @@ -15,8 +15,8 @@ % Spline interpolation onto the output wind = -5:5; -mid1 = round(numel(sxy1)/2); -mid2 = round(numel(sxy2)/2); +mid1 = round(numel(sxy1)/2)-10; +mid2 = round(numel(sxy2)/2)-10; ssx1 = sxy1(mid1+wind); try diff --git a/Background-Subtraction-Codes/Rolling Ball/mex_CUDA_win64.xml b/Background-Subtraction-Codes/Rolling Ball/mex_CUDA_win64.xml index ca654c2..31587c4 100644 --- a/Background-Subtraction-Codes/Rolling Ball/mex_CUDA_win64.xml +++ b/Background-Subtraction-Codes/Rolling Ball/mex_CUDA_win64.xml @@ -37,7 +37,7 @@ COMPILER="nvcc" COMPFLAGS="--compiler-options=/Zp8,/GR,/W3,/EHs,/nologo,/MD $ARCHFLAGS" - ARCHFLAGS="-gencode=arch=compute_20,code=sm_20 -gencode=arch=compute_30,code=sm_30 -gencode=arch=compute_50,code=\"sm_50,compute_50\"" + ARCHFLAGS="" COMPDEFINES="--compiler-options=/D_CRT_SECURE_NO_DEPRECATE,/D_SCL_SECURE_NO_DEPRECATE,/D_SECURE_SCL=0,$MATLABMEX" MATLABMEX="/DMATLAB_MEX_FILE" OPTIMFLAGS="--compiler-options=/O2,/Oy-,/DNDEBUG" diff --git a/DBScan/DB_scan.m b/DBScan/DB_scan.m index 9ade6d5..00fffb8 100644 --- a/DBScan/DB_scan.m +++ b/DBScan/DB_scan.m @@ -2,6 +2,17 @@ % This function applies the DB scan algorithm C = 1; +<<<<<<< HEAD +clust = (1:numel(data(:,1)))*0; + +% scatter(data(:,1),data(:,2),2,[0 0 0],'Filled') +count = 1; +for i = 1:numel(data(:,1)) + + if clust(i) == 0 % If we haven't visited the point yet, perform analysis + N = func_range_scan(data(i,:),data,eps); % Find mols in neighborhood + N = N(:); % N is an array of molecules in neighborhood +======= visit = (1:numel(data(:,1)))*0; clust = visit; scatter(data(:,1),data(:,2),2,[0 0 0],'Filled') @@ -11,10 +22,44 @@ visit(i) = 1; % mark as visited N = func_range_scan(data(i,:),data,eps); % Find mols in neighborhood N = N(:); +>>>>>>> master if numel(N) >= minp % noise criteria clust(i) = C; C = C + 1; +<<<<<<< HEAD + + + test = clust(N) == 0; % see if anything in neighborhood is unvisited + while sum(test) > 0 % while there are unvisited points + ind = find(clust(N) == 0); + ind2 = clust(N) == -1; % Incorporate all 'noise' variables into the cluster + clust(ind2) = C; + ID = N(ind(1)); + + ns = func_range_scan(data(ID,:),data,eps); + if clust(ID) == 0 % visit criteria + if numel(ns) >= minp % noise criteria + N = unique([N;ns(:)]); + clust(ID) = C; % mark point as apart of cluster + else + clust(ID) = -1; % mark point as apart of noise + end + end + + test = clust(N) == 0; + end + else % Minimum number of points not satified + clust(i) = -1; % Set cluster to 'noise' cluster + end + end +end +% movie2gif(M,'C:\Users\AJN Lab\Dropbox\Lab Meetings AJN\July 22 2019\DB_example.gif','DelayTime', 0.03,'LoopCount',Inf); +end + +function ids = func_range_scan(me, them, eps) +% Algorithm to find molecules in neighborhood +======= hold on scatter(data(i,1),data(i,2),5,clust(i),'Filled'); scatter(data(N,1),data(N,2),3,ones(numel(N),1)*clust(i),'Filled'); @@ -46,6 +91,7 @@ end function ids = func_range_scan(me, them, eps) +>>>>>>> master dists = them(:,1)*0; for i = 1:numel(them(1,:)) dists = dists + (them(:,i) - me(i)).^2; diff --git a/DBScan/DB_scan_old.m b/DBScan/DB_scan_old.m new file mode 100644 index 0000000..040a2ca --- /dev/null +++ b/DBScan/DB_scan_old.m @@ -0,0 +1,55 @@ +function [clust] = DB_scan(data, eps, minp) +% This function applies the DB scan algorithm + +C = 1; +visit = (1:numel(data(:,1)))*0; +clust = visit; +% scatter(data(:,1),data(:,2),2,[0 0 0],'Filled') +count = 1; +for i = 1:numel(data(:,1)) + + if visit(i) == 0 % If we haven't visited the point yet, perform analysis + visit(i) = 1; % mark as visited + N = func_range_scan(data(i,:),data,eps); % Find mols in neighborhood + N = N(:); + + if numel(N) >= minp % noise criteria + clust(i) = C; + C = C + 1; +% hold on +% scatter(data(i,1),data(i,2),5,clust(i),'Filled'); +% scatter(data(N,1),data(N,2),3,ones(numel(N),1)*clust(i),'Filled'); + + test = visit(N) == 0; % see if anything in neighborhood is unvisited + while sum(test) > 0 % while there are unvisited points + ind = find(visit(N) == 0); + ID = N(ind(1)); + visit(ID) = 1; % mark point as visited + ns = func_range_scan(data(ID,:),data,eps); +% scatter(data(ID,1),data(ID,2),3,[1 0 0]) +% scatter(data(ns,1),data(ns,2),3,ones(numel(ns),1)*clust(i),'Filled'); +% M(count) = getframe(gcf); + count = count +1; +% drawnow + if numel(ns) >= minp % noise criteria + N = unique([N;ns(:)]) + + end + if clust(ID) == 0 + clust(ID) = clust(i); + end + test = visit(N) == 0; + end + end + end +end +% movie2gif(M,'C:\Users\AJN Lab\Dropbox\Lab Meetings AJN\July 22 2019\DB_example.gif','DelayTime', 0.03,'LoopCount',Inf); +end + +function ids = func_range_scan(me, them, eps) +dists = them(:,1)*0; +for i = 1:numel(them(1,:)) + dists = dists + (them(:,i) - me(i)).^2; +end +ids = find(dists <= eps^2); +end \ No newline at end of file diff --git a/Functions/Overhead Codes/cuda_compile.m b/Functions/Overhead Codes/cuda_compile.m index 2f1fb00..46353c6 100644 --- a/Functions/Overhead Codes/cuda_compile.m +++ b/Functions/Overhead Codes/cuda_compile.m @@ -1,4 +1,7 @@ function cuda_compile(fname) clc -ipaths = {['-I' fullfile('C:\ProgramData\NVIDIA Corporation\CUDA Samples\v8.0\common\inc')]}; +ipaths = {['-I' fullfile('C:\ProgramData\NVIDIA Corporation\CUDA Samples\v10.0\common\inc')]}; +% if exist('mex_CUDA_win64.xml') ==2 +% copyfile(['C:\Users\andre\Desktop\Code Development\Matlab Testing Folder\Neural Quhzx\Source Codes\Open this after you read the text file\mex_CUDA_win64.xml']); +% end mex(ipaths{:},fname) \ No newline at end of file diff --git a/Hurricane/Hurricane Tolerance/H2_tol.m b/Hurricane/Hurricane Tolerance/H2_tol.m index 8dd374b..98a0e5f 100644 --- a/Hurricane/Hurricane Tolerance/H2_tol.m +++ b/Hurricane/Hurricane Tolerance/H2_tol.m @@ -49,7 +49,14 @@ if tol.r.flims(2) == -1 tol.r.flims(2) = max(cdata.red.framenumber); +<<<<<<< HEAD + try tol.o.flims(2) = max(cdata.orange.framenumber); + catch + end +======= + tol.o.flims(2) = max(cdata.orange.framenumber); +>>>>>>> master end save('Tolfile.mat','tol'); @@ -61,16 +68,28 @@ for i = 1:numel(fnames) cdata.red.(fnames{i})(ind,:) = []; end +<<<<<<< HEAD +try + cdata.orange.xf = cdata.orange.xf(:); + cdata.orange.yf = cdata.orange.yf(:); + ind = cdata.orange.xf < 0 | cdata.orange.xf > 200 | cdata.orange.yf < 0 | cdata.orange.yf > 200; +======= cdata.orange.xf = cdata.orange.xf(:); cdata.orange.yf = cdata.orange.yf(:); ind = cdata.orange.xf < 0 | cdata.orange.xf > 200 | cdata.orange.yf < 0 | cdata.orange.yf > 200; +>>>>>>> master fnames = fieldnames(cdata.orange); for i = 1:numel(fnames) cdata.orange.(fnames{i})(ind,:) = []; end +<<<<<<< HEAD +catch +end +======= +>>>>>>> master %Define fractionals cdata.red.fr_N = cdata.red.crlbs(:,3).^0.5./cdata.red.fits(:,3); @@ -80,7 +99,11 @@ cdata.red.ilv = cdata.red.llv(:)./cdata.red.fits(:,3); cdata.red.eps = abs(cdata.red.fits(:,4)./cdata.red.fits(:,5)); cdata.red.snr =(cdata.red.fits(:,3)./(cdata.red.fits(:,3)+(2*6+1)^2*cdata.red.fits(:,6)).^0.5); +<<<<<<< HEAD +try +======= +>>>>>>> master cdata.orange.fr_N = cdata.orange.crlbs(:,3).^0.5./cdata.orange.fits(:,3); cdata.orange.fr_sx = cdata.orange.crlbs(:,4).^0.5./cdata.orange.fits(:,4); cdata.orange.fr_sy = cdata.orange.crlbs(:,5).^0.5./cdata.orange.fits(:,5); @@ -88,7 +111,12 @@ cdata.orange.ilv = cdata.orange.llv(:)./cdata.orange.fits(:,3); cdata.orange.eps = abs(cdata.orange.fits(:,4)./cdata.orange.fits(:,5)); cdata.orange.snr =(cdata.orange.fits(:,3)./(cdata.orange.fits(:,3)+(2*6+1)^2*cdata.orange.fits(:,6)).^0.5); +<<<<<<< HEAD +catch +end +======= +>>>>>>> master @@ -109,7 +137,11 @@ rind = ind; % Index represents all molecules that PASS tolerances % Scream here for sanity +<<<<<<< HEAD +try +======= +>>>>>>> master % Build Orange ind = cdata.orange.fits(:,3) > tol.o.N_tol(1) & cdata.orange.fits(:,3) < tol.o.N_tol(2); % Photon Tolerance ind = ind & cdata.orange.snr >= tol.o.minsnr; % Signal to Noise tolerances @@ -122,7 +154,12 @@ ind = ind & cdata.orange.fr_N < tol.o.frac_lim & abs(cdata.orange.fr_o) < tol.o.off_frac; % Fraction photon tolerance ind = ind & cdata.orange.fr_sx < tol.o.frac_lim & cdata.orange.fr_sy < tol.o.frac_lim; % Fraction width tolerance oind = ind; % Index represents all molecules that PASS tolerances +<<<<<<< HEAD +catch +end +======= +>>>>>>> master % save('Tolfile.mat','flims','minsnr','zlims','lat_max','N_tol','s_tol','iln','frac_lim','off_frac','offlim'); notind = logical(1-ind); @@ -130,7 +167,14 @@ r = str2num(fname(strfind(fname,'_r')+3)); r=0; zfr = func_shift_correct(cdata.red.zf*q,cdata.red.framenumber,r); +<<<<<<< HEAD +try +zfo = func_shift_correct(cdata.orange.zf*q,cdata.orange.framenumber,r); +catch +end +======= zfo = func_shift_correct(cdata.orange.zf*q,cdata.orange.framenumber,r); +>>>>>>> master % zf = ncoords(:,3)*q; f = figure; @@ -141,7 +185,14 @@ ax = axes(t21); plot3(ax, cdata.red.xf*q,cdata.red.yf*q,cdata.red.zf*q,'.r','MarkerSize',mwidth); hold on +<<<<<<< HEAD +try +plot3(ax, cdata.orange.xf*q,cdata.orange.yf*q,cdata.orange.zf*q,'.b','MarkerSize',mwidth); +catch +end +======= plot3(ax, cdata.orange.xf*q,cdata.orange.yf*q,cdata.orange.zf*q,'.b','MarkerSize',mwidth); +>>>>>>> master legend('Red','Orange') hold off % s.MarkerFaceColor = s.MarkerEdgeColor; @@ -153,7 +204,14 @@ ax = axes(t2); plot(ax, cdata.red.xf(rind),cdata.red.yf(rind),'.r','MarkerSize',mwidth); hold on +<<<<<<< HEAD +try +plot(ax, cdata.orange.xf(oind),cdata.orange.yf(oind),'.b','MarkerSize',mwidth); +catch +end +======= plot(ax, cdata.orange.xf(oind),cdata.orange.yf(oind),'.b','MarkerSize',mwidth); +>>>>>>> master legend('Red','Orange') hold off xlabel(ax,'microns'); @@ -163,7 +221,14 @@ ax = axes(t3); plot3(ax, cdata.red.xf(rind)*q,cdata.red.yf(rind)*q,zfr(rind),'.r','MarkerSize',mwidth); hold on +<<<<<<< HEAD +try plot3(ax, cdata.orange.xf(oind)*q,cdata.orange.yf(oind)*q,zfo(oind),'.b','MarkerSize',mwidth); +catch +end +======= +plot3(ax, cdata.orange.xf(oind)*q,cdata.orange.yf(oind)*q,zfo(oind),'.b','MarkerSize',mwidth); +>>>>>>> master legend('Red','Orange') hold off xlabel(ax,'microns'); @@ -180,7 +245,14 @@ ax = axes(t21); histogram(ax,cdata.red.zf(rind)*q); hold on +<<<<<<< HEAD +try histogram(ax,cdata.orange.zf(oind)*q); +catch +end +======= +histogram(ax,cdata.orange.zf(oind)*q); +>>>>>>> master legend('Red','Orange') xlabel('Z-position(um)') ylabel('Frequency') @@ -190,7 +262,14 @@ ax = axes(t22); histogram(ax,cdata.red.fits(rind,3)); hold on +<<<<<<< HEAD +try histogram(ax,cdata.orange.fits(oind,3)); +catch +end +======= +histogram(ax,cdata.orange.fits(oind,3)); +>>>>>>> master hold off legend('Red','Orange') xlabel('N (photons)') @@ -202,8 +281,16 @@ histogram(ax,cdata.red.fits(rind,4)); hold on histogram(ax,cdata.red.fits(rind,5)); +<<<<<<< HEAD +try +histogram(ax,cdata.orange.fits(oind,4)); +histogram(ax,cdata.orange.fits(oind,5)); +catch +end +======= histogram(ax,cdata.orange.fits(oind,4)); histogram(ax,cdata.orange.fits(oind,5)); +>>>>>>> master hold off legend(ax,'Red \sigma_x','Red \sigma_y', 'Orange \sigma_x','Orange \sigma_y') xlabel(ax,'\sigma (pixels)') diff --git a/Hurricane/Hurricane Tolerance/H2_tol.mlx b/Hurricane/Hurricane Tolerance/H2_tol.mlx new file mode 100644 index 0000000..f265714 Binary files /dev/null and b/Hurricane/Hurricane Tolerance/H2_tol.mlx differ diff --git a/Hurricane/Hurricane.m b/Hurricane/Hurricane.m index a00f8d3..5ef694d 100644 --- a/Hurricane/Hurricane.m +++ b/Hurricane/Hurricane.m @@ -19,7 +19,7 @@ an_dir = 'Analysis'; % name of analysis directory angle = 0; %approximate astig rotation in degrees sv_im = 'n'; % set to y/Y to save image of localizations -thresh = 2; +thresh = 10; %% Optionals % This section is dedicated to a list of variables for the user to select @@ -36,11 +36,11 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -p = mfilename('fullpath'); -[fpath, fname, fext] = fileparts([p, '.m']); -addpath(fpath); -addpath([fpath,'\da_c']); -addpath([fpath,'\da_functions']); +% p = mfilename('fullpath'); +% [fpath, fname, fext] = fileparts([p, '.m']); +% addpath(fpath); +% addpath([fpath,'\da_c']); +% addpath([fpath,'\da_functions']); % [dname, dpath] = uigetfile('*.tif'); % cd(dpath); % we change to the path we are going to use @@ -67,13 +67,15 @@ elseif varys(4) == 1 mkdir('Rolling_Ball'); end -for i = 1:numel(files) - tic - func_da_storm(files(i).name, dpath, an_dir, q, pix2pho, pixw,thresh, angle, sv_im, mi1, varys); - clc; - disp(['File number ' , num2str(i) , ' out of ', num2str(numel(files))]); - t(i) = toc; - ajn_wait(t,i,numel(files)); +for i = 3:numel(files) +% if isempty(strfind(files(i).name,'scan')) + tic + func_da_storm(files(i).name, dpath, an_dir, q, pix2pho, pixw,thresh, angle, sv_im, mi1, varys); + clc; + % disp(['File number ' , num2str(i) , ' out of ', num2str(numel(files))]); + t(i) = toc; +% end +% ajn_wait(t,i,numel(files)); close all end diff --git a/Hurricane/Hurricane_live_script.mlx b/Hurricane/Hurricane_live_script.mlx new file mode 100644 index 0000000..c034b82 Binary files /dev/null and b/Hurricane/Hurricane_live_script.mlx differ diff --git a/Hurricane/hurricane_functions/func_da_storm.m b/Hurricane/hurricane_functions/func_da_storm.m index 1574b5f..9e2c478 100644 --- a/Hurricane/hurricane_functions/func_da_storm.m +++ b/Hurricane/hurricane_functions/func_da_storm.m @@ -4,7 +4,7 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ % pix2pho = single(pix2pho); q = double(q); try - load('C:\Users\AJN Lab\Documents\GitHub\Matlab-Master\Hurricane\hurricane_functions\z_calib.mat') + load('C:\Users\andre\Documents\GitHub\Matlab-Master\Hurricane\hurricane_functions\z_calib.mat') % if exist([data_d, 'z_calib.mat']) % cal = load([data_d 'z_calib.mat']); % else @@ -45,8 +45,14 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ % If we're doing 2 color, block out frame we're not interested in if choices(5) == 1 - load('C:\Users\AJN Lab\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); ifind = func_image_block(ifind,split); + elseif choices(5) == 2 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block_orange(ifind,split); + elseif choices(5) == 2 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block_orange(ifind,split); end % ifind = func_image_block2(ifind,split); % thrsh = thresh/100*mean(max(max(iprod))); @@ -99,9 +105,9 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'pixw','q','ncoords','fits','crlbs','llv','framenumber','cal'); else % load('C:\Users\AJN Lab\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); - id = cents(:,1) < 180; % Identify localizations below the split + id = cents(:,1) < split; % Identify localizations below the split %% First fit is all red, so those can be immediately - if sum(id)>0 + if sum(id)>0&& choices(5) ~= 3 [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.red.ang); % As everywhere in the equations used sigma is squared, we can without @@ -110,11 +116,13 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ fits(:,5) = abs(fits(:,5)); % Put data into cdata structure - cdata.red.fits = fits; - cdata.red.crlbs = crlbs; + for i = 1:6 + cdata.red.fits(:,i) = fits(:,i); + cdata.red.crlbs(:,i) = crlbs(:,i); + end cdata.red.llv = llv; cdata.red.framenumber = framenumber; - save + % Z calculations zf = getdz(cdata.red.fits(:,4),cdata.red.fits(:,5),cal.red.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere ncoords = astig_tilt([cdata.red.fits(:,1:2),zf],cal.red); % corrections due to astigmatism @@ -129,7 +137,7 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ end %% Repeat above for orange id = logical(1-id); % Changes 0 -> 1 and 1 -> 0 flipping the ID so now we can fit orange - if sum(id) >0 + if sum(id) >0 && choices(5) ~= 2 [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.orange.ang); % As everywhere in the equations used sigma is squared, we can without % loss of generality make these fits positive definite @@ -137,8 +145,10 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ fits(:,5) = abs(fits(:,5)); % Put data into cdata structure - cdata.orange.fits = fits; - cdata.orange.crlbs = crlbs; + for i = 1:6 + cdata.orange.fits(:,i) = fits(:,i); + cdata.orange.crlbs(:,i) = crlbs(:,i); + end cdata.orange.llv = llv; cdata.orange.framenumber = framenumber; @@ -149,15 +159,43 @@ function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_ x = o2rx.'*vec.'; y = o2ry.'*vec.'; % Assign fixed coordinates - cdata.orange.xf = x.'; - cdata.orange.yf = y.'; + cdata.orange.xf = ncoords(:,1); + cdata.orange.yf = ncoords(:,2); cdata.orange.zf = ncoords(:,3); + cal.o2rx = o2rx; + cal.o2ry = o2ry; end + % Remove problem entries + orange_index = cdata.orange.xf < 0 | cdata.orange.yf <0 ; + red_index = cdata.red.xf < 0 | cdata.red.yf <0 ; + + + field_names = fieldnames(cdata.red); + + + if sum(red_index)>0 + cdata.red.fits(red_index,:) = []; + cdata.red.crlbs(red_index,:) = []; + end + if sum(orange_index)>0 + cdata.orange.fits(orange_index,:) = []; + cdata.orange.crlbs(orange_index,:) = []; + end + for k=3:numel(field_names) + if sum(red_index) > 0 + cdata.red.(field_names{k})(red_index) = []; + end + if sum(orange_index) > 0 + cdata.orange.(field_names{k})(orange_index) = []; + end + end + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cdata', 'pixw','q','cal'); end catch lsterr disp(lsterr) + waitforbuttonpress end % save('results_of_bump.mat','fnum','q','iloc','cal','cents'); diff --git a/Hurricane/hurricane_functions/func_da_storm_with_working_xfo.m b/Hurricane/hurricane_functions/func_da_storm_with_working_xfo.m new file mode 100644 index 0000000..1df750a --- /dev/null +++ b/Hurricane/hurricane_functions/func_da_storm_with_working_xfo.m @@ -0,0 +1,205 @@ +function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_im, mi1, choices) + +% Convert Variabls +% pix2pho = single(pix2pho); +q = double(q); +try + load('C:\Users\andre\Documents\GitHub\Matlab-Master\Hurricane\hurricane_functions\z_calib.mat') + % if exist([data_d, 'z_calib.mat']) + % cal = load([data_d 'z_calib.mat']); + % else + % cal = load('z_calib.mat'); + % end + % mi1 = 0 + % Load file and dark current background subtraction + i1 = (readtiff(fname) - mi1); + % i1 = sum(i1,3); + % i1 = i1.*(i1>0); + [m,n,o] = size(i1); + % i1(1:30,:,:) = 0; + % i1(m-30:m,:,:) = 0; + % i1(:,1:30,:) = 0; + % i1(:,n-30:n,:) = 0; + % Rolling Ball Background Subtract + % iprod = rollingball(i1); + iprod = gpu_rball(i1); + if choices(4) == 1 + writetiff(iprod,[data_d,'\Rolling_Ball\',fname(1:end-4),'_rb.tif']); + end + % iprod = bp_subtract(i1); + % iprod = imgaussfilt(i1,0.8947); + % iprod = i1; + % Peak Detection + + + % thrsh = 300/pix2pho; + % diprod = diff(iprod,3); + % for i = 1:o + % ifind = denoise_psf(iprod,2); + iwaves = gpu_waves(iprod); + se = strel('Disk',1); + ifind = imerode(iwaves,se); + if choices(1) == 1 + writetiff(ifind,[data_d,'\Waves\',fname(1:end-4),'_waves.tif']); + end + + % If we're doing 2 color, block out frame we're not interested in + if choices(5) == 1 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block(ifind,split); + elseif choices(5) == 2 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block_orange(ifind,split); + elseif choices(5) == 2 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block_orange(ifind,split); + end +% ifind = func_image_block2(ifind,split); + % thrsh = thresh/100*mean(max(max(iprod))); + % tic + % thrsh = 3*std(iprod(:)) + mean(iprod(:)); + % toc + % thrsh = thresh/100*mean(max(max(diprod))); + % surf(max(ifind,[],3)); + % thrsh = input('What should the threshold be? '); + % thrsh = min(iprod(:))*thresh/100; + % excerpt = 1; + dps = cpu_peaks(ifind,thresh,pixw); + if choices(2) == 1 + in_d_eye(iprod, dps, pixw); + end + + clear ip ipf i1 + + % divide up the data + + [iloc, fnum, cents] = divide_up(iprod, pixw, dps); + [m,n,o] = size(iloc); + % remove duplicate data + [ind] = find_fm_dupes(cents,fnum,pixw*1.5); + iloc(:,:,ind) = []; + cents(ind,:) = []; + fnum(ind) = []; + + + + % Localize the Data + % [xf_all,xf_crlb, yf_all,yf_crlb,sigx_all, sigx_crlb, sigy_all, sigy_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, y, inloc, xin, yin] = da_locs_sigs(iloc, fnum, cents, angle); + % zf_all = getdz(sigx_all,sigy_all)/q; + % [xf_all,xf_crlb, yf_all,yf_crlb,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, y, inloc, xin, yin] = da_locs(iloc, fnum, cents, angle);zf_all = zf_all/q; % This is to handle Z informtation uncomment once calibration is fixed + % [xf_all,xf_crlb, yf_all,yf_crlb,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, iters, cex, cey] = da_splines(iloc, fnum, cents, cal, pixw); + % [~,~, ~,~,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, iters, cex, cey] = da_splines(iloc, fnum, cents, cal, pixw); + % i2 = reshape(iloc,m*n,o); + % save('thisbit.mat','iloc','cents','fnum','cal'); + if choices(3) == 1 + writetiff(iloc,[data_d,'\psfs\',fname(1:end-4),'_psfs.tif']); + end + if choices(5) == 0 + [fits, crlbs, llv, framenumber] = slim_locs(iloc, fnum, cents, cal.red.ang); + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + zf = getdz(abs(fits(:,4)),abs(fits(:,5)),cal.red.z_cal,2)/q; + coords = [fits(:,1:2),zf]; + [ncoords] = astig_tilt(coords,cal.red); + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'pixw','q','ncoords','fits','crlbs','llv','framenumber','cal'); + else + % load('C:\Users\AJN Lab\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + id = cents(:,1) < split; % Identify localizations below the split + %% First fit is all red, so those can be immediately + if sum(id)>0&& choices(5) ~= 3 + [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.red.ang); + + % As everywhere in the equations used sigma is squared, we can without + % loss of generality make these fits positive definite + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + % Put data into cdata structure + cdata.red.fits = fits; + cdata.red.crlbs = crlbs; + cdata.red.llv = llv; + cdata.red.framenumber = framenumber; + + % Z calculations + zf = getdz(cdata.red.fits(:,4),cdata.red.fits(:,5),cal.red.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere + ncoords = astig_tilt([cdata.red.fits(:,1:2),zf],cal.red); % corrections due to astigmatism + + % Assign fixed coordinates + cdata.red.xf = ncoords(:,1); + cdata.red.yf = ncoords(:,2); + cdata.red.zf = ncoords(:,3); + + + clear fits crlbs llv framenumber + end + %% Repeat above for orange + id = logical(1-id); % Changes 0 -> 1 and 1 -> 0 flipping the ID so now we can fit orange + if sum(id) >0 && choices(5) ~= 2 + [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.orange.ang); + % As everywhere in the equations used sigma is squared, we can without + % loss of generality make these fits positive definite + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + % Put data into cdata structure + cdata.orange.fits = fits; + cdata.orange.crlbs = crlbs; + cdata.orange.llv = llv; + cdata.orange.framenumber = framenumber; + + % Z calculations + zf = getdz(cdata.orange.fits(:,4),cdata.orange.fits(:,5),cal.orange.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere + ncoords = astig_tilt([cdata.orange.fits(:,1:2),zf],cal.orange); % corrections due to astigmatism + vec = xy_feature(ncoords(:,1),ncoords(:,2)); + x = o2rx.'*vec.'; + y = o2ry.'*vec.'; + % Assign fixed coordinates + cdata.orange.xf = x.'; + cdata.orange.yf = y.'; + cdata.orange.xfo = ncoords(:,1); + cdata.orange.yfo = ncoords(:,2); + cdata.orange.zf = ncoords(:,3); + end + % Remove problem entries + orange_index = cdata.orange.xf < 0 | cdata.orange.yf <0 ; + red_index = cdata.red.xf < 0 | cdata.red.yf <0 ; + + + field_names = fieldnames(cdata.red); + for k=1:numel(field_names) + cdata.red.(field_names{k}) (red_index) = []; + cdata.orange.(field_names{k}) (orange_index) = []; + end + cdata.orange.xfo(orange_index) = []; + cdata.orange.yfo(orange_index) = []; + + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cdata', 'pixw','q','cal'); + end + +catch lsterr + disp(lsterr) + waitforbuttonpress +end +% save('results_of_bump.mat','fnum','q','iloc','cal','cents'); + +% save('for_trial.mat','iloc' + +% Save the Analysis +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'zf_all','sigx_all' ,'sigy_all','sigx_crlb','sigy_crlb','y','iloc','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','pixw','q','pix2pho'); +% if strcmp(sv_im,'Y') || strcmp(sv_im,'y') +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cents','zf_all','zf_crlb','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','iters','pixw','q','pix2pho','ilocs'); +% else +% figure +% imagesc(mean(iprod,3)); +% hold on +% plot(fits(:,1),fits(:,2),'rx') +% hold off +% colormap('gray'); + + +% end +% catch lsterr +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'zf_all','sigx_all' ,'sigy_all','sigx_crlb','sigy_crlb','y','iloc','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','pixw','q','pix2pho'); +end + diff --git a/Hurricane/hurricane_functions/func_da_storm_with_xfo.m b/Hurricane/hurricane_functions/func_da_storm_with_xfo.m new file mode 100644 index 0000000..548a40f --- /dev/null +++ b/Hurricane/hurricane_functions/func_da_storm_with_xfo.m @@ -0,0 +1,208 @@ +function func_da_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_im, mi1, choices) + +% Convert Variabls +% pix2pho = single(pix2pho); +q = double(q); +try + load('C:\Users\andre\Documents\GitHub\Matlab-Master\Hurricane\hurricane_functions\z_calib.mat') + % if exist([data_d, 'z_calib.mat']) + % cal = load([data_d 'z_calib.mat']); + % else + % cal = load('z_calib.mat'); + % end + % mi1 = 0 + % Load file and dark current background subtraction + i1 = (readtiff(fname) - mi1); + % i1 = sum(i1,3); + % i1 = i1.*(i1>0); + [m,n,o] = size(i1); + % i1(1:30,:,:) = 0; + % i1(m-30:m,:,:) = 0; + % i1(:,1:30,:) = 0; + % i1(:,n-30:n,:) = 0; + % Rolling Ball Background Subtract + % iprod = rollingball(i1); + iprod = gpu_rball(i1); + if choices(4) == 1 + writetiff(iprod,[data_d,'\Rolling_Ball\',fname(1:end-4),'_rb.tif']); + end + % iprod = bp_subtract(i1); + % iprod = imgaussfilt(i1,0.8947); + % iprod = i1; + % Peak Detection + + + % thrsh = 300/pix2pho; + % diprod = diff(iprod,3); + % for i = 1:o + % ifind = denoise_psf(iprod,2); + iwaves = gpu_waves(iprod); + se = strel('Disk',1); + ifind = imerode(iwaves,se); + if choices(1) == 1 + writetiff(ifind,[data_d,'\Waves\',fname(1:end-4),'_waves.tif']); + end + + % If we're doing 2 color, block out frame we're not interested in + if choices(5) == 1 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block(ifind,split); + elseif choices(5) == 2 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block_orange(ifind,split); + elseif choices(5) == 2 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block_orange(ifind,split); + end +% ifind = func_image_block2(ifind,split); + % thrsh = thresh/100*mean(max(max(iprod))); + % tic + % thrsh = 3*std(iprod(:)) + mean(iprod(:)); + % toc + % thrsh = thresh/100*mean(max(max(diprod))); + % surf(max(ifind,[],3)); + % thrsh = input('What should the threshold be? '); + % thrsh = min(iprod(:))*thresh/100; + % excerpt = 1; + dps = cpu_peaks(ifind,thresh,pixw); + if choices(2) == 1 + in_d_eye(iprod, dps, pixw); + end + + clear ip ipf i1 + + % divide up the data + + [iloc, fnum, cents] = divide_up(iprod, pixw, dps); + [m,n,o] = size(iloc); + % remove duplicate data + [ind] = find_fm_dupes(cents,fnum,pixw*1.5); + iloc(:,:,ind) = []; + cents(ind,:) = []; + fnum(ind) = []; + + + + % Localize the Data + % [xf_all,xf_crlb, yf_all,yf_crlb,sigx_all, sigx_crlb, sigy_all, sigy_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, y, inloc, xin, yin] = da_locs_sigs(iloc, fnum, cents, angle); + % zf_all = getdz(sigx_all,sigy_all)/q; + % [xf_all,xf_crlb, yf_all,yf_crlb,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, y, inloc, xin, yin] = da_locs(iloc, fnum, cents, angle);zf_all = zf_all/q; % This is to handle Z informtation uncomment once calibration is fixed + % [xf_all,xf_crlb, yf_all,yf_crlb,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, iters, cex, cey] = da_splines(iloc, fnum, cents, cal, pixw); + % [~,~, ~,~,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, iters, cex, cey] = da_splines(iloc, fnum, cents, cal, pixw); + % i2 = reshape(iloc,m*n,o); + % save('thisbit.mat','iloc','cents','fnum','cal'); + if choices(3) == 1 + writetiff(iloc,[data_d,'\psfs\',fname(1:end-4),'_psfs.tif']); + end + if choices(5) == 0 + [fits, crlbs, llv, framenumber] = slim_locs(iloc, fnum, cents, cal.red.ang); + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + zf = getdz(abs(fits(:,4)),abs(fits(:,5)),cal.red.z_cal,2)/q; + coords = [fits(:,1:2),zf]; + [ncoords] = astig_tilt(coords,cal.red); + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'pixw','q','ncoords','fits','crlbs','llv','framenumber','cal'); + else + % load('C:\Users\AJN Lab\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + id = cents(:,1) < split; % Identify localizations below the split + %% First fit is all red, so those can be immediately + if sum(id)>0&& choices(5) ~= 3 + [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.red.ang); + + % As everywhere in the equations used sigma is squared, we can without + % loss of generality make these fits positive definite + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + % Put data into cdata structure + for i = 1:6 + cdata.red.fits(:,i) = fits(:,i); + cdata.red.crlbs(:,i) = crlbs(:,i); + end + cdata.red.llv = llv; + cdata.red.framenumber = framenumber; + + % Z calculations + zf = getdz(cdata.red.fits(:,4),cdata.red.fits(:,5),cal.red.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere + ncoords = astig_tilt([cdata.red.fits(:,1:2),zf],cal.red); % corrections due to astigmatism + + % Assign fixed coordinates + cdata.red.xf = ncoords(:,1); + cdata.red.yf = ncoords(:,2); + cdata.red.zf = ncoords(:,3); + + + clear fits crlbs llv framenumber + end + %% Repeat above for orange + id = logical(1-id); % Changes 0 -> 1 and 1 -> 0 flipping the ID so now we can fit orange + if sum(id) >0 && choices(5) ~= 2 + [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.orange.ang); + % As everywhere in the equations used sigma is squared, we can without + % loss of generality make these fits positive definite + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + % Put data into cdata structure + for i = 1:6 + cdata.orange.fits(:,i) = fits(:,i); + cdata.orange.crlbs(:,i) = crlbs(:,i); + end + cdata.orange.llv = llv; + cdata.orange.framenumber = framenumber; + + % Z calculations + zf = getdz(cdata.orange.fits(:,4),cdata.orange.fits(:,5),cal.orange.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere + ncoords = astig_tilt([cdata.orange.fits(:,1:2),zf],cal.orange); % corrections due to astigmatism + vec = xy_feature(ncoords(:,1),ncoords(:,2)); + x = o2rx.'*vec.'; + y = o2ry.'*vec.'; + % Assign fixed coordinates + cdata.orange.xf = x.'; + cdata.orange.yf = y.'; + cal.o2rx = o2rx; + cal.o2ry = o2ry; + end + % Remove problem entries + orange_index = cdata.orange.xf < 0 | cdata.orange.yf <0 ; + red_index = cdata.red.xf < 0 | cdata.red.yf <0 ; + + + field_names = fieldnames(cdata.red); + for k=1:numel(field_names) + cdata.red.(field_names{k}) (red_index) = []; + cdata.orange.(field_names{k}) (orange_index) = []; + end + cdata.orange.xfo(orange_index) = []; + cdata.orange.yfo(orange_index) = []; + + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cdata', 'pixw','q','cal'); + end + +catch lsterr + disp(lsterr) + waitforbuttonpress +end +% save('results_of_bump.mat','fnum','q','iloc','cal','cents'); + +% save('for_trial.mat','iloc' + +% Save the Analysis +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'zf_all','sigx_all' ,'sigy_all','sigx_crlb','sigy_crlb','y','iloc','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','pixw','q','pix2pho'); +% if strcmp(sv_im,'Y') || strcmp(sv_im,'y') +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cents','zf_all','zf_crlb','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','iters','pixw','q','pix2pho','ilocs'); +% else +% figure +% imagesc(mean(iprod,3)); +% hold on +% plot(fits(:,1),fits(:,2),'rx') +% hold off +% colormap('gray'); + + +% end +% catch lsterr +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'zf_all','sigx_all' ,'sigy_all','sigx_crlb','sigy_crlb','y','iloc','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','pixw','q','pix2pho'); +end + diff --git a/Hurricane/hurricane_functions/func_image_block_orange.m b/Hurricane/hurricane_functions/func_image_block_orange.m new file mode 100644 index 0000000..29eac49 --- /dev/null +++ b/Hurricane/hurricane_functions/func_image_block_orange.m @@ -0,0 +1,12 @@ +function i2 = func_image_block_orange(i1,split,frame) +[m,n,o] = size(i1); +i2 = i1*0; +oblock = ones(m,n); +rblock = oblock; + +oblock(:,split:end) = 0; +rblock(:,1:split) = 0; + +for i = 1:o + i2(:,:,i) = i1(:,:,i).*oblock; +end \ No newline at end of file diff --git a/Hurricane/hurricane_functions/func_image_block_red.m b/Hurricane/hurricane_functions/func_image_block_red.m new file mode 100644 index 0000000..e94c555 --- /dev/null +++ b/Hurricane/hurricane_functions/func_image_block_red.m @@ -0,0 +1,12 @@ +function i2 = func_image_block_red(i1,split,frame) +[m,n,o] = size(i1); +i2 = i1*0; +oblock = ones(m,n); +rblock = oblock; + +oblock(:,split:end) = 0; +rblock(:,1:split) = 0; + +for i = 1:o + i2(:,:,i) = i1(:,:,i).*rblock; +end \ No newline at end of file diff --git a/Hurricane/hurricane_functions/z_calib.mat b/Hurricane/hurricane_functions/z_calib.mat index e54b151..5a7eada 100644 Binary files a/Hurricane/hurricane_functions/z_calib.mat and b/Hurricane/hurricane_functions/z_calib.mat differ diff --git a/Hurricane/hurricane_functions/z_calib_jan_to_jul_2020.mat b/Hurricane/hurricane_functions/z_calib_jan_to_jul_2020.mat new file mode 100644 index 0000000..e54b151 Binary files /dev/null and b/Hurricane/hurricane_functions/z_calib_jan_to_jul_2020.mat differ diff --git a/Minutia/Neural_Hurricane.m b/Minutia/Neural_Hurricane.m new file mode 100644 index 0000000..d36e959 --- /dev/null +++ b/Minutia/Neural_Hurricane.m @@ -0,0 +1,13 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Neural Hurricane +% +% This is a script that will be based off Neural Learning, but incorporate +% Matlab functions instead of homemade ones to improve upon learning and +% Reliability +% AJN 8-2-20 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +close all +clearvars; +clc + +%% User Controlled Variables diff --git a/Minutia/Neural_Learning.m b/Minutia/Neural_Learning.m new file mode 100644 index 0000000..ada6953 --- /dev/null +++ b/Minutia/Neural_Learning.m @@ -0,0 +1,159 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Neural Learning +% +% This script will generate a random set of neural thetas, analyze an +% image, check the quality of the data, use that to relearn neural thetas, +% and repeat the analysis. This will repeat until a currently undetermined +% point. +% +% +% +% AJN 3-12-16 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +close all +clear all +clc + +%% User Controlled Variables +sp = 50; +lambda = 0.3; +trainper = 0.9; +tossper = 0.05; +epsint = 0.12; +wvlngth = 575; % peak wavelength of fluorophore emission +NA = 1.4; +q = 0.128; % Pixel Size +n_start = 1; % Starting frame +n_end = 0; % End 0 for end frame set to 0 for final frame +n_bkgn = 100; % Frame to use for background measurement +pix2pho = 12.78; % Pixel to photon ratio (for sCMOS set to 1 after pre-processing) +n_start = 1; +n_end = 0; +bkgn = 1; + +%% Tolerance Variables +N_tol_min = 10; % minimum number of photons +N_tol_max = 5000; % maximum number of photons +off_min = -10; % minimum number of offset photons +off_max = 12; % maximum number of offset photons +N_crlb_min = 0; % minimum variance in number of photons +N_crlb_max = 15000; % maximum variance in number of photons +sigma_min = 1; % minimum sigma value +sigma_max = 10; % maximum sigma value +xf_crlb_min = 0; % minimum error in number of x-position +xf_crlb_max = 0.91; % maximum error in number of x-position +yf_crlb_min = 0; % minimum error in number of y-position +yf_crlb_max = 0.91; % maximum error in number of y-position +off_crlb_min = 0; % minimum error in number of offset photons +off_crlb_max = 2; % maximum error in number of offset photons +sig_crlb_min = 0; +sig_crlb_max = 2; +fr_unc_N = 0.5; % fractional uncertainty in N +fr_unc_off = 0.5; % Fractional uncertainty in offset +fr_unc_sig = 0.5; % Fractional uncertainty in width +x=[]; +y=[]; + +% Code Begins +% build random thetas +theta1 = rand(30, 50)*2*epsint - epsint; +theta2 = rand(1, 31)*2*epsint - epsint; + + + +% load('C:\Users\AJN Lab\Desktop\5-30-17 munc13-halo\Subfolder\Neural_thetas_it_19.mat'); +% theta1 = Theta1; +% theta2 = Theta2; +%% In principle we would choose a random file but for now +[fname, fpath] = uigetfile('*.tif'); +cd(fpath); +finfo = dir('*.tif'); +% save these thetas +% save('neural_thetas_it_0.mat','theta1','theta2'); +data_dir = fpath; +an_dir = fpath; +it = 1; +h1 = figure('Units','Normalized','Outerposition',[0 0 1 1]); +x1 = []; +y1 = []; +county = 1; + +% amem = 5.9367*10^9; % This is for the Geforce 1060 +amem = 956174336; % this is for the geforce 550 Ti +while true + % analyze file and return fitting variables + if it < 5 + n_end = 10; + sp = 2; + elseif it < 15 && it >=5 + n_end = 50; + sp = 10; + elseif it < 25 && it >=15 + n_end = 500; + sp = 50; + else + n_end = 1000; + sp = 50; + end + + % figure('Units','Normalized','OuterPosition',[0 0 1 1]); + [savename] = func_Quhzx_02_7_for_learning(data_dir,finfo(randi(numel(finfo))).name,an_dir,bkgn,theta1, theta2,q,pix2pho,n_start,n_end, wvlngth, NA, amem, sp); + % saveas(gcf,['figure_it_', num2str(it),'.tif']); + % p1 = pca; + % M(county) = getframe(h1); + % county = county +1; + % Tolerancing file + try + + + % if it < 10 + % [x1, y1] = app_gpu_tol_all_color_learning_easy(data_dir, savename, sigma_min, sigma_max, sig_crlb_min, sig_crlb_max, fr_unc_sig, N_crlb_min, N_crlb_max, xf_crlb_min, xf_crlb_max, yf_crlb_min, yf_crlb_max, off_crlb_min, off_crlb_max, N_tol_min,N_tol_max, off_min, off_max, fr_unc_N, fr_unc_off); + % elseif it >= 10 + [x1, y1] = app_gpu_tol_all_color_learning(it + 20,data_dir, savename, sigma_min, sigma_max, sig_crlb_min, sig_crlb_max, fr_unc_sig, N_crlb_min, N_crlb_max, xf_crlb_min, xf_crlb_max, yf_crlb_min, yf_crlb_max, off_crlb_min, off_crlb_max, N_tol_min,N_tol_max, off_min, off_max, fr_unc_N, fr_unc_off); + % end + %% add tole + x = [x;x1]; + y = [y;y1]; + catch lsterr + end + sum(y) + % if sum(y) < 40 && it == 1 % if there are no good fits found, put in images of gaussians + % [x,y] = func_build_gauss(x,y, N_tol_min, N_tol_max, sigma_min, sigma_max); + % end + if sum(y) > 40 || it > 1 + + %% At this point x contains a group of all images and y is a key of whether the image passed the tolerances or not + % disp(['Iteration number: ', num2str(it-1), ' and ', num2str(100*sum(y1)/numel(y1)),' % good regions found and ', num2str(sum(y1)), '# of regions found']); + progress(it,1) = it -1; + progress(it,2) = sum(y1)/numel(y1); + progress(it,3) = sum(y1); + + + % now the training set should contain positive and negative examples + if sum(y)/numel(y) < 1.99 % if neural net picked out less than 70% good regions train new set + [theta1, theta2] = func_neural_teach( theta1, theta2, x, y, lambda, it, trainper); + it = it +1; + + + tindex = find(y == 1); + tosspoint = round(tossper*numel(tindex)); + index = randperm(numel(tindex)); + if it -1 >10 + x(index(1:tosspoint),:) = []; + y(index(1:tosspoint)) = []; + end + else + break + end + % subplot(1,2,2); + % plot(progress(:,1),progress(:,2)*100,'.b','MarkerSize',30); + % xlabel('Iteration Number') + % ylabel('Percentage of Good Regions') + else + theta1 = rand(30, 82)*2*epsint - epsint; + theta2 = rand(1, 31)*2*epsint - epsint; +% [x,y] = func_build_gauss(x,y, N_tol_min, N_tol_max, sigma_min, sigma_max); + end + + % M(it-1) = getframe(gcf); +end \ No newline at end of file diff --git a/Minutia/func_neural_storm.m b/Minutia/func_neural_storm.m new file mode 100644 index 0000000..5484c67 --- /dev/null +++ b/Minutia/func_neural_storm.m @@ -0,0 +1,185 @@ +function func_neural_storm(fname,data_d, an_dir, q, pix2pho, pixw,thresh, angle, sv_im, mi1, choices) + +% Convert Variabls +% pix2pho = single(pix2pho); +q = double(q); +try + load('C:\Users\andre\Documents\GitHub\Matlab-Master\Hurricane\hurricane_functions\z_calib.mat') + % if exist([data_d, 'z_calib.mat']) + % cal = load([data_d 'z_calib.mat']); + % else + % cal = load('z_calib.mat'); + % end + % mi1 = 0 + % Load file and dark current background subtraction + i1 = (readtiff(fname) - mi1); + % i1 = sum(i1,3); + % i1 = i1.*(i1>0); + [m,n,o] = size(i1); + % i1(1:30,:,:) = 0; + % i1(m-30:m,:,:) = 0; + % i1(:,1:30,:) = 0; + % i1(:,n-30:n,:) = 0; + % Rolling Ball Background Subtract + % iprod = rollingball(i1); + iprod = gpu_rball(i1); + if choices(4) == 1 + writetiff(iprod,[data_d,'\Rolling_Ball\',fname(1:end-4),'_rb.tif']); + end + % iprod = bp_subtract(i1); + % iprod = imgaussfilt(i1,0.8947); + % iprod = i1; + % Peak Detection + + + % thrsh = 300/pix2pho; + % diprod = diff(iprod,3); + % for i = 1:o + % ifind = denoise_psf(iprod,2); + iwaves = gpu_waves(iprod); + se = strel('Disk',1); + ifind = imerode(iwaves,se); + if choices(1) == 1 + writetiff(ifind,[data_d,'\Waves\',fname(1:end-4),'_waves.tif']); + end + + % If we're doing 2 color, block out frame we're not interested in + if choices(5) == 1 + load('C:\Users\andre\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + ifind = func_image_block(ifind,split); + end +% ifind = func_image_block2(ifind,split); + % thrsh = thresh/100*mean(max(max(iprod))); + % tic + % thrsh = 3*std(iprod(:)) + mean(iprod(:)); + % toc + % thrsh = thresh/100*mean(max(max(diprod))); + % surf(max(ifind,[],3)); + % thrsh = input('What should the threshold be? '); + % thrsh = min(iprod(:))*thresh/100; + % excerpt = 1; + dps = cpu_peaks(ifind,thresh,pixw); + dps = neural_peaks(ifind, theta1, theta2); + if choices(2) == 1 + in_d_eye(iprod, dps, pixw); + end + + clear ip ipf i1 + + % divide up the data + + [iloc, fnum, cents] = divide_up(iprod, pixw, dps); + [m,n,o] = size(iloc); + % remove duplicate data + [ind] = find_fm_dupes(cents,fnum,pixw*1.5); + iloc(:,:,ind) = []; + cents(ind,:) = []; + fnum(ind) = []; + + + + % Localize the Data + % [xf_all,xf_crlb, yf_all,yf_crlb,sigx_all, sigx_crlb, sigy_all, sigy_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, y, inloc, xin, yin] = da_locs_sigs(iloc, fnum, cents, angle); + % zf_all = getdz(sigx_all,sigy_all)/q; + % [xf_all,xf_crlb, yf_all,yf_crlb,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, y, inloc, xin, yin] = da_locs(iloc, fnum, cents, angle);zf_all = zf_all/q; % This is to handle Z informtation uncomment once calibration is fixed + % [xf_all,xf_crlb, yf_all,yf_crlb,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, iters, cex, cey] = da_splines(iloc, fnum, cents, cal, pixw); + % [~,~, ~,~,zf_all, zf_crlb, N, N_crlb,off_all, off_crlb, framenum_all, llv, iters, cex, cey] = da_splines(iloc, fnum, cents, cal, pixw); + % i2 = reshape(iloc,m*n,o); + % save('thisbit.mat','iloc','cents','fnum','cal'); + if choices(3) == 1 + writetiff(iloc,[data_d,'\psfs\',fname(1:end-4),'_psfs.tif']); + end + if choices(5) == 0 + [fits, crlbs, llv, framenumber] = slim_locs(iloc, fnum, cents, cal.red.ang); + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + zf = getdz(abs(fits(:,4)),abs(fits(:,5)),cal.red.z_cal,2)/q; + coords = [fits(:,1:2),zf]; + [ncoords] = astig_tilt(coords,cal.red); + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'pixw','q','ncoords','fits','crlbs','llv','framenumber','cal'); + else + % load('C:\Users\AJN Lab\Documents\GitHub\Matlab-Master\2-Channel Codes\2_color_calibration.mat', 'split', 'o2rx','o2ry'); + id = cents(:,1) < 180; % Identify localizations below the split + %% First fit is all red, so those can be immediately + if sum(id)>0 + [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.red.ang); + + % As everywhere in the equations used sigma is squared, we can without + % loss of generality make these fits positive definite + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + % Put data into cdata structure + cdata.red.fits = fits; + cdata.red.crlbs = crlbs; + cdata.red.llv = llv; + cdata.red.framenumber = framenumber; + + % Z calculations + zf = getdz(cdata.red.fits(:,4),cdata.red.fits(:,5),cal.red.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere + ncoords = astig_tilt([cdata.red.fits(:,1:2),zf],cal.red); % corrections due to astigmatism + + % Assign fixed coordinates + cdata.red.xf = ncoords(:,1); + cdata.red.yf = ncoords(:,2); + cdata.red.zf = ncoords(:,3); + + + clear fits crlbs llv framenumber + end + %% Repeat above for orange + id = logical(1-id); % Changes 0 -> 1 and 1 -> 0 flipping the ID so now we can fit orange + if sum(id) >0 + [fits, crlbs, llv, framenumber] = slim_locs(iloc(:,:,id), fnum(id), cents(id,:), cal.orange.ang); + % As everywhere in the equations used sigma is squared, we can without + % loss of generality make these fits positive definite + fits(:,4) = abs(fits(:,4)); + fits(:,5) = abs(fits(:,5)); + + % Put data into cdata structure + cdata.orange.fits = fits; + cdata.orange.crlbs = crlbs; + cdata.orange.llv = llv; + cdata.orange.framenumber = framenumber; + + % Z calculations + zf = getdz(cdata.orange.fits(:,4),cdata.orange.fits(:,5),cal.orange.z_cal,2)/q; % Z assignment, this variable will be updated and stored elsewhere + ncoords = astig_tilt([cdata.orange.fits(:,1:2),zf],cal.orange); % corrections due to astigmatism + vec = xy_feature(ncoords(:,1),ncoords(:,2)); + x = o2rx.'*vec.'; + y = o2ry.'*vec.'; + % Assign fixed coordinates + cdata.orange.xf = x.'; + cdata.orange.yf = y.'; + cdata.orange.zf = ncoords(:,3); + end + save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cdata', 'pixw','q','cal'); + end + +catch lsterr + disp(lsterr) + waitforbuttonpress +end +% save('results_of_bump.mat','fnum','q','iloc','cal','cents'); + +% save('for_trial.mat','iloc' + +% Save the Analysis +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'zf_all','sigx_all' ,'sigy_all','sigx_crlb','sigy_crlb','y','iloc','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','pixw','q','pix2pho'); +% if strcmp(sv_im,'Y') || strcmp(sv_im,'y') +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'cents','zf_all','zf_crlb','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','iters','pixw','q','pix2pho','ilocs'); +% else +% figure +% imagesc(mean(iprod,3)); +% hold on +% plot(fits(:,1),fits(:,2),'rx') +% hold off +% colormap('gray'); + + +% end +% catch lsterr +% save([an_dir,'\', fname(1:end-4),'_dast.mat'], 'zf_all','sigx_all' ,'sigy_all','sigx_crlb','sigy_crlb','y','iloc','xf_all' , 'xf_crlb' , 'yf_all' , 'yf_crlb' , 'N' , 'N_crlb' ,'off_all' , 'off_crlb', 'framenum_all', 'llv','pixw','q','pix2pho'); +end + diff --git a/Minutia/image_neural_3.mexw64 b/Minutia/image_neural_3.mexw64 new file mode 100644 index 0000000..3e7eb9e Binary files /dev/null and b/Minutia/image_neural_3.mexw64 differ diff --git a/Minutia/image_neural_4.cu b/Minutia/image_neural_4.cu new file mode 100644 index 0000000..71fbf5f --- /dev/null +++ b/Minutia/image_neural_4.cu @@ -0,0 +1,334 @@ +/* +* psf_neural_discover.cu is a program to take matlab images and process them through a neural network to determine activation values for pixels in a 7x7 region centered on each pixel +* We are assuming a single output for the neural network, and a single hidden layer of 100 nodes this will perform the necessary calculation on each pixel in the image and return +* an image of 0 or 1 to be used later in the computation +* V 1.0 +* we expect a format of [im_activate] = image_neural [i1, theta1, threta2, numoframes]; +* AJN 11/2/15 +*/ + +#include "mex.h" +#include "cuda_runtime.h" +#include "device_launch_parameters.h" +#include +#include +#define PI 3.14159265358979323846 +#define O_TILE_WIDTH 20 // variable to determine how many output tiles will be considered in a block +# define BLOCK_WIDTH (O_TILE_WIDTH + (9-1)) // block width needs to be output tiles + mask_width - 1 to ensure enough pixels are covered for calculation + + + + +/* +* Device code +* +* +*/ + + + +void __global__ activate(float *d_iall, // the gaussian is a separable filter and be treated as such + float *d_theta1, // makes these elements eligible for constant caching + float *d_theta2, + float *d_ifin, + int irow, + int icol, + int numi) +{ + // Declare variables + //float d_i3[(50)] = { 0 }; + __shared__ float d_i2[(BLOCK_WIDTH)][(BLOCK_WIDTH)]; // preallocate space for shared image + /*__shared__ float d_th1[(50)][(100)]; + __shared__ float d_th2[(100)]; +*/ + // Coordinate building + int tx = threadIdx.x; // local x coord + int ty = threadIdx.y; // local y coord + int tz = threadIdx.z; + // location of output pixel being analyzed + int row_output = blockIdx.x*O_TILE_WIDTH + tx; // gives y coordinate as a function of tile width **these lose meaning for (ty || tx) >= O_TILE_WIDTH and the same is true for ** + int col_output = blockIdx.y*O_TILE_WIDTH + ty; // gives x coordinate as a function of tile width + int imnum = blockIdx.z; + if (imnum < numi){ + // initialize location of apron this forces the first pixel to take care of both the first output pixel, and loading the first input pixel + // BLOCK_WIDTH is larger than O_TILE_WIDTH so there are more threads being used than output pixels being calculated + int row_input = row_output - 4; // EACH thread should load 1 input tile to the shared image as there are [BLOCK_WIDTH]x[BLOCK_WIDTH] threads in a block + int col_input = col_output - 4; // and BLOCK_WIDTH = O_TILE_WIDTH + MASK_WIDTH-1 +/* + // Buffer data into block + for (int tcol = 0; tcol < 100; tcol++){ // buffer theta1 matrix into shared block space + for (int trow = 0; trow < 50; trow++){ + d_th1[trow][tcol] = d_theta1[trow + 50 * tcol]; + } + } + + // Buffer theta2 into shared block space + for (int t2row = 0; t2row < 101; t2row++){ + d_th2[t2row] = d_theta2[t2row]; + } */ + // buffer shared image into d_i2 + // row/col_input represents the row/col of the input pixel being considered by + // thread [blockIdx.y*BLOCK_WIDTH+ty][blockIdx.x*BLOCK_WIDTH+tx] + if ((row_input >= 0) && (row_input < irow) && (col_input >= 0) && (col_input < icol)){ // if statement checks the row/col indices to ensure they fall onto the input image + d_i2[ty][tx] = d_iall[row_input + col_input*irow + imnum*irow*icol]; // if true, the value of the image is written to the shared array at location d_i2[ty][tx] and stored locally + } // on the block + else{ + d_i2[ty][tx] = 0; // If row/col do not satisfy boundary condtions then assign a 0 to the value to build and apron of + } // of pixels that will not contribute to the calculation + + __syncthreads(); // each thread uploads to a shared array later accessed by all threads, it is imperative to synch threads here + //d_i3[0] = 1.0; + // convolution calculation + float z1[(30)] = { 0 }; + float a1[(31)] = { 0 }; + float a = 0.0; + float z2 = 0.0; + a1[0] = 1; + if (ty < O_TILE_WIDTH && tx < O_TILE_WIDTH) { // check that the local thread should be apart of the calcualtion + /*for (int rowcount = 0; rowcount < 7; rowcount++){ + for (int colcount = 0; colcount < 7; colcount++){ + d_i3[rowcount + 7 * colcount + 1] = d_i2[rowcount + ty][colcount + tx]; // linearize region to prep for matrix math ensure coloumn major setup for proper neural net behavior + } // end coloumn image for loop + }*/ + // At this point d_i3 should be linearized as matlab would do so + // perform matrix calculation here + for (int th1count = 0; th1count < 30; th1count++){ + z1[th1count] = d_theta1[82*th1count]; + for (int rowcount = 0; rowcount < 9; rowcount++){ + for (int colcount = 0; colcount < 9; colcount++){ + z1[th1count] += d_i2[rowcount+ty][colcount + tx] * d_theta1[rowcount + 9*colcount + 1 + 82*th1count]; + } + } + a1[th1count + 1] = powf(1.0 + exp(-z1[th1count]), -1.0); + } // this completes the first half of the calculation + + for (int th2count = 0; th2count<31; th2count++){ + z2 += a1[th2count] * d_theta2[th2count]; + } + a = powf(1 + exp(-z2), -1.0); // activation value + //a = z2; + + + if (row_output < irow && col_output < icol && imnum < numi){ + + d_ifin[row_output + col_output*irow + imnum*irow*icol] = tx; // assign to output variable THIS SECTION WILL CORRECTLY WRITE TO d_ifin + + } + /*else{ + __syncthreads(); + d_ifin[row_output + col_output*irow + imnum*irow*icol] = d_ifin[row_output + col_output*irow]; // assign to output variable THIS SECTION WILL CORRECTLY WRITE TO d_ifin + __syncthreads(); + } // end if else statement tyo decide what to write to final image + */ + } // end if statement to decide whether to calculate output + } +} // end gpu void activate + + + + +/* + +THIS IS THE SECTION FOR IDENTIFICATION + +*/ + + + + + +/* +* Host code +* +* +*/ + + +void mexFunction(int nlhs, mxArray *plhs[], + int nrhs, mxArray const *prhs[]) +{ + /* Declare all variables.*/ + float *iall; // the pointer to the array of all images to be analyzed + float *theta1; // pointer to theta1 matrix + float *theta2; // pointer to theta2 matrix + float *d_iall; // Pointer to image array on gpu + float *d_theta1; // Pointer to d_theta1 on gpu + float *d_theta2; // Pointer to d_theta2 on gpu + float *d_ifin; // pointer to d_ifin on gpu + int irow; // number of pixels in a row which should also be the number in a coloumn + int icol; // n + int numi; // number of images imported + const size_t *idims, *th1dims, *th2dims; + + + /* Throw an error if the input does not match expectations. */ +/* if (nrhs != 4) { + printf("Must have 4 inputs ( i1, theta1, theta2, numthreads) line: %d\n", __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + if (!mxIsfloat(prhs[0]) || mxIsComplex(prhs[0])){ + printf("i1 must be a n x m x numel(i1(1,1,:)) float array\n"); + mexErrMsgTxt("See Error above!\n"); + + } + if (!mxIsfloat(prhs[1]) || mxIsComplex(prhs[1])){ + printf("Theta1 must be a n +1 x l float array\n"); + mexErrMsgTxt("See Error above!\n"); + } + if (!mxIsfloat(prhs[2]) || mxIsComplex(prhs[2])){ + printf("Theta2 must be a l x 1 float array\n"); + mexErrMsgTxt("See Error above!\n"); + } + if (!mxIsfloat(prhs[3]) || mxIsComplex(prhs[3])){ + printf("number of threads per block must be an integer between 1 and 1024\n"); + mexErrMsgTxt("See Error above!\n"); + } +*/ + + // get pointer to input arguments + iall = (float *)mxGetPr(prhs[0]); // matlab linearizes in a coloumn major format which affects indexing (Writing MAtlab C/MEX Code - Research Gate) + idims = mxGetDimensions(prhs[0]); // get dimensions of image array + icol = (int)idims[1]; + irow = (int)idims[0]; + numi = mxGetScalar(prhs[3]); // get number of images perblock from matlab + + + // get theta1 dims + theta1 = (float *)mxGetPr(prhs[1]); + th1dims = mxGetDimensions(prhs[1]); + int th1row = (int)th1dims[0]; // number of rows in theta1 + int th1col = (int)th1dims[1]; // number of coloumns in theta1 + + // get theta2 dims + theta2 = (float *)mxGetPr(prhs[2]); + th2dims = mxGetDimensions(prhs[2]); + int th2row = (int)th2dims[0]; // number of rows in theta2 + int th2col = 1; // number of coloumns in theta2 + +/* + // EVERYONE LOVES SOME GOOD VARIABLE CHECKING!!!!!!!!!!!!!! + if (th1row != 50){ + printf("Theta1 must have 50 rows for what you want to do\n"); + mexErrMsgTxt("See Above Error!\n"); + } + + if (th1col != 30 || th2row != 31){ + printf("Theta2 must have oone more row than Theta1 has coloumns and that number should be 100\n"); + mexErrMsgTxt("See Above Error!\n"); + } + + if (th2col != 1){ + printf("Theta2 must have 1 coloumn\n"); + mexErrMsgTxt("See Above Error!\n"); + } + + // Did the User declare an output? + if (nlhs != 1){ + printf("Declare an output variable [im_activate] = image_neural(i1, Theta1.', Theta2.')\n"); // oh user...... TEACH THEM A LESSON!!!! + mexErrMsgTxt("See Error above!\n"); + }*/ + cudaDeviceReset(); + // allocate memory on the gpu device + cudaError_t err1 = cudaMalloc((void**)&d_iall, irow*icol*(numi)*sizeof(float)); // allocate image memory + if (err1 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err1), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + cudaError_t err2 = cudaMalloc((void**)&d_theta1, th1row*th1col*sizeof(float)); // allocate theta1 memory + if (err2 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err2), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + cudaError_t err3 = cudaMalloc((void**)&d_theta2, th2row*sizeof(float)); // allocate theta2 memory + if (err3 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err3), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + cudaError_t err4 = cudaMalloc((void**)&d_ifin, irow*icol*(numi)*sizeof(float)); // allocate completed activation image memory this will be a float for convience + if (err4 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err4), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + // copy data from host to device + cudaError_t err6 = cudaMemcpy(d_iall, iall, irow*icol*(numi)*sizeof(float), cudaMemcpyHostToDevice); // copy image data to gpu + if (err6 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err6), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + cudaError_t err7 = cudaMemcpy(d_theta1, theta1, th1row*th1col*sizeof(float), cudaMemcpyHostToDevice); // copy theta1 data to gpu + if (err7 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err7), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + cudaError_t err8 = cudaMemcpy(d_theta2, theta2, th2row*sizeof(float), cudaMemcpyHostToDevice); // copy theta2 data to gpu + if (err8 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err8), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } +/* + cudaError_t err10 = cudaMemcpy(d_ifin, iall, irow*icol*numi*sizeof(float), cudaMemcpyHostToDevice); // copy image data to gpu + if (err10 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err10), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } +*/ + /* Run GPU kernel*/ + dim3 dimBlock(BLOCK_WIDTH, BLOCK_WIDTH); // run 2-D gpu kernel to help with indexing + dim3 dimGrid((irow - 1) / O_TILE_WIDTH + 1, (icol - 1) / O_TILE_WIDTH + 1, numi ); + + + + + //printf("numi = %d, irow = %d, icol = %d, th1row = %d, th1col = %d, th2row = %d, th2col = %d\n", numi, irow, icol, th1row, th1col, th2row, th2col); + activate << < dimGrid, dimBlock >> >(d_iall, d_theta1, d_theta2, d_ifin, irow, icol, (numi)); + + //ident << < dimGrid, dimBlock >> >(d_ifin, d_iout, irow, icol, numi); + + /* copy data back to mxarray pointers for output + * + * + * Duplicate the input array of equal size to the output array + * Send the pointer to a variable + * copy data to place pointer points to, which is output + */ + +/* + cudaError_t errk1 = cudaPeekAtLastError(); + if (errk1 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(errk1), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + */ + cudaError_t errk2 = cudaThreadSynchronize(); + if (errk2 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(errk2), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + + + plhs[0] = mxDuplicateArray(prhs[0]); + float *iout = (float *)mxGetPr(plhs[0]); + //printf("%d\n",numi); + cudaError_t err9 = cudaMemcpy(iout, d_ifin, irow*icol*(numi)*sizeof(float), cudaMemcpyDeviceToHost); // copy xf_all data + if (err9 != cudaSuccess){ + printf("%s in %s at line %d\n", cudaGetErrorString(err9), __FILE__, __LINE__); + mexErrMsgTxt("See Error above!\n"); + } + // printf("irow %f, icol %f, numi %f, line %d\n", ifin[0], ifin[1], ifin[2], __LINE__); + + cudaDeviceReset(); +/* + cudaFree(d_iall); + cudaFree(d_theta1); + cudaFree(d_theta2); + cudaFree(d_ifin); +*/ + return; +} \ No newline at end of file diff --git a/Minutia/image_neural_4.mexw64 b/Minutia/image_neural_4.mexw64 new file mode 100644 index 0000000..16799dc Binary files /dev/null and b/Minutia/image_neural_4.mexw64 differ diff --git a/Minutia/mex_CUDA_win64.xml b/Minutia/mex_CUDA_win64.xml new file mode 100644 index 0000000..ca654c2 --- /dev/null +++ b/Minutia/mex_CUDA_win64.xml @@ -0,0 +1,166 @@ + + + + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Minutia/neural_net_live_Script.mlx b/Minutia/neural_net_live_Script.mlx new file mode 100644 index 0000000..7f24d02 Binary files /dev/null and b/Minutia/neural_net_live_Script.mlx differ