Code covered by the BSD License  

Highlights from
Droste Effect Tool

image thumbnail

Droste Effect Tool

by

 

13 May 2009 (Updated )

Apply the Droste Effect to an image using a GUI or function call.

droste_effect(img, reg, varargin)
function img_out = droste_effect(img, reg, varargin)
% Apply the Droste Effect to an image.
%
% The Droste Effect is when an image recursively includes itself. The basic
% method is to insert a scaled copy of the image, but more interesting
% results can be achieved via conformal mapping.
%
% img_out = droste_effect(img, reg) will apply the Droste Effect to matrix
% img within the region specified by reg. The img matrix can be 2-D or 3-D,
% corresponding to a grayscale or RGB image. The reg vector specifies
% [left, top, width, height] in pixels.
%
% img_out = droste_effect(img, reg, 'p1', v1, 'p2', v2, ...) specifies
% optional parameters (pn) and values (vn) associated with the Droste
% Effect. The table below describes these parameters.
%
%    Param  | Definition
% ==========|==============================================================
%      zoom | Scale. Default = 1.0x.
%    hshift | Horizontal shift, normalized [0...1]. Default = 0.
%    vshift | Vertical shift, normalized [0...1]. Default = 0.
%    rotate | Rotate. Default = 0 degrees.
%   nspiral | Number of spirals. Default = 1.
%   ncopies | Number of copies of image per spiral rotation. Default = 1.
%  maxrecur | Max recursion allowed to compute output. Default = 10.
%     width | Width of output in pixels. Default = same resolution as input.
%  aafactor | Anti-aliasing factor. Default = 1, higher = smoother. (0 = disable)
%
% Example:
%     load penny
%     reg = [45 45 40 40];
%     img_out = droste_effect(P, reg, 'width', 400);
%     figure
%     imshow(round(img_out),copper(256))
%
% Running out of memory? Reduce the size of the input image before calling
% this function.
%
% by Steve Hoelzer
%
% 2009-05-13
%
% Inspired by Jon McLoone:
%    http://blog.wolfram.com/2009/04/24/droste-effect-with-mathematica/

% crop image to make things simpler later
img = crop_image(img,reg);

% default params
zoom = 1;
hshift = 0;
vshift = 0;
rotate = 0;
nspiral = 1;
ncopies = 1;
maxrecur = 10;
width = size(img,2);
aafactor = 1;

% process inputs an overwrite defaults
nargchk(2,20,nargin);
valid_param_names = {'zoom', 'hshift', 'vshift', 'rotate', 'nspiral', ...
    'ncopies', 'maxrecur', 'width', 'aafactor'};
for n = 1:2:length(varargin)
    switch varargin{n}
        case valid_param_names
            eval(sprintf('%s = %f;',varargin{n},varargin{n+1}))
        otherwise
            error('Invalid input parameter ''%s''',varargin{n})
    end
end

% init some vars
[r,c,d] = size(img); %#ok<NASGU>
aspect = r/c; % height/width aspect ratio

% scaling factor
xscale = reg(3) / c;
yscale = reg(4) / r;
dscale = max(xscale, yscale);

% generate zero-centered normalized coordinates [-1...1] for output image
ci = width;
ri = round(ci*aspect);
[xi,yi] = meshgrid(1:ci,1:ri);
xi = (xi/ci)*2 - 1;
yi = (yi/ri)*2 - 1;

% apply transform to coordinates
[xi, yi] = xform(xi, yi, zoom, hshift, vshift, ...
    rotate, ncopies, nspiral, dscale);

% tweak coordinates to get self-replicating effect
[xi, yi] = replicate(xi, yi, dscale, aspect, maxrecur);

% change to non-zero-centered normalized coordinates [0...1]
xi = (xi + 1) / 2;
yi = (yi + 1) / 2;

% convert img type to double for the math below
img = double(img);

