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.

createSurface(varargin)
function obj = createSurface(varargin)
% createPlane - create a structure which describes a plane in space or a
% rough surface. The rough surface is computed with respect to the object's
% coordinate system.

%% parsing inputs

p = inputParser;
p.addParamValue('type','surface', @ischar);
p.addParamValue('subType', 'plane', @(x) strcmp(x,'plane') || strcmp(x,'function') || strcmp(x,'rough'));
p.addParamValue('name', '', @ischar);
p.addParamValue('fileName', '', @ischar);
p.addParamValue('folder', 'objects', @ischar);

% CUDE cufftplan ifft2() fails for certain matrix sizes, CPU FFTW version
% does not (see file generateSurfaceRough, line 74)
p.addParamValue('useGpu', false, @isboolean);
% split the source plane in N parts to speed up the calculation
p.addParamValue('stride', [], @(x)isnumeric(x));
% add uniformly distributed x/y jitter with boundaries +/-[wx, wy]
p.addParamValue('xyJitter', [], @(x)isnumeric(x) && numel(x) == 2);
% filename to save surface data only: X, Y, Z, ra, surfaceNGrid
p.addParamValue('surfaceFileName', '', @ischar);
% overall size of the grid, [x y] im meters
p.addParamValue('size', [], @(x)isnumeric(x) && numel(x) == 2);
% resolution [dx dy]
p.addParamValue('resolution', [], @(x)isnumeric(x) && numel(x) == 2);
% number of points per dimension of the grid, [x y] 
p.addParamValue('nGrid', [], @(x)isnumeric(x) && numel(x) == 2);
% grid locations [x1 x2 ...; y1 y2 ...; z1 z2 ...] in meters
p.addParamValue('grid', [], @(x) isnumeric(x) && size(x,1) == 3);
% x/y grid in reference position without translation/rotation
p.addParamValue('referenceGrid', [], @(x) isnumeric(x) && size(x,1) == 3);
% aperture criterion, a function of (x,y) in reference position
% should return logicals, i.e. a matrix of true/false
p.addParamValue('apertureFunction', [], @(x)isa(x,'function_handle'));
% indices of points in a rectangular grid within the aperture
p.addParamValue('aperture', [], @isnumeric);
% electric field on grid points
p.addParamValue('E', [], @(x) isnumeric(x) && size(x,1) == 1);
% the wavelength of light illuminating the surface
p.addParamValue('wavelength', [], @isnumeric);
% the refractive index of the material, will be used if this surface is the
% input E-field
p.addParamValue('refractiveIndex',1,@isnumeric);
% number of points per dimension of the surface grid, [x y] 
p.addParamValue('surfaceNGrid', [], @(x)isnumeric(x) && numel(x) == 2);
% relative surface matrix [x1 x2 ...; y1 y2 ...; z1 z2 ...] in meters, no
% translations/rotations
p.addParamValue('surface', [], @(x) isnumeric(x) && size(x,1) == 3);
% a function describing a surface
p.addParamValue('function',@(x,y,P1,P2) 0,@(x)isa(x,'function_handle'));
% parameters P1 of the surface function
p.addParamValue('functionParameters', [], @isnumeric);
% parameters P2 of the surface function to be scaled with the image size
p.addParamValue('scaledFunctionParameters', [.1; .1], @(x) isnumeric(x) && size(x,1) == 2);
% surface roughness Ra, arithmetic average of absolute values
p.addParamValue('ra', 0, @isnumeric);
% low-pass filter in the frequency domain to smoothen the surface
% p.addParamValue('filter', @(x,y,P1,P2) P1(1).*exp(-((x-P1(2)).^2./(2.*P2(1).^2)+(y-P1(3)).^2./(2.*P2(2).^2))), @(x)isa(x,'function_handle'));
p.addParamValue('filter', @(x,y,P1,P2)P1(1).*exp(-max(x.^2./(P2(1)).^2+y.^2./(P2(2)).^2-1, zeros(size(a)))), @(x)isa(x,'function_handle'));
% parameters P1 of the filter function (here [amplitude])
p.addParamValue('filterParameters', 1, @isnumeric);
% parameters P2 of the filter function to be scaled with the image size;
p.addParamValue('scaledFilterParameters', [.2; .1], @(x) isnumeric(x) && size(x,1) == 2);
% the shortest distance from the global origin O to the plane
p.addParamValue('normalDistance', [], @isnumeric);
% the surface normal
p.addParamValue('normal', [], @(x)isnumeric(x) && numel(x)==3);
% rotated coordinate system, columnwise unit vectors: [xx yx zx; xy yy zy; xz yz zz]
p.addParamValue('coordinateSystem', [], @(x)isnumeric(x) && size(x,1)==3 && size(x,2) == 3);
% surface normal in each grid point, calculated by its neighbours
p.addParamValue('gridNormals', [], @(x)isnumeric(x));
% vector of ideal reflection of a lightsource in every grid point
p.addParamValue('mainReflectionVectors', [], @(x)isnumeric(x));
% rotations of the base plane around its origin, order x -> y -> z
p.addParamValue('rotations', [], @(x)isnumeric(x) && numel(x)==3);
% translation of the rotated plane along the surface normal 
p.addParamValue('translations', [], @(x)isnumeric(x) && numel(x)==3);
% camera saturation function, applied to the intensity I = E.*conj(E)
p.addParamValue('saturation', [], @(x)isa(x, 'function_handle'));
% strip object of unused data, e.g. surface
p.addParamValue('strip', false, @isboolean);


