function [y,quincunx_filters] = perform_quincunx_wavelet_transform(x,Jmin,dir,options)

% perform_quincunx_wavelet_transform - compute quincunx transform
%
% Forward transform
%   [MW,options.quincunx_filters] = perform_quincunx_wavelet_transform(M,Jmin,+1,options);
% Backward transform
%   M = perform_quincunx_wavelet_transform(MW,Jmin,-1,options);
%
% This is a simple wrapper to the code of
%   Dimitri Van De Ville
%   Biomedical Imaging Group, BIO-E/STI
%   Swiss Federal Institute of Technology Lausanne
%   CH-1015 Lausanne EPFL, Switzerland
%
% For more information, see
%   Isotropic Polyharmonic B-Splines: Scaling Functions and Wavelets
%   D. Van De Ville, T. Blu, M. Unser
%   IEEE Transactions on Image Processing, vol. 14, no. 11, pp. 1798-1813, November 2005.
%
% options.type is the type of transform:
%    - McClellan: O/B/D,
%    - Polyharmonic Rabut: Po/Pb/Pd
%    - Polyharmonic isotropic: PO/PB/PD
% options.gamma is the degree of the wavelet transform

options.null = 0;


% Setup parameters of wavelet transform
% -------------------------------------
% 1. Type of wavelet transform
%    - McClellan: O/B/D,
%    - Polyharmonic Rabut: Po/Pb/Pd
%    - Polyharmonic isotropic: PO/PB/PD
if isfield(options, 'type')
    type = options.type;
else
    type='PO';
end

% 2. Degree of the wavelet transform
%    (order gamma polyharmonic=degree alpha+2)
if isfield(options, 'gamma')
    gamma = options.gamma;
else
    gamma = 4;
end
alpha=gamma-2;

% 3. Number of iterations
n = size(x,1);
Jmax = log2(n)-1;
J = 2*(Jmax-Jmin+1);


% Precalculate filters
% --------------------
if isfield(options, 'quincunx_filters')
    quincunx_filters = options.quincunx_filters;
    FA = quincunx_filters{1};
    FS = quincunx_filters{2};
else
    [FA,FS]=FFTquincunxfilter2D([size(x)],alpha,type);
    quincunx_filters{1} = FA;
    quincunx_filters{2} = FS;
end

% Perform wavelet analysis
% ------------------------
if dir==1
    y = FFTquin2D_analysis(x,J,FA); y=real(y);
else
    y = FFTquin2D_synthesis(x,J,FS);
end


function A0=autocorr2d(H)

% AUTOCORR2D
%       Iterative frequency domain calculation of the 2D autocorrelation
%       function. A=autocorr2d(H) computes the frequency response of
%       the autocorrelation filter A(exp(2*i*Pi*nu)) corresponding to the
%       scaling function with refinement filter H.
%       Please note that the 2D grid of the refinement filter (nu),
%       corresponds to a uniformly sampled grid, twice as coarse
%       (in each dimension) as the given filter H.
%
% Please treat this source code as confidential!
% Its content is part of work in preparation to be submitted:
%
% T. Blu, D. Van De Ville, M. Unser, ''Numerical methods for the computation
% of wavelet correlation sequences,'' SIAM Numerical Analysis.
%
% Biomedical Imaging Group, EPFL, Lausanne, Switzerland.


len1=size(H,1)/2;
len2=size(H,2)/2;

% stop criterion
crit=1e-4;
improvement=1.0;

% initial "guess"
A0=ones(len1,len2);

% calculation loop

