Skip to content

Commit

Permalink
Drift-correction improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
ezrabru committed Aug 25, 2023
1 parent 38465e6 commit 5baddba
Show file tree
Hide file tree
Showing 6 changed files with 144 additions and 22 deletions.
Binary file modified POLCAM-SR/POLCAM-SR.mlappinstall
Binary file not shown.
6 changes: 3 additions & 3 deletions POLCAM-SR/POLCAM-SR.prj
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
<deployment-project plugin="plugin.apptool" plugin-version="1.0">
<configuration build-checksum="593887742" file="/Users/ezrabruggeman/Documents/GitHub/POLCAM-SR/POLCAM-SR/POLCAM-SR.prj" location="/Users/ezrabruggeman/Documents/GitHub/POLCAM-SR/POLCAM-SR" name="POLCAM-SR" target="target.mlapps" target-name="Package App">
<configuration build-checksum="1671076407" file="/Users/ezrabruggeman/Documents/GitHub/POLCAM-SR/POLCAM-SR/POLCAM-SR.prj" location="/Users/ezrabruggeman/Documents/GitHub/POLCAM-SR/POLCAM-SR" name="POLCAM-SR" target="target.mlapps" target-name="Package App">
<param.appname>POLCAM-SR</param.appname>
<param.authnamewatermark>Ezra Bruggeman</param.authnamewatermark>
<param.email>eb758@cam.ac.uk</param.email>
Expand All @@ -12,8 +12,8 @@
</param.icons>
<param.summary>Application to process single-molecule and diffraction-limited polarisation camera image data</param.summary>
<param.description>Process single-molecule and diffraction-limited polarisation camera data</param.description>
<param.screenshot>/private/var/folders/6j/9pmd1ms95nn1t7d0_vmbff240000gn/T/tp5fcb071b_fce0_4f70_936f_97e7dfa348c7.png</param.screenshot>
<param.version>1.1.1</param.version>
<param.screenshot>/private/var/folders/6j/9pmd1ms95nn1t7d0_vmbff240000gn/T/tpbac4b73a_95f9_4dae_8f03_79497799f656.png</param.screenshot>
<param.version>1.1.2</param.version>
<param.products.name />
<param.products.id />
<param.products.version />
Expand Down
Binary file modified POLCAM-SR/POLCAM_SR.mlapp
Binary file not shown.
Binary file modified POLCAM-SR/dialog apps/DialogAppDriftCorrection.mlapp
Binary file not shown.
32 changes: 13 additions & 19 deletions POLCAM-SR/lib/+Drift/correctDriftRCC.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [x,y,dx,dy] = correctDriftRCC(frame,x,y,filter,segmentation,pixelsize_hist)
function [x,y,dx_spline,dy_spline,dx,dy] = correctDriftRCC(frame,x,y,filter,segmentation,pixelsize_hist)

% select only filtered localisations for drift correction
frame_filter = frame(filter);
Expand All @@ -9,6 +9,8 @@
numFrames = max(frame_filter);
numSegments = floor(numFrames/segmentation);

%%

% Get edges for 2d histogram rendering
xEdges = min(x_filter)-pixelsize_hist:pixelsize_hist:max(x_filter)+pixelsize_hist;
yEdges = min(y_filter)-pixelsize_hist:pixelsize_hist:max(y_filter)+pixelsize_hist;
Expand Down Expand Up @@ -65,10 +67,6 @@
x_peak_subpix = x_centroid - (w+1);
y_peak_subpix = y_centroid - (w+1);

% % pixel-resolution
% corr_offset = [(xpeak-size(segment2,2))
% (ypeak-size(segment2,1))];

% sub-pixel resolution
corr_offset = [((xpeak + x_peak_subpix) - size(segment2,2))
((ypeak + y_peak_subpix) - size(segment2,1))];
Expand All @@ -82,23 +80,19 @@
dx = -cumsum(dx);
dy = -cumsum(dy);

% interpolate drift
t_segment = (0:numSegments-1)'*segmentation;
t_interp = 0:(numSegments*segmentation-1);
dx_spline = spline(t_segment,dx,t_interp);
dy_spline = spline(t_segment,dy,t_interp);

% Undo calculated drift
x_driftCorrected = zeros(size(x));
y_driftCorrected = zeros(size(y));
for i=1:numSegments
if i < numSegments
lb = (i-1)*segmentation; % first frame substack
ub = i*segmentation; % last frame substack
keep = logical((frame > lb).*(frame < ub));
x_driftCorrected(keep) = x(keep) - dx(i);
y_driftCorrected(keep) = y(keep) - dy(i);
else
lb = (i-1)*segmentation; % first frame substack
keep = logical((frame > lb));
x_driftCorrected(keep) = x(keep) - dx(i);
y_driftCorrected(keep) = y(keep) - dy(i);
end

