% opticalLowpassInterpolation Filter and interpolate microscopy image while avoiding artifacts.
% Images acquired with lenses (microscope, camera) can possess fine
% features only upto the spatial frequency cut-off of the optics. Any fine
% features beyond that are due to noise. A simple and effective noise
% removal strategy is to remove the above-cutoff spatial frequencies.
% opticalLowpassInterpolation implements filtering of the imaging data with
% super-gaussian filter in such a way that filtering artifacts are
% The image can be optionally interpolated in frequency domain while
% maintaining physical variation in intensity.
% USAGE: [out, outpix] =
% out - filtered and interpolated output image.
% outpix - pixel-size in the output image.
% in - raw image
% inpix - pixel-size in the input image (numerical
% value in your chosen units of distance).
% fcut - spatial-frequency cutoff
% (numerical value in the inverse units of distance).
% IntFac - Interpolation factor.
% Author and Copyright: Shalin Mehta (www.mshalin.com)
% License: BSD
% Version history: April 2012, initial implementation with gaussian filter.
% August 2012, use super-gaussian filter.
% Feb 10, 2013, added functionality to interpolate in
% frequency domain.
% Feb 27, 2013, Order of super-gaussian is intelligently estimated from the sampling.
% Resolved bug that caused the center of
% intensity to shift (phase-shift in space)
% when interpolation is used.
%% Pad the input image to avoid edge artifacts.
%% Obtain the spectrum with DC at center of image.
%% Pad the spectrum to achieve spatial interpolation.
%% Generate frequency grid for padded spectrum.
mcut=1/(2*outscale); %Cut-off of frequency grid.
[ylen, xlen]=size(specpad); % The first return value is the height and the second is the width.
%% Estimate the sharpness of supergaussian and generate filter over above grid.
% Equation of supergaussian is a=exp(-(f/fo)^(1/2n)).
% I choose the transition of the filter response from 0.99 to 0.01 to
% be sampled over 5 frequency bins.
% The filter response is equal to a, when
% f/fo = -Log[a]^(1/2n).
% The super-gaussian filter is = 1/e when f=fo.
% The transition region normalized by fcut is thus given by,
% mtrans/fcut = -Log[0.01])^(1/2n) + Log[0.99])^(1/2n).
% Using Mathematica, above is simplified to
% Above value is computed for various integer values of n and tabulated
% here. To estimate n from given mtrans/fcut, we just find the index
% where required mtrans/fcut is closest to the entry in the table.
% Size of frequency step in radial direction.
dm=sqrt( (mx(2)-mx(1))^2 + (my(2)-my(1))^2 );
% We want transition period to be atleast 3 frequency bins and transition region should start at fcut.
% transition period normalized by fcut.
% The suitable order of super-gaussian.
% Use super-gaussian as our frequency filter.
%% Filter the padded spectrum and obtain interpolated spatial image.
% Compute filtered, interpolated, zero-padded output.
% Crop the center of the output image.