diff --git a/data_examples/nonlinear_drift_correction_synthetic_dataset_for_testing.mat b/data_examples/nonlinear_drift_correction_synthetic_dataset_for_testing.mat new file mode 100644 index 0000000..6b379c6 Binary files /dev/null and b/data_examples/nonlinear_drift_correction_synthetic_dataset_for_testing.mat differ diff --git a/matlab/SPmakeImage.m b/matlab/SPmakeImage.m new file mode 100644 index 0000000..331358b --- /dev/null +++ b/matlab/SPmakeImage.m @@ -0,0 +1,84 @@ +function [sMerge] = SPmakeImage(sMerge,indImage,indLines) + +% This function generates a resampled scanning probe image with dimensions +% of imageSize, from a an array of N scan lines given in scaneLines, +% (lines specified as image rows), from an array of Nx2 origins in scanOr. +% scanDir is a 2 element vector specifying the direction of the scan. +% All arrays are stored inside struct sMerge. ind specified update index. +% indLines is a vector of binary values specifying which lines to include. + +if nargin == 2 + % indLines = 1:size(sMerge.scanLines,2); + indLines = true(1,size(sMerge.scanLines,2)); +end + +% Expand coordinates +t = repmat(1:size(sMerge.scanLines,2),[sum(indLines) 1]); +x0 = repmat(sMerge.scanOr(indLines,1,indImage),[1 size(sMerge.scanLines,2)]); +y0 = repmat(sMerge.scanOr(indLines,2,indImage),[1 size(sMerge.scanLines,2)]); +xInd = x0(:) + t(:)*sMerge.scanDir(indImage,1); +yInd = y0(:) + t(:)*sMerge.scanDir(indImage,2); + +% Prevent pixels from leaving image boundaries +xInd = max(min(xInd,sMerge.imageSize(1)-1),1); +yInd = max(min(yInd,sMerge.imageSize(2)-1),1); + +% Convert to bilinear interpolants and weights +xIndF = floor(xInd); +yIndF = floor(yInd); +xAll = [xIndF xIndF+1 xIndF xIndF+1]; +yAll = [yIndF yIndF yIndF+1 yIndF+1]; +dx = xInd-xIndF; +dy = yInd-yIndF; +w = [(1-dx).*(1-dy) dx.*(1-dy) (1-dx).*dy dx.*dy]; +indAll = sub2ind(sMerge.imageSize,xAll,yAll); + +% Generate image +sL = sMerge.scanLines(indLines,:,indImage); +sig = reshape(accumarray(indAll(:),[ ... + w(:,1).*sL(:); + w(:,2).*sL(:); + w(:,3).*sL(:); + w(:,4).*sL(:)],... + [prod(sMerge.imageSize) 1]),sMerge.imageSize); +count = reshape(accumarray(indAll(:),[ ... + w(:,1);w(:,2);w(:,3);w(:,4)],... + [prod(sMerge.imageSize) 1]),sMerge.imageSize); + +% Apply KDE +r = max(ceil(sMerge.KDEsigma*3),5); +sm = fspecial('gaussian',2*r+1,sMerge.KDEsigma); +sm = sm / sum(sm(:)); +sig = conv2(sig,sm,'same'); +count = conv2(count,sm,'same'); +sub = count > 0; +sig(sub) = sig(sub) ./ count(sub); +sMerge.imageTransform(:,:,indImage) = sig; + +% Estimate sampling density +bound = count == 0; +bound([1 end],:) = true; +bound(:,[1 end]) = true; +sMerge.imageDensity(:,:,indImage) = ... + sin(min(bwdist(bound)/sMerge.edgeWidth,1)*pi/2).^2; + + +% % Plot testing +% % dens = sMerge.imageDensity(:,:,indImage); +% figure(1) +% clf +% imagesc(sig) +% % imagesc(sig.*dens + (1-dens)*mean(sig(dens>.5))) +% % imagesc(sMerge.imageTransform(:,:,indImage)) +% axis equal off +% set(gca,'position',[0 0 1 1]) +% % imagesc(t) +% % hold on +% % ind = 17:3213:numel(xInd); +% % scatter(yInd(ind),xInd(ind),'r.') +% % % hold off +% % % axis equal off +% % % colormap(gray(256)) + + +end \ No newline at end of file diff --git a/matlab/SPmerge01linear.m b/matlab/SPmerge01linear.m new file mode 100644 index 0000000..1bff947 --- /dev/null +++ b/matlab/SPmerge01linear.m @@ -0,0 +1,279 @@ +function [sMerge] = SPmerge01linear(scanAngles,varargin) + +% Colin Ophus, National Center for Electron Microscopy, Molecular Foundry, +% Lawrence Berkeley National Laboratory, Berkeley, CA, USA. (Feb 2016). + +% New version of SPmerge01.m - This script now searches over linear drift +% vectors, aligned to first two images. This search is performed twice. + +% Merge multiple scanning probe images. Assume scan origin is upper left +% corner of the image, and that the scan direction for zero degrees is +% horizontal (along MATLAB columns). All input images must have fast scan +% direction along the array rows (horizontal direction). Original data is +% stored in 3D arrray sMerge.scanLines + +% Inputs: +% scanAngles - Vector containg scan angles in degrees. +% images - 2D image arrays, order in scanAngles, all same size. +% sigmaLP = 32; % - sigma value in pixels for initial alignment LP filter. +flagReportProgress = 1; % Set to true to see updates on console.v = +paddingScale = (1+1/4);%1.5; % - padding amount for scaling of the output. +sMerge.KDEsigma = 1/2; % - Smoothing between pixels for KDE. +sMerge.edgeWidth = 1/128; % - size of edge blending relative to input images. +sMerge.linearSearch = linspace(-0.02,0.02,1+2*2); % Initial linear search vector, relative to image size. +% sMerge.linearSearch = linspace(-0.04,0.04,1+2*4); % Initial linear search vector, relative to image size. + +% flagSkipInitialAlignment = 0; % Set to true to skip initial phase correlation. +% flagCrossCorrelation = 0+1; % Set to true to use cross correlation. +% sigmaAlignLP = 64; % low pass filtering for initial alignment (pixels). + +% Initialize struct containing all data +sMerge.imageSize = round([size(varargin{1},1) size(varargin{1},2)]*paddingScale/4)*4; +sMerge.scanAngles = scanAngles; +sMerge.numImages = length(scanAngles); +sMerge.scanLines = zeros(size(varargin{1},1),size(varargin{1},2),sMerge.numImages); +sMerge.scanOr = zeros(size(varargin{1},1),2,sMerge.numImages); +sMerge.scanDir = zeros(sMerge.numImages,2); +sMerge.imageTransform = zeros(sMerge.imageSize(1),sMerge.imageSize(2),sMerge.numImages); +sMerge.imageDensity = zeros(sMerge.imageSize(1),sMerge.imageSize(2),sMerge.numImages); + +% Generate all initial variables +for a0 = 1:sMerge.numImages + % Raw images + if nargin == (sMerge.numImages+1) + % Sequential images + sMerge.scanLines(:,:,a0) = varargin{a0}; + elseif nargin == 2 + % 3D image stack + sMerge.scanLines(:,:,a0) = varargin{1}(:,:,a0); + else + + end + + % Generate pixel origins + if nargin == (sMerge.numImages+1) + xy = [(1:size(varargin{a0},1))' ones(size(varargin{a0},1),1)]; + xy(:,1) = xy(:,1) - size(varargin{a0},1)/2; + xy(:,2) = xy(:,2) - size(varargin{a0},2)/2; + elseif nargin == 2 + xy = [(1:size(varargin{1}(:,:,a0),1))' ones(size(varargin{1}(:,:,a0),1),1)]; + xy(:,1) = xy(:,1) - size(varargin{1}(:,:,a0),1)/2; + xy(:,2) = xy(:,2) - size(varargin{1}(:,:,a0),2)/2; + end + xy = [xy(:,1)*cos(scanAngles(a0)*pi/180) ... + - xy(:,2)*sin(scanAngles(a0)*pi/180) ... + xy(:,2)*cos(scanAngles(a0)*pi/180) ... + + xy(:,1)*sin(scanAngles(a0)*pi/180)]; + xy(:,1) = xy(:,1) + sMerge.imageSize(1)/2; + xy(:,2) = xy(:,2) + sMerge.imageSize(2)/2; + xy(:,1) = xy(:,1) - mod(xy(1,1),1); + xy(:,2) = xy(:,2) - mod(xy(1,2),1); + sMerge.scanOr(:,:,a0) = xy; + + % Scan direction + sMerge.scanDir(a0,:) = [cos(scanAngles(a0)*pi/180 + pi/2) ... + sin(scanAngles(a0)*pi/180 + pi/2)]; +end + +if flagReportProgress == true + reverseStr = ''; % initialize console message piece +end + + +% First linear alignment, search over possible linear drift vectors. +sMerge.linearSearch = sMerge.linearSearch * size(sMerge.scanLines,1); +[yDrift,xDrift] = meshgrid(sMerge.linearSearch); +sMerge.linearSearchScore1 = zeros(length(sMerge.linearSearch)); +inds = linspace(-0.5,0.5,size(sMerge.scanLines,1))'; +N = size(sMerge.scanLines); +w2 = circshift(padarray(hanningLocal(N(1))*hanningLocal(N(2))',... + sMerge.imageSize-N(1:2),0,'post'),round((sMerge.imageSize-N(1:2))/2)); +% IcorrNorm = ifft2(abs(fft2(w2)).^2,'symmetric'); +for a0 = 1:length(sMerge.linearSearch) + for a1 = 1:length(sMerge.linearSearch) + % Calculate time dependent lienar drift + xyShift = [inds*xDrift(a0,a1) inds*yDrift(a0,a1)]; + + % Apply linear drift to first two images + sMerge.scanOr(:,:,1:2) = sMerge.scanOr(:,:,1:2) ... + + repmat(xyShift,[1 1 2]); + % Generate trial images + sMerge = SPmakeImage(sMerge,1); + sMerge = SPmakeImage(sMerge,2); + % Measure alignment score with hybrid correlation + m = fft2(w2.*sMerge.imageTransform(:,:,1)) ... + .* conj(fft2(w2.*sMerge.imageTransform(:,:,2))); + Icorr = ifft2(sqrt(abs(m)).*exp(1i*angle(m)),'symmetric'); + sMerge.linearSearchScore1(a0,a1) = max(Icorr(:)); + % Remove linear drift from first two images + sMerge.scanOr(:,:,1:2) = sMerge.scanOr(:,:,1:2) ... + - repmat(xyShift,[1 1 2]); + if flagReportProgress == true + comp = (a1 / length(sMerge.linearSearch) ... + + a0 - 1) / length(sMerge.linearSearch); + msg = sprintf(['Initial Linear Drift Search = ' ... + sprintf('%.02f',100*comp) ' percent complete']); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); + end + end +end + + +% Second linear alignment, refine possible linear drift vectors. +[~,ind] = max(sMerge.linearSearchScore1(:)); +[xInd,yInd] = ind2sub(size(sMerge.linearSearchScore1),ind); +% if xInd == 1 +% xInds = 1:2; +% elseif xInd == size(sMerge.linearSearchScore1,1) +% xInds = (-1:0) + size(sMerge.linearSearchScore1,1); +% else +% if sMerge.linearSearchScore1(xInd-1,yInd) ... +% > sMerge.linearSearchScore1(xInd+1,yInd) +% xInds = xInd + (-1:0); +% else +% xInds = xInd + (0:1); +% end +% end +% if yInd == 1 +% yInds = 1:2; +% elseif yInd == size(sMerge.linearSearchScore1,2) +% yInds = (-1:0) + size(sMerge.linearSearchScore1,2); +% else +% if sMerge.linearSearchScore1(xInd,yInd-1) ... +% > sMerge.linearSearchScore1(xInd,yInd+1) +% yInds = yInd + (-1:0); +% else +% yInds = yInd + (0:1); +% end +% end +step = sMerge.linearSearch(2) - sMerge.linearSearch(1); +xRefine = sMerge.linearSearch(xInd) ... + + linspace(-0.5,0.5,length(sMerge.linearSearch))*step; +yRefine = sMerge.linearSearch(yInd) ... + + linspace(-0.5,0.5,length(sMerge.linearSearch))*step; +[yDrift,xDrift] = meshgrid(yRefine,xRefine); +sMerge.linearSearchScore2 = zeros(length(sMerge.linearSearch)); +for a0 = 1:length(sMerge.linearSearch) + for a1 = 1:length(sMerge.linearSearch) + % Calculate time dependent lienar drift + xyShift = [inds*xDrift(a0,a1) inds*yDrift(a0,a1)]; + + % Apply linear drift to first two images + sMerge.scanOr(:,:,1:2) = sMerge.scanOr(:,:,1:2) ... + + repmat(xyShift,[1 1 2]); + % Generate trial images + sMerge = SPmakeImage(sMerge,1); + sMerge = SPmakeImage(sMerge,2); + % Measure alignment score with hybrid correlation + m = fft2(w2.*sMerge.imageTransform(:,:,1)) ... + .* conj(fft2(w2.*sMerge.imageTransform(:,:,2))); + Icorr = ifft2(sqrt(abs(m)).*exp(1i*angle(m)),'symmetric'); + sMerge.linearSearchScore2(a0,a1) = max(Icorr(:)); + % Remove linear drift from first two images + sMerge.scanOr(:,:,1:2) = sMerge.scanOr(:,:,1:2) ... + - repmat(xyShift,[1 1 2]); + if flagReportProgress == true + comp = (a1 / length(sMerge.linearSearch) ... + + a0 - 1) / length(sMerge.linearSearch); + msg = sprintf(['Linear Drift Refinement = ' ... + sprintf('%.02f',100*comp) ' percent complete']); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); + end + end +end +[~,ind] = max(sMerge.linearSearchScore2(:)); +[xInd,yInd] = ind2sub(size(sMerge.linearSearchScore2),ind); +sMerge.xyLinearDrift = [xDrift(xInd) yDrift(yInd)]; + + +% Apply linear drift to all images +xyShift = [inds*sMerge.xyLinearDrift(1) inds*sMerge.xyLinearDrift(2)]; +for a0 = 1:sMerge.numImages + sMerge.scanOr(:,:,a0) = sMerge.scanOr(:,:,a0) + xyShift; + sMerge = SPmakeImage(sMerge,a0); +end + +% Estimate initial alignment +dxy = zeros(sMerge.numImages,2); +G1 = fft2(w2.*sMerge.imageTransform(:,:,1)); +for a0 = 2:sMerge.numImages + G2 = fft2(w2.*sMerge.imageTransform(:,:,a0)); + m = G1.*conj(G2); + Icorr = ifft2(sqrt(abs(m)).*exp(1i*angle(m)),'symmetric'); + [~,ind] = max(Icorr(:)); + [dx,dy] = ind2sub(size(Icorr),ind); + dx = mod(dx - 1 + size(Icorr,1)/2,size(Icorr,1)) - size(Icorr,1)/2; + dy = mod(dy - 1 + size(Icorr,2)/2,size(Icorr,2)) - size(Icorr,2)/2; + dxy(a0,:) = dxy(a0-1,:) + [dx dy]; + G1 = G2; +end +dxy(:,1) = dxy(:,1) - mean(dxy(:,1)); +dxy(:,2) = dxy(:,2) - mean(dxy(:,2)); +% Apply alignments and regenerate images +for a0 = 1:sMerge.numImages + sMerge.scanOr(:,1,a0) = sMerge.scanOr(:,1,a0) + dxy(a0,1); + sMerge.scanOr(:,2,a0) = sMerge.scanOr(:,2,a0) + dxy(a0,2); + sMerge = SPmakeImage(sMerge,a0); +end +% Set reference point +sMerge.ref = round(sMerge.imageSize/2); + + +% Plot results, image with scanline origins overlaid +imagePlot = mean(sMerge.imageTransform,3); +dens = prod(sMerge.imageDensity,3); +% Scale intensity of image +mask = dens>0.5; +imagePlot = imagePlot - mean(imagePlot(mask)); +imagePlot = imagePlot / sqrt(mean(imagePlot(mask).^2)); + +figure(1) +clf +imagesc(imagePlot) +hold on +cvals = [1 0 0; + 0 .7 0; + 0 .6 1; + 1 .7 0; + 1 0 1; + 0 0 1]; +for a0 = 1:sMerge.numImages + scatter(sMerge.scanOr(:,2,a0),sMerge.scanOr(:,1,a0),'marker','.',... + 'sizedata',25,'markeredgecolor',cvals(mod(a0-1,size(cvals,1))+1,:)) +end +scatter(sMerge.ref(2),sMerge.ref(1),... + 'marker','+','sizedata',500,... + 'markeredgecolor',[1 1 0],'linewidth',6) +scatter(sMerge.ref(2),sMerge.ref(1),... + 'marker','+','sizedata',500,... + 'markeredgecolor','k','linewidth',2) +hold off +axis equal off +colormap(gray(256)) +set(gca,'position',[0 0 1 1]) +caxis([-3 3]) % units of image RMS + +if flagReportProgress == true + fprintf([reverseStr ' ']); +end +end + +% function [Iscale] = scaleImage(I,intRange,sigmaLP) +% Iscale = I - mean(I(:)); +% Iscale = Iscale / sqrt(mean(I(:).^2)); +% Iscale = (Iscale - intRange(1))/(intRange(2)-intRange(1)); +% if sigmaLP > 0 +% sm = fspecial('gaussian',2*ceil(3*sigmaLP)+1,sigmaLP); +% Iscale = conv2(Iscale,sm,'same') ... +% ./ conv2(ones(size(Iscale)),sm,'same'); +% end +% Iscale(Iscale<0) = 0; +% Iscale(Iscale>1) = 1; +% end + +function [w] = hanningLocal(N) +% Replacement for simplest 1D hanning function to remove dependency. +w = sin(pi*((1:N)'/(N+1))).^2; +end \ No newline at end of file diff --git a/matlab/SPmerge02.m b/matlab/SPmerge02.m new file mode 100644 index 0000000..5340437 --- /dev/null +++ b/matlab/SPmerge02.m @@ -0,0 +1,574 @@ +function [sMerge] = SPmerge02(sMerge,refineMaxSteps,initialRefineSteps) +% tic + +% Colin Ophus, National Center for Electron Microscopy, Molecular Foundry, +% Lawrence Berkeley National Laboratory, Berkeley, CA, USA. (Mar 2015). + +% Refinement function for scanning probe image merge. Requires struct +% input from SPmerge01. This script will perform both the initial +% alignment if it has not been performed, and refinement steps. + +% Inputs: +% sMerge - struct containing data for STEM alignment from SPmerge01 +% refineMaxSteps - Maximum number of refinement steps (optional). If this +% value is not specified, it will be set to 32 iterations. +% Set to 0 to perform only initial refinement. +flagPlot = 1; +flagReportProgress = 1; % Set to true to see updates on console. +densityCutoff = 0.8; % density cutoff for image boundaries (norm. to 1). +distStart = mean([size(sMerge.scanLines,1) ... + size(sMerge.scanLines,2)])/16; % Radius of # of scanlines used for initial alignment. +% initialRefineSteps = 8/2; % Number of initial refinement iterations. +initialShiftMaximum = 1/4; % Max number of pixels shifted per line for the +% initial alignment step. This value should +% have a maximum of 1, but can be set lower +% to stabilize initial alignment. +refineInitialStep = 1/2; % Initial step size for main refinement, in pixels. +pixelsMovedThreshold = 0.1;% If number of pixels shifted (per image) is +% below this value, refinement will be halted. +stepSizeReduce = 1/2; % When a scanline origin does not move, +% step size will be reduced by this factor. +flagPointOrder = 1; % Use this flag to force origins to be ordered, i.e. +% disallow points from changing their order. +flagGlobalShift = 1; % If this flag is true, a global phase correlation +% performed each primary iteration (This is meant to +% fix unit cell shifts and similar artifacts). +% This option is highly recommended! +flagGlobalShiftIncrease = 0; % If this option is true, the global scoring +% function is allowed to increase after global +% phase correlation step. (false is more stable) +minGlobalShift = 4; % Global shifts only if shifts > this value (pixels) +densityDist = mean([size(sMerge.scanLines,1) ... + size(sMerge.scanLines,2)])/32; % density mask edge threshold +% To generate a moving average along the scanline origins +% (make scanline steps more linear), use the settings below: +originWindowAverage = 16; % Window sigma in px for smoothing scanline origins. +% Set this value to zero to not use window avg. +% This window is relative to linear steps. +originInitialAverage = mean([size(sMerge.scanLines,1) ... + size(sMerge.scanLines,2)])/4; % Window sigma in px for initial smoothing. +resetInitialAlignment = 0; % Set this value to true to redo initial alignment. + + +% Outputs: +% sMerge - struct containing data for STEM alignment + + +% Default number of iterations if user does not provide values +if nargin == 1 + refineMaxSteps = 32; + initialRefineSteps = 4*0; +elseif nargin == 2 + initialRefineSteps = 4*0; +end + + +% Make kernel for moving average of origins +if originInitialAverage > 0 % || originLinearFraction > 0 + if originInitialAverage > 0 + r = ceil(3*originInitialAverage); + v = (-r:r)'; + KDEorigin = exp(-v.^2/(2*originInitialAverage^2)); + else + KDEorigin = 1; + end + KDEnorm = 1./convn(ones(size(sMerge.scanOr)),KDEorigin,'same'); + % indsOr = repmat((1:size(sMerge.scanOr,1))',... + % [2 sMerge.numImages]); + basisOr = [ones(size(sMerge.scanLines,1),1) ... + (1:size(sMerge.scanLines,1))']; + scanOrLinear = zeros(size(sMerge.scanOr)); +end + + +% If initial alignment has not been performed, do it. +if flagReportProgress == true + reverseStr = ''; % initialize console message piece +end +if (~isfield(sMerge,'scanActive') || resetInitialAlignment == true) ... + || ((initialRefineSteps > 0) && (nargin == 3)) + for aInit = 1:initialRefineSteps + + sMerge.scanActive = false(size(sMerge.scanLines,1),... + sMerge.numImages); + + % Get starting scanlines for initial alignment + indStart = zeros(sMerge.numImages,1); + for a0 = 1:sMerge.numImages + % Scan line direction and origins + v = [-sMerge.scanDir(a0,2) sMerge.scanDir(a0,1)]; + or = sMerge.scanOr(:,1:2,a0); + + % Determine closest scanline origin from point-line distance + c = -sum(sMerge.ref.*v); + dist = abs(v(1)*or(:,1)+v(2)*or(:,2) + c) / norm(v); + [~,indStart(a0)] = min(dist); + sub = dist < distStart; + sMerge.scanActive(sub,a0) = true; + end + + % Rough initial alignment of scanline origins, to nearest pixel + inds = 1:size(sMerge.scanLines,1); + dxy = [0 0; + 1 0; + -1 0; + 0 1; + 0 -1]; + score = zeros(size(dxy,1),1); + for a0 = 1:sMerge.numImages + % Determine which image to align to, based on orthogonality + [~,indAlign] = min(abs(sum(repmat(sMerge.scanDir(a0,:), ... + [sMerge.numImages 1]).* sMerge.scanDir,2))); + + % Generate alignment image + if ~isfield(sMerge,'imageRef') + sMerge = SPmakeImage(sMerge,indAlign,sMerge.scanActive(:,indAlign)); + imageAlign = sMerge.imageTransform(:,:,indAlign) ... + .* (sMerge.imageDensity(:,:,indAlign)>densityCutoff); + else + imageAlign = sMerge.imageRef; + end + + % sub = sMerge.imageDensity(:,:,indAlign)>densityCutoff; + % imageAlign = sMerge.imageTransform(:,:,indAlign).*sub ... + % + (1-sub)*intensityMedian; + + + % figure(111) + % clf + % imagesc(imageAlign) + % axis equal off + % colormap(violetFire) + % drawnow + + + % align origins + xyStep = mean(sMerge.scanOr(2:end,:,a0)-sMerge.scanOr(1:(end-1),:,a0),1); + indAligned = false(size(sMerge.scanLines,1),1); + indAligned(indStart(a0)) = true; + while all(indAligned) == 0 + % Determine scanline indices to check next + v = bwmorph(indAligned,'dilate',1); + v(indAligned) = false; + indMove = inds(v); + + % currently active scanlines + indsActive = inds(indAligned); + +% % % If needed, get linear estimates +% % if weightInitialLinear > 0 +% % A = [ones(length(indsActive),1) indsActive']; +% % beta = A \ sMerge.scanOr(indAligned,1:2,a0); +% % end + + % Align selected scanlines + for a1 = 1:length(indMove) + % determine starting point from neighboring scanline + [~,minDistInd] = min(abs(indMove(a1)-indsActive)); + + % Step perpendicular to scanDir + xyOr = sMerge.scanOr(indsActive(minDistInd),1:2,a0) ... + + xyStep * (indMove(a1)-indsActive(minDistInd)); + + % Refine score by moving origin of this scanline + xInd = round(xyOr(1) + inds*sMerge.scanDir(a0,1)); + yInd = round(xyOr(2) + inds*sMerge.scanDir(a0,2)); + % % Prevent pixels from leaving image boundaries + % xInd = max(min(xInd,sMerge.imageSize(1)-1),1); + % yInd = max(min(yInd,sMerge.imageSize(2)-1),1); + for a2 = 1:size(dxy,1) + score(a2) = sum(abs(imageAlign(sub2ind(sMerge.imageSize,... + xInd + dxy(a2,1),yInd + dxy(a2,2))) ... + - sMerge.scanLines(indMove(a1),:,a0))); + end + [~,ind] = min(score); + sMerge.scanOr(indMove(a1),1:2,a0) = xyOr ... + + dxy(ind,1:2) * initialShiftMaximum; + indAligned(indMove(a1)) = true; + end + + % Report progess if flag is set to true + if flagReportProgress == true + comp = sum(indAligned) / numel(indAligned); + msg = sprintf(['Initial refinement ' ... + num2str(aInit) '/' num2str(initialRefineSteps) ... + ' of image ' ... + num2str(a0) '/' num2str(sMerge.numImages) ' is ' ... + num2str(round(100*comp)) ' percent complete']); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); + end + end + end + + % If required, compute moving average of origins using KDE. + if originInitialAverage > 0 % || originLinearFraction > 0 + % Linear fit to scanlines + for a0 = 1:sMerge.numImages + ppx = basisOr \ sMerge.scanOr(:,1,a0); + ppy = basisOr \ sMerge.scanOr(:,2,a0); + scanOrLinear(:,1,a0) = basisOr*ppx; + scanOrLinear(:,2,a0) = basisOr*ppy; + end + % Subtract linear fit + sMerge.scanOr = sMerge.scanOr - scanOrLinear; + % Moving average of scanlines using KDE + sMerge.scanOr = convn(sMerge.scanOr,KDEorigin,'same') ... + .* KDEnorm; + % Add linear fit back into to origins, and/or linear weighting + % sMerge.scanOr = sMerge.scanOr *(1-originLinearFraction) ... + % + scanOrLinear; + sMerge.scanOr = sMerge.scanOr + scanOrLinear; + end + end +end + + + +% Make kernel for moving average of origins +if originWindowAverage > 0 % || originLinearFraction > 0 + if originWindowAverage > 0 + r = ceil(3*originWindowAverage); + v = (-r:r)'; + KDEorigin = exp(-v.^2/(2*originWindowAverage^2)); + else + KDEorigin = 1; + end + KDEnorm = 1./convn(ones(size(sMerge.scanOr)),KDEorigin,'same'); + % indsOr = repmat((1:size(sMerge.scanOr,1))',... + % [2 sMerge.numImages]); + basisOr = [ones(size(sMerge.scanLines,1),1) ... + (1:size(sMerge.scanLines,1))']; + scanOrLinear = zeros(size(sMerge.scanOr)); +end + + + +% Main alignment steps +if flagReportProgress == true + msg = sprintf('Beginning primary refinement ...'); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); +end +scanOrStep = ones(size(sMerge.scanOr,1),size(sMerge.scanOr,3)) ... + *refineInitialStep; +inds = 1:size(sMerge.scanLines,2); +dxy = [0 0; + 1 0; + -1 0; + 0 1; + 0 -1]; +score = zeros(size(dxy,1),1); +alignStep = 1; +sMerge.stats = zeros(refineMaxSteps+1,2); +indsLoop = 1:sMerge.numImages; +while alignStep <= refineMaxSteps + % Reset pixels moved count + pixelsMoved = 0; + + % Compute all images from current origins + for a0 = indsLoop + sMerge = SPmakeImage(sMerge,a0); + end + + % Get mean absolute difference as a fraction of the mean scanline intensity. + Idiff = mean(abs(sMerge.imageTransform ... + - repmat(mean(sMerge.imageTransform,3),[1 1 sMerge.numImages])),3); + meanAbsDiff = mean(Idiff(min(sMerge.imageDensity,[],3)>densityCutoff)) ... + / mean(abs(sMerge.scanLines(:))); + sMerge.stats(alignStep,1:2) = [alignStep-1 meanAbsDiff]; + + % If required, check for global alignment of images + if flagGlobalShift == true + if flagReportProgress == true + msg = sprintf('Checking global alignment ...'); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); + end + + % save current origins, step size and score + scanOrCurrent = sMerge.scanOr; + scanOrStepCurrent = scanOrStep; + meanAbsDiffCurrent = meanAbsDiff; + + % Align to windowed image 1 or imageRef + intensityMedian = median(sMerge.scanLines(:)); + densityMask = sin((pi/2)*min(bwdist( ... + sMerge.imageDensity(:,:,1) minGlobalShift + % apply global origin shift, if possible + xNew = sMerge.scanOr(:,1,a0) + dx; + yNew = sMerge.scanOr(:,2,a0) + dy; + + % Verify shifts are within image boundaries + if min(xNew) >= 1 ... + && max(xNew) < sMerge.imageSize(1)-1 ... + && min(yNew) >= 1 ... + && max(yNew) < sMerge.imageSize(2)-1 + sMerge.scanOr(:,1,a0) = xNew; + sMerge.scanOr(:,2,a0) = yNew; + + % Recompute image with new origins + sMerge = SPmakeImage(sMerge,a0); + + % Reset search values for this image + scanOrStep(:,a0) = refineInitialStep; + end + end + end + + if flagGlobalShiftIncrease == false + % Verify global shift did not make mean abs. diff. increase. + Idiff = mean(abs(sMerge.imageTransform ... + - repmat(mean(sMerge.imageTransform,3),[1 1 sMerge.numImages])),3); + meanAbsDiffNew = mean(Idiff(min(sMerge.imageDensity,[],3)>densityCutoff)) ... + / mean(abs(sMerge.scanLines(:))); + + if meanAbsDiffNew < meanAbsDiffCurrent + % If global shift decreased mean absolute different, keep. + sMerge.stats(alignStep,1:2) = [alignStep-1 meanAbsDiff]; + else + % If global shift incresed mean abs. diff., return origins + % and step sizes to previous values. + sMerge.scanOr = scanOrCurrent; + scanOrStep = scanOrStepCurrent; + end + + end + end + + % Refine each image in turn, against the sum of all other images + for a0 = indsLoop + % Generate alignment image, mean of all other scanline datasets, + % unless user has specified a reference image. + if ~isfield(sMerge,'imageRef') + indsAlign = 1:sMerge.numImages; + indsAlign(a0) = []; + imageAlign = sum(sMerge.imageTransform(:,:,indsAlign) ... + .* (sMerge.imageDensity(:,:,indsAlign)>densityCutoff),3); + dens = sum(sMerge.imageDensity(:,:,indsAlign)>densityCutoff,3); + sub = dens > 0; + imageAlign(sub) = imageAlign(sub) ./ dens(sub); + imageAlign(~sub) = mean(imageAlign(sub)); % Replace zeros with mean + else + imageAlign = sMerge.imageRef; + end + + % If ordering is used as a condition, determine parametric positions + if flagPointOrder == true + % Use vector perpendicular to scan direction (negative 90 deg) + n = [sMerge.scanDir(a0,2) -sMerge.scanDir(a0,1)]; + vParam = n(1)*sMerge.scanOr(:,1,a0) + n(2)*sMerge.scanOr(:,2,a0); + end + + % Loop through each scanline and perform alignment + for a1 = 1:size(sMerge.scanLines,1) + % Refine score by moving the origin of this scanline + orTest = repmat(sMerge.scanOr(a1,1:2,a0),[size(dxy,1) 1]) ... + + dxy * scanOrStep(a1,a0); + + % If required, force ordering of points + if flagPointOrder == true + vTest = n(1)*orTest(:,1) + n(2)*orTest(:,2); + if a1 == 1 + vBound = [-Inf vParam(a1+1)]; + elseif a1 == size(sMerge.scanLines,1) + vBound = [vParam(a1-1) Inf]; + else + vBound = [vParam(a1-1) vParam(a1+1)]; + end + for a2 = 1:size(dxy,1) + if vTest(a2) < vBound(1) + orTest(a2,:) = orTest(a2,:) + n*(vBound(1)-vTest(a2)); + elseif vTest(a2) > vBound(2) + orTest(a2,:) = orTest(a2,:) + n*(vBound(2)-vTest(a2)); + end + end + end + + % Loop through origin tests + for a2 = 1:size(dxy,1) + xInd = orTest(a2,1) + inds*sMerge.scanDir(a0,1); + yInd = orTest(a2,2) + inds*sMerge.scanDir(a0,2); + % Prevent pixels from leaving image boundaries + xInd = max(min(xInd,sMerge.imageSize(1)-1),1); + yInd = max(min(yInd,sMerge.imageSize(2)-1),1); + % Bilinear coordinates + xF = floor(xInd); + yF = floor(yInd); + score(a2) = calcScore(imageAlign,xF,yF,xInd-xF,yInd-yF,... + sMerge.scanLines(a1,:,a0)); + end + % Note that if moving origin does not change score, dxy = (0,0) + % will be selected (ind = 1). + [~,ind] = min(score); + if ind == 1 + % Reduce the step size for this origin + scanOrStep(a1,a0) = scanOrStep(a1,a0) * stepSizeReduce; + else + pixelsMoved = pixelsMoved ... + + norm(orTest(ind,:)-sMerge.scanOr(a1,1:2,a0)); + sMerge.scanOr(a1,1:2,a0) = orTest(ind,:); + end + + % Report progress if requested + if flagReportProgress == true && mod(a1,16) == 0 + comp = (a1 / size(sMerge.scanLines,1) ... + + a0 - 1) / sMerge.numImages; + msg = sprintf([ ... + 'Mean Abs. Diff. = ' sprintf('%.04f',100*meanAbsDiff) ... + ' percent, ' ... + 'iter. ' num2str(alignStep) ... + '/' num2str(refineMaxSteps) ', ' ... + num2str(round(100*comp)) ' percent complete, ' ... + num2str(round(pixelsMoved)) ' px moved']); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); + end + end + end + + % If required, compute moving average of origins using KDE. + if originWindowAverage > 0 % || originLinearFraction > 0 + % Linear fit to scanlines + for a0 = 1:sMerge.numImages + ppx = basisOr \ sMerge.scanOr(:,1,a0); + ppy = basisOr \ sMerge.scanOr(:,2,a0); + scanOrLinear(:,1,a0) = basisOr*ppx; + scanOrLinear(:,2,a0) = basisOr*ppy; + end + % Subtract linear fit + sMerge.scanOr = sMerge.scanOr - scanOrLinear; + % Moving average of scanlines using KDE + sMerge.scanOr = convn(sMerge.scanOr,KDEorigin,'same') ... + .* KDEnorm; + % Add linear fit back into to origins, and/or linear weighting + % sMerge.scanOr = sMerge.scanOr *(1-originLinearFraction) ... + % + scanOrLinear; + sMerge.scanOr = sMerge.scanOr + scanOrLinear; + end + + % If pixels moved is below threshold, halt refinement + if (pixelsMoved/sMerge.numImages) < pixelsMovedThreshold + alignStep = refineMaxSteps + 1; + else + alignStep = alignStep + 1; + end +end + + + +% Remake images for plotting +if flagReportProgress == true + msg = sprintf('Recomputing images and plotting ...'); + fprintf([reverseStr, msg]); + reverseStr = repmat(sprintf('\b'),1,length(msg)); +end +for a0 = indsLoop + sMerge = SPmakeImage(sMerge,a0); +end + + +if flagPlot == 1 + imagePlot = sum(sMerge.imageTransform.*sMerge.imageDensity,3); + dens = sum(sMerge.imageDensity,3); + mask = dens>0; + imagePlot(mask) = imagePlot(mask) ./ dens(mask); + % Scale intensity of image + mask = dens>0.5; + imagePlot = imagePlot - mean(imagePlot(mask)); + imagePlot = imagePlot / sqrt(mean(imagePlot(mask).^2)); + + + % Plot results, image with scanline origins overlaid + figure(1) + clf + imagesc(imagePlot) + hold on + cvals = [1 0 0; + 0 .7 0; + 0 .6 1; + 1 .7 0; + 1 0 1; + 0 0 1]; + % plot(squeeze(sMerge.scanOr(1,2,:)),... + % squeeze(sMerge.scanOr(1,1,:)),... + % 'linewidth',4,'color',[1 1 1]) + for a0 = 1:sMerge.numImages + scatter(sMerge.scanOr(:,2,a0),sMerge.scanOr(:,1,a0),'marker','.',... + 'sizedata',25,'markeredgecolor',cvals(mod(a0-1,size(cvals,1))+1,:)) + end + hold off + axis equal off + colormap(gray(256)) + set(gca,'position',[0 0 1 1]) + caxis([-3 3]) % units of image RMS + + % Get final stats + % Idiff = mean(abs(sMerge.imageTransform ... + % - repmat(mean(sMerge.imageTransform,3),[1 1 sMerge.numImages])),3); + % meanAbsDiff = mean(Idiff(min(sMerge.imageDensity,[],3)>densityCutoff)) ... + % / mean(sMerge.scanLines(:)); + % sMerge.stats(end,1:2) = [refineMaxSteps meanAbsDiff]; + Idiff = mean(abs(sMerge.imageTransform ... + - repmat(mean(sMerge.imageTransform,3),[1 1 sMerge.numImages])),3); + meanAbsDiff = mean(Idiff(min(sMerge.imageDensity,[],3)>densityCutoff)) ... + / mean(abs(sMerge.scanLines(:))); + sMerge.stats(alignStep,1:2) = [alignStep-1 meanAbsDiff]; + + % Plot statistics + if size(sMerge.stats,1) > 1 + figure(2) + clf + plot(sMerge.stats(:,1),sMerge.stats(:,2)*100,'linewidth',2,'color','r') + xlabel('Iteration [Step Number]') + ylabel('Mean Absolute Difference [%]') + end +end + +if flagReportProgress == true + fprintf([reverseStr ' ']); +end +% toc +end + +function [score] = calcScore(image,xF,yF,dx,dy,intMeas) +% imageSample = interp2(image,yF+dy,xF+dx,'linear'); % Interp2 is too slow +imageSample = ... + image(sub2ind(size(image),xF, yF)) .*(1-dx).*(1-dy) ... + + image(sub2ind(size(image),xF+1,yF)) .*dx .*(1-dy) ... + + image(sub2ind(size(image),xF, yF+1)) .*(1-dx).*dy ... + + image(sub2ind(size(image),xF+1,yF+1)) .*dx .*dy; +score = sum(abs(imageSample-intMeas)); +end diff --git a/matlab/SPmerge03.m b/matlab/SPmerge03.m new file mode 100644 index 0000000..cd6f5fe --- /dev/null +++ b/matlab/SPmerge03.m @@ -0,0 +1,192 @@ +function [imageFinal,signalArray,densityArray] = SPmerge03(sMerge) +tic + +% Colin Ophus, National Center for Electron Microscopy, Molecular Foundry, +% Lawrence Berkeley National Laboratory, Berkeley, CA, USA. (Mar 2015). + +% Final scanning probe merge script. This script uses KDE and Fourier +% filtering to produce a combined image from all component scans. +% The Fourier weighting used (if flag is enabled) is cos(theta)^2, where +% theta is the angle from the scan direction. This essentially zeros out +% the slow scan direction, while weighting the fast scan direction at 100%. + +% Inputs: +% sMerge -struct containing data for STEM alignment from SPmerge02 +KDEsigma = 0.5 * 2; % -Gaussian sigma value used in kernel density estimator. (pixels) +% KDEcutoff = 4; % -Cutoff limit of the kernel in units of sigma. +upsampleFactor = 2; % -upsampling factor used in image generation (integer). +% Using a large upsampleFactor can be very slow. +sigmaDensity = 8; % -Smoothing sigma value for density estimation. (pixels) +boundary = 8; % -Thickness of windowed boundary. (pixels) +flagFourierWeighting = 1; % Set to true to enable cos(theta)^2 Fourier weights. +flagDownsampleOutput = 1; % Set to true to downsample output to the original +% resolution, as opposed to that of "upsampleFactor." + +% Outputs: +% imageCombine - final combined image. +% signalArray - image stack containing estimated image from each scan. +% densityArray - image stack containing estimated density of each scan. + + +% Initialize arrays +signalArray = zeros([sMerge.imageSize*upsampleFactor sMerge.numImages]); +densityArray = zeros([sMerge.imageSize*upsampleFactor sMerge.numImages]); +imageFinal = zeros(sMerge.imageSize*upsampleFactor); + + +% kernel generation in upsampled coordinates +x = makeFourierCoords(upsampleFactor*sMerge.imageSize(1),... + 1/sMerge.imageSize(1))'; +y = makeFourierCoords(upsampleFactor*sMerge.imageSize(2),... + 1/sMerge.imageSize(2)); +kernel = fft2(exp(-x.^2/(2*KDEsigma^2)) * exp(-y.^2/(2*KDEsigma^2))); +smoothDensityEstimate = fft2(exp(-x.^2/(2*sigmaDensity^2)) ... + *exp(-y.^2/(2*sigmaDensity^2)) ... + / (2*pi*sigmaDensity^2*upsampleFactor^2)); + + +% Loop over scans and create images / densities +t = repmat(1:size(sMerge.scanLines,2),[size(sMerge.scanLines,1) 1]); +for a0 = 1:sMerge.numImages + % Expand coordinates + x0 = repmat(sMerge.scanOr(:,1,a0),[1 size(sMerge.scanLines,2)]); + y0 = repmat(sMerge.scanOr(:,2,a0),[1 size(sMerge.scanLines,2)]); + xInd = x0(:)*upsampleFactor - (upsampleFactor-1)/2 ... + + (t(:)*sMerge.scanDir(a0,1))*upsampleFactor; + yInd = y0(:)*upsampleFactor - (upsampleFactor-1)/2 ... + + (t(:)*sMerge.scanDir(a0,2))*upsampleFactor; + xInd(:) = min(max(xInd,1),sMerge.imageSize(1)*upsampleFactor-1); + yInd(:) = min(max(yInd,1),sMerge.imageSize(2)*upsampleFactor-1); + + % Create bilinear coordinates + xIndF = floor(xInd); + yIndF = floor(yInd); + xAll = [xIndF xIndF+1 xIndF xIndF+1]; + yAll = [yIndF yIndF yIndF+1 yIndF+1]; + dx = xInd-xIndF; + dy = yInd-yIndF; + w = [(1-dx).*(1-dy) dx.*(1-dy) (1-dx).*dy dx.*dy]; + + % Generate image and density + image = sMerge.scanLines(:,:,a0); + signalArray(:,:,a0) = accumarray([xAll(:) yAll(:)],... + [image(:).*w(:,1); + image(:).*w(:,2); + image(:).*w(:,3); + image(:).*w(:,4)],... + sMerge.imageSize*upsampleFactor); + densityArray(:,:,a0) = accumarray([xAll(:) yAll(:)],... + [w(:,1); w(:,2); w(:,3); w(:,4)],... + sMerge.imageSize*upsampleFactor); +end + + +% Apply KDE to both arrays +for a0 = 1:sMerge.numImages + signalArray(:,:,a0) = real(ifft2(fft2(signalArray(:,:,a0)).*kernel)); + densityArray(:,:,a0) = real(ifft2(fft2(densityArray(:,:,a0)).*kernel)); +end +% Normalize image intensity by sampling density +sub = densityArray > 1e-8; +signalArray(sub) = signalArray(sub) ./ densityArray(sub); + + +% Calculate smooth density estimate, set max value to 2 to reduce edge +% effects, apply to images. +intensityMedian = median(sMerge.scanLines(:)); +for a0 = 1:sMerge.numImages + densityArray(:,:,a0) = real(ifft2(fft2(min(densityArray(:,:,a0),2)) ... + .* smoothDensityEstimate)); + densityArray(:,:,a0) = sin((pi/2)*min((bwdist( ... + densityArray(:,:,a0) < 0.5) ... + / (boundary*upsampleFactor)),1)).^2; + % Apply mask to each image + signalArray(:,:,a0) = signalArray(:,:,a0).*densityArray(:,:,a0) ... + + (1-densityArray(:,:,a0))*intensityMedian; +end + + + +% Combine scans to produce final image +if flagFourierWeighting == true + % Make Fourier coordinates + qx = makeFourierCoords(size(imageFinal,1),1); + qy = makeFourierCoords(size(imageFinal,2),1); + [qya,qxa] = meshgrid(qy,qx); + qTheta = atan2(qya,qxa); + + % Generate Fourier weighted final image + weightTotal = zeros(size(imageFinal)); + for a0 = 1:sMerge.numImages + % Filter array + thetaScan = atan2(sMerge.scanDir(a0,2),sMerge.scanDir(a0,1)); + qWeight = cos(qTheta-thetaScan).^2; +% qWeight = abs(cos(qTheta-thetaScan)); + + qWeight(1,1) = 1; + + % Filtered image + imageFinal = imageFinal ... + + fft2(signalArray(:,:,a0)).*qWeight; + weightTotal = weightTotal + qWeight; + end + imageFinal = real(ifft2(imageFinal ./ weightTotal)); + + % apply global density mask + density = prod(densityArray,3); + % density = max(sum(densityArray,3),1); % Temporary line for huge drifts + imageFinal = imageFinal.*density ... + + (1-density)*intensityMedian; +else + % Density weighted average + density = prod(densityArray,3); + % density = max(sum(densityArray,3),1); % Temporary line for huge drifts + imageFinal(:) = mean(signalArray,3).*density ... + + (1-density)*median(sMerge.scanLines(:)); + imageFinal = imageFinal / upsampleFactor^2; % Keep intensity constant +end + + +% Downsample outputs if required +if (flagDownsampleOutput == true) && upsampleFactor > 1 + % Fourier subsets for downsampling + xVec = [(1:(sMerge.imageSize(1)/2)) ... + (((-sMerge.imageSize(1)/2+1):0) ... + + sMerge.imageSize(1)*upsampleFactor)]; + yVec = [(1:(sMerge.imageSize(2)/2)) ... + (((-sMerge.imageSize(2)/2+1):0) ... + + sMerge.imageSize(2)*upsampleFactor)]; + + % Downsample output image + imageFinal = fft2(imageFinal); + imageFinal = real(ifft2(imageFinal(xVec,yVec))) / upsampleFactor^2; + if nargout > 1 + signalArray = fft2(signalArray); + signalArray = real(ifft2(signalArray(xVec,yVec,:))); + end + if nargout > 2 + densityArray = fft2(densityArray); + densityArray = real(ifft2(densityArray(xVec,yVec,:))); + end +end + + +figure(1) +clf +imagesc(imageFinal) +axis equal off +colormap(gray(256)) +set(gca,'position',[0 0 1 1]) + +toc +end + + +function [q] = makeFourierCoords(N,pSize) +% This function generates Fourier coordinates +if mod(N,2) == 0 + q = circshift(((-N/2):(N/2-1))/(N*pSize),[0 -N/2]); +else + q = circshift(((-N/2+.5):(N/2-.5))/((N-1)*pSize),[0 -N/2+.5]); +end +end \ No newline at end of file diff --git a/matlab/outdated_SPmerge_instructions.txt b/matlab/outdated_SPmerge_instructions.txt new file mode 100644 index 0000000..1965e46 --- /dev/null +++ b/matlab/outdated_SPmerge_instructions.txt @@ -0,0 +1,68 @@ +USAGE GUIDE +Scanning probe nonlinear drift correction code. + +Included MATLAB scripts used and their purposes: +SPmerge01.m -Initalize data structures, initial alignment, find reference position. +SPmerge02.m -Initial and primary refinement of scanline origins. +SPmerge03.m -Generate final corrected image. +SPmakeImage.m -Script for generating KDE images from new scanline orgins. +These scripts primarily work with a matlab struct with default name "sMerge." + +Author of all included scripts. +Colin Ophus, National Center for Electron Microscopy, Molecular Foundry, +Lawrence Berkeley National Laboratory, Berkeley, CA, USA. (Mar 2015). +cophus@gmail.com + +Usage steps for included synthetic data: image00deg and image90deg + +1. First verify alignment of images. All images should be passed into + SPmerge01.m with horizontal fast scan directions and without any rotation. + SPmerge01.m will rotate and align the data. Image alignment can be verifed + by plotting the data and ensuring horizontal scanlines for all images. + +Running these lines: +>> imagesc(image00deg); axis equal off; colormap(gray(256)) +>> imagesc(image90deg); axis equal off; colormap(gray(256)) +will show that both images have horizontal scanline directions. + +To ensure that the scanline orientation directions match those expected by +SPmerge01.m, we can plot both images at the same time, rotating the second +one by 90degrees: +>> imagesc([image00deg rot90(image90deg)]); axis equal off; colormap(gray(256)) + +Remember that the 90 degree rotation is just for visualization, and that the +images should not be rotated when passing them into SPmerge01.m. If the second +image were flipped 180 degrees in the two-image plot, we would need to use +a scanline direction of either -90 or 270 degrees in SPmerge01.m. + +2. Next, initialize the correction struct: +>> sMerge = SPmerge01([0 90],image00deg,image90deg); + +The initial alignment appears quite poor, due to the nonlinear drift. This +is not a problem, as the drift and alignment will be corrected in subsequent +steps. The yellow "+" mark indicates the starting position for the rough +alignment. If this point is inside a region of the images with no contrast, +you must manually adjust it by changing the coordinates of the struct in +your workspace, sMerge.ref, to the position in the image where the initial +alignment is best. + +3. Perform the drift correction: +>> sMerge = SPmerge02(sMerge); + +Running this script will perform first an initial rough drift correction, and +then the full nonlinear drift correction described in the manuscript. +After the alignment is completed, the script will make two plots: +-An alignment plot showing the current state of the correction. +-A plot showing the mean absolute difference at each iteration. This plot +is used to show convergence. In this synthetic dataset, the large drop after +iteration 1 indicates where a global phase shift took place. The error does +not decrease quite to zero, due to the undersampling effect of nonlinear drift. + +4. Generate a final reconstructed output image: +>> imageFinal = SPmerge03(sMerge); + +The resulting image should be fully drift corrected. To compare this image +to an image taken with perfect sampling and no drift, plot imageIdeal: +>> figure(1000); imagesc(imageIdeal); axis equal off; colormap(gray(256)) + +We see that the reconstruction has succeeded.