%% parse the results, calculate missing parameters, return object
p.parse(varargin{:});
obj=p.Results;

displayMessage(1, ['START create ' obj.subType ' ' obj.type]);

% check basic inputs like the surface normal, the normal distance, 
% translations and rotations
obj = checkBasics(obj);

% compute size, number of grid points and resolution
obj = checkGrid(obj);

% compute aperture function
if ~isempty(obj.apertureFunction)
	displayMessage(2, [obj.name ': aperture function']);
	[X, Y] = generateMesh(obj);
	obj.aperture = obj.apertureFunction(X,Y);
end

% compute grid locations (in the x/y plane -> z-displacement 0)
if isempty(obj.grid)
	displayMessage(2, [obj.name ': grid locations']);
	[obj.grid obj.referenceGrid]= generateGridLocations(obj, 0);
else
	displayMessage(2, [obj.name ': skip grid locations']);
end

% save memory
if isempty(obj.xyJitter)
    obj.referenceGrid = [];
end

% compute surface roughness
if strcmp(obj.subType, 'rough')
	if isempty(obj.ra)
		displayMessage(-1, [obj.name ': ERROR no surface roughness defined'])
	end

	if isempty(obj.surface)                     % does it exist?
		displayMessage(2, [obj.name ': compute rough surface']);
		[obj.surface obj.surfaceNGrid] = generateSurfaceRough(obj);
	else										% is it properly scaled
		displayMessage(2, [obj.name ': skip rough surface generation, rescale existing']);
		if isempty(obj.surfaceNGrid)
			displayMessage(-1, [obj.name ': ERROR missing parameter ''surfaceNGrid''']);
		end
		obj.surface = scaleSurface(obj.surface, obj.surfaceNGrid, obj.ra);
	end

	displayMessage(2, [obj.name ': resample and add surface']);
	obj.grid = addSurface(obj);
	
elseif strcmp(obj.subType, 'function')
	if isempty(obj.surface)
		displayMessage(2, [obj.name ': compute surface function']);
		[obj.surface obj.surfaceNGrid] = generateSurfaceFunction(obj);
	else
		displayMessage(2, [obj.name ': skip surface function']);
		if isempty(obj.surfaceNGrid)
			displayMessage(-1, [obj.name ': ERROR missing parameter ''surfaceNGrid''']);
		end	
	end
	
	displayMessage(2, [obj.name ': resample and add surface']);
	obj.grid = addSurface(obj);
end

% save surface data if a filename is given
if ~strcmp(obj.surfaceFileName, '')
	X = obj.surface(1,:);
	Y = obj.surface(2,:);
	Z = obj.surface(3,:);
	ra = obj.ra;
	gridSize = obj.surfaceNGrid;
	name = obj.name;
	save(obj.surfaceFileName, 'X', 'Y', 'Z', 'ra', 'gridSize', 'name');
	clear X Y Z ra gridSize name;
end

% clear surface data to save memory
if obj.strip==true
	obj.surface = [];
end

% compute grid normals
if isempty(obj.gridNormals)
	displayMessage(2, [obj.name ': get grid normals']);
	% recover meshgrids
	X = reshape(obj.grid(1,:), obj.nGrid(1),[]);
	Y = reshape(obj.grid(2,:), obj.nGrid(1),[]);
	Z = reshape(obj.grid(3,:), obj.nGrid(1),[]);
	obj.gridNormals = getGridNormals(X, Y, Z);
else
	displayMessage(2, [obj.name ': skip grid normals']);
end

displayMessage(1, ['DONE create ' obj.subType ' ' obj.type]);

Contact us