adding some minor helper functions

compute_dead_leaves_image.m

@@ -0,0 +1,100 @@
+function M = compute_dead_leaves_image(n,sigma,options)
+
+% compute_dead_leaves_image - compute a random image using the dead-leaves model
+%
+%   M = compute_dead_leaves_image(n,sigma,options);
+%
+%   n is the size of the image.
+%   sigma>0 control the repartition of the size of the basic shape:
+%       sigma --> 0  gives more uniform repartition of shape
+%       sigma=3 gives a nearly scale invariant image.
+%   options.nbr_iter put a bound on the number of iteration.
+%   options.shape can be 'disk' or 'square'
+%
+% References :
+%   Dead leaves correct simulation :
+%    http://www.warwick.ac.uk/statsdept/staff/WSK/dead.html
+%
+%   Mathematical analysis
+%    The dead leaves model : general results and limits at small scales
+%    Yann Gousseau, Fran¸cois Roueff, Preprint 2003
+%
+%   Scale invariance in the case sigma=3
+%    Occlusion Models for Natural Images: A Statistical Study of a Scale-Invariant Dead Leaves Model
+%     Lee, Mumford and Huang, IJCV 2003.
+%
+%   Copyright (c) 2005 Gabriel Peyré
+
+options.null = 0;
+
+if nargin<2
+    sigma = 3;
+end
+if isfield(options,'rmin')
+    rmin = options.rmin;
+else
+    rmin = 0.01;    % maximum proba for rmin, shoult be >0
+end
+if isfield(options,'rmax')
+    rmax = options.rmax;
+else
+    rmax = 1;    % maximum proba for rmin, shoult be >0
+end
+if isfield(options,'nbr_iter')
+    nbr_iter = options.nbr_iter;
+else
+    nbr_iter = 5000;
+end
+if isfield(options,'shape')
+    shape = options.shape;
+else
+    shape = 'disk';
+end
+
+M = zeros(n)+Inf;
+
+x = linspace(0,1,n);
+[Y,X] = meshgrid(x,x);
+
+
+% compute radius distribution
+k = 200;        % sampling rate of the distrib
+r_list = linspace(rmin,rmax,k);
+r_dist = 1./r_list.^sigma;
+if sigma>0
+    r_dist = r_dist - 1/rmax^sigma;
+end
+r_dist = rescale( cumsum(r_dist) ); % in 0-1
+
+m = n^2;
+
+for i=1:nbr_iter
+
+
+	% compute scaling using inverse mapping
+    r = rand(1);
+	[~,I] = min( abs(r-r_dist) );
+    r = r_list(I);
+
+    x = rand(1);    % position
+    y = rand(1);
+    a = rand(1);    % albedo
+
+    if strcmp(shape, 'disk')
+        I = find(isinf(M) & (X-x).^2 + (Y-y).^2 < r^2);
+    else
+        I = find(isinf(M) & abs(X-x)<r & abs(Y-y)<r );
+    end
+
+    m = m - length(I);
+    M(I) = a;
+
+    if m==0
+        % the image is covered
+        break;
+    end
+end
+
+% remove remaining background
+I = find(isinf(M));
+M(I) = 0;

gabor.m

@@ -15,6 +15,7 @@
 %
 %   Parameters (Default):
 %       theta (2*pi*rand): orientation of the gabor patch
+%       phi   (2*pi*rand): phase of the sinusoid
 %       lambda (20): spatial wavelength
 %       Sigma (10): standard deviation of gaussian window
 %       width (256): width of generated image
@@ -25,7 +26,7 @@
 % EXAMPLE
 %   [x y F] = gabor;
 %   pcolor(x,y,F); axis image;
-%   shading('interp'); colormap copper;
+%   shading('interp'); colormap gray;
 %
 % VERSION 1.1, Thu Jan 10 17:47:35 2013     Added ability to modify the phase of the sinusoid
 % VERSION 1.0, Thu Jul 12 09:47:52 2012     Initial version

gaborbank.m

@@ -0,0 +1,52 @@
+%Generates a bank of gabor filters
+%
+% AUTHOR: Niru Maheswaranathan
+%         [email protected]
+
+width = 64;
+height = 64;
+
+thetas = 0;
+phis = linspace(0,2*pi,11);
+lambdas = round(logspace(1, 2, 6));
+Sigmas = 10;
+pxs = 0.5;
+pys = 0.5;
+%pxs = 0.1:0.2:0.9;
+%pys = 0.1:0.2:0.9;
+
+numGabors = length(thetas)*length(phis)*length(lambdas)*length(Sigmas)*length(pxs)*length(pys);
+gabors = zeros(numGabors, width*height);
+gidx = 1;
+
+for tidx = 1:length(thetas)
+    for pidx = 1:length(phis)
+        for lidx = 1:length(lambdas)
+            for sidx = 1:length(Sigmas)
+                for xidx = 1:length(pxs)
+                    for yidx = 1:length(pys)
+
+                        % generate gabor
+                        [x y F] = gabor('theta',  thetas(tidx),  ...
+                                        'phi',    phis(pidx),    ...
+                                        'lambda', lambdas(lidx), ...
+                                        'Sigma',  Sigmas(sidx),  ...
+                                        'px',     pxs(xidx),     ...
+                                        'py',     pys(yidx),     ...
+                                        'width',  width,         ...
+                                        'height', height         ...
+                                        );
+
+                        % store the function
+                        gabors(gidx,:) = F(:)';
+                        gidx = gidx + 1;
+
+                        % update progress
+                        progressbar(gidx, numGabors, 50);
+
+                    end
+                end
+            end
+        end
+    end
+end

helper_functions/makeWindow.m

