69 lines
2.6 KiB
Mathematica
69 lines
2.6 KiB
Mathematica
|
function [optrefimages] = fringeremoval(absimages, refimages, bgmask)
|
||
|
|
% FRINGEREMOVAL - Fringe removal and noise reduction from absorption images.
|
||
|
|
% Creates an optimal reference image for each absorption image in a set as
|
||
|
|
% a linear combination of reference images, with coefficients chosen to
|
||
|
|
% minimize the least-squares residuals between each absorption image and
|
||
|
|
% the optimal reference image. The coefficients are obtained by solving a
|
||
|
|
% linear set of equations using matrix inverse by LU decomposition.
|
||
|
|
%
|
||
|
|
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
|
||
|
|
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
|
||
|
|
%
|
||
|
|
% Syntax:
|
||
|
|
% [optrefimages] = fringeremoval(absimages,refimages,bgmask);
|
||
|
|
%
|
||
|
|
% Required inputs:
|
||
|
|
% absimages - Absorption image data,
|
||
|
|
% typically 16 bit grayscale images
|
||
|
|
% refimages - Raw reference image data
|
||
|
|
% absimages and refimages are both cell arrays containing
|
||
|
|
% 2D array data. The number of refimages can differ from the
|
||
|
|
% number of absimages.
|
||
|
|
%
|
||
|
|
% Optional inputs:
|
||
|
|
% bgmask - Array specifying background region used,
|
||
|
|
% 1=background, 0=data. Defaults to all ones.
|
||
|
|
% Outputs:
|
||
|
|
% optrefimages - Cell array of optimal reference images,
|
||
|
|
% equal in size to absimages.
|
||
|
|
%
|
||
|
|
|
||
|
|
% Dependencies: none
|
||
|
|
%
|
||
|
|
% Authors: Shannon Whitlock, Caspar Ockeloen
|
||
|
|
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
|
||
|
|
% S. Whitlock, Improved detection of small atom numbers through
|
||
|
|
% image processing, arXiv:1007.2136
|
||
|
|
% Email:
|
||
|
|
% May 2009; Last revision: 11 August 2010
|
||
|
|
|
||
|
|
% Process inputs
|
||
|
|
|
||
|
|
% Set variables, and flatten absorption and reference images
|
||
|
|
nimgs = size(absimages,3);
|
||
|
|
nimgsR = size(refimages,3);
|
||
|
|
xdim = size(absimages(:,:,1),2);
|
||
|
|
ydim = size(absimages(:,:,1),1);
|
||
|
|
|
||
|
|
R = single(reshape(refimages,xdim*ydim,nimgsR));
|
||
|
|
A = single(reshape(absimages,xdim*ydim,nimgs));
|
||
|
|
optrefimages=zeros(size(absimages)); % preallocate
|
||
|
|
|
||
|
|
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
|
||
|
|
k = find(bgmask(:)==1); % Index k specifying background region
|
||
|
|
|
||
|
|
% Ensure there are no duplicate reference images
|
||
|
|
%R=unique(R','rows')'; % comment this line if you run out of memory
|
||
|
|
|
||
|
|
% Decompose B = R*R' using singular value or LU decomposition
|
||
|
|
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
|
||
|
|
|
||
|
|
for j=1:nimgs
|
||
|
|
b=R(k,:)'*A(k,j);
|
||
|
|
% Obtain coefficients c which minimise least-square residuals
|
||
|
|
lower.LT = true; upper.UT = true;
|
||
|
|
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
|
||
|
|
|
||
|
|
% Compute optimised reference image
|
||
|
|
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
|
||
|
|
end
|