while improvement>crit,

    % sinc interpolation
    if 1,
        Af=fftshift(ifftn(A0));
        Af(1,:)=Af(1,:)/2;
        Af(:,1)=Af(:,1)/2;
        Af(len1+1,1:len2)=Af(1,len2:-1:1);
        Af(1:len1,len2+1)=Af(len1:-1:1,1);
        Ai=zeros(2*len1,2*len2);
        Ai(len1-len1/2+1:len1+len1/2+1,len2-len2/2+1:len2+len2/2+1)=Af;
        Ai=fftn(ifftshift(Ai));
    else
        Ai=interp2(A0,1,'nearest'); Ai(:,end+1)=Ai(:,end); Ai(end+1,:)=Ai(end,:);
    end;

    % recursion
    A1=Ai(1:len1,1:len2).*H(1:len1,1:len2)+Ai(len1+1:2*len1,1:len2).*H(len1+1:2*len1,1:len2)+Ai(1:len1,len2+1:2*len2).*H(1:len1,len2+1:2*len2)+Ai(len1+1:2*len1,len2+1:2*len2).*H(len1+1:2*len1,len2+1:2*len2);

    improvement=mean(abs(A1(:)-A0(:)));

    A0=A1;

end;


function y2=FFTquin2D_analysis(x,J,FA);

% FFTQUIN2D_ANALYSIS Perform 2D quincunx wavelet analysis
% 	y2 = FFTquin2D_analysis(x,J,FA)
%
% 	Input:
%       x  = input data
%       J  = number of iterations
%       FA = analysis filters
%
%       Output:
%       y2 = result (folded)
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland

% Dimensions of data
[m,n]=size(x);

% Dimensions of remaining folded lowpass subband
cm=m; cn=n;

% Check dimensions of input data
for iter=1:2,
    tmp=size(x,iter)/2^floor((J+2-iter)/2);
    if round(tmp)~=tmp,
        disp(' ')
        disp('The size of the input signal must be a sufficient power of two!')
        disp(' ')
        y2=[];
        return;
    end;
end;

% Analysis filters:
%  H1,   G1   : filters for iterations with mod(j,2)=1
%  H1D,  G1D  : filters for iterations with mod(j,2)=0
% (premultiply analysis filters by the downsampling factor; i.e., det(D)=2)
H1=FA(:,:,1)/2;
G1=FA(:,:,2)/2;
H1D=FA(:,:,3)/2;
G1D=FA(:,:,4)/2;

% Compute the FFT of the input data
X=fftn(x); clear x;

% Initialize cube for folded wavelet coefficients
y2=zeros(cm,cn);

for j=1:J,

    mj=mod(j,2);

    % Select filter
    switch mj,
        case 1,
            H=H1;
            G=G1;
        case 0,
            H=H1D;
            G=G1D;
    end;

    % Filtering
    Y1=H(1:size(X,1),1:size(X,2)).*X;
    Y2=G(1:size(X,1),1:size(X,2)).*X;

    % Downsampling & upsampling
    switch mj,
        case 1,
            if j~=J,
                % This operates like: Y1=Y1+fftshift(Y1); Y1=Y1(1:cm,1:cn/2);
                Y1(1:cm/2,1:cn/2)=Y1(1:cm/2,1:cn/2)+Y1(cm/2+1:cm,cn/2+1:cn);
                Y1(cm/2+1:cm,1:cn/2)=Y1(cm/2+1:cm,1:cn/2)+Y1(1:cm/2,cn/2+1:cn);
                Y1=Y1(1:cm,1:cn/2);
            else
                Y1=Y1+fftshift(Y1);
            end;
            Y2=Y2 + fftshift(Y2);
            %! y2(1:cm,cn/2+1:cn)=fold2D(real(ifftn(Y2)),mj,m,n);
            y2(1:cm,cn/2+1:cn)=fold2D((ifftn(Y2)),mj,m,n);
            cn=cn/2;
        case 0,
            m2=m/2;n2=n/2;
            Y1=Y1(1:m2,1:n2)+Y1(m2+1:2*m2,1:n2);
            Y2=Y2(1:m2,1:n2)+Y2(m2+1:2*m2,1:n2);
            y2(cm/2+1:cm,1:cn)=real(ifftn(Y2));
            cm=cm/2;

            % Preparing filters for next three iterations
            lm=1:2:m; ln=1:2:n;
            H1=H1(lm,ln);
            G1=G1(lm,ln);
            H1D=H1D(lm,ln);
            G1D=G1D(lm,ln);
            m=m2; n=n2;
    end;

    X=Y1;
