% CLOUD9 - Cloud movie of 1/f noise.%% Function to create a movie of noise images having 1/f amplitude spectum% properties.%% Usage: CloudMovie = cloud9(size, factor, nturns, velfactor, nframes)%%        size   - [rows cols] size of image to produce%        factor - controls spectrum = 1/(f^factor)%        nturns - No of 2pi cycles phase can change over the whole sequence%        lowvel  - phase velocity at 0 frequency %        velfactor - phase velocity = freq^velfactor%        nframes - no of frames in movie%%        factor = 0   - raw Gaussian noise image%               = 1   - gives the 1/f `standard' drop-off for `natural' images%               = 1.5 - seems to give the most intersting `cloud patterns'%               = 2 or greater - produces `blobby' images%% Favourite parameters:%                     m = cloud9([480 640], 1.5, 4, 1, .1, 100);%% Copyright (c) 2000 Peter Kovesi% School of Computer Science & Software Engineering% The University of Western Australia% http://www.csse.uwa.edu.au/% % Permission is hereby granted, free of charge, to any person obtaining a copy% of this software and associated documentation files (the "Software"), to deal% in the Software without restriction, subject to the following conditions:% % The above copyright notice and this permission notice shall be included in % all copies or substantial portions of the Software.%% The Software is provided "as is", without warranty of any kind.%  April 2000function CloudMovie = cloud9(sze, factor, nturns, lowvel, velfactor, nframes)rows = sze(1);cols = sze(2);phase = i*random('Uniform',0,2*pi,rows,cols);  % Random uniform distribution 0 - 2pi% Create two matrices, x and y. All elements of x have a value equal to its % x coordinate relative to the centre, elements of y have values equal to % their y coordinate relative to the centre.  From these two matrices produce% a radius matrix that gives distances from the middlex = ones(rows,1) * (-cols/2 : (cols/2 - 1)); % x = x/(cols/2);y = (-rows/2 : (rows/2 - 1))' * ones(1,cols);% y = y/(rows/2);radius = sqrt(x.^2 + y.^2);         % Matrix values contain radius from centre.radius(rows/2+1,cols/2+1) = 1;      % .. avoid division by zero.amp = 1./(radius.^factor);          % Construct the amplitude spectrumamp = fftshift(amp);phasemod = round(fftshift(radius.^velfactor + lowvel));phasechange = 2*pi*((random('unid',nturns+1,rows,cols) -1 - nturns/2) .* phasemod );maxturns = max(max(phasechange/(2*pi)))maxturns = min(min(phasechange/(2*pi)))minturns = min(min(abs(phasechange)/(2*pi)))dphase = i*phasechange/(nframes-1);   % premultiply by i to save time in th eloop% Construct fft of noise image with the specified amplitude spectrumfig = figure(1), warning off, imagesc(zeros(rows,cols)), axis('off'), truesize(1)set(fig,'DoubleBuffer','on');%set(gca,'xlim',[-80 80],'ylim',[-80 80],...%    	   'NextPlot','replace','Visible','off')mov = avifile('cloud')a = 0.7;  % Set up colormapmap = a*bone + (1-a)*gray;for n = 1:nframes  fprintf('frame %d/%d \r',n, nframes);  phase = phase + dphase;  newfft =  amp .* exp(phase);  im = real(ifft2(newfft));            % Invert to obtain final noise image  imagesc(im), colormap(bone),axis('equal'), axis('off'), truesize(1)   %if n==1    % CloudMovie = moviein(nframes);    % initialise movie storage    %end   F = getframe(gca);   mov = addframe(mov,F);  %CloudMovie(:,n) = getframe;endfprintf('\n');warning on mov = close(mov);%movie(CloudMovie,5,30);%save('CloudMovie','CloudMovie');