Code covered by the BSD License  

Highlights from
finterp

image thumbnail

finterp

by

Michael Hawks

 

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