Simulation of speckle patterns

Matlab code

% Program to generate movie of scattering pattern of increasing number of

% particles but at a roughly constant particle density whereby

% interference between scattering from individual particles generates

% speckle, whose characteristic separation is inversely proportional to

% the linear size of the envelope of the ensemble of the particles. Three

% options are possible: identical circles, variable-radius circles, or

% variable size polygons.

clear; close all

disp('Choose one of the following ensembles: Identical circles (I),');

prompt = 'variable-radius circles (V), Random polygons (P): ';

str = input(prompt,'s');

if (str == 'I') || (str == 'i')

vid = VideoWriter('identicalCircleSpeckle.mp4','MPEG-4');

elseif (str == 'V') || (str == 'v')

vid = VideoWriter('variableCircleSpeckle.mp4','MPEG-4');

else

vid = VideoWriter('polygonSpeckle.mp4','MPEG-4');

end

vid.Quality = 100;

vid.FrameRate = 1; % 1 frame/s

open(vid);

figure('units','pixels','position',[0 0 4*4096 4*4096],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);

% Plotting areas of two components of video

pos1 = [0 0 0.5 1]; % Ensemble of particles

pos2 = [0.5 0 0.5 1]; % FT of ensemble

ra = 0.04; % Radius of particles

R = ra*2.5; % Initial size of radius of ensemble

xmax = R*2^4.52;

coords = zeros(2,1024); % Initialize radii and angles of circle origins

for jj = 1:9

% Subplot of ensemble of particles

pic1 = subplot('position',pos1);

% Increase area by two to accommodate doubling of number of particles

R = R*sqrt(2);

if (jj == 1) || (jj == 2)

startval = jj;

else

startval = 2^(jj-2)+1;

end

endval = 2^(jj-1);

for ii = startval:endval % Double number of particles

if (jj == 1)

theta = 0;

delr = 0;

else

theta = 2*pi*rand;

delr = R*(1 - 1/sqrt(2))*rand;

end

% Place particles within annular ring with area equal to circle

% within it

coords(1,ii) = (R/sqrt(2)+delr)*cos(theta);

coords(2,ii) = (R/sqrt(2)+delr)*sin(theta);

end

hold off

for kk = 1:endval

hold on

if (str == 'P') || (str == 'p') % Polygons and circles

polyVal = 0.2 + rand*0.8; % Determines number of sides of polygon

polyVarSize = 0.5+rand; % Determines size of polygon

if (polyVal >= 0.3) % Not a circle

sides = round(10*polyVal);

th = 0:2*pi/sides:2*pi;

randAng = 2*pi*rand; % Determines angular orientation of polygon

x_circle = ra*polyVarSize*cos(th+randAng) + coords(1,kk);

y_circle = ra*polyVarSize*sin(th+randAng) + coords(2,kk);

plot(x_circle, y_circle);

fill(x_circle, y_circle, 'k')

else % Circle

th = 0:2*pi/50:2*pi;

x_circle = ra*polyVarSize*cos(th) + coords(1,kk);

y_circle = ra*polyVarSize*sin(th) + coords(2,kk);

plot(x_circle, y_circle);

fill(x_circle, y_circle, 'k')

end

elseif (str == 'V') || (str == 'v') % Variable-radius circles

th = 0:2*pi/50:2*pi;

VarSize = 0.5+rand; % Radius between 0.5 and 1.5 of average

x_circle = ra*VarSize*cos(th) + coords(1,kk);

y_circle = ra*VarSize*sin(th) + coords(2,kk);

plot(x_circle, y_circle);

fill(x_circle, y_circle, 'k')

else % Identical circles

th = 0:2*pi/50:2*pi;

x_circle = ra * cos(th) + coords(1,kk);

y_circle = ra * sin(th) + coords(2,kk);

plot(x_circle, y_circle);

fill(x_circle, y_circle, 'k')

end

end

% Add white circles at corners to stop saved image shrinking to size of

% present ensemble

for aa = -1:2:1

for bb = -1:2:1

x_circle = ra * cos(th) + aa*xmax;

y_circle = ra * sin(th) + bb*xmax;

plot(x_circle, y_circle,'w*','Linewidth',0.0001);

pp = fill(x_circle, y_circle, 'w');

pp.LineWidth = 0.0001;

pp.EdgeColor = [0.99 0.99 0.99];

end

end

axis square

xlim([-xmax xmax])

ylim([-xmax xmax])

axis off

filename = 'circles.png';

myAxes=findobj(pic1,'Type','Axes');

exportgraphics(myAxes,filename,'Resolution',300);

% Subplot of FT of ensemble

subplot('position',pos2);

hold off

I1 = double(I);

I1=imcomplement(rgb2gray(I1));

F = fft2(I1);

F = fftshift(F);

minInt = min(min((abs(F)).^0.4)); % Power 0.4 compresses intensity scale

maxInt = max(max((abs(F)).^0.4));

xsize = round(0.1*size(F,1)); % 10% of linear dimension of FT image

F = F(4*xsize:6*xsize,4*xsize:6*xsize); % Central 20% x 20% of FT

imshow(((abs(F)).^0.4),[minInt,maxInt/(1+jj/5)]); % Scale varies with jj to keep nice intensity range

greenBlack = customcolormap([0 1], {'#00ff00','#000000'}); % Similar to laser green speckle

colormap(greenBlack);

frame = getframe(gcf);

writeVideo(vid,frame);

end

close(vid)