Commit 6a3ce55e authored by Simon Maretzke's avatar Simon Maretzke
Browse files

Clean up: consistently use spaces instead of tabs (Matlab standard), remove...

Clean up: consistently use spaces instead of tabs (Matlab standard), remove author and last-modified-notes, formatting fixes
parent a0879d12
Pipeline #100437 passed with stage
in 1 minute and 38 seconds
......@@ -59,29 +59,26 @@ function varargout = centeredGrid(N, dx, computeMeshgrid)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Last modified on March 19 2019 by Simon Maretzke
% Default: unit spacing
if nargin < 2
dx = 1;
dx = 1;
end
dx = dx(:).' .* ones([1,numel(N)]);
ndim = numel(N);
% Grid in Fourier space
for jj = 1:ndim
x{jj} = ( -(0.5*(N(jj)-1)) : (0.5*(N(jj)-1)) ) * dx(jj);
x{jj} = ( -(0.5*(N(jj)-1)) : (0.5*(N(jj)-1)) ) * dx(jj);
end
% Optionally assemble ndgrid (occupies more memory!)
if nargin == 3 && computeMeshgrid
[x{:}] = ndgrid(x{:});
[x{:}] = ndgrid(x{:});
elseif ndim > 1
for jj = 1:ndim
x{jj} = reshape( x{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] );
end
for jj = 1:ndim
x{jj} = reshape( x{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] );
end
end
% Unpack cell-array
......
......@@ -68,23 +68,17 @@ function combinedImage = combineInFourierspace(highFrequencyImage,lowFrequencyIm
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Original author: Martin Krenkel
% Last modified on March 22 2019 by Mareike Toepperwien
% default settings
defaults.combineCutoff = 0.1;
defaults.combineSigma = 0.01;
if (nargin == 0)
combinedImage = defaults;
return
end
if (nargin < 3)
settings = struct;
end
settings = completeStruct(settings, defaults);
......@@ -106,6 +100,5 @@ fftLowFrec = fft2(lowFrequencyImage);
% result
combinedImage = ifft2((1-weight).*fftHighFrec+weight.*fftLowFrec);
end
......@@ -48,20 +48,17 @@ function sNew = completeStruct(s, sRef)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Last modified on March 16 2019 by Simon Maretzke
sNew = s;
fnames = fieldnames(sRef);
for jj = 1:length(fnames);
fname = fnames{jj};
x = getfield(sRef, fname);
if ~isfield(sNew, fname)
sNew = setfield(sNew, fname, x);
else if isstruct(x) % recursively apply the completion to sub-structs
sNew = setfield(sNew, fname, completeStruct(getfield(sNew, fname), x));
end
fname = fnames{jj};
x = getfield(sRef, fname);
if ~isfield(sNew, fname)
sNew = setfield(sNew, fname, x);
else if isstruct(x) % recursively apply the completion to sub-structs
sNew = setfield(sNew, fname, completeStruct(getfield(sNew, fname), x));
end
end
end
......@@ -50,13 +50,10 @@ function arrayCropped = croparray(array, cropPre, cropPost)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Last modified on March 12 2019 by Simon Maretzke
% Default: symmetric cropping
if nargin < 3
cropPost = cropPre;
cropPost = cropPre;
end
N = size(array);
......@@ -70,11 +67,11 @@ cropPost = [cropPost(:).', zeros([1,ndim-numel(cropPost)])];
% Crop array dimension by dimension
arrayCropped = array;
for dim = 1:ndim
if cropPre(dim) > 0 || cropPost(dim) > 0
c1 = repmat({':'}, [dim-1,1]);
c2 = repmat({':'}, [ndim-dim,1]);
arrayCropped = arrayCropped(c1{:}, cropPre(dim)+1:N(dim)-cropPost(dim), c2{:});
end
if cropPre(dim) > 0 || cropPost(dim) > 0
c1 = repmat({':'}, [dim-1,1]);
c2 = repmat({':'}, [ndim-dim,1]);
arrayCropped = arrayCropped(c1{:}, cropPre(dim)+1:N(dim)-cropPost(dim), c2{:});
end
end
end
......@@ -55,29 +55,26 @@ function varargout = fftfreq(N, dx, computeMeshgrid)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Last modified on March 19 2019 by Simon Maretzke
% Default: unit spacing
if nargin < 2
dx = 1;
dx = 1;
end
dx = dx(:).' .* ones([1,numel(N)]);
ndim = numel(N);
% Grid in Fourier space
for jj = 1:ndim
xi{jj} = ifftshift( ( -floor(0.5*N(jj)) : floor(0.5*(N(jj)-1)) ).' * ( 2*pi / (N(jj)*dx(jj)) ) );
xi{jj} = ifftshift( ( -floor(0.5*N(jj)) : floor(0.5*(N(jj)-1)) ).' * ( 2*pi / (N(jj)*dx(jj)) ) );
end
% Optionally assemble ndgrid (occupies more memory!)
if nargin == 3 && computeMeshgrid
[xi{:}] = ndgrid(xi{:});
[xi{:}] = ndgrid(xi{:});
elseif ndim > 1
for jj = 1:ndim
xi{jj} = reshape( xi{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] );
end
for jj = 1:ndim
xi{jj} = reshape( xi{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] );
end
end
% Unpack cell-array
......
......@@ -40,8 +40,6 @@ function XiSq = fftfreqNormSq(N, dx)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Last modified on May 15 2019 by Simon Maretzke
% Default: Standard grid spacing assigned by function fftfreq.
......
......@@ -48,8 +48,6 @@ function arraySlice = getSlice(array, sliceIdx, dim)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Last modified on March 13 2019 by Simon Maretzke
ndim = ndims(array);
c1 = repmat({':'}, [dim-1,1]);
......
......@@ -24,10 +24,6 @@ function gpuIsAvailable = checkGPUSupport()
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Author: Simon Maretzke
% Last checked and modified: April 24 2019
try
gpuDevice();
......@@ -35,4 +31,5 @@ try
catch ME
gpuIsAvailable = false;
end
end
......@@ -31,12 +31,9 @@ function array = gatherIf(array, condition)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Author: Simon Maretzke
% Last checked and modified: April 24 2019
if condition
array = gather(array);
end
end
......@@ -29,10 +29,7 @@ function arrayOnHost = gatherIfOnGPU(array)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Author: Simon Maretzke
% Last checked and modified: April 24 2019
arrayOnHost = gatherIf(array, isOnGPU(array));
end
......@@ -32,13 +32,10 @@ function arrayOnGPUIf = gpuArrayIf(array, condition)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Author: Simon Maretzke
% Last checked and modified: April 24 2019
arrayOnGPUIf = array;
if condition
arrayOnGPUIf = gpuArray(arrayOnGPUIf);
end
end
......@@ -29,10 +29,7 @@ function arrayIsOnGPU = isOnGPU(array)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Author: Simon Maretzke
% Last checked and modified: April 24 2019
arrayIsOnGPU = strcmp(class(array), 'gpuArray');
end
......@@ -21,7 +21,7 @@ function [imPolar, R, Phi] = imageToPolarCoords(im, settings)
% If numAngles == [], a suitable number is determined based on the size of the input image.
% interpMethod : interpMethod = 'linear'
% The interpolation method used in the image-transformation. The admissible choices are
% the same as for the MATLAB-functions interp1, interp2, etc.
% the same as for the MATLAB-functions interp1, interp2, etc.
%
% Returns
% -------
......@@ -64,9 +64,6 @@ function [imPolar, R, Phi] = imageToPolarCoords(im, settings)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
%
% Last modified on March 13 2019 by Simon Maretzke
% Default settings
defaults.numAngles = [];
......@@ -74,12 +71,12 @@ defaults.numRadii = [];
defaults.interpMethod = 'linear';
if (nargin == 0)
imPolar = defaults;
return
imPolar = defaults;
return
end
if (nargin < 2)
settings = struct;
settings = struct;
end
settings = completeStruct(settings, defaults);
......@@ -88,10 +85,10 @@ N = size(im);
% If radial and/or angular resolutions are not assigned, choose automatically
if isempty(settings.numRadii)
numRadii = 0.5*min(N);
numRadii = 0.5*min(N);
end
if isempty(settings.numAngles)
numAngles = 2*numRadii;
numAngles = 2*numRadii;
end
% Construct polar-equispaced coordinate grid
......
......@@ -30,14 +30,14 @@ function [] = initOctave()
% along with this program. If not, see <http://www.gnu.org/licenses/>.
if isOctave
% Load additional (dummy-)functions that are needed but not natively implemented in OCTAVE
addpath(fullfile(fileparts(mfilename('fullpath')), '..', '..', 'octave'));
% Load additional (dummy-)functions that are needed but not natively implemented in OCTAVE
addpath(fullfile(fileparts(mfilename('fullpath')), '..', '..', 'octave'));
% Load image-processing toolbox
pkg load image;
% Load image-processing toolbox
pkg load image;
else
error(['This function performs the required initialization to run the toolbox in OCTAVE ', ...
'and is not meant to be called from within MATLAB.']);
error(['This function performs the required initialization to run the toolbox in OCTAVE ', ...
'and is not meant to be called from within MATLAB.']);
end
end
......@@ -28,6 +28,6 @@ function calledFromOctave = isOctave()
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
calledFromOctave = exist('OCTAVE_VERSION');
calledFromOctave = exist('OCTAVE_VERSION');
end
......@@ -49,9 +49,6 @@ function window = kaiserBesselWindow(N, beta)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Original author: Klaus Giewekemeyer
% Last modified on March 11 2019 by Simon Maretzke
if nargin < 2
beta = 8;
......
......@@ -23,7 +23,7 @@ function N_FFTfriendly = makeSizeFFTfriendly(N)
%
% .. code-block:: matlab
%
% N = [2053, 1031]; % these sizes are primes
% N = [2053, 1031]; % these sizes are primes
% N_FFTfriendly = makeSizeFFTfriendly(N);
% disp(N_FFTfriendly);
%
......@@ -54,11 +54,11 @@ function N_FFTfriendly = makeSizeFFTfriendly(N)
N_FFTfriendly = N;
for jj = 1:numel(N)
n = N(jj);
while(max(factor(n)) > 7)
n = n + 1;
end
N_FFTfriendly(jj) = n;
n = N(jj);
while(max(factor(n)) > 7)
n = n + 1;
end
N_FFTfriendly(jj) = n;
end
end
......@@ -45,11 +45,9 @@ function scaledImage = scale01(image)
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Original authors: Martin Krenkel and Mareike Toepperwien
% Last modified on March 26 2019 by Mareike Toepperwien
scaledImage = image-min(image(:));
scaledImage = scaledImage/max(scaledImage(:));
scaledImage = image-min(image(:));
scaledImage = scaledImage/max(scaledImage(:));
end
......@@ -92,205 +92,207 @@ classdef FresnelPropagator
% along with this program. If not, see <http://www.gnu.org/licenses/>.
properties (SetAccess = 'private')
settings;
sizeIn;
fresnelNumbers
numDistances;
numDims;
padPre;
padPost;
cropPre;
cropPost;
propKernel;
end % end properties
properties (SetAccess = 'private')
settings;
sizeIn;
fresnelNumbers
numDistances;
numDims;
padPre;
padPost;
cropPre;
cropPost;
propKernel;
methods
%% Constructor
function obj = FresnelPropagator(sizeIn, fresnelNumbers, settings)
end % end properties
if nargin < 3
settings = struct;
end
obj.settings = completeStruct(settings, FresnelPropagator.defaultSettings());
obj.sizeIn = sizeIn;
obj.fresnelNumbers = fresnelNumbers;
obj.numDims = numel(obj.sizeIn);
obj.numDistances = size(fresnelNumbers,2);
% Check consistency of input-parameters
%
% Input-, output- and padded array sizes
if isempty(obj.settings.sizePad)
obj.settings.sizePad = sizeIn;
elseif numel(obj.settings.sizePad) ~= numel(sizeIn) || any(obj.settings.sizePad < sizeIn)
error('settings.sizePad must be of the same length as sizeIn and settings.sizePad >= sizeIn.');
end
if isempty(obj.settings.sizeOut)
obj.settings.sizeOut = obj.settings.sizePad;
elseif numel(obj.settings.sizeOut) ~= numel(sizeIn) || any(obj.settings.sizeOut > obj.settings.sizePad)
error('settings.sizeOut must be of the same length as sizeIn and settings.sizeOut <= settings.sizePad.');
end
%
% Check sanity of Fresnel numbers
if size(obj.fresnelNumbers,1) > 1 && size(obj.fresnelNumbers,1) ~= numel(sizeIn)
error(['size(fresnelNumbers,1) must be either 1 (no astigmatism) or equal', ...
' to the number of dimensions numel(sizeIn) (possibly astigmatic setting).']);
end
%
% Check whether GPU-computations are possible
if obj.settings.gpuFlag && ~checkGPUSupport();
warning(['Unable to initialize GPU. Fresnel-propagation will be performed on CPU instead. ', ...
'Set settings.gpuFlag = false to suppress this warning.']);
obj.settings.gpuFlag = false;
end
% Ensure that Fresnel-propagation will be performed in the correct precision (single/double)
% by casting Fresnel numbers to single/double
if obj.settings.singlePrecision
obj.fresnelNumbers = single(obj.fresnelNumbers);
else
obj.fresnelNumbers = double(obj.fresnelNumbers);
end
% Determine padding- and cropping-amounts (and whether to pad/crop at all)
obj.padPre = ceil((obj.settings.sizePad-obj.sizeIn)/2);
obj.padPost = floor((obj.settings.sizePad-obj.sizeIn)/2);
obj.cropPre = ceil((obj.settings.sizePad-obj.settings.sizeOut)/2);
obj.cropPost = floor((obj.settings.sizePad-obj.settings.sizeOut)/2);
methods
%% Constructor
function obj = FresnelPropagator(sizeIn, fresnelNumbers, settings)
% Construct propagation kernel to be applied in Fourier-space
for idxFresnel = 1:obj.numDistances
obj.propKernel{idxFresnel} = fresnelPropagationKernel(obj.settings.sizePad, obj.fresnelNumbers(:,idxFresnel), obj.settings);
end
obj.propKernel = cat(obj.numDims+1, obj.propKernel{:});
end %end constructor
if nargin < 3
settings = struct;
end
obj.settings = completeStruct(settings, FresnelPropagator.defaultSettings());
obj.sizeIn = sizeIn;
obj.fresnelNumbers = fresnelNumbers;
obj.numDims = numel(obj.sizeIn);
obj.numDistances = size(fresnelNumbers,2);
% Check consistency of input-parameters
%
% Input-, output- and padded array sizes
if isempty(obj.settings.sizePad)
obj.settings.sizePad = sizeIn;
elseif numel(obj.settings.sizePad) ~= numel(sizeIn) || any(obj.settings.sizePad < sizeIn)
error('settings.sizePad must be of the same length as sizeIn and settings.sizePad >= sizeIn.');
end
if isempty(obj.settings.sizeOut)
obj.settings.sizeOut = obj.settings.sizePad;
elseif numel(obj.settings.sizeOut) ~= numel(sizeIn) || any(obj.settings.sizeOut > obj.settings.sizePad)
error('settings.sizeOut must be of the same length as sizeIn and settings.sizeOut <= settings.sizePad.');
end
%
% Check sanity of Fresnel numbers
if size(obj.fresnelNumbers,1) > 1 && size(obj.fresnelNumbers,1) ~= numel(sizeIn)
error(['size(fresnelNumbers,1) must be either 1 (no astigmatism) or equal', ...
' to the number of dimensions numel(sizeIn) (possibly astigmatic setting).']);
end
%
% Check whether GPU-computations are possible
if obj.settings.gpuFlag && ~checkGPUSupport();
warning(['Unable to initialize GPU. Fresnel-propagation will be performed on CPU instead. ', ...
'Set settings.gpuFlag = false to suppress this warning.']);
obj.settings.gpuFlag = false;
end
function imProp = propagate(obj, im)
% PROPAGATE Fresnel-propagates a given image according to the internal parameters
% (Fresnel-numbers, padding, etc) of the Fresnel-propagator object.
%
% ``imProp = propagate(im)``
% Transfer to GPU if required
imProp = gpuArrayIf(im, obj.settings.gpuFlag && ~isOnGPU(im));
% Optional padding
if any(obj.settings.sizePad ~= obj.sizeIn)
imProp = padarray(padarray(imProp, obj.padPre, obj.settings.padVal, 'pre'), ...
obj.padPost, obj.settings.padVal, 'post');
end
% Propagation step
imProp = obj.iFT(obj.propKernel .* fftn(imProp));
% Optional cropping
if any(obj.settings.sizePad ~= obj.settings.sizeOut)
imProp = croparray(imProp, obj.cropPre, obj.cropPost);
end
% Transfer back from GPU to host if required
imProp = gatherIf(imProp, obj.settings.gpuFlag && ~isOnGPU(im));
% Ensure that Fresnel-propagation will be performed in the correct precision (single/double)
% by casting Fresnel numbers to single/double
if obj.settings.singlePrecision
obj.fresnelNumbers = single(obj.fresnelNumbers);
else
obj.fresnelNumbers = double(obj.fresnelNumbers);
end
function imBack = backPropagate(obj, im)
% BACKPROPAGATE Fresnel-back-propagates a given image, i.e. applies the inverse
% operation to propagate. In the case of Fresnel-propagation to more than one
% plane, the result is the sum of the back-propagated images from each of the planes.
%
% ``imBack = backPropagate(im)``
% Determine padding- and cropping-amounts (and whether to pad/crop at all)
obj.padPre = ceil((obj.settings.sizePad-obj.sizeIn)/2);
obj.padPost = floor((obj.settings.sizePad-obj.sizeIn)/2);
obj.cropPre = ceil((obj.settings.sizePad-obj.settings.sizeOut)/2);
obj.cropPost = floor((obj.settings.sizePad-obj.settings.sizeOut)/2);
% Transfer to GPU if required
imBack = gpuArrayIf(im, obj.settings.gpuFlag && ~isOnGPU(im));
% Padding from size obj.settings.sizeOut to obj.settings.sizePad
if any(obj.settings.sizePad ~= obj.settings.sizeOut)
imBack = padarray(padarray(imBack, obj.cropPre, obj.settings.padVal, 'pre'), ...
obj.cropPost, obj.settings.padVal, 'post');
end
% Back-propagation step (plus summation in the case of multiple Fresnel-numbers)
imBack = conj(obj.propKernel) .* obj.FT(imBack);
if obj.numDistances > 1
imBack = sum(imBack, obj.numDims+1);
end
imBack = ifftn(imBack);
% Cropping from size obj.sizePad to obj.sizeIn
if any(obj.settings.sizePad ~= obj.sizeIn)
imBack = croparray(imBack, obj.padPre, obj.padPost);
end