% cid.m
% dxjones (at) gmail.com
% Copyright (c) 2008 David Jones.  All Rights Reserved.

function cid
    
    figure(1);
    fig = 1;
    pic = [];
    
    intensity = 0.90;
    for complexity = linspace(0.05,0.95,4)
        row = [];
        for disorder = linspace(0.05,0.95,4)

            rand('twister',7);
            img = instruction(complexity,intensity,disorder);
            
            subplot(4,4,fig), fig = fig + 1;
            imagesc(img), colormap(gray(1024)), axis image, axis off
            title(sprintf('C = %.2f, D = %.2f', complexity, disorder));
            drawnow
            row = [row img];
        end
        pic = [pic; row];
    end
    figure(2);
    imshow(pic);
end

%%
function pic = instruction(complexity,intensity,disorder)
    N = [270 480];
    
    smin = 1/32;
    smax = 1/3;
    sigma = (smin + (1-disorder)^2 * (smax-smin))* sqrt(N(1)*N(2));
    
    lambda = sqrt(N(1)*N(2)) / 8;
    [img1a img1b] = blendnoise(N,1,lambda,complexity,intensity,sigma);
    [img2a img2b] = blendnoise(N,2,lambda/1.5,complexity,intensity,sigma);
    
    pic1 = zeros(N(1),N(2),3);
    pic1(:,:,1) = img1a;
    pic1(:,:,2) = img1b;
    pic1(:,:,3) = 0.5 * img1a + 0.5 * img1b;
    
    pic2 = zeros(N(1),N(2),3);
    pic2(:,:,1) = img2a;
    pic2(:,:,2) = img2a;
    pic2(:,:,3) = img2b;
    
    [mask1 mask2] = mask(N(1),N(2),sigma);
    pic = zeros(N(1),N(2),3);
    for i = 1:3
        pic(:,:,i) = mask1 .* pic1(:,:,i) + mask2 .* pic2(:,:,i);
    end
    pic = max(0,min(1,pic));
end

%%
function [a b] = blendnoise(N,dir,lambda,complexity,intensity,sigma)
    x = filtered_noise(N,dir,lambda,complexity,sigma);
    a = filtered_noise(N,dir,lambda,complexity,sigma);
    b = filtered_noise(N,dir,lambda,complexity,sigma);
    a = (1-intensity) * x + intensity * a;
    b = (1-intensity) * x + intensity * b;
    a = 0.5 + 0.5 * a;
    b = 0.5 + 0.5 * b;
end

%%
function img = filtered_noise(N,dir,lambda,complexity,sigma)
    sigma = 0.75 * sigma;
    g = gausskernel(sigma);    
    sigma = 1 + (1-complexity).^3 * 1.5 * lambda;
    f = gaborkernel(lambda,sigma);
    if dir == 1
        img = wrapconv(wrapconv(rand(N)-0.5,g'),f);
    else
        img = wrapconv(wrapconv(rand(N)-0.5,g),f');
    end
    img = (0.99/prctile(abs(img(:)),99)) * img;
end

%%
function img = wrapconv(img,kernel)
    [NR NC] = size(img);
    if size(kernel,1) > size(kernel,2)
        img = conv2([img;img;img],kernel,'same');
        img = img(NR+(1:NR),:);
    else
        img = conv2([img,img,img],kernel,'same');
        img = img(:,NC+(1:NC));
    end
end

%%
function [mask1 mask2] = mask(NR,NC,sigma)
    blur = gausskernel(sigma);
    mask = wrapconv(wrapconv(rand(NR,NC) - 0.5,blur),blur');
    cutoff = prctile(abs(mask(:)),75);
    mask = (1/cutoff) * min(cutoff, max(-cutoff, mask));
    mask1 = 0.5 + 0.5 * mask;
    mask2 = 1 - mask1;
    if sum(mask1(:)) < sum(mask2(:))
        mask1 = mask2;
        mask2 = 1 - mask1;
    end
end

%%
function v = normalize(v)
    v = (1/sum(v)) * v;
end

%%
function f = gaborkernel(lambda,sigma)
	f = normalize(gabor(0,1,1/lambda,sigma,-4*sigma : 4*sigma));
end

%%
function g = gausskernel(sigma)
    g = normalize(gauss(0,sigma,-4*sigma:4*sigma));
end
    
%% gabor.m
function [ filter ] = gabor(order,amp,freq,sigma,Xdata)
    if mod(order,2) == 0
        filter = amp * cos(2*pi*freq*Xdata) .* gauss(0,sigma,Xdata);
    else
        filter = amp * sin(2*pi*freq*Xdata) .* gauss(0,sigma,Xdata);
    end
end

%% gauss.m
function [ g ] = gauss(order,sigma,Xdata)
    Xdata = (1/sigma) * Xdata;
    g=(1/(sigma*sqrt(2*pi))) * exp(-(Xdata.*Xdata)/2);
    if order == 0
        return;
    end
    g=((-1/sigma)^order) * g;

    if order == 1
        g = g.*Xdata;
    elseif order == 2
        g = g.*(Xdata.*Xdata - 1);
    elseif order == 3
        g = g.*(Xdata.*(Xdata.*Xdata - 3));
    elseif order == 4
        g = g.*(Xdata.*Xdata.*(Xdata.*Xdata-6)+3);
    elseif order == 5
        g = g.*(Xdata.*(Xdata.*Xdata.*(Xdata.*Xdata-10)+15));
    elseif order == 6
        g = g.*(Xdata.*Xdata.*(Xdata.*Xdata.*(Xdata.*Xdata-15)+45)-15);
    elseif order == 7
        g = g.*(Xdata.*(Xdata.*Xdata.*(Xdata.*Xdata.*(Xdata.*Xdata-21)+105)-105));
    else
        fprintf('Gauss order %d > 7 is too high\n', order);
    end
end