end;

% Insert folded lowpass subband
y2(1:cm,1:cn)=fold2D(real(ifftn(Y1)),mj,m,n);
return;


% ----------------------------------------------------------------------
% Fold 2D subband
% ----------------------------------------------------------------------
function f=fold2D(u,j,m,n)
switch j,
    case 0,
        f=u;
    case 1,
        f=u(:,1:2:n)+u(:,2:2:n);
end;


function x=FFTquin2D_synthesis(y2,J,FS);

% FFTQUIN2D_SYNTHESIS Perform 2D quincunx wavelet synthesis
% 	x = FFTquin2D_synthesis(y2,J,FS)
%
% 	Input:
%       y2 = data (folded)
%       J = number of interations
%       FS = synthesis filters
%
%       Output:
%       x = result (unfolded)
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland

% Dimensions of input data
[m,n]=size(y2);

% Synthesis filters:
%  H1,   G1   : filters for iterations with mod(j,2)=1
%  H1D,  G1D  : filters for iterations with mod(j,2)=0
H1=FS(:,:,1);
G1=FS(:,:,2);
H1D=FS(:,:,3);
G1D=FS(:,:,4);

% How many "double-iterations"?
l=2^((J-mod(J,2))/2);
M=m/l; N=n/l;

% Dimensions of full data
[M0,N0,P0]=size(H1);

% Dimensions of folded lowpass subband
cm=M0; cn=N0;
cn=cn/2^floor((J+1)/2);
cm=cm/2^floor((J)/2);

% Extract folded lowpass subband and unfold
y1=y2(1:cm,1:cn);
y1=unfold2D(y1,mod(J,2),M,N);
Y1=fftn(y1); clear y1;

j=J;
while j>0;
    mj=mod(j,2);

    % unfold highpass subband and select filters
    switch mj,
        case 1,
            f=unfold2D(y2(1:M,N/2+1:N),mj,M,N);
            H=H1(1:l:M0,1:l:N0);
            G=G1(1:l:M0,1:l:N0);
        case 0,
            f=unfold2D(y2(M+1:2*M,1:N),mj,M,N);
            l=l/2;
            Y1=[Y1 Y1;Y1 Y1];
            H=H1D(1:l:M0,1:l:N0);
            G=G1D(1:l:M0,1:l:N0);
    end;

    Y2=fftn(f); clear f;
    if mj==0,
        Y2=[Y2 Y2;Y2 Y2];
        M=2*M; N=2*N;
    end;

    % filter data
    Y1=Y1.*H + Y2.*G;

    j=j-1;
end;
x=real(ifftn(Y1));
return;

% ----------------------------------------------------------------------
% Unfold 2D subband
% ----------------------------------------------------------------------
function u=unfold2D(y,j,M,N)
switch j,
    case 0,
        u=y;
    case 1,
        u=zeros(M,N);
        u(1:2:M,1:2:N)=y(1:2:M,:);
        u(2:2:M,2:2:N)=y(2:2:M,:);
end;
return;


function [FA,FS]=FFTquincunxfilter2D(dim,alpha,type);

% FFTQUINCUNXFILTER2D Supplies filters for the 2D quincunx transform.
% 	[FFTanalysisfilters,FFTsynthesisfilters]=FFTquincunxfilter2D(dim,alpha,type)
% 	computes the frequency response of low- and high-pass filters.
%
% 	Input:
% 	dim = dimensions of the input signal
% 	alpha = degree of the filters, any real number >=0
%   type =
%     class I: McClellan-based filters
%       o = orthogonal; b = discrete linear B-spline (synthesis); d = dual
%     class II: Polyharmonic B-spline wavelets (Rabut-style localization)
%       Po = orthogonal; Pb = B-spline; Pd = dual
%     class III: Polyharmonic B-spline wavelets (isotropic-style
%                localization)
%       PO = orthogonal {=Po}; PB = B-spline; PD = dual
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland

m=dim(1);
n=dim(2);

gamma=alpha+2;

% Setup downsampling matrix
D=[1 1; 1 -1]';

% Coordinates in the Fourier domain
[xo,yo]=ndgrid(2*pi*([1:m]-1)/m,2*pi*([1:n]-1)/n);

for iter=0:1,

    if iter==0, % First filter: no downsampling
        x=xo;
        y=yo;
    end;
    if iter>0,  % Coordinates after downsampling
        x=D(1,1)*xo+D(1,2)*yo;
        y=D(2,1)*xo+D(2,2)*yo;
        % D=D*D;    % prepare filter for next iteration
    end;

    if lower(type(1)) == 'p', % Isotropic polyharmonic splines

        if iter==0,
            [ac,acD,loc,B] = FFTquincunxpolyfilter(x,y,gamma,type(2));
        else
            %ac = fftshift(interp2(xo,yo,ac,mod(x,2*pi),mod(y,2*pi),'*nearest'),1);
            %acD= fftshift(interp2(xo,yo,acD,mod(x,2*pi),mod(y,2*pi),'*nearest'));
            %loc= fftshift(interp2(xo,yo,loc,mod(x,2*pi),mod(y,2*pi),'*nearest'),1);
            %B  = fftshift(interp2(xo,yo,B,mod(x,2*pi),mod(y,2*pi),'*nearest'),1);
            ac = interp2(xo,yo,ac0,mod(x,2*pi),mod(y,2*pi),'*nearest');
            acD= interp2(xo,yo,acD0,mod(x,2*pi),mod(y,2*pi),'*nearest');
            loc= interp2(xo,yo,loc0,mod(x,2*pi),mod(y,2*pi),'*nearest');
            B  = interp2(xo,yo,B0,mod(x,2*pi),mod(y,2*pi),'*nearest');
            %[ac,acD,loc,B] = FFTquincunxpolyfilter(x,y,gamma,type(2));
        end;
        ortho = acD./ac;

    else

        [ortho,B] = FFTquincunxmcclellan(x,y,alpha);

    end;

    % Lowpass filter
    [H,H1] = FFTquincunxfilter2D_lowpass(x,y,type,B,ortho);

    % Reversed filter
    B0 = B;
    if iter==0,
        B=fftshift(B); ortho=fftshift(ortho);
    else
        B=fftshift(B,1); ortho=fftshift(ortho,1);
    end;

    if lower(type(1)) == 'p',
        ac0=ac; acD0=acD; loc0=loc;
        if iter==0,
            ac=fftshift(ac);
            acD=fftshift(acD);
            loc=fftshift(loc);
        else
            ac=fftshift(ac,1);
            acD=fftshift(acD);
            loc=fftshift(loc,1);
        end;
    else
        ac=0; ac0=0; acD=0; acD0=0; loc=0; loc0=0;
    end;

    [Hd,H1d] = FFTquincunxfilter2D_highpass(x-pi,y-pi,type,B,ortho,ac,ac0,acD,acD0,loc,loc0);

    % Highpass filters: modulation
    G  = exp(i*x).*H1d;
    G1 = exp(-i*x).*Hd;

    % Fill up array with analysis and synthesis filters
    FA(:,:,iter*2+1) = H1;
    FA(:,:,iter*2+2) = G1;
    FS(:,:,iter*2+1) = H;
    FS(:,:,iter*2+2) = G;
end;



function [H,H1]=FFTquincunxfilter2D_highpass(x,y,type,B,ortho,ac,ac0,acD,acD0,loc,loc0);