@@ -0,0 +1,23 @@
+%Window generator for tapered FFTs
+%
+% Generates a 2D tapered windowing function. Used by multiplying with an image
+% before taking a 2D Fourier transform to reduce artifacts due to boundary effects
+% Thu Jan 10 13:57:26 2013
+%
+%
+% USAGE:
+%   win = makeWindow(M, N)
+%
+% PARAMETERS:
+%   [M N] is the size of the image
+%
+% VERSION 1.0, Fri Jan 11 15:35:15 2013     Initial version
+%
+% AUTHOR: Niru Maheswaranathan
+%         [email protected]
+
+function win = makeWindow(M, N)
+
+    w1 = cos(linspace(-pi/2, pi/2, M));
+    w2 = cos(linspace(-pi/2, pi/2, N));
+    win = w1' * w2;

helper_functions/normalize.m

@@ -0,0 +1,22 @@
+%Normalize (z-score) an image patch
+%
+% Z-scores the pixels in an image patch. Useful pre-processing step.
+%
+% USAGE:
+%   zimg = normalize(img)
+%
+% PARAMETERS:
+%   img: the image to be normalized
+%
+% RETURNS:
+%   zimg: the same image with zero mean and unit variance
+%
+%
+% VERSION 1.0, Thu Jan 10 14:20:31 2013         Initial version
+%
+% AUTHOR: Niru Maheswaranathan
+%         [email protected]
+
+function zimg = normalize(img)
+
+    zimg = (img - mean(img(:)))./std(img(:));

helper_functions/removeDC.m

@@ -0,0 +1,20 @@
+%removeDC: Removes DC component from a signal
+%
+% Centers the *rows* of a signal x (vector or matrix)
+%
+% USAGE:
+%   [y mu] = removeDC(x);
+%
+% RETURNS:
+%   y: the centered signal(s)
+%  mu: the mean that was removed
+%
+% VERSION 1.0, Tue Nov  6 10:40:24 2012     Initial version
+%
+% AUTHOR: Niru Maheswaranathan
+%         [email protected]
+
+function [y mu] = removeDC(x)
+
+    mu = mean(x,2);
+    y = x - mu*ones(1,size(x,2));

helper_functions/rescale.m

@@ -0,0 +1,25 @@
+function y = rescale(x,a,b)
+
+% rescale - rescale data in [a,b]
+%
+%   y = rescale(x,a,b);
+%
+%   Copyright (c) 2004 Gabriel Peyr?
+
+if nargin<2
+    a = 0;
+end
+if nargin<3
+    b = 1;
+end
+
+m = min(x(:));
+M = max(x(:));
+
+if M-m<eps
+    y = x;
+else
+    y = (b-a) * (x-m)/(M-m) + a;
+end
+
+

helper_functions/spike.m

@@ -0,0 +1,44 @@
+% creates a spike image
+% Thu Dec  6 11:41:44 2012
+% NOTE: This function is from Lawrence Cormack
+
+function s=spike(alpha, r, c)
+%
+% SPIKE -   creates a spike image, an array in which the
+%           value is inversely proportional to the distance
+%           from the center.
+%
+% usage:
+%     S = SPIKE(ALPHA, R, C) produces an RxC array, the value
+%                           of which falls off as 1/(d^ALPHA) where d is
+%                           the euclidean distance from the array center.
+%
+%     Useful for making 1/f^alpha noise images.
+%
+% See also:  oneoverf(), spike3d(), oneoverf3d()
+%
+% Lawrence K. Cormack
+
+% history:
+% ??/??/2000    lkc     wrote it
+% 11/12/2000    lkc     updated for any rxc size
+% 11/12/2000    lkc     added odd/even dependent shift to compensate for the
+%                       location of the DC after fftshift()
+% 3/31/2002     lkc     added the exponent, alpha, to this function (where it belongs)
+%                       rather than dealing with it in oneoverf()
+% 4/1/2002      lkc     changed the distance computations to a more intuitive and
+%                       straight-forward method using meshgrid().
+
+if nargin == 2 c=r; end
+
+s = ones(r, c);
+
+y = linspace(-r./2, r./2, r);
+x = linspace(-c./2, c./2, c);
+
+[X, Y] = meshgrid(x, y);
+
+dists = sqrt(X.^2 + Y.^2);
+s = s./dists.^alpha;
+
+

makeWindow.m

@@ -0,0 +1,23 @@
+%Window generator for tapered FFTs
+%
+% Generates a 2D tapered windowing function. Used by multiplying with an image
+% before taking a 2D Fourier transform to reduce artifacts due to boundary effects
+% Thu Jan 10 13:57:26 2013
+%
+%
+% USAGE:
+%   win = makeWindow(M, N)
+%
+% PARAMETERS:
+%   [M N] is the size of the image
+%
+% VERSION 1.0, Fri Jan 11 15:35:15 2013     Initial version
+%
+% AUTHOR: Niru Maheswaranathan
+%         [email protected]
+
+function win = makeWindow(M, N)
+
+    w1 = cos(linspace(-pi/2, pi/2, M));
+    w2 = cos(linspace(-pi/2, pi/2, N));
+    win = w1' * w2;

normalize.m

@@ -0,0 +1,22 @@
+%Normalize (z-score) an image patch
+%
+% Z-scores the pixels in an image patch. Useful pre-processing step.
+%
+% USAGE:
+%   zimg = normalize(img)
+%
+% PARAMETERS:
+%   img: the image to be normalized
+%
+% RETURNS:
+%   zimg: the same image with zero mean and unit variance
+%
+%
+% VERSION 1.0, Thu Jan 10 14:20:31 2013         Initial version
+%
+% AUTHOR: Niru Maheswaranathan
+%         [email protected]
+
+function zimg = normalize(img)
+
+    zimg = (img - mean(img(:)))./std(img(:));