Code covered by the BSD License  

Highlights from
Huygens-Fresnel integral approximation, free-form apertures and rough surfaces

image thumbnail

Huygens-Fresnel integral approximation, free-form apertures and rough surfaces

by

 

Simulate wave propagation through free-form apertures, or off rough surfaces.

fig08_11_speckle_images.m
% fig08_11_speckle_images.m
%-----------------------------
% Compute objective speckle images of rough surfaces. Plot Histograms and
% the resulting images using different wavelengths, observation distances
% and roughnesses.
%
% http://dx.doi.org/10.1016/j.measurement.2013.05.003
% Section 5, Fig.  8 "Histograms"
% Section 5, Fig.  9 "Wavelength variation"
% Section 5, Fig. 10 "Distance variation"
% Section 5, Fig. 11 "Roughness variation"
%
%
% NEEDS GPUmat (http://sourceforge.net/projects/gpumat/)

% prepare workspace
clear; clc;									% clear up

addpath('gpuSpeckle');						% add function folder to PATH variable
GPUstart();									% start GPUmat

%% define optical components in the path of light

options = createStruct('options', ...
	'display', 500, ...			% display a message every x frames
	'save', 100000, ...			% save a snapshot every y frames
    'title','fig08_11_speckle tests', ...
	'description', 'objective speckles', ...	% description text
    'overwrite', false, ...     % overwrite existing data
	'useGpu', true, ...			% use GPU acceleration
	'singlePrecisionIsError', true, ...	% throw an error for older/cheaper graphics cards
	'dryRun', false ...			% discard computing, test object/folder creation
	);
% useGpu = true: 500 px in ~13 s, false: 500 px in ~5 min
% Intel Core i7-2600 @ 3.4 GHz, 8 GB RAM, Win7-x64, Matlab 2011b
% Nvidia GeForce GTX 560 Ti, 1024 MB RAM

% light source
lightsource = createStruct('source', ...
	'subtype', 'planeWave', ...
	'name', 'light source', ...
	'wavelength', 632.8e-9, ...				% wavelength
	'irradiance', 1, ...                    % power in watts/m^2
	'normalDistance', 0.2, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180 ...		% rotations of aperture with respect to the origin [alpha beta gamma]
	);

% master surface
surface = createStruct('surface', ...
	'subtype', 'rough', ...
	'name', 'rough surface', ...
    'size', [3 3] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [.75 .75] * 1e-6, ...		% source grid spacing [dx dy]
	'rotations', [0 0 0]*pi/180, ...		% rotations of plane with respect to the origin [alpha beta gamma]
	'translations', [0 0 0]', ...			% translation of the plane with respect to the origin [x y z]
	'ra', 1e-6, ...							% surface roughness					
	'filter', @(x,y,P1,P2) exp(-(max(x.^2./P2(1).^2 + y.^2./P2(2).^2 - 1,0))), ...
    'filterParameters', [0 0], ...
	'scaledFilterParameters', [.1; .1] ...
	);

% save it, to be loaded separately for each other testcase
masterSurface = surface.surface;
masterSurfaceNGrid = surface.surfaceNGrid;
save('masterSurface.mat', 'surface', 'masterSurface', 'masterSurfaceNGrid');


%% define camera chip and calculate path of light

% case 1: camera distance 100 mm, Ra = 0.1 m, lambda = 632.8 nm
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u1_distance100_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.10, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180 ...		% rotations with respect to the origin [alpha beta gamma]
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E1_roughness0u1_distance0100_lambda6328.mat', 'E');


% case 2: camera distance 150 mm, Ra = 0.1 m, lambda = 632.8 nm
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u1_distance150_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.15, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180, ...		% rotations with respect to the origin [alpha beta gamma]
	'apertureFunction',@(x,y) (x <= 0) & (y >= 0) ... % 2nd quadrant
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E2_roughness0u1_distance0150_lambda6328.mat', 'E');


% case 3: camera distance 200 mm, Ra = 0.1 m, lambda = 632.8 nm
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u1_distance200_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.20, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180 ...		% rotations with respect to the origin [alpha beta gamma]
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E3_roughness0u1_distance0200_lambda6328.mat', 'E');


% case 4: camera distance 250 mm, Ra = 0.1 m, lambda = 632.8 nm
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u1_distance250_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.25, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180, ...		% rotations with respect to the origin [alpha beta gamma]
	'apertureFunction',@(x,y) (x >= 0) & (y <= 0) ... % 4th quadrant
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E4_roughness0u1_distance0250_lambda6328.mat', 'E');


% case 5: camera distance 200 mm, Ra = 0.1 m, lambda = 640 nm
lightsource.wavelength = 640e-9;
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u1_distance200_lambda640', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.2, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180 ...		% rotations with respect to the origin [alpha beta gamma]
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E5_roughness0u1_distance0200_lambda640.mat', 'E');


% case 6: camera distance 100 mm, Ra = 0.075 m, lambda = 632.8 nm
% use master surface from above, scale it to a new Ra
surface = createStruct('surface', ...
	'subtype', 'rough', ...
	'name', 'rough surface', ...
    'size', [3 3] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [.75 .75] * 1e-6, ...		% source grid spacing [dx dy]
	'rotations', [0 0 0]*pi/180, ...		% rotations of plane with respect to the origin [alpha beta gamma]
	'translations', [0 0 0]', ...			% translation of the plane with respect to the origin [x y z]
	'ra', .075e-6, ...                      % surface roughness					
	'surface', masterSurface, ...
	'surfaceNGrid', masterSurfaceNGrid, ...
	'filter', @(x,y,P1,P2) exp(-(max(x.^2./P2(1).^2 + y.^2./P2(2).^2 - 1,0))), ...
    'functionParameters', [0 0], ...
	'scaledFunctionParameters', [.1; .1] ...
	);
lightsource.wavelength = 632.8e-9;			% wavelength to be used for this simulation
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u075_distance100_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.10, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180,...		% rotations with respect to the origin [alpha beta gamma]
	'apertureFunction',@(x,y) (x <= 0) & (y >= 0) ... % 2nd quadrant
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E6_roughness0u075_distance0100_lambda6328.mat', 'E');


% case 7: camera distance 100 mm, Ra = 0.05 m, lambda = 632.8 nm
surface = createStruct('surface', ...
	'subtype', 'rough', ...
	'name', 'rough surface', ...
    'size', [3 3] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [.75 .75] * 1e-6, ...		% source grid spacing [dx dy]
	'rotations', [0 0 0]*pi/180, ...		% rotations of plane with respect to the origin [alpha beta gamma]
	'translations', [0 0 0]', ...			% translation of the plane with respect to the origin [x y z]
	'ra', .05e-6, ...						% surface roughness					
	'surface', masterSurface, ...
	'surfaceNGrid', masterSurfaceNGrid, ...
	'filter', @(x,y,P1,P2) exp(-(max(x.^2./P2(1).^2 + y.^2./P2(2).^2 - 1,0))), ...
    'functionParameters', [0 0], ...
	'scaledFunctionParameters', [.1; .1] ...
	);
lightsource.wavelength = 632.8e-9;			% wavelength to be used for this simulation
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u05_distance100_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.10, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180, ...		% rotations with respect to the origin [alpha beta gamma]
	'apertureFunction',@(x,y) (x >= 0) & (y >= 0) ... % 1st quadrant
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E7_roughness0u05_distance0100_lambda6328.mat', 'E');


% case 8: camera distance 100 mm, Ra = 0.01 m, lambda = 632.8 nm
surface = createStruct('surface', ...
	'subtype', 'rough', ...
	'name', 'rough surface', ...
    'size', [3 3] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [.75 .75] * 1e-6, ...		% source grid spacing [dx dy]
	'rotations', [0 0 0]*pi/180, ...		% rotations of plane with respect to the origin [alpha beta gamma]
	'translations', [0 0 0]', ...			% translation of the plane with respect to the origin [x y z]
	'ra', .01e-6, ...						% surface roughness					
	'surface', masterSurface, ...
	'surfaceNGrid', masterSurfaceNGrid, ...
	'filter', @(x,y,P1,P2) exp(-(max(x.^2./P2(1).^2 + y.^2./P2(2).^2 - 1,0))), ...
    'functionParameters', [0 0], ...
	'scaledFunctionParameters', [.1; .1] ...
	);
lightsource.wavelength = 632.8e-9;			% wavelength to be used for this simulation
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'ra0u01_distance100_lambda6328', ...
	'nGrid', [512 512], ...					% number imaging grid points [m n]
	'resolution', [4.65 4.65] * 1e-6, ...	% imaging grid spacing [dx dy]
	'normalDistance', 0.10, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180, ...		% rotations with respect to the origin [alpha beta gamma]
	'apertureFunction',@(x,y) (x >= 0) & (y <= 0) ... % 4th quadrant
	);
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);
E = reshape(camera.E,camera.nGrid(1),[]);
save('E8_roughness0u01_distance0100_lambda6328.mat', 'E');


