% function [cloud,im1,im2] = imcloud(sz,im)
%
% This function creates a stimulus looking like a cloud. Optionally,
% another black&white image can be added on top to create the
% Anderson-Winaver illusion, where only the white part of the image is
% used (if a grayscale image is provided, the 90th percentile intensity is
% picked to binarize the image [most of the time not pretty] and will be
% wrong when background is white!]
%
% Input:
%       sz: width and height of image (default 200, if an image
%           is provided, the size will be adapted to the input image size)
%       im (optional): 2D image loaded into Matlab or path to image as string
%
% Output:
%       cloud: original cloud
%       im1: cloud with image on top (object has same luminance as im2!)
%       im2: dark version with image on top
%
% Examples:
%
%   cloud = imcloud(250);
%   figure, imshow(cloud)
%
%   im = imread('circles.png');
%   [cloud,im1,im2] = imcloud([],im);
%   figure, subplot(1,2,1), imshow(im1), subplot(1,2,2), imshow(im2)
%
% by Martin Hebart (2009)

function [cloud,im1,im2] = imcloud(sz,im)

if nargin > 1 % if image is provided
    
    if ischar(im)
        im = imread(im);
    end
    
    if isempty(sz)
        sz = max(size(im));
    elseif sz < max(size(im))
        error('Input sz must be equal or larger than the maximum size of the image (here %.0f). Set sz = [] to scale it automatically to the image size!',max(size(im)))
    end
    
    im = mean(im,3); % merge color channels
    
end


[nx,ny] = meshgrid(linspace(-1,1,sz));
d_sq = nx.^2 + ny.^2 + 1e-10; % creates a circle around the center (without square-root)

af = 1./d_sq .* exp(2j * pi * rand(sz)); % function from Anderson & Winaver, with random phase

ax=ifft2(af); % inverse fourier

aa= real(ax .* conj(ax)); % conj (i.e. real - i*imag) is important for proper reconstruction (additional real for rounding errors)
cloud = (aa-min(aa(:))) / (max(aa(:))-min(aa(:))); % rescale image


if nargin > 1 % if an image is provided
    
    % resize image
    cloud = cloud(1:size(im,1),1:size(im,2));
    
    if length(unique(im(:)))>2 % if image is grayscale
        tmp = sort(im(:));
        mask = im>=tmp(round(0.9*length(tmp))); % 90th percentile
    else % if image is binary
        mask = im==max(im(:));
    end   
    
    im1 = zeros(size(mask));
    im2 = zeros(size(mask));
    
    im1(~mask) = cloud(~mask)+0.3;
    im1(mask) = cloud(mask);
    im1 = min(im1,1);
    
    im2(~mask) = 0.7*cloud(~mask);
    im2(mask) = cloud(mask);
    im2 = min(im2,1);
    
end