for i=1:numSegments*segmentation
idx = logical(frame == i);
x_driftCorrected(idx) = x(idx) - dx_spline(i);
y_driftCorrected(idx) = y(idx) - dy_spline(i);
end

x = x_driftCorrected;
Expand Down
128 changes: 128 additions & 0 deletions POLCAM-SR/lib/+Drift/testDriftCorection_RCC_interpolation.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
%% Test drift correction on simulated dataset
%
% Script to simulate an SMLM localization file with drift to test
% drift-correction approaches on.
%
% The 'sample' is a circle and has two fiducials in it. To be somewhat
% realistic, the fiducial is not detected in all frames.
%
% Author: Ezra Bruggeman, Lee Lab, University of Cambridge
% eb758@cam.ac.uk
%
% Last modified: 10 Jan 2021


clear all
close all
clc

segmentation = 1000;
linear_drift_x = 5/segmentation; % nm/frame
linear_drift_y = -10/segmentation; % nm/frame

numFrames = 10000; % total number of frames simulated
sigma_sample = 20; % nm, standard deviation gaussian noise on localisation error sample
sigma_fiducial = 5; % nm, standard deviation gaussian noise on localisation error fiducial
pixelsize = 57.5; % nm, pixel size (in sample space)

rate_sample = 10; % localisations per frame from sample (should be > 1)
rate_fiducial = 1; % localisations per frame from fiducial (should be <= 1)

radius_circle = 400; % nm, radius of sample, which is a circle

% can be multiple
x_fiducial = [1000,1100];
y_fiducial = [1000,600];



%% Simulate dataset

[frame,x,y] = Drift.simulateDriftDataset(linear_drift_x,linear_drift_y,...
numFrames,radius_circle,x_fiducial,y_fiducial,sigma_sample,sigma_fiducial,...
rate_sample,rate_fiducial);

figure;
scatter(x,y,20,frame,'.')
colormap jet; axis equal; grid on
xlabel('x (nm)')
ylabel('y (nm)')
title('Before drift correction')

filter = true(size(x));


%% RCC

pixelsize_hist = pixelsize/10;

[x_corr,y_corr,dx_spline,dy_spline,dx,dy] = Drift.correctDriftRCC(frame,x,y,filter,segmentation,pixelsize_hist);

%%

figure('Name','RCC drift correction');
set(0,'DefaultAxesTitleFontWeight','normal');
subplot(2,2,1)
t = 1:numFrames;
plot(t,linear_drift_x*(t-1),'b'); hold on
plot(t,dx_spline,'--b');
plot(t,linear_drift_y*(t-1),'r');
plot(t,dy_spline,'--r');
legend('dx_{true}','dx_{est}','dy_{true}','dy_{est}','Location','northwest');
xlabel('Frame'); ylabel('Drift (nm)')
xlim([-inf inf]);

subplot(2,2,2)
plot(dx_spline,dy_spline,'Color',0.5*[1,1,1]); hold on
scatter(dx_spline,dy_spline,10,t); axis equal
xlabel('x (nm)'); ylabel('y (nm)'); grid on
title('Drift trajectory')
xlim([-inf inf]); ylim([-inf inf])

subplot(2,2,3); scatter(x,y,20,frame,'.')
colormap parula; grid on; axis equal
xlabel('x (nm)')
ylabel('y (nm)')
title('Before drift correction')

subplot(2,2,4); scatter(x_corr,y_corr,20,frame,'.')
colormap parula; grid on; axis equal
xlabel('x (nm)')
ylabel('y (nm)')
title('After drift correction')


%% Fiducial-based drift correction

% [x_corr,y_corr,dx,dy] = Drift.correctDriftFiducial(frame,x,y,pixelsize);
%
% figure('Name','Fiducial-based drift correction');
% set(0,'DefaultAxesTitleFontWeight','normal');
%
% subplot(2,2,1)
% t = linspace(0,numFrames,floor(numFrames/segmentation));
% plot(dx,'b'); hold on
% plot(dy,'r');
% legend('dx_{est}','dy_{est}','Location','northwest');
% xlabel('Frame'); ylabel('Drift (nm)')
% xlim([-inf inf]);
%
% subplot(2,2,2)
% plot(dx,dy,'Color',0.5*[1,1,1]); hold on
% scatter(dx,dy,20,1:length(dx))
% xlabel('x (nm)'); ylabel('y (nm)'); grid on
% title('Drift trajectory')
% xlim([-inf inf]); ylim([-inf inf])
%
% subplot(2,2,3); scatter(x,y,20,frame,'.')
% colormap parula; grid on
% xlabel('x (nm)')
% ylabel('y (nm)')
% title('Before drift correction')
% subplot(2,2,4); scatter(x_corr,y_corr,20,frame,'.')
% colormap parula; grid on
% xlabel('x (nm)')
% ylabel('y (nm)')
% title('After drift correction')


0 comments on commit 5baddba

Please sign in to comment.