%% predefine necessary functions
intensity = @(x) x.*conj(x);
norm = @(x) x./max(x(:));

% images size
M = camera.nGrid(1);
N = camera.nGrid(2);

% image coordinate
X = reshape(camera.grid(1,:), camera.nGrid(1),[]);
Y = reshape(camera.grid(2,:), camera.nGrid(1),[]);
x = X(:,1);
y = Y(1,:);

% logical indexing of quadrants
q1 = (X >= 0) & (Y >= 0);					% upper right quarter image
q2 = (X < 0) & (Y >= 0);					% upper left quarter image
q3 = (X < 0) & (Y < 0);						% lower left quarter image
q4 = (X >= 0) & (Y < 0);					% lower right quarter image

%% Section 5, Fig. 8: (a) Normalized histogram of the amplitudes A of a
% simulated speckle pattern at z = 100 mm, superimposed with its
% theoretical Rayleigh distribution. (b) Histogram of the intensities I
% with its theoretical exponential distribution.

% load image, approx. intensity
screen = load('E1_roughness0u1_distance0100_lambda6328.mat');
E1 = screen.E;
I1 = intensity(E1);

A = 0:1/1023:1;

% amplitudes: histogram, sum over quadratic errors, search for a parameter
% set with minimum error
h1 = norm(hist(abs(E1(:)),1024));
f1 = @(x) sum((h1-(x(2)*A./x(1)^2.*exp(-(A.^2)./(2*x(1)^2)))).^2);
[x1 fval1] = fminsearch(f1, [1 1], optimset('TolX',1e-14,'TolFun',1e-14,'MaxIter',1e12));

