function out=distort(inp,distortion,base)
%DISTORT distort image.
%   DISTORT distorts an image by the specified amount.
%
%   OUT = DISTORT( INP , DISTORTION) returns an image which is the same
%   size as the input image INP, but which has been distorted by the
%   amount DISTORTION.  If DISTORTION=0, the output image is the same
%   as the input.  If DISTORTION<0 then  the center of the image is magnified
%   and the sides bow outward.  For DISTORTION>0, the center of the image
%   shrinks and the sides bow inward.
%
%   OUT = DISTORT( INP, DISTORTION, BASE) zooms the final image by
%   the factor BASE.
%
%   If no output is specified, the distorted image is displayed.
%
%   For example, DISTORT( INP, 0.2) displays the image with 20% distortion.
%
%   Class Support
%   -------------
%   The input image can be of class uint8 or double. The output
%   image is of the same class as the input image.

%   John Loomis  25 Feb 2000

error(nargchk(2,3,nargin));      % Check for 2 input arguments.
[vlen,hlen,colors]=size(inp);    % Size of input image.

% Define output image as same type as input (double or uint8).
output=zeros(vlen,hlen,colors);
if isa(inp,'uint8')
   output=im2uint8(output);
end;

xcen=hlen/2.0;                 % Horizontal center, in pixels.
ycen=vlen/2.0;                 % Vertical center, in pixels.
ds  = 4.0/(hlen+vlen);         % intermediate circle (inscribed in square)
%ds = 2.0/min(hlen,vlen);      % inscribed circle
%ds = 2.0/sqrt(hlen^2+vlen^2); % circumscribed circle

% Base magnification factor.
if (nargin<3)
   if (distortion<0.0)
      rbase = min(hlen,vlen)*ds/2;
   else
      rbase = sqrt(hlen^2+vlen^2)*ds/2;
	end;
	base = 1.0/(1.0+distortion*rbase^2);
end;


% x and y are arrays whose elements are x and y positions relative to the
% center, normalized to 1 or -1 at the edges.
[x,y]=meshgrid(linspace(0.5-xcen,hlen-0.5-xcen,hlen)*ds,linspace(0.5-ycen,vlen-0.5-ycen,vlen)*ds);

rf = sqrt(x.^2+y.^2);           % Normalized dist. from center.
m = base*mag(rf/base,distortion);    % Magnification as function of final radius
nx=round(xcen+x./(ds*m));       % Array of distorted x indices.
ny=round(ycen+y./(ds*m));       % Array of distorted y indices.

% Look for out-of-range indices for input image.
bad = find(nx<1 | nx>hlen | ny<1 | ny>vlen);
% Set to an arbitrary in-range value.
nx(bad) = 1;
ny(bad) = 1;

% reshape arrays as one-dimensional vectors
inp = reshape(inp,vlen*hlen,colors);
output = reshape(output,vlen*hlen,colors);
idx = reshape(ny+(nx-1)*vlen,1,vlen*hlen);

% set output image to indexed values from input image
output = inp(idx,:);
% set out-of-range values to black
output(bad,:)=0;
% restore output to two-dimensional matrix
output = reshape(output,vlen,hlen,colors);


% If no output argument was given, display result.
if nargout==0
   imshow(output);
else
   out=output;
end;

function m = mag(rfinal,distortion)
%   MAG calculates the magnification at a specified radius in the
%   output image.
%
%         rfinal = mag * rinput
%         where mag = 1.0 + distortion * rinput^2
%
%   Given rfinal, find rinput and then calculate and return the magnification.
%
%   This is a cubic equation of the form rfinal = rinput + distortion * rinput^3
%   which can be solved analytically. However, Newton's method converges to the
%   solution in about three iterations and the resulting calculation avoids the
%   use of cube roots.
%

%
%  no valid inverse solution exists beyond rmax (r' vs r has horizontal slope)
%  so we limit the range of rfinal and later set the magnification to near zero.
%
if (distortion<0.0)
	rmax = 2/sqrt(-27.0*distortion);
	idx = find(rfinal>rmax);
	rfinal(idx)=rmax;
end;

r = rfinal; % initial value
for i=1:3
	rsq = r.^2;
	r = (rfinal + 2*distortion*rsq.*r)./(1+3*distortion*rsq);
end;

m = 1.0+distortion*r.^2;

if (distortion<0.0)
	m(idx)=0.001;
end

Published with MATLAB® R2017a