Code covered by the BSD License  

Highlights from
finterp

image thumbnail

finterp

by

 

Use Fourier interpolation to increase sampling of an image. Also includes an optional filter.

finterp(f,newsize,apod)
function f1 = finterp(f,newsize,apod)
% NEWIMAGE = finterp(IMAGE,NEWSIZE,APOD)
% Resizes the image IMAGE using Fourier interpolation.  
% 
% NEWSIZE is the desired size of the interpolated image.  This can be a
% 2-element vector for non-square images.  If NEWSIZE is a scalar, the 
% output image will be square (i.e. NEWSIZE=256 will create a 256x256 image). 
%
% APOD is an optional parameter to apodize the image (by a Hanning filter)
% to reduce ringing in the output image
%   APOD = 0 or blank (DEFAULT) -- no apodization
%   0 < APOD < 1 -- Apodize
%   
%   The Hanning filter is = 1 for x<xL
%                         = 0 for x>xH
%           between xL and xH, the filter transmission is a half-cycle of cos^2
%   For this function, xH is assumed to be the highest frequency in f, and
%   xL = APOD * xH
%
% Michael Hawks, Department of Engineering Physics, Air Force Institute of
% Technology
% 5 July 2013

if nargin<3, apod=0; 
elseif apod>=1, error('ERROR: APOD must be between 0 and 1');
end
if ndims(f)~=2, error('ERROR: FINTERP only defined for 2-D arrays'); beep; end
newsize=[1,1].*newsize;   % ensure newsize is 2D -- if input is one number, this makes a square array
if any(size(f)>newsize), f1=f; return; end;

maxI = max(max(f)); 
minI = min(min(f));

pad = newsize - size(f); % number of elements to add
pad1 = floor(pad/2);     % portion to add to right/top
pad2 = pad - pad1;       % remainder add to left/bottom

F = fft2(f);
Fpad = padarray(fftshift(F),pad1,'pre');
Fpad = padarray(Fpad,pad2,'post');

if apod>0
    [Nx,Ny]=size(Fpad);  
    H = 0.5.*size(F);  % use this as the high-freq cutoff (from -H to +H)
    if isodd(Nx);   x = -floor(Nx/2):floor(Nx/2);
    else            x = -((Nx/2)-1):(Nx/2);
    end
    if isodd(Ny);   y = -floor(Ny/2):floor(Ny/2);
    else            y = -((Ny/2)-1):(Ny/2);
    end
    filt = hanning(x,apod*H(1),H(1))' * hanning(y,apod*H(2),H(2)); 
    Fpad = Fpad.*filt;
end

ft = ifft2(ifftshift(Fpad));
f1 = (ft.*conj(ft)).^0.5; % assume we want a real-valued image out, but asymmetric padding could add an
                          % imaginary componenet so we approximate by SQRT(F F*)

% rescale/normalize the image.  This is kind of a hack, but it should take
% care of everything
f1 = f1 - min(min(f1));  f1 = f1./max(max(f1)); 
f1 = f1 .* maxI + minI;


% ...............
function i = isodd(x)
i = ( floor(x/2) ~= (x/2) );


% ...............
function f = hanning(x,xL,xH)

if ndims(squeeze(x))==1, x=1:x; dx=1;
else dx=abs(x(2)-x(1)); 
end

f=zeros(1,length(x));

if x(1) >= 0;  % one-sided filter
    jL = find(abs(x-xL)<dx,1); if isnan(jL); jL=1; end
    jH = find(abs(x-xH)<dx,1); if isnan(jH); jH=length(x); end
    f(1:jL)=1;
    f(jH:end)=0;
    k = jL:jH;      f(k)=0.5 + 0.5*cos(pi*(k-jL)/(jH-jL));
else
    j1 = find(abs(x+xL)<dx,1);  if isnan(j1); j1=1; end
    j2 = find(abs(x-xL)<dx,1);  if isnan(j2); j2=length(x); end
    j3 = find(abs(x+xH)<dx,1);  if isnan(j3); j3=1; end
    j4 = find(abs(x-xH)<dx,1);  if isnan(j4); j4=length(x); end

    f(1:j3)=0;
    f(j4:end)=0;
    f(j1:j2)=1;
    k = j2:j4;      f(k)=0.5 + 0.5*cos(pi*(k-j2)/(j4-j2)); 
    k = j3:j1;      f(k)=0.5 - 0.5*cos(pi*(k-j3)/(j3-j1));
end

Contact us