% FFTQUINCUNXFILTER2D_HIGHPASS Supplies highpass filters for 2D quincunx transform.
% 	[Hsynthesis,Hanalysis]=FFTquincunxfilter2D_highpass(x,y,alpha,type)
% 	computes the frequency response of highpass filters.
%
% 	Input:
% 	x,y = coordinates in the frequency domain
% 	alpha = degree of the filters, any real number >=0
%       type = type of filter (ortho, dual, bspline)
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland


if lower(type(1))=='p', % Polyharmonic splines

    type = [ type(2:length(type)) ];

    switch lower(type(1)),

        % Orthogonal filters
        case 'o',
            H = B*sqrt(2) ./ sqrt(ortho);
            H1  = conj(H);

            % Dual filters (B-spline at the analysis side)
        case 'd',
            H = sqrt(2)*conj(B).*ac;
            H1  = sqrt(2)*B./acD;

            % B-spline filters (B-spline at the synthesis side)
        case 'b',
            H1 = sqrt(2)*conj(B).*ac;
            H  = sqrt(2)*B./acD;

        case 'i',
            H = sqrt(2)*loc0./ac0;
            H1 = sqrt(2)*(B.^2).*ac./acD.*ac0./loc0;
            H1(find(isinf(H1)))=0;

    end;

else  % Simple fractional iterate with McClellan transform
    [H,H1]=FFTquincunxfilter2D_lowpass(x,y,type,B,ortho);
end;



function [H,H1]=FFTquincunxfilter2D_lowpass(x,y,type,B,ortho);

% FFTQUINCUNXFILTER2D_LOWPASS Supplies lowpass filters for 2D quincunx transform.
% 	[Hsynthesis,Hanalysis]=FFTquincunxfilter2D_lowpass(x,y,alpha,type)
% 	computes the frequency response of lowpass filters.
%
% 	Input:
% 	x,y = coordinates in the frequency domain
% 	alpha = degree of the filters, any real number >=0
%       type = type of filter (ortho, dual, bspline)
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland

H=B;

if lower(type(1)) == 'p',
    type = [ type(2:length(type)) ];
end;

switch lower(type(1)),

    % Orthogonal filters
    case 'o',
        H = H*sqrt(2) ./ sqrt(ortho);
        H1  = conj(H);

        % Dual filters (B-spline at the analysis side)
    case 'd',
        H1 = sqrt(2)*conj(H);
        H = conj(H1./(ortho));

        % B-spline filters (B-spline at the synthesis side)
    case 'b',
        H = sqrt(2)*H;
        H1 = conj(H./(ortho));

        % Isotropic wavelet filter
    case 'i',
        H1 = sqrt(2)*conj(H);
        H = conj(H1./(ortho));

end;

function [ortho,H] = FFTquincunxmcclellan(x,y,alpha)

% FFTQUINCUNXMCCLELLAN Supplies orthonormalizing part of the 2D McClellan
%       filters for 2D quincunx transform.
%
% 	v=FFTquincunmcclellan(x,y,alpha)
%
% 	Input:
% 	x,y = coordinates in the frequency domain
% 	alpha = degree of the filters, any real number >=1
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland


% McClellan transform filters
H  = (2*ones(size(x)) + (cos(x) + cos(y))).^((alpha+1)/2);
H  = H/4^((alpha+1)/2);
Hc = (2*ones(size(x)) - (cos(x) + cos(y))).^((alpha+1)/2);
Hc = Hc/4^((alpha+1)/2);

% Orthonormalizing denominator (l_2 dual)
ortho = H.*conj(H) + Hc.*conj(Hc);

function [ac,acD,loc,B] = FFTquincunxpolyfilter(x,y,gamma,type)