% intensities: histogram, sum over quadratic errors, search for a parameter
% set with minimum error
h2 = norm(hist(abs(I1(:)),1024));
f2 = @(x) sum((h2-(x(2)./(2*x(1)^2).*exp(-A./(2*x(1)^2)))).^2);
[x2 fval] = fminsearch(f2,[1 1],optimset('TolX',1e-8,'TolFun',1e-8,'MaxIter',1e6));

% plot and annotate both histograms
set(0,'defaulttextinterpreter','none');
h = figure('units','centimeters', 'Position',[10 10 10 10/4*3]); % position: [llx lly width height]
ah = plot(A, h1, A, x1(2)*A./x1(1)^2.*exp(-A.^2./(2*x1(1)^2)), A, h2, A, x2(2)./(2*x2(1)^2).*exp(-A./(2*x2(1)^2)),'LineWidth', 1.5);
grid on; hold on;
set(ah(2),'Color',[0 0 0],'LineWidth',.8);
set(ah(1),'Color',[0.7 0.7 0.7],'LineWidth',1.5);
set(ah(4),'Color',[0 0 0],'LineWidth',.8);
set(ah(3),'Color',[0.7 0.7 0.7],'LineWidth',1.5);
ylim([0,1])
xlabel('$I/I_{\text{max}}$ \qquad\qquad $A/A_{\text{max}}$');
ylabel('$\propto$ p$\left(A/A_{\text{max}}\right)$ \qquad\qquad $\propto$ p$\left(I/I_{\text{max}}\right)$')
text(0.225,0.70,'a)')
text(0.225,0.12,'b)')
set(gca,'looseinset',get(gca,'tightinset'))

% mymlf2pdf(h,'08histograms', 'fontsize', 9, 'packages',{{'', 'amsmath'}})


%% Section 5, Fig. 9: Sections of speckle images for wavelengths
% lambda_1 = 632.8 nm and lambda_2 = 640 nm (upper row), the left being
% slightly scaled by lambda_2/Lambda_1 with respect to the optical axis
% through the image center. The sum and the magnitude of the difference
% of both intensity images are shown in the lower row, where the scaling
% relative to the center becomes apparent, as the differences grow with
% increasing radius.

% load images, approx. normalized intensities
screen = load('E3_roughness0u1_distance0200_lambda6328.mat');
I2 = norm(intensity(screen.E));

screen = load('E5_roughness0u1_distance0200_lambda640.mat');
I3 = norm(intensity(screen.E));

% build sum and difference
I2p3 = norm(abs(I2+I3));
I2m3 = norm(abs(I2-I3));

% compose image quarters into a single image,
% separate them by a line of width 2 * 'b'
% camera is flipped by beta=180 about the y-axis -> quadrants 1/2 and 3/4
% change place
b=1;
img = I2;
img(q1)	= I3(q1);
img(q2)	= I2m3(q2);
img(q3) = I2p3(q3);
img(M/2-b:M/2+b,:) = 1;
img(:,N/2-b:N/2+b) = 1;