% apply anti-aliasing filter to source image
aafactor = aafactor / zoom; % adjust by zoom level
img = filter_image(img, ri, aafactor);

% interpolate input image with transformed coordinates to get output image
img_out = interp_image(img, xi, yi);


% -------------------------------------------------------------------------
function img = crop_image(img,reg)
% crop input image so droste region is centered and aspect ratio is the
% same as droste region aspect ratio

% center of droste region
xc = reg(1) + reg(3)/2;
yc = reg(2) + reg(4)/2;

% calculate max distance from center to edge
[r,c,d] = size(img); %#ok<NASGU>
xd = min(xc-1, c-xc);
yd = min(yc-1, r-yc);

% make aspect ratio the same as droste region aspect ratio
daspect = reg(4)/reg(3);
xd = min(xd, yd/daspect);
yd = min(yd, xd*daspect);

% crop image
img = img(round(yc-yd):round(yc+yd), round(xc-xd):round(xc+xd), :);


% -------------------------------------------------------------------------
function [x, y] = xform(x, y, zoom, xshift, yshift, ...
    rotate, ncopies, nspiral, dscale)

% apply transform
base = (x + i*y) + (-1*xshift + i*yshift);
expo = ncopies + i*(nspiral * log(dscale) / (2*pi));
p = (1/zoom) * exp(i*rotate*pi/180) * base.^expo;
p(base==0) = 0; % 0^i gives NaN, so set to zero

% separate x and y for output
x = real(p);
y = imag(p);


% -------------------------------------------------------------------------
function [x, y] = replicate(x, y, dscale, aspect, maxrecur)

% reached recursion limit?
if maxrecur == 0
    return
end

% if coordinates outside image, scale down
dndx = find(abs(x) > 1 | abs(y) > 1);
x(dndx) = x(dndx) * dscale;
y(dndx) = y(dndx) * dscale;

% if coordinates inside droste region, scale up
undx = find(abs(x) < dscale & abs(y) < dscale);
x(undx) = x(undx) / dscale;
y(undx) = y(undx) / dscale;

% recheck coordinates that were just scaled
if length(undx) > 1 || length(dndx) > 1
    ndx = unique([undx; dndx]);
    [x(ndx), y(ndx)] = replicate(x(ndx), y(ndx), dscale, aspect, maxrecur-1);
end


% -------------------------------------------------------------------------
function img = filter_image(img, ri, aafactor)

% check if anti-aliasing filter is disabled
if aafactor == 0
    return
end

% calculate cutoff freq for anti-aliasing filter
[r,c,d] = size(img);
cutoff = ri / r; % due to output/input size ratio
cutoff = cutoff / aafactor; % adjust by aafactor

% check if cutoff is valid
if cutoff >= 1
    return
end

% apply anti-aliasing filter to source image
n = round(5/cutoff) * 2;
h = fir1(n,cutoff);
img = [img, zeros(r,n/2,d); zeros(n/2,c+n/2,d)]; % pad with zeros for group delay
img = filter(h,1,img,[],1); % filter cols
img = filter(h,1,img,[],2); % filter rows
img = img(n/2+1:end,n/2+1:end,:); % trim to account for group delay


% -------------------------------------------------------------------------
function img_out = interp_image(img, xi, yi)

% generate normalized coordinates [0...1] for source image
[r,c,d] = size(img);
% Note: Coordinate grid expanded from [1...c, 1...r] to [0...c+1, 0...r+1]
% so that output coordinates always fit inside this grid. Otherwise,
% interpolation will not work and some pixels will be filled with NaNs.
[x,y] = meshgrid(linspace(0,c+1,c),linspace(0,r+1,r));
x = x / c;
y = y / r;

% init the output image
[ri,ci] = size(xi);
img_out = zeros(ri,ci,d);

% interpolate pixels
for i = 1:d
    img_out(:,:,i) = interp2(x, y, img(:,:,i), xi, yi, 'linear');
end

Contact us