% FFTQUINCUNXPOLYFILTER Supplies orthonormalizing part of the 2D polyharmonic
%       spline filters for 2D quincunx transform.
%
% 	v=FFTquincunxpolyfilter(x,y,gamma,type)
%
% 	Input:
% 	x,y = coordinates in the frequency domain
% 	gamma = order of the filters, any real number >=0
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland


% Get dimensions
[m n]=size(x);
m=2*m; n=2*n;

% Construct fine grid [0,2pi[ x [0,2pi[
[xo,yo]=ndgrid(2*pi*([1:m]-1)/m,2*pi*([1:n]-1)/n);

warning off MATLAB:divideByZero;

% Refinement filter for [2 0; 0 2]
if upper(type(1)) == type(1),
    loc = 8/3*(sin(xo/2).^2+sin(yo/2).^2)+2/3*(sin((xo+yo)/2).^2+sin((xo-yo)/2).^2);
    H  = 2^(-gamma)*((8/3*(sin(xo).^2+sin(yo).^2)+2/3*(sin(xo+yo).^2+sin(xo-yo).^2))./loc).^(gamma/2);
else
    loc = sin(xo/2).^2+sin(yo/2).^2;
    H  = 2^(-gamma)*((sin(xo).^2+sin(yo).^2)./loc).^(gamma/2);
end;
H(find(isnan(H))) = 1;

% Autocorrelation function
ac = autocorr2d(H.^2);

% Construct fine grid [0,2pi] x [0,2pi]
[xo2,yo2]=ndgrid(2*pi*([1:m/2+1]-1)/(m/2),2*pi*([1:n/2+1]-1)/(n/2));

% Extend autocorrelation function
ac(m/2+1,:)=ac(1,:);
ac(:,n/2+1)=ac(:,1);

% Compute autocorrelation function on subsampled grid; D=[1 1; 1 -1]
x2=mod(xo2+yo2,2*pi);
y2=mod(xo2-yo2,2*pi);

% Compute values on grid (x,y)
%------------------------------
% Autocorrelation filter
acD= interp2(xo2,yo2,ac,mod(x+y,2*pi),mod(x-y,2*pi),'*nearest');
ac = interp2(xo2,yo2,ac,mod(x,2*pi),mod(y,2*pi),'*nearest');

% Localization operator
loc = interp2(xo,yo,loc,mod(x,2*pi),mod(y,2*pi),'*nearest');
loc = loc.^(gamma/2);

% Scaling filter for quincunx lattice
if upper(type(1)) == type(1),
    B = 2^(-gamma/2)*((8/3*(sin((x+y)/2).^2+sin((x-y)/2).^2)+2/3*(sin(x).^2+sin(y).^2))./(8/3*(sin(x/2).^2+sin(y/2).^2)+2/3*(sin((x+y)/2).^2+sin((x-y)/2).^2))).^(gamma/2);
else
    B = 2^(-gamma/2)*((sin((x+y)/2).^2+sin((x-y)/2).^2)./(sin(x/2).^2+sin(y/2).^2)).^(gamma/2);
end;
B(find(isnan(B))) = 1;

function [v,fv] = poly(alpha,type)

% POLY - compute the polyharmonic B-spline
%
% 	Input:
% 	alpha = degree of the polyharmonic spline
%       type = type of B-spline (b,d,o,B,D,O)
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland


warning off MATLAB:divideByZero;

N=16; % spatial support: [-N/4:N/4[
Z=16; % spatial zoom factor
%N=64;
%Z=2;

step=2*pi/N;
[x,y]=meshgrid(-Z*pi:step:Z*pi-step);
w1=x; w2=y;

% compute frequency response
gamma=alpha+2;

x = x/2; y = y/2;
loc = sin(x).^2+sin(y).^2;

% alternative localization
if upper(type(1)) == type(1),
    loc = ( 8/3*loc + 2/3*(sin(x+y).^2+sin(x-y).^2) ) / 4;
end;

pow=(x.^2+y.^2);

loc=loc.^(gamma/2);
pow=pow.^(gamma/2);

fv  = ( loc ./ pow );
fv(find(isnan(fv))) = 1;

% autocorrelation function
[xo,yo]=meshgrid(0:step/2:2*pi-step/2);

% refinement filter for [2 0; 0 2]
if upper(type(1)) == type(1),
    H  = 2^(-gamma)*((8/3*(sin(xo).^2+sin(yo).^2)+2/3*(sin(xo+yo).^2+sin(xo-yo).^2))./(8/3*(sin(xo/2).^2+sin(yo/2).^2)+2/3*(sin((xo+yo)/2).^2+sin((xo-yo)/2).^2))).^(gamma/2);
else
    H  = 2^(-gamma)*((sin(xo).^2+sin(yo).^2)./(sin(xo/2).^2+sin(yo/2).^2)).^(gamma/2);
end;
H(find(isnan(H))) = 1;

% autocorrelation function
ac0 = autocorr2d(H.^2);

sx = size(ac0,1); sy = size(ac0,2);
ac = zeros(size(fv));
for iterx=1:Z,
    for itery=1:Z,
        bx=(iterx-1)*sx+1;
        by=(itery-1)*sy+1;
        ac(bx:bx+sx-1,by:by+sy-1)=ac0;
    end;
end;
ac = fftshift(ac);

if lower(type(1)) == 'o',
    fv = fv ./ sqrt(ac);
end;

if lower(type(1)) == 'd',
    fv = fv ./ ac;
end;

% compute spatial version
fv = ifftshift(fv);              % shift to [0:2pi[
v = real((ifft2(fv)))*Z^2;       % compensate for zooming

% only keep the half [-N/4:N/4[ instead of [-N/2:N/2[ (to avoid aliasing)
[x,y]=meshgrid(-N/4:1/Z:N/4-1/Z);
v=fftshift(v);
v=v((N*Z)/4:(N*Z)*3/4-1,(N*Z)/4:(N*Z)*3/4-1);
v=ifftshift(v);

% optional: normalize energy (for movies)
if length(type)>1 & type(2) == 'n',
    s=sum(v(:).^2)/(Z*Z);
    v=v*sqrt(1/s);
end;

% surface plot with lighting
surfl(x,y,fftshift(v)); shading flat; view(-26,44);
xlabel('x_1');
ylabel('x_2');

function [v,fv] = polyw(alpha,type)

% POLYW - polyharmonic B-spline wavelet
%
% 	Input:
%
%       type
%         lowercase - Rabut style
%         uppercase - isotropic brand
%
% 	alpha = degree of the polyharmonic spline
%
% Dimitri Van De Ville
%      Biomedical Imaging Group, BIO-E/STI
%      Swiss Federal Institute of Technology Lausanne
%      CH-1015 Lausanne EPFL, Switzerland


gamma=alpha+2;

warning off MATLAB:divideByZero;

N=16;
Z=16;
step=2*pi/N;
[w1,w2]=meshgrid(-Z*pi:step:Z*pi-step);

% downsampled grid
w1D=(w1+w2)/2;
w2D=(w1-w2)/2;

% scaling function
w1D  = w1D/2; w2D = w2D/2;

loc  = 4*sin(w1D).^2+4*sin(w2D).^2;

if upper(type(1)) == type(1),
    loc = loc - 8/3*sin(w1D).^2.*sin(w2D).^2;
end;

pow = (2*w1D).^2+(2*w2D).^2;

loc=loc.^(gamma/2);
pow=pow.^(gamma/2);

fv1  = loc ./  pow;
fv1(find(isnan(fv1))) = 1;
w1D  = 2*w1D; w2D = 2*w2D;

% autocorrelation function
[w1o,w2o]=meshgrid(0:step/4:2*pi-step/4);

% refinement filter for [2 0; 0 2]
if upper(type(1)) == type(1),
    H  = 2^(-gamma)*((8/3*(sin(w1o).^2+sin(w2o).^2)+2/3*(sin(w1o+w2o).^2+sin(w1o-w2o).^2))./(8/3*(sin(w1o/2).^2+sin(w2o/2).^2)+2/3*(sin((w1o+w2o)/2).^2+sin((w1o-w2o)/2).^2))).^(gamma/2);
else
    H  = 2^(-gamma)*((sin(w1o).^2+sin(w2o).^2)./(sin(w1o/2).^2+sin(w2o/2).^2)).^(gamma/2);
end;
H(find(isnan(H))) = 1;

% autocorrelation function
ac0 = autocorr2d(H.^2);
ac0(size(ac0,1)+1,:)=ac0(1,:);
ac0(:,size(ac0,2)+1)=ac0(:,1);

[w1o,w2o]=meshgrid(0:step/2:2*pi);
ac = interp2(w1o,w2o,ac0,mod(w1D,2*pi),mod(w2D,2*pi),'*linear');

if lower(type(1)) == 'o' || lower(type(1)) == 'j',
    fv1 = fv1 ./ sqrt(ac);
end;
if lower(type(1)) == 'd',
    fv1 = fv1 ./ ac;
end;

% wavelet filter
w1D = -w1D-pi; w2D = -w2D-pi;

acD = interp2(w1o,w2o,ac0,mod(w1D,2*pi),mod(w2D,2*pi),'*linear');

w1D = w1D/2;  w2D = w2D/2;

% Quincunx scaling filter
if upper(type(1)) == type(1),
    t1 = 8/3*(sin(w1D+w2D).^2 + sin(w1D-w2D).^2) + 2/3*(sin(2*w1D).^2 + sin(2*w2D).^2);
    t2 = 8/3*(sin(w1D).^2 + sin(w2D).^2) + 2/3*(sin(w1D+w2D).^2 + sin(w1D-w2D).^2);
else
    t1 = sin(w1D+w2D).^2 + sin(w1D-w2D).^2;
    t2 = sin(w1D).^2 + sin(w2D).^2;
end;

% Dyadic scaling filter
if length(type)>1 & type(2)=='D',
    t1 = sin(2*w1D).^2 + sin(2*w2D).^2;
    t2 = sin(w1D).^2 + sin(w2D).^2;
end;

t=( 0.5 * t1./t2 ).^(gamma/2);
fv2 = conj(t);
fv2(find(isnan(fv2))) = 1;
w1D = 2*w1D; w2D = 2*w2D;
w1D = -w1D-pi; w2D = -w2D-pi;

if lower(type(1)) == 'i',
    fv2 = loc;
end;
if lower(type(1)) == 'j',
    fv2 = loc;
    %  fv2(find(isinf(abs(fv2)))) = 0;
end;

% shift
fv2 = fv2 .* exp(-i*w1D);

acD2 = interp2(w1o,w2o,ac0,mod(w1,2*pi),mod(w2,2*pi),'*linear');

if lower(type(1)) == 'o',
    fv3 = sqrt(acD./acD2);
end;
if lower(type(1)) == 'b',
    fv3 = acD;
end;
if lower(type(1)) == 'd',
    fv3 = 1./acD2;
end;
if lower(type(1)) == 'i',
    fv3 = 1./ac;
end;

% scale relation
fv = 2^(gamma/2-1)*fv1.*fv2.*fv3;
%fv=fv1;
fv = ifftshift(fv);

% compute spatial version
v = ifft2(fv)*Z^2;

[x1,x2]=meshgrid(-N/4:1/Z:N/4-1/Z);
v=fftshift(v);
v=v((N*Z)/4:(N*Z)*3/4-1,(N*Z)/4:(N*Z)*3/4-1);
v=ifftshift(v);

surfl(x1,x2,fftshift(real(v))); shading flat; view(-26,44);
xlabel('x_1');
ylabel('x_2');