%
%function alias_demo()
%
% Aliasing demo
%
function alias_demo()

    close all

    imsize = 500;

    vessel_fig = figure('name','vessel','position', ...
                        [250,1000,2*imsize,imsize],...
                        'visible','off') ;

    brain_fig = figure('name','brain','position', ...
                       [250+2*imsize,1000,2*imsize,imsize],...
                       'visible','off') ;


    % ----- User input --------- %

    fprintf('--------- User input ---------\n');

    % Vessel frequency

    vpow_def = 2;
    istr = sprintf('Vessel frequency exponent (power of 2), as function of phase encoding): [2^%g] ',vpow_def);
    vpow = input(istr);
    if isempty(vpow),
        vpow = vpow_def;
    end
    vfreq = 2^vpow;

    fprintf('Vessel frequency = %g\n',vfreq);

    n = 128;
    brain = phantom('Modified Shepp-Logan',n);

    %vessel = zeros(n);
    %vessel(n/2,n/2) = 1;

    vessel = disk2D(n,.05);

    % Vessel n/2 ghost 
    fb = fft2d(brain);
    fv = fft2d(vessel);
    fv(1:vfreq:end) = 2*fv(1:vfreq:end);

    % add brain and vessel data
    fbv = fb+fv;

    % create brain and vessel image
    ifbv = ifft2d(fbv);

    figure(vessel_fig)
    subplot(121)
    imagesc(ifbv)
    colormap gray
    axis equal tight
    xlabel('x')
    ylabel('x')
    clearticklabels
    title('Image[brain+vessel]')

    subplot(122)
    imagesc(log(abs(ifbv-brain)))
    colormap gray
    axis equal tight
    xlabel('x')
    ylabel('x')
    clearticklabels
    title('log(abs(Image[brain+vessel] - Image[brain]))')

    % Brain n/2 ghost 

    figure(brain_fig)
    subplot(121)
    imagesc(brain)
    colormap gray
    axis equal tight
    xlabel('x')
    ylabel('x')
    clearticklabels
    title('Image[brain]')

    fb = fft2d(brain);
    fb(1:vfreq:end) = 2*fb(1:vfreq:end);
    ifb = ifft2d(fb);

    subplot(122)
    imagesc(ifb)
    colormap gray
    axis equal tight
    xlabel('x')
    ylabel('x')
    clearticklabels
    title('Image[brain] - even/odd problems')

    %keyboard
    return

%-------- Functions ---------%

function disk = disk2D(ngrid,radpct)

    x = symspace(1,ngrid);
    y = x;
    [X,Y] = meshgrid(x,y);

    rad = radpct*max(abs(x));

    disk = zeros(ngrid);
    for i=1:ngrid
        for j=1:ngrid
            xij = X(i,j);
            yij = Y(i,j);
            r = sqrt(xij.^2+yij.^2);
            if r<rad
                disk(i,j) = 1;
            else
                disk(i,j) = 0;
            end
        end
    end

    return

function fx = fft2d(x)
    fx = fftshift(fft2(fftshift(x)));
    return

function ifx = ifft2d(x)
    ifx = fftshift(ifft2(fftshift(x)));
    return
    %function q = clearlabels()
    %
    % Clear tick labels marks from all axes of plots
    %

function q = clearticklabels()
    set(gca,'xticklabel',[])
    set(gca,'yticklabel',[])
    set(gca,'zticklabel',[])
    return