% plot and annotate image
set(0,'defaulttextinterpreter','none');
h = figure('units','centimeters', 'Position',[10 10 10 10/4*3]); % position: [llx lly width height]
imagesc(img); colormap gray;
set(gca,'XTick',[], 'YTick',[]);
text(0,-25,'$\lambda=632.8$\,nm', 'FontSize',9);
text(N,-25,'$\lambda=640$\,nm', 'FontSize',9,'HorizontalAlignment','right');
text(0,M+25,'$|\Sigma(632.8$\,nm, $640$\,nm$)|$', 'FontSize',9);
text(N,M+25,'$|\Delta(632.8$\,nm, $640$\,nm$)|$', 'FontSize',9,'HorizontalAlignment','right');

mlf2pdf(h,'09wavelength', 'fontsize', 9)


%% Section 5, Fig. 10: Sections of fully developed speckle patterns,
% simulated for increasing observation distances
% z = {100, 150, 200, 250} mm in clockwise arrangement, show greater
% speckle size for greater distances.

% load images, approx. normalized intensities
screen = load('E1_roughness0u1_distance0100_lambda6328.mat');
I4 = norm(intensity(screen.E));
screen = load('E2_roughness0u1_distance0150_lambda6328.mat');
I5 = norm(intensity(screen.E));
screen = load('E3_roughness0u1_distance0200_lambda6328.mat');
I6 = norm(intensity(screen.E));
screen = load('E4_roughness0u1_distance0250_lambda6328.mat');
I7 = norm(intensity(screen.E));

% compose image quarters into a single image,
% separate them by a line of width 2 * 'b'
% camera is flipped by beta=180 about the y-axis -> quadrants 1/2 and 3/4
% change place
b=1;
img = I4;
img(q1) = I5(q1);
img(q2) = I6(q2);
img(q3) = I7(q3);
img(M/2-b:M/2+b,:)=1;
img(:,N/2-b:N/2+b)=1;

% plot and annotate image (flipud img, since x-axis should point upwards)
set(0,'defaulttextinterpreter','none');
h = figure('units','centimeters', 'Position',[10 10 10 10/4*3]); % position: [llx lly width height]
imagesc(img); colormap gray;
set(gca,'XTick',[], 'YTick',[]);
text(0,-25,'$z=100$\,mm', 'FontSize',9);
text(N,-25,'$z=150$\,mm', 'FontSize',9,'HorizontalAlignment','right');
text(0,M+25,'$z=250$\,mm', 'FontSize',9);
text(N,M+25,'$z=200$\,mm', 'FontSize',9,'HorizontalAlignment','right');

% mymlf2pdf(h,'10distance', 'fontsize', 9)


%% Section 5, Fig. 11: The top left section shows a fully developed
% speckle pattern of a scattering surface. As the surface roughness
% Ra = {100, 75, 50, 10} nm decreases, shown in adjacent sections in
% clockwise order, the speckle pattern is less developed and converges to
% the diffraction pattern of the scattering aperture.

% load images, approx. normalized intensities
screen = load('E1_roughness0u1_distance0100_lambda6328.mat');
I8 = norm(intensity(screen.E));
screen = load('E6_roughness0u075_distance0100_lambda6328.mat');
I9 = norm(intensity(screen.E));
screen = load('E7_roughness0u05_distance0100_lambda6328.mat');
I10 = norm(intensity(screen.E));
screen = load('E8_roughness0u01_distance0100_lambda6328.mat');
I11 = norm(intensity(screen.E));

% compose image quarters into a single image,
% separate them by a line of width 2 * 'b'
% camera is flipped by beta=180 about the y-axis -> quadrants 1/2 and 3/4
% change place
b=1;
img = I8;
img(q1) = I9(q1);
img(q2) = I10(q2);
img(q3) = I11(q3);
img(M/2-b:M/2+b,:) = 1;
img(:,N/2-b:N/2+b) = 1;

% plot and annotate image (flipud img, since x-axis should point upwards)
set(0,'defaulttextinterpreter','none');
h=figure('units','centimeters', 'Position',[10 10 10 10/4*3]); % position: [llx lly width height]
imagesc(img); colormap gray;
set(gca,'XTick',[], 'YTick',[]);
text(0,-25,'$R_a=100$\,nm', 'FontSize',9);
text(N,-25,'$R_a=75$\,nm', 'FontSize',9,'HorizontalAlignment','right');
text(0,M+25,'$R_a=10$\,nm', 'FontSize',9);
text(N,M+25,'$R_a=50$\,nm', 'FontSize',9,'HorizontalAlignment','right');

% mymlf2pdf(h,'11roughness', 'fontsize', 9)

Contact us