diff --git a/Cpufit/python/pycpufit/cpufit.py b/Cpufit/python/pycpufit/cpufit.py index 8e79bc01..771a753d 100755 --- a/Cpufit/python/pycpufit/cpufit.py +++ b/Cpufit/python/pycpufit/cpufit.py @@ -29,9 +29,9 @@ # %% cpufit_func.restype = c_int -cpufit_func.argtypes = [c_size_t, c_size_t, POINTER(c_float), POINTER(c_float), c_int, POINTER(c_float), - c_float, c_int, POINTER(c_int), c_int, c_size_t, - POINTER(c_char), POINTER(c_float), POINTER(c_int), POINTER(c_float), POINTER(c_int)] +cpufit_func.argtypes = [c_size_t, c_size_t, POINTER(c_float), POINTER(c_float), c_int, POINTER(c_float), + c_float, c_int, POINTER(c_int), c_int, c_size_t, + POINTER(c_char), POINTER(c_float), POINTER(c_int), POINTER(c_float), POINTER(c_int)] # int cpufit # ( @@ -254,7 +254,7 @@ def fit_constrained(data, weights, model_id, initial_parameters, constraints=Non else: constraints_p = None if user_info is not None: - user_info_p = user_info.ctypes.data_as(cpufit_func.argtypes[13]) + user_info_p = user_info.ctypes.data_as(cpufit_func.argtypes[11]) else: user_info_p = None diff --git a/Gpufit/matlab/matlab/gauss1d.m b/Gpufit/matlab/matlab/gauss1d.m deleted file mode 100644 index 53dd2492..00000000 --- a/Gpufit/matlab/matlab/gauss1d.m +++ /dev/null @@ -1,51 +0,0 @@ -function gauss1d() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% 1D fitting performance estimation -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); -rng(0); - -%% parameter -n_points = 20; % number of data points -true_parameters = single([30, 9, 8 / 2.35, 1].'); % true model parameters -tolerance = 1e-4; -max_n_iterations = 10; -model_id = ModelID.GAUSS_1D; -estimator_id = EstimatorID.LSE; -n_fits = 5e7; - -%% generate data and weights -x = (0 : n_points-1).'; -p = true_parameters; -m = p(1) * exp(-(x-p(2)).^2/(2*p(3)^2))+p(4); - -data = repmat(m, 1, n_fits); -data = data + randn(size(data)) .* sqrt(data); -weights = 1 ./ max(1, data); - -%% generate initial parameters -initial_parameters = repmat(true_parameters, 1, n_fits); -initial_parameters(1, :) = initial_parameters(1, :) + rand(1, n_fits) * 3 - 1.5; -initial_parameters(2, :) = initial_parameters(2, :) + rand(1, n_fits) * 3 - 1.5; -initial_parameters(3, :) = initial_parameters(3, :) + rand(1, n_fits) * 1 - 0.5; -initial_parameters(4, :) = initial_parameters(4, :) + rand(1, n_fits) * 1 - 0.5; - -%% run Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit(data, weights, ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id); - -for i = 1 : 4 - fprintf('p%d: %.2f est: %.2f (%.2f)\n', i, p(i), mean(parameters(i, :)), std(parameters(i, :), [], 2)); -end - -fprintf('evaluation of %d million fits took %.1fs (Mfits/s = %.2f)\n', n_fits / 1e6, time, n_fits / 1e6 / time); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/gauss2d.m b/Gpufit/matlab/matlab/gauss2d.m deleted file mode 100644 index f525eae4..00000000 --- a/Gpufit/matlab/matlab/gauss2d.m +++ /dev/null @@ -1,188 +0,0 @@ -function gauss2d() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% Multiple fits of a 2D Gaussian peak function with Poisson distributed noise -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); - -% perform some 2D Gaussian peak fits with a symmetrical Gaussian peak -fit_gauss2d(); - -% perform some 2D Gaussian peak fits with an asymmetrical, rotated Gaussian peak -fit_gauss2d_rotated(); - -end -function fit_gauss2d() - -%% number of fits and fit points -number_fits = 1e4; -size_x = 20; -number_parameters = 5; - -%% set input arguments - -% true parameters -true_parameters = single([20, 9.5, 9.5, 3, 10]); - -% initialize random number generator -rng(0); - -% initial parameters (randomized) -initial_parameters = repmat(single(true_parameters'), [1, number_fits]); -% randomize relative to width for positions -initial_parameters([2,3], :) = initial_parameters([2,3], :) + true_parameters(4) * (-0.2 + 0.4 * rand(2, number_fits)); -% randomize relative for other parameters -initial_parameters([1,4,5], :) = initial_parameters([1,4,5], :) .* (0.8 + 0.4 * rand(3, number_fits)); - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% generate data with Poisson noise -data = gaussian_2d(x, y, true_parameters); -data = repmat(data(:), [1, number_fits]); -data = poissrnd(data); - -% tolerance -tolerance = 1e-3; - -% maximum number of iterations -max_n_iterations = 20; - -% estimator id -estimator_id = EstimatorID.MLE; - -% model ID -model_id = ModelID.GAUSS_2D; - -%% run Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, []); - -%% displaying results -display_results('2D Gaussian peak', model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations); - -end - -function fit_gauss2d_rotated() - -%% number of fits and fit points -number_fits = 1e4; -size_x = 20; -number_parameters = 7; - -%% set input arguments - -% true parameters -true_parameters = single([200, 9.5, 9.5, 3, 4, 10, 0.5]); - -% initialize random number generator -rng(0); - -% initial parameters (randomized) -initial_parameters = repmat(single(true_parameters'), [1, number_fits]); -% randomize relative to width for positions -initial_parameters(2, :) = initial_parameters(2, :) + true_parameters(4) * (-0.2 + 0.4 * rand(1, number_fits)); -initial_parameters(3, :) = initial_parameters(3, :) + true_parameters(5) * (-0.2 + 0.4 * rand(1, number_fits)); -% randomize relative for other parameters -initial_parameters([1,4,5,6,7], :) = initial_parameters([1,4,5,6,7], :) .* (0.8 + 0.4 * rand(5, number_fits)); - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% generate data with Poisson noise -data = gaussian_2d_rotated(x, y, true_parameters); -data = repmat(data(:), [1, number_fits]); -data = poissrnd(data); - -% tolerance -tolerance = 1e-3; - -% maximum number of iterations -max_n_iterations = 20; - -% estimator id -estimator_id = EstimatorID.MLE; - -% model ID -model_id = ModelID.GAUSS_2D_ROTATED; - -%% run Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, []); - -%% displaying results -display_results('2D rotated Gaussian peak', model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations); - - -end - -function g = gaussian_2d(x, y, p) -% Generates a 2D Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#gauss-2d -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -g = p(1) * exp(-((x - p(2)).^2 + (y - p(3)).^2) / (2 * p(4)^2)) + p(5); - -end - -function g = gaussian_2d_rotated(x, y, p) -% Generates a 2D rotated elliptic Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#d-rotated-elliptic-gaussian-peak -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -% cosine and sine of rotation angle -cp = cos(p(7)); -sp = sin(p(7)); - -% Gaussian peak with two axes -arga = (x - p(2)) .* cp - (y - p(3)) .* sp; -argb = (x - p(2)) .* sp + (y - p(3)) .* cp; -ex = exp(-0.5 .* (((arga / p(4)) .* (arga / p(4))) + ((argb / p(5)) .* (argb / p(5))))); -g = p(1) .* ex + p(6); - -end - -function display_results(name, model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations) - -%% displaying results -converged = states == 0; -fprintf('\nGpufit of %s\n', name); - -% print summary -fprintf('\nmodel ID: %d\n', model_id); -fprintf('number of fits: %d\n', number_fits); -fprintf('fit size: %d x %d\n', size_x, size_x); -fprintf('mean chi-square: %6.2f\n', mean(chi_squares(converged))); -fprintf('mean iterations: %6.2f\n', mean(n_iterations(converged))); -fprintf('time: %6.2f s\n', time); - -% get fit states -number_converged = sum(converged); -fprintf('\nratio converged %6.2f %%\n', number_converged / number_fits * 100); -fprintf('ratio max it. exceeded %6.2f %%\n', sum(states == 1) / number_fits * 100); -fprintf('ratio singular hessian %6.2f %%\n', sum(states == 2) / number_fits * 100); -fprintf('ratio neg curvature MLE %6.2f %%\n', sum(states == 3) / number_fits * 100); - -% mean and std of fitted parameters -converged_parameters = parameters(:, converged); -converged_parameters_mean = mean(converged_parameters, 2); -converged_parameters_std = std(converged_parameters, [], 2); -fprintf('\nparameters of %s\n', name); -for i = 1 : number_parameters - fprintf('p%d true %6.2f mean %6.2f std %6.2f\n', i, true_parameters(i), converged_parameters_mean(i), converged_parameters_std(i)); -end - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/gauss2d_comparison.m b/Gpufit/matlab/matlab/gauss2d_comparison.m deleted file mode 100644 index 9618bfc0..00000000 --- a/Gpufit/matlab/matlab/gauss2d_comparison.m +++ /dev/null @@ -1,212 +0,0 @@ -function gauss2d_comparison() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% Multiple fits of a 2D Gaussian peak function with Poisson distributed noise -% compared to a generic Matlab implementation using fminunc and supplying -% the gradient by the user (uses quasi-newton as algorithm) -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); - -%% number of fits and fit points -number_fits = 1e3; -size_x = 20; -number_parameters = 5; - -%% set input arguments - -% true parameters -true_parameters = single([10, 9.5, 9.5, 3, 10]); - -% initialize random number generator -rng(0); - -% initial parameters (randomized) -initial_parameters = repmat(single(true_parameters'), [1, number_fits]); -% randomize relative to width for positions -initial_parameters([2,3], :) = initial_parameters([2,3], :) + true_parameters(4) * (-0.2 + 0.4 * rand(2, number_fits)); -% randomize relative for other parameters -initial_parameters([1,4,5], :) = initial_parameters([1,4,5], :) .* (0.8 + 0.4 * rand(3, number_fits)); - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% generate data with Poisson noise -data = gaussian_2d(x, y, true_parameters); -data = repmat(data(:), [1, number_fits]); -data = poissrnd(data); - -% tolerance -tolerance = 1e-4; - -% maximum number of iterations -max_n_iterations = 20; - -% estimator id -estimator_id = EstimatorID.MLE; - -% model ID -model_id = ModelID.GAUSS_2D; % Gaussian peak in 2D - -%% run Gpufit -fprintf('run Gpufit\n'); -[gf_parameters, gf_states, gf_chi_squares, gf_n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, []); - -% display results -display_results('Gpufit', gf_parameters, gf_states, gf_chi_squares, gf_n_iterations, time, true_parameters); - -% store parameters - -%% run Matlab - -% convert data and initial_parameters to double (otherwise causes an error -% in fminunc) -data = double(data); -initial_parameters = double(initial_parameters); -xi = double(x(:)'); -yi = double(y(:)'); - -% set fit options -options = optimoptions(@fminunc,'Display', 'off', 'MaxIter', max_n_iterations, 'Algorithm', 'quasi-newton', 'TolFun', tolerance, 'GradObj', 'on', 'DerivativeCheck', 'off', 'Diagnostics', 'off'); - -% initialize output arrays -m_parameters = zeros(number_parameters, number_fits); -m_states = zeros(1, number_fits); -m_chi_squares = zeros(1, number_fits); -m_n_iterations = zeros(1, number_fits); - -% loop over each fit -fprintf('\n') -progress = 0; -L = 50; % length of progressbar -tic; -for i = 1 : number_fits - - % get data and initial_parameters - d = data(:, i)'; - p0 = initial_parameters(:, i); - - % define minimizer function (give grid and data as implicit parameters) - fun = @(p) minimizer(p, xi, yi, d); - - % call to fminunc - [p, fval, exitflag, output] = fminunc(fun, p0, options); - - % copy to output - m_parameters(:, i) = p; - m_chi_squares(i) = fval; - m_states(i) = exitflag - 1; - m_n_iterations(i) = output.iterations; - - progress = progress + 1; - if progress >= number_fits / L - progress = 0; - fprintf('|'); - end -end -time = toc; -fprintf(repmat('\b', [1, L])); - -% display results -display_results('Matlab (one CPU kernel)', m_parameters, m_states, m_chi_squares, m_n_iterations, time, true_parameters); - -end - -function [f, g] = minimizer(p, xi, yi, d) -% calls the model with the current parameters, then the likelihood function -% and returns value and derivatives of the likelihood function -% -% p - current parameters -% xi, yi - grid positions -% d - current data - -if nargout > 1 - [m, mg] = gaussian_2d_with_gradient(xi, yi, p); - [f, g] = poisson_likelihood(m, mg, d); -else - m = gaussian_2d(xi, yi, p); - f = poisson_likelihood(m, [], d); -end - -end - -function [f, g] = poisson_likelihood(m, mg, d) -% Calculates value and derivatives of the poisson likelihood function for -% given model and model derivatives - -h = d > 0; -f = 2 * (sum(m-d) - sum(d(h) .* log(m(h) ./ d(h)))); - -if nargout > 1 % gradient required - h = 2 * (1 - d ./ max(m, 1e-6)); - h = repmat(h, [size(mg, 1), 1]); - g = h .* mg; - g = sum(g, 2); -end - -end - - -function display_results(name, parameters, ~, chi_squares, n_iterations, time, true_parameters) -% displaying results - -fprintf('*%s*\n', name); -number_parameters = size(parameters, 1); -number_fits = size(parameters, 2); - -% print summary -fprintf('\nnumber of fits: %d\n', number_fits); -fprintf('mean chi-square: %6.2f\n', mean(chi_squares)); -fprintf('mean iterations: %6.2f\n', mean(n_iterations)); -fprintf('time: %6.2f s\n', time); -fprintf('fits per second: %.0f\n', number_fits / time); - -% mean and std of fitted parameters -parameters_mean = mean(parameters, 2); -parameters_std = std(parameters, [], 2); -fprintf('\nparameters of 2D Gaussian peak\n'); -for i = 1 : number_parameters - fprintf('p%d true %6.2f mean %6.2f std %6.2f\n', i, true_parameters(i), parameters_mean(i), parameters_std(i)); -end - -end - -function f = gaussian_2d(x, y, p) -% Generates a 2D Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#gauss-2d -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -f = p(1) * exp(-((x - p(2)).^2 + (y - p(3)).^2) / (2 * p(4)^2)) + p(5); - -end - -function [f, g] = gaussian_2d_with_gradient(x, y, p) -% Computes the gradient for a 2D Gaussian peak with respect to parameters. - -dx = x - p(2); -dy = y - p(3); -p42 = p(4)^2; -arg = (dx.^2 + dy.^2) / p42; -exp_f = exp(-0.5 * arg); -p1_exp_f = p(1) * exp_f; - -f = p1_exp_f + p(5); - -g1 = exp_f; -g2 = p1_exp_f .* dx / p42; -g3 = p1_exp_f .* dy / p42; -g4 = p1_exp_f .* arg / p(4); -g5 = ones(size(x)); -g = [g1; g2; g3; g4; g5]; - -end diff --git a/Gpufit/matlab/matlab/gauss2d_constrained.m b/Gpufit/matlab/matlab/gauss2d_constrained.m deleted file mode 100644 index 263ae47c..00000000 --- a/Gpufit/matlab/matlab/gauss2d_constrained.m +++ /dev/null @@ -1,141 +0,0 @@ -function gauss2d_constrained() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% Multiple fits of a 2D symmetric Gaussian peak function -% Comparison of unconstrained vs. constrained fit -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); - -%% number of fits and fit points -number_fits = 2e5; -size_x = 20; -number_parameters = 5; - -%% set input arguments - -% true parameters (amplitude, center, center, sigma, background) -true_parameters = single([2, 9.5, 9.5, 3, 10]); - -% initialize random number generator -rng(0); - -% initial parameters (randomized) -initial_parameters = repmat(single(true_parameters'), [1, number_fits]); -% randomize relative to width for positions -initial_parameters([2,3], :) = initial_parameters([2,3], :) + true_parameters(4) * (-0.2 + 0.4 * rand(2, number_fits)); -% randomize relative for other parameters -initial_parameters([1,4,5], :) = initial_parameters([1,4,5], :) .* (0.8 + 0.4 * rand(3, number_fits)); - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% generate data with Poisson noise -data = gaussian_2d(x, y, true_parameters); -data = repmat(data(:), [1, number_fits]); -data = data + randn(size(data)).* sqrt(data); - -% tolerance -tolerance = 1e-3; - -% maximum number of iterations -max_n_iterations = 20; - -% estimator id -estimator_id = EstimatorID.LSE; - -% model ID -model_id = ModelID.GAUSS_2D; - -%% run unconstrained Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, []); - -%% displaying results -display_results('unconstrained 2D Gaussian peak', model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations); - -%% set constraints -constraints = zeros([2*number_parameters, number_fits], 'single'); -constraints(7, :) = 2.9; -constraints(8, :) = 3.1; -constraint_types = int32([ConstraintType.LOWER, ConstraintType.FREE, ConstraintType.FREE, ConstraintType.LOWER_UPPER, ConstraintType.LOWER]); - -%% run constrained Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit_constrained(data, [], ... - model_id, initial_parameters, constraints, constraint_types, tolerance, max_n_iterations, [], estimator_id, []); - -%% displaying results -display_results('constrained 2D Gaussian peak', model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations); - -% gist: if I know the width of a peak really well before-hand, I can estimate its position better - -end - -function g = gaussian_2d(x, y, p) -% Generates a 2D Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#gauss-2d -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -g = p(1) * exp(-((x - p(2)).^2 + (y - p(3)).^2) / (2 * p(4)^2)) + p(5); - -end - -function g = gaussian_2d_rotated(x, y, p) -% Generates a 2D rotated elliptic Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#d-rotated-elliptic-gaussian-peak -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -% cosine and sine of rotation angle -cp = cos(p(7)); -sp = sin(p(7)); - -% Gaussian peak with two axes -arga = (x - p(2)) .* cp - (y - p(3)) .* sp; -argb = (x - p(2)) .* sp + (y - p(3)) .* cp; -ex = exp(-0.5 .* (((arga / p(4)) .* (arga / p(4))) + ((argb / p(5)) .* (argb / p(5))))); -g = p(1) .* ex + p(6); - -end - -function display_results(name, model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations) - -%% displaying results -converged = states == 0; -fprintf('\nGpufit of %s\n', name); - -% print summary -fprintf('\nmodel ID: %d\n', model_id); -fprintf('number of fits: %d\n', number_fits); -fprintf('fit size: %d x %d\n', size_x, size_x); -fprintf('mean chi-square: %6.2f\n', mean(chi_squares(converged))); -fprintf('mean iterations: %6.2f\n', mean(n_iterations(converged))); -fprintf('time: %6.2f s\n', time); - -% get fit states -number_converged = sum(converged); -fprintf('\nratio converged %6.2f %%\n', number_converged / number_fits * 100); -fprintf('ratio max it. exceeded %6.2f %%\n', sum(states == 1) / number_fits * 100); -fprintf('ratio singular hessian %6.2f %%\n', sum(states == 2) / number_fits * 100); -fprintf('ratio neg curvature MLE %6.2f %%\n', sum(states == 3) / number_fits * 100); - -% mean and std of fitted parameters -converged_parameters = parameters(:, converged); -converged_parameters_mean = mean(converged_parameters, 2); -converged_parameters_std = std(converged_parameters, [], 2); -fprintf('\nparameters of %s\n', name); -for i = 1 : number_parameters - fprintf('p%d true %6.2f mean %6.2f std %6.2f\n', i, true_parameters(i), converged_parameters_mean(i), converged_parameters_std(i)); -end - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/gauss2d_fixed_parameters.m b/Gpufit/matlab/matlab/gauss2d_fixed_parameters.m deleted file mode 100644 index 55e18f03..00000000 --- a/Gpufit/matlab/matlab/gauss2d_fixed_parameters.m +++ /dev/null @@ -1,197 +0,0 @@ -function gauss2d_fixed_parameters() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% Multiple fits of a 2D Gaussian peak function ones without fixed -% parameters and once with fixed parameters, measures performance -% Longest times is spent in preparing the data -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); - -% perform some 2D Gaussian peak fits with an assymetrical Gaussian peak -fit_gauss2d_elliptic(); - -% perform some 2D Gaussian peak fits with an asymmetrical, rotated Gaussian -% peak (but with fixed rotation angle) -fit_gauss2d_rotated_fixed(); - -end -function fit_gauss2d_elliptic() - -%% number of fits and fit points -number_fits = 1e4; -size_x = 20; -number_parameters = 6; - -%% set input arguments - -% true parameters -true_parameters = single([20, 9.5, 9.5, 3, 4, 10]); - -% initialize random number generator -rng(0); - -% initial parameters (randomized) -initial_parameters = repmat(single(true_parameters'), [1, number_fits]); -% randomize relative to width for positions -initial_parameters(2, :) = initial_parameters(2, :) + true_parameters(4) * (-0.2 + 0.4 * rand(1, number_fits)); -initial_parameters(3, :) = initial_parameters(3, :) + true_parameters(5) * (-0.2 + 0.4 * rand(1, number_fits)); -% randomize relative for other parameters -initial_parameters([1,4,5,6], :) = initial_parameters([1,4,5,6], :) .* (0.8 + 0.4 * rand(4, number_fits)); - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% generate data with Poisson noise -data = gaussian_2d_rotated(x, y, [true_parameters, 0]); -data = repmat(data(:), [1, number_fits]); -data = poissrnd(data); - -% tolerance -tolerance = 1e-3; - -% maximum number of iterations -max_n_iterations = 20; - -% estimator id -estimator_id = EstimatorID.MLE; - -% model ID -model_id = ModelID.GAUSS_2D_ELLIPTIC; - -%% run Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, []); - -%% displaying results -display_results('2D elliptic Gaussian peak', model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations); - -end - -function fit_gauss2d_rotated_fixed() -% with the rotation as fixed parameter - -%% number of fits and fit points -number_fits = 1e4; -size_x = 20; -number_parameters = 7; - -%% set input arguments - -% true parameters -true_parameters = single([20, 9.5, 9.5, 3, 4, 10, 0]); - -% initialize random number generator -rng(0); - -% initial parameters (randomized) -initial_parameters = repmat(single(true_parameters'), [1, number_fits]); -% randomize relative to width for positions -initial_parameters(2, :) = initial_parameters(2, :) + true_parameters(4) * (-0.2 + 0.4 * rand(1, number_fits)); -initial_parameters(3, :) = initial_parameters(3, :) + true_parameters(5) * (-0.2 + 0.4 * rand(1, number_fits)); -% randomize relative for other parameters -initial_parameters([1,4,5,6], :) = initial_parameters([1,4,5,6], :) .* (0.8 + 0.4 * rand(4, number_fits)); - -% set fixed parameters -parameters_to_fit = ones([number_parameters, 1], 'int32'); -parameters_to_fit(7) = 0; - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% generate data with Poisson noise -data = gaussian_2d_rotated(x, y, true_parameters); -data = repmat(data(:), [1, number_fits]); -data = poissrnd(data); - -% tolerance -tolerance = 1e-3; - -% maximum number of iterations -max_n_iterations = 20; - -% estimator id -estimator_id = EstimatorID.MLE; - -% model ID -model_id = ModelID.GAUSS_2D_ROTATED; - -%% run Gpufit -[parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, parameters_to_fit, estimator_id, []); - -%% displaying results -display_results('2D rotated Gaussian peak (fixed parameter)', model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations); - - -end - -function g = gaussian_2d(x, y, p) -% Generates a 2D Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#gauss-2d -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -g = p(1) * exp(-((x - p(2)).^2 + (y - p(3)).^2) / (2 * p(4)^2)) + p(5); - -end - -function g = gaussian_2d_rotated(x, y, p) -% Generates a 2D rotated elliptic Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/api.html#d-rotated-elliptic-gaussian-peak -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -% cosine and sine of rotation angle -cp = cos(p(7)); -sp = sin(p(7)); - -% Gaussian peak with two axes -arga = (x - p(2)) .* cp - (y - p(3)) .* sp; -argb = (x - p(2)) .* sp + (y - p(3)) .* cp; -ex = exp(-0.5 .* (((arga / p(4)) .* (arga / p(4))) + ((argb / p(5)) .* (argb / p(5))))); -g = p(1) .* ex + p(6); - -end - -function display_results(name, model_id, number_fits, number_parameters, size_x, time, true_parameters, parameters, states, chi_squares, n_iterations) - -%% displaying results -converged = states == 0; -fprintf('\nGpufit of %s\n', name); - -% print summary -fprintf('\nmodel ID: %d\n', model_id); -fprintf('number of fits: %d\n', number_fits); -fprintf('fit size: %d x %d\n', size_x, size_x); -fprintf('mean chi-square: %6.2f\n', mean(chi_squares(converged))); -fprintf('mean iterations: %6.2f\n', mean(n_iterations(converged))); -fprintf('time: %6.2f s\n', time); - -% get fit states -number_converged = sum(converged); -fprintf('\nratio converged %6.2f %%\n', number_converged / number_fits * 100); -fprintf('ratio max it. exceeded %6.2f %%\n', sum(states == 1) / number_fits * 100); -fprintf('ratio singular hessian %6.2f %%\n', sum(states == 2) / number_fits * 100); -fprintf('ratio neg curvature MLE %6.2f %%\n', sum(states == 3) / number_fits * 100); - -% mean and std of fitted parameters -converged_parameters = parameters(:, converged); -converged_parameters_mean = mean(converged_parameters, 2); -converged_parameters_std = std(converged_parameters, [], 2); -fprintf('\nparameters of %s\n', name); -for i = 1 : number_parameters - fprintf('p%d true %6.2f mean %6.2f std %6.2f\n', i, true_parameters(i), converged_parameters_mean(i), converged_parameters_std(i)); -end - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/gauss2d_plot.m b/Gpufit/matlab/matlab/gauss2d_plot.m deleted file mode 100644 index df7a0f8d..00000000 --- a/Gpufit/matlab/matlab/gauss2d_plot.m +++ /dev/null @@ -1,123 +0,0 @@ -function gauss2d_plot() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% Multiple fits of a 2D Gaussian peak function with Poisson distributed noise -% repeated for a different total number of fits each time and plotting the -% results -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); - -%% number of fit points -size_x = 5; -n_points = size_x * size_x; - -%% set input arguments - -% mean true parameters -mean_true_parameters = single([100, 3, 3, 1, 10]); - -% average noise level -average_noise_level = 10; - -% initialize random number generator -rng(0); - -% tolerance -tolerance = 1e-4; - -% max number of itetations -max_n_iterations = 10; - -% model id -model_id = ModelID.GAUSS_2D; - -%% loop over different number of fits -n_fits_all = round(logspace(2, 6, 20)); - -% generate x and y values -g = single(0 : size_x - 1); -[x, y] = ndgrid(g, g); - -% loop -speed = zeros(length(n_fits_all), 1); -for i = 1:length(n_fits_all) - n_fits = n_fits_all(i); - - % vary positions of 2D Gaussians peaks slightly - test_parameters = repmat(mean_true_parameters', [1, n_fits]); - test_parameters([2,3], :) = test_parameters([2,3], :) + mean_true_parameters(4) * (-0.2 + 0.4 * rand(2, n_fits)); - - % generate data - data = gaussians_2d(x, y, test_parameters); - data = reshape(data, [n_points, n_fits]); - - % add noise - data = data + average_noise_level * randn(size(data), 'single'); - - % initial parameters (randomized) - initial_parameters = repmat(mean_true_parameters', [1, n_fits]); - % randomize relative to width for positions - initial_parameters([2,3], :) = initial_parameters([2,3], :) + mean_true_parameters(4) * (-0.2 + 0.4 * rand(2, n_fits)); - % randomize relative for other parameters - initial_parameters([1,4,5], :) = initial_parameters([1,4,5], :) .* (0.8 + 0.4 * rand(3, n_fits)); - - % run Gpufit - [parameters, states, chi_squares, n_iterations, time] = gpufit(data, [], ... - model_id, initial_parameters, tolerance, max_n_iterations); - - % analyze result - converged = states == 0; - speed(i) = n_fits / time; - precision_x0 = std(parameters(2, converged) - test_parameters(2, converged)); - - % display result - fprintf(' iterations: %.2f | time: %.3f s | speed: %8.0f fits/s\n', ... - mean(n_iterations(converged)), time, speed(i)); -end - -%% plot -figure(); -semilogx(n_fits_all, speed, 'bo-') -xlabel('number of fits per function call') -ylabel('fits per second') -legend('Gpufit', 'Location', 'NorthWest') -grid on; -xlim(n_fits_all([1,end])); - -end - -function g = gaussians_2d(x, y, p) -% Generates many 2D Gaussians peaks for a given set of parameters - -n_fits = size(p, 2); -msg = sprintf('generating %d fits ', n_fits); -fprintf(msg); - -g = zeros([size(x), n_fits], 'single'); - -progress = 0; -L = 50; % length of progressbar -l = 0; -for i = 1 : n_fits - - pi = p(:, i); - g(:, :, i) = pi(1) * exp(-((x - pi(2)).^2 + (y - pi(3)).^2) / (2 * pi(4)^2)) + pi(5); - - progress = progress + 1; - if progress >= n_fits / L - progress = 0; - fprintf('|'); - l = l + 1; - end -end -fprintf(repmat('\b', [1, length(msg) + l])); -fprintf('%7d fits', n_fits); - -end diff --git a/Gpufit/matlab/matlab/gaussian_peak_1d.m b/Gpufit/matlab/matlab/gaussian_peak_1d.m deleted file mode 100644 index 2b7c7a6f..00000000 --- a/Gpufit/matlab/matlab/gaussian_peak_1d.m +++ /dev/null @@ -1,12 +0,0 @@ -function f = gaussian_peak_1d(x, p) -% Generates a 1D Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/fit_model_functions.html#d-gaussian-function -% -% x - x grid position values -% p - parameters (amplitude, center position, width, offset) - -assert(nargin == 2); - -f = p(1) * exp(-(x - p(2)).^2 / (2 * p(3)^2)) + p(4); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/gaussian_peak_2d.m b/Gpufit/matlab/matlab/gaussian_peak_2d.m deleted file mode 100644 index fe2f1304..00000000 --- a/Gpufit/matlab/matlab/gaussian_peak_2d.m +++ /dev/null @@ -1,12 +0,0 @@ -function f = gaussian_peak_2d(x, y, p) -% Generates a 2D Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/fit_model_functions.html#d-gaussian-function-cylindrical-symmetry -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, width, offset) - -assert(nargin == 3); - -f = p(1) * exp(-((x - p(2)).^2 + (y - p(3)).^2) / (2 * p(4)^2)) + p(5); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/gaussian_peak_2d_rotated.m b/Gpufit/matlab/matlab/gaussian_peak_2d_rotated.m deleted file mode 100644 index 7f25ce6b..00000000 --- a/Gpufit/matlab/matlab/gaussian_peak_2d_rotated.m +++ /dev/null @@ -1,18 +0,0 @@ -function g = gaussian_peak_2d_rotated(x, y, p) -% Generates a 2D rotated elliptic Gaussian peak. -% http://gpufit.readthedocs.io/en/latest/fit_model_functions.html#d-gaussian-function-elliptical-rotated -% -% x,y - x and y grid position values -% p - parameters (amplitude, x,y center position, x,y-width, offset, rotation angle) - -% cosine and sine of rotation angle -cp = cos(p(7)); -sp = sin(p(7)); - -% Gaussian peak with two axes -arga = (x - p(2)) .* cp - (y - p(3)) .* sp; -argb = (x - p(2)) .* sp + (y - p(3)) .* cp; -ex = exp(-0.5 .* (((arga / p(4)) .* (arga / p(4))) + ((argb / p(5)) .* (argb / p(5))))); -g = p(1) .* ex + p(6); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/linear_regression.m b/Gpufit/matlab/matlab/linear_regression.m deleted file mode 100644 index bc939cbb..00000000 --- a/Gpufit/matlab/matlab/linear_regression.m +++ /dev/null @@ -1,38 +0,0 @@ -function linear_regression() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% 1D linear regression with custom x values given as user information -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -% x values for the 1D linear regression model must be single -x = single(0:10); -parameters = single([0;1]); -y = single(parameters(1)+parameters(2)*x'); - -% fit parameter -tolerance = 1e-9; -max_n_iterations = 1000; -estimator_id = EstimatorID.LSE; -model_id = ModelID.LINEAR_1D; -initial_parameters = single([0;1]); % should result in correct parameters with only 1 iteration - -% with user info -[parameters, states, chi_squares, n_iterations, time] = gpufit(y, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, x); -fprintf('first fit with user info: fitted parameters = [%.2f, %.2f]\n', parameters); - -% without user info -[parameters, states, chi_squares, n_iterations, time] = gpufit(y, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id); -fprintf('first fit without user info: fitted parameters = [%.2f, %.2f]\n', parameters); - -% now more meaningful -x = single([1,3,4,5,5.5]); -y = single(parameters(1)+parameters(2)*x'); -[parameters, states, chi_squares, n_iterations, time] = gpufit(y, [], ... - model_id, initial_parameters, tolerance, max_n_iterations, [], estimator_id, x); -fprintf('second fit with user info: fitted parameters = [%.2f, %.2f]\n', parameters); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/simple.m b/Gpufit/matlab/matlab/simple.m deleted file mode 100644 index 82e4cf1d..00000000 --- a/Gpufit/matlab/matlab/simple.m +++ /dev/null @@ -1,32 +0,0 @@ -function simple() -% Example of the Matlab binding of the Gpufit library implementing -% Levenberg Marquardt curve fitting in CUDA -% https://github.com/gpufit/Gpufit -% -% Simple example demonstrating a minimal call of all needed parameters for the Matlab interface -% http://gpufit.readthedocs.io/en/latest/bindings.html#matlab - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end - -assert(gpufit_cuda_available(), 'CUDA not available'); - -% number of fits, number of points per fit -number_fits = 10; -number_points = 10; - -% model ID and number of parameter -model_id = ModelID.GAUSS_1D; -number_parameter = 4; - -% initial parameters -initial_parameters = zeros(number_parameter, number_fits, 'single'); - -% data -data = zeros(number_points, number_fits, 'single'); - -% run Gpufit -[parameters, states, chi_squares, number_iterations, execution_time] = gpufit(data, [], model_id, initial_parameters); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/splinefit_1d.m b/Gpufit/matlab/matlab/splinefit_1d.m deleted file mode 100644 index aafa987f..00000000 --- a/Gpufit/matlab/matlab/splinefit_1d.m +++ /dev/null @@ -1,148 +0,0 @@ -function splinefit_1d() -% Spline fit 1D example -% -% -% -% Requires Gpuspline (https://github.com/gpufit/Gpuspline) additionally - -% there are two coordinate systems -% the spline coordinate system: 0:1:size-1 -% the user coordinate system: whatever -% be careful to not mix them up (fit coordinates must be in the spline coordinate system!) - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end -if isempty(which('spline_coefficients.m')) - error('Gpuspline library not found in Matlab path.'); -end - -% initialize random number generator -rng(0); - -% data size -size_x = 25; - -% tolerances -tolerance = 1e-30; -max_n_iterations = 100; -estimator_id = EstimatorID.LSE; -model_id = ModelID.SPLINE_1D; - -% derived values -x = (0 : size_x - 1)'; -SF = 2; % scaling factor -x_spline = (0 : SF * (size_x - 1))'; -x2 = single(x_spline / SF); % 0, 0.5, .., size_x -1 - -%% generate PSF (two Gaussians) -psf_parameters = single([100, (size_x-1)/2, 1.5, 10]); -psf = calculate_non_gaussian_psf(x2, psf_parameters); -psf_normalized = (psf - psf_parameters(4)) / psf_parameters(1); -psf_normalized = psf_normalized(1:end-5); % splines do not cover all the data - -%% calculate spline coefficients of the PSF template -coefficients = spline_coefficients(psf_normalized); -n_intervals = size(psf_normalized, 1) - 1; - -%% add noise to PSF data (no shift) -snr = 10; -amplitude = psf_parameters(1); -noise_std_dev = amplitude ./ (snr * log(10.0)); -noise = noise_std_dev * randn(size(psf)); -noisy_psf = psf + noise; - -%% set user info -user_info = [n_intervals, reshape(coefficients,1,numel(coefficients))]; - -%% true fit parameters -true_fit_parameters = zeros(3, 1, 'single'); -true_fit_parameters(1) = psf_parameters(1); % amplitude -true_fit_parameters(2) = 0; % center shift -true_fit_parameters(3) = psf_parameters(4); % offset - -%% set initial fit parameters -pos_shift = -1.2 * SF; -amp_shift = 30; -off_shift = 20; - -fit_initial_parameters = true_fit_parameters + [amp_shift; pos_shift; off_shift]; % small deviation - -gauss_fit_initial_parameters = psf_parameters' + [amp_shift, pos_shift, 0, off_shift]'; - -%% call to gpufit with spline fit -[parameters_gpufit_spline, states_gpufit_spline, chi_squares_gpufit_spline, n_iterations_gpufit_spline, time_gpufit_spline]... - = gpufit(noisy_psf, [], model_id, fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, user_info); -assert(all(states_gpufit_spline == 0)); - -%% call to cpufit with spline fit -[parameters_cpufit_spline, states_cpufit_spline, chi_squares_cpufit_spline, n_iterations_cpufit_spline, time_cpufit_spline]... - = cpufit(noisy_psf, [], model_id, fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, user_info); -assert(all(states_cpufit_spline == 0)); - -%% call to gpufit with gauss1d fit -[parameters_gpufit_gauss, states_gpufit_gauss, chi_squares_gpufit_gauss, n_iterations_gpufit_gauss, time_gpufit_gauss]... - = gpufit(noisy_psf, [], ModelID.GAUSS_1D, gauss_fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, x2); -assert(all(states_gpufit_gauss == 0)); - -%% get data to plot -a = true_fit_parameters(1); -x = x_spline-true_fit_parameters(2); -b = true_fit_parameters(3); -spline_model = a * spline_values(coefficients, x) + b; -a = fit_initial_parameters(1); -x = x_spline-fit_initial_parameters(2); -b = fit_initial_parameters(3); -initial_spline_fit = a * spline_values(coefficients, x) + b; -a = parameters_gpufit_spline(1); -x = x_spline-parameters_gpufit_spline(2); -b = parameters_gpufit_spline(3); -final_fit_gpufit = a * spline_values(coefficients, x) + b; -a = parameters_cpufit_spline(1); -x = x_spline-parameters_cpufit_spline(2); -b = parameters_cpufit_spline(3); -final_fit_cpufit = a * spline_values(coefficients, x) + b; -gauss_final_fit = gaussian_peak_1d(x2, parameters_gpufit_gauss(:, 1)); - -%% make a figure of function values -figure(1); -hold off; -plot(x2, gauss_final_fit, '--.y', 'MarkerSize', 8, 'LineWidth', 2); -hold on; -plot(x2, noisy_psf(:, 1), 'ks', 'MarkerSize', 8, 'LineWidth', 2); -plot(x2, initial_spline_fit,'--sg', 'MarkerSize', 8, 'LineWidth', 2); -plot(x2, final_fit_cpufit,'-xc', 'MarkerSize', 8, 'LineWidth', 2); -plot(x2, final_fit_gpufit,'--+b', 'MarkerSize', 8, 'LineWidth', 2); -plot(x2, spline_model, ':r', 'MarkerSize', 8, 'LineWidth', 1.5); -ylim([0, max(initial_spline_fit)]); -legend(... - 'final gauss fit',... - 'noisy data',... - 'initial spline fit',... - 'final spline fit cpu',... - 'final spline fit gpu',... - 'true parameters spline'); - -end - -function y = calculate_non_gaussian_psf(x, p) -% Test PSF consists of sum of two Gaussians with different centers and -% width -% TODO define center of PSF (for now we assume that center is at middle) - -% p(1) - amplitude first Gaussian -% p(2) - center (both Gaussians shifted a bit towards left and right -% p(3) - Standard deviation (second Gaussian is a bit wider) -% p(4) - constant background - -assert(nargin == 2); - -distance = x(end) / 12; -arg1 = ((x - p(2)-distance).^2) / (2 * p(3)^2); -arg2 = ((x - p(2)+distance).^2) / (2 * p(3)^2); -ex = exp(-arg1) + 0.5 * exp(-arg2); - -y = ex / max(ex); % normalized to [0,1] -y = p(1) * y + p(4); % scale with amplitude and background - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/splinefit_2d.m b/Gpufit/matlab/matlab/splinefit_2d.m deleted file mode 100644 index 87ee0a71..00000000 --- a/Gpufit/matlab/matlab/splinefit_2d.m +++ /dev/null @@ -1,166 +0,0 @@ -function splinefit_2d() -% spline fit 2D rectangular -% -% Requires Gpuspline (https://github.com/gpufit/Gpuspline) additionally - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end -if isempty(which('spline_coefficients.m')) - error('Gpuspline library not found in Matlab path.'); -end - -% initialize random number generator -rng(0); - -% data size -size_x = 12; -size_y = 14; - -% tolerances -tolerance = 1e-30; -max_n_iterations = 100; -estimator_id = EstimatorID.LSE; - -% derived values -x = single(0 : size_x - 1)'; -y = single(0 : size_y - 1); - -SF = 2; % scaling factor -x_spline = single(0 : SF * (size_x - 1))'; -y_spline = single(0 : SF * (size_y - 1)); - -x2 = x_spline / SF; -y2 = y_spline / SF; - -%% generate PSF -psf_parameters = single([100, (size_x-1)/2, (size_y-1)/2, 1, 10]); - -%% calculate PSF template -% calculate PSF on fine grid -psf = calculate_psf(x2, y2, psf_parameters); -% PSF template (normalized to minimum = 0 and maximum = 1) -psf_normalized = (psf - psf_parameters(5)) / psf_parameters(1); - -%% calculate spline coefficients of the PSF template -coefficients = spline_coefficients(psf_normalized); -n_intervals = size(psf_normalized) - 1; - -%% add noise to PSF data -snr = 50; -amplitude = psf_parameters(1); -noise_std_dev = amplitude ./ (snr * log(10.0)); -noise = noise_std_dev * randn(size(psf)); -noisy_psf = psf + noise; - -%% set user info -user_info = [... - numel(x2); numel(y2);... - n_intervals';... - coefficients(:)]; - -%% true fit parameters -true_fit_parameters(1,:) = psf_parameters(1); % amplitude -true_fit_parameters(2,:) = 0; % center x shift -true_fit_parameters(3,:) = 0; % center y shift -true_fit_parameters(4,:) = psf_parameters(5); % offset - -%% set initial fit parameters -pos_shift_x = 1.; -pos_shift_y = -2.1; -amp_shift = 20; -off_shift = 13; - -spline_fit_initial_parameters = true_fit_parameters + [amp_shift; pos_shift_x * SF; pos_shift_y * SF; off_shift]; - -gauss_fit_initial_parameters(1,:) = (psf_parameters(1) + amp_shift); -gauss_fit_initial_parameters(2,:) = (psf_parameters(2) + pos_shift_x) * SF; -gauss_fit_initial_parameters(3,:) = (psf_parameters(3) + pos_shift_y) * SF; -gauss_fit_initial_parameters(4,:) = (psf_parameters(4) + 0) * SF; -gauss_fit_initial_parameters(5,:) = (psf_parameters(5) + off_shift); - -%% reshape data -linear_noisy_psf = reshape(noisy_psf, numel(noisy_psf), 1); -linear_noisy_psf_gauss = reshape(noisy_psf(:, 3:end-2), 23^2, 1); % this will only work if size_x,y aren't changed above (try to cut off data symmetrical to make a square data array for gauss fit) - -%% call to cpufit with spline fit -[parameters_cpufit_spline, states_cpufit_spline, chi_squares_cpufit_spline, n_iterations_cpufit_spline, time_cpufit_spline]... - = cpufit(linear_noisy_psf, [], ModelID.SPLINE_2D, spline_fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, user_info); -assert(all(states_cpufit_spline == 0)); - -%% call to gpufit with spline fit -[parameters_gpufit_spline, states_gpufit_spline, chi_squares_gpufit_spline, n_iterations_gpufit_spline, time_gpufit_spline]... - = gpufit(linear_noisy_psf, [], ModelID.SPLINE_2D, spline_fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, user_info); -assert(all(states_gpufit_spline == 0)); - -%% call to gpufit with gauss1d fit -[parameters_gpufit_gauss, states_gpufit_gauss, chi_squares_gpufit_gauss, n_iterations_gpufit_gauss, time_gpufit_gauss]... - = gpufit(linear_noisy_psf_gauss, [], ModelID.GAUSS_2D, gauss_fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id); -assert(all(states_gpufit_gauss == 0)); - -%% get data to plot -a = spline_fit_initial_parameters(1); -xx = x_spline - spline_fit_initial_parameters(2); -yy = y_spline - spline_fit_initial_parameters(3); -b = spline_fit_initial_parameters(4); -initial_spline_fit = a * spline_values(coefficients, xx, yy) + b; -a = parameters_cpufit_spline(1); -xx = x_spline - parameters_cpufit_spline(2); -yy = y_spline - parameters_cpufit_spline(3); -b = parameters_cpufit_spline(4); -final_spline_fit = a * spline_values(coefficients, xx, yy) + b; -initial_gauss_fit = gaussian_peak(x_spline, y_spline, gauss_fit_initial_parameters); -final_gauss_fit = gaussian_peak(x_spline, y_spline, parameters_gpufit_gauss); - -%% make a figure of psf, psf template, initial and final gauss fit, initial and final spline fit -figure(1); -min_noisy_psf = min(min(noisy_psf)); -max_noisy_psf = max(max(noisy_psf)); -min_temp = min(min([initial_gauss_fit initial_spline_fit psf_normalized final_gauss_fit final_spline_fit])); -max_temp = max(max([initial_gauss_fit initial_spline_fit psf_normalized final_gauss_fit final_spline_fit])); -min_value = min(min_noisy_psf, min_temp); -max_value = max(max_noisy_psf, max_temp); -clims = [min_value max_value]; -subplot(231); imagesc(x, y, noisy_psf, clims); colorbar; title('noisy PSF'); axis square; -subplot(232); imagesc(x2, y2, initial_gauss_fit, clims); colorbar; title('initial Gaussian fit'); axis square; -subplot(233); imagesc(x2, y2, initial_spline_fit, clims); colorbar; title('initial spline fit'); axis square; -subplot(234); imagesc(x2, y2, psf_normalized); colorbar; title('PSF template'); axis square; -subplot(235); imagesc(x2, y2, final_gauss_fit, clims); colorbar; title('final Gaussian fit'); axis square; -subplot(236); imagesc(x2, y2, final_spline_fit, clims); colorbar; title('final spline fit'); axis square; -colormap('hot'); - -end - -function z = calculate_psf(x, y, p) -% Test PSF consists of an elliptic 2D Gaussian - -% p(1) - amplitude -% p(2) - center x -% p(3) - center y -% p(4) - Standard deviation -% p(5) - constant background -assert(nargin == 3); - -sx = p(4) + 0.1; -sy = p(4) - 0.2; - -arg_ex = exp(-1/2*((x-p(2))/sx).^2-1/2*((y-p(3))/sy).^2); - -z = p(1) .* arg_ex + p(5); % scale with amplitude and background - -end - -function f = gaussian_peak(x, y, p) - -% x,y - x and y grid position values -% p(1) - amplitude -% p(2) - center x -% p(3) - center y -% p(4) - Standard deviation -% p(5) - constant background - -assert(nargin == 3); - -f = p(1) * exp(-((x - p(2)).^2 + (y - p(3)).^2) / (2 * p(4)^2)) + p(5); - -end \ No newline at end of file diff --git a/Gpufit/matlab/matlab/splinefit_3d.m b/Gpufit/matlab/matlab/splinefit_3d.m deleted file mode 100644 index cd2916e9..00000000 --- a/Gpufit/matlab/matlab/splinefit_3d.m +++ /dev/null @@ -1,153 +0,0 @@ -function splinefit_3d() -% spline fit 3D -% -% Requires Gpuspline (https://github.com/gpufit/Gpuspline) additionally - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end -if isempty(which('spline_coefficients.m')) - error('Gpuspline library not found in Matlab path.'); -end - -% initialize random number generator -rng(0); - -% data size -size_x = 18; -size_y = 13; -size_z = 100; - -% tolerances -tolerance = 1e-30; -max_n_iterations = 100; -estimator_id = EstimatorID.LSE; - -% derived values -x = single(0 : size_x - 1); -y = single(0 : size_y - 1); -z = single(0 : size_z - 1); - -%% generate PSF -psf_parameters = single([100, (size_x-1)/2, (size_y-1)/2, (size_z-1)/2+1, 1, 10]); -psf = calculate_psf(x, y, z, psf_parameters); -z_slice_index = 61; - -%% add noise -snr = 10; -amplitude = psf_parameters(1); -noise_std_dev = amplitude ./ (snr * log(10.0)); -noise = noise_std_dev * randn(size(psf)); -noisy_psf = psf + noise; - -%% calculate PSF template -psf_template = (psf - psf_parameters(6)) / psf_parameters(1); - -%% calculate spline coefficients of the PSF template -coefficients = spline_coefficients(psf_template); -n_intervals = size(psf_template) - 1; -coefficients = reshape(coefficients, 64, n_intervals(1), n_intervals(2), n_intervals(3)); - -%% set user info -user_info = [... - size_x; size_y; 1;... - n_intervals';... - coefficients(:)]; - -%% true fit parameters -true_fit_parameters(1,:) = psf_parameters(1); % amplitude -true_fit_parameters(2,:) = 0; % center x shift -true_fit_parameters(3,:) = 0; % center y shift -true_fit_parameters(4,:) = 1 - z_slice_index; % center z shift -true_fit_parameters(5,:) = psf_parameters(6); % offset - -%% set initial fit parameters -pos_shift_x = 1.2; -pos_shift_y = -0.3; -pos_shift_z = -20; -amp_shift = -15; -off_shift = 2; - -spline_fit_initial_parameters = true_fit_parameters + [amp_shift; pos_shift_x; pos_shift_y; pos_shift_z; off_shift]; - -%% reshape data -linear_psf = reshape(noisy_psf(:,:,z_slice_index), numel(noisy_psf(:,:,z_slice_index)), 1); - -%% call to gpufit with spline fit -[parameters_gpufit_spline, states_gpufit_spline, chi_squares_gpufit_spline, n_iterations_gpufit_spline, time_gpufit_spline]... - = gpufit(linear_psf, [], ModelID.SPLINE_3D, spline_fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, user_info); -assert(all(states_gpufit_spline == 0)); - -%% get data to plot -% spline with true parameters -a_true = true_fit_parameters(1); -x_true = x-true_fit_parameters(2); -y_true = y-true_fit_parameters(3); -z_true = -true_fit_parameters(4); -b_true = true_fit_parameters(5); -true_spline_fit = a_true * spline_values(coefficients, x_true, y_true, z_true) + b_true; - -% spline with initial fit parameters -a_init = spline_fit_initial_parameters(1); -x_init = x-spline_fit_initial_parameters(2); -y_init = y-spline_fit_initial_parameters(3); -z_init = -spline_fit_initial_parameters(4); -b_init = spline_fit_initial_parameters(5); -initial_spline_fit = a_init * spline_values(coefficients, x_init, y_init, z_init) + b_init; - -% spline with fit parameters -a_fit = parameters_gpufit_spline(1); -x_fit = x-parameters_gpufit_spline(2); -y_fit = y-parameters_gpufit_spline(3); -z_fit = -parameters_gpufit_spline(4); -b_fit = parameters_gpufit_spline(5); -final_spline_gpufit = a_fit * spline_values(coefficients, x_fit, y_fit, z_fit) + b_fit; - -%% figure -figure(2); -current_slice = noisy_psf(:,:,z_slice_index); -min_noisy_psf = min(current_slice(:)); -max_noisy_psf = max(current_slice(:)); -min_temp = min(min([initial_spline_fit final_spline_gpufit])); -max_temp = max(max([initial_spline_fit final_spline_gpufit])); -min_value = min(min_noisy_psf, min_temp); -max_value = max(max_noisy_psf, max_temp); -clims = [min_value max_value]; -subplot(2,2,1); imagesc(x, y', true_spline_fit.', clims); title(sprintf('true spline z=%.2f', true_fit_parameters(4))); axis image; -subplot(2,2,2); imagesc(x, y', current_slice.', clims); title(sprintf('noisy psf z=%.2f', true_fit_parameters(4))); axis image; -subplot(2,2,3); imagesc(x, y', initial_spline_fit.', clims); title(sprintf('initial spline fit z=%.2f', spline_fit_initial_parameters(4))); axis image; -subplot(2,2,4); imagesc(x, y', final_spline_gpufit.', clims); title(sprintf('final spline gpufit z=%.2f', parameters_gpufit_spline(4))); axis image; -colormap('hot'); - -end - -function f = calculate_psf(x, y, z, p) - -size_x = numel(x); -size_y = numel(y); -size_z = numel(z); - -s_max = p(5) * 5; -s_min = p(5) / 5; - -sx = linspace(s_max, s_min, numel(z)); -sy = linspace(s_min, s_max, numel(z)); -sz = p(5) * 10; - -f = zeros(size_x, size_y, size_z, 'single'); - -for x = 0 : size_x-1 - for y = 0 : size_y-1 - for z = 0 : size_z-1 - - arg_x = exp(-1/2*((x+1-p(2))/sx(z+1)).^2); - arg_y = exp(-1/2*((y-p(3))/sy(z+1)).^2); - arg_z = exp(-1/2*((z-p(4))/sz ).^2); - - f(x+1,y+1,z+1) = p(1) * arg_x * arg_y * arg_z + p(6); - - end - end -end - -end diff --git a/Gpufit/matlab/matlab/splinefit_3d_multichannel.m b/Gpufit/matlab/matlab/splinefit_3d_multichannel.m deleted file mode 100644 index 4c564de3..00000000 --- a/Gpufit/matlab/matlab/splinefit_3d_multichannel.m +++ /dev/null @@ -1,164 +0,0 @@ -function splinefit_3d_multichannel() -% Spline fit 3D 4 channel example. Can for example be used in 4Pi-Storm microscopy. -% -% Requires Gpuspline (https://github.com/gpufit/Gpuspline) additionally - -if isempty(which('gpufit.m')) - error('Gpufit library not found in Matlab path.'); -end -if isempty(which('spline_coefficients.m')) - error('Gpuspline library not found in Matlab path.'); -end - -% initialize random number generator -rng(0); - -% data size -size_x = 19; -size_y = 25; -size_z = 50; -n_channels = 4; - -% tolerances -tolerance = 1e-30; -max_n_iterations = 200; -estimator_id = EstimatorID.LSE; - -% derived values -x = single(0 : size_x - 1); -y = single(0 : size_y - 1); -z = single(0 : size_z - 1); - -%% generate PSF -psf_parameters = single([100, (size_x-1)/2+1, (size_y-1)/2-1, (size_z-1)/2, 1, 10]); -psf = calculate_psf(size_x, size_y, size_z, psf_parameters, n_channels); -z_slice_index = 25; - -%% add noise -rng(0); -snr = 5; -amplitude = psf_parameters(1); -noise_std_dev = amplitude ./ (snr * log(10.0)); -noise = noise_std_dev * randn(size(psf)); -noisy_psf = psf + noise; - -%% calculate PSF template -psf_normalized = (psf - psf_parameters(6)) / psf_parameters(1); - -%% calculate spline coefficients of the PSF template -coefficients = zeros([64 size(psf(:,:,:,1))-1 n_channels],'single'); -for ch=1:n_channels; coefficients(:,:,:,:,ch) = spline_coefficients(psf_normalized(:,:,:,ch)); end -n_intervals = size(psf_normalized) - 1; -n_intervals = n_intervals(1:3); - -%% set user info -user_info = [... - n_channels,... - size_x, size_y, 1,... - n_intervals,... - reshape(coefficients, 1, numel(coefficients))]; - -%% true fit parameters -true_fit_parameters(1,:) = psf_parameters(1); % amplitude -true_fit_parameters(2,:) = 0; % center x shift -true_fit_parameters(3,:) = 0; % center y shift -true_fit_parameters(4,:) = 1 - z_slice_index; % z index -true_fit_parameters(5,:) = psf_parameters(6); % offset - -%% set initial fit parameters -pos_shift_x = -0.9; -pos_shift_y = 0.7; -pos_shift_z = 10; -amp_shift = -20; -off_shift = 5; - -spline_fit_initial_parameters = true_fit_parameters + [amp_shift; pos_shift_x; pos_shift_y; pos_shift_z; off_shift]; - -%% reshape data -linear_psf = reshape(noisy_psf(:,:,z_slice_index,:), size_x * size_y * n_channels, 1); - -%% call to gpufit with spline fit -[parameters_gpufit_spline, states_gpufit_spline, chi_squares_gpufit_spline, n_iterations_gpufit_spline, time_gpufit_spline]... - = gpufit(linear_psf, [], ModelID.SPLINE_3D_MULTICHANNEL, spline_fit_initial_parameters, tolerance, max_n_iterations, [], estimator_id, user_info); - -% check if gpufit succeeded -assert(all(states_gpufit_spline == 0)); - -%% get data to plot -merge_channels = @(x) [x(:, :, 1), x(:, :, 2); x(:, :, 3), x(:, :, 4)]; - -psf_4ch = merge_channels(psf(:,:,z_slice_index,:)); - -noisy_psf_4ch = merge_channels(noisy_psf(:,:,z_slice_index,:)); - -initial_spline_fit = spline_values_3d_multichannel(coefficients, n_intervals, n_channels, spline_fit_initial_parameters); -initial_spline_fit_4ch = merge_channels(initial_spline_fit); - -final_spline_gpufit = spline_values_3d_multichannel(coefficients, n_intervals, n_channels, parameters_gpufit_spline); -final_spline_gpufit_4ch = merge_channels(final_spline_gpufit); - -%% plot -figure(5); -min_value = min([psf(:); noisy_psf(:); initial_spline_fit(:); final_spline_gpufit(:)]); -max_value = max([psf(:); noisy_psf(:); initial_spline_fit(:); final_spline_gpufit(:)]); -clims = [min_value max_value]; - -subplot(2,2,1); imagesc(x, y', psf_4ch, clims); title(sprintf('psf z=%.2f', true_fit_parameters(4))); axis image; -subplot(2,2,2); imagesc(x, y', noisy_psf_4ch, clims); title(sprintf('noisy psf z=%.2f', true_fit_parameters(4))); axis image; -subplot(2,2,3); imagesc(x, y', initial_spline_fit_4ch, clims); title(sprintf('initial spline fit z=%.2f', spline_fit_initial_parameters(4))); axis image; -subplot(2,2,4); imagesc(x, y', final_spline_gpufit_4ch, clims); title(sprintf('final gpufit z=%.2f', parameters_gpufit_spline(4))); axis image; -colormap('hot'); - -end - - -%% calculate psf -function f = calculate_psf(size_x, size_y, size_z, p, n_channels) - -s_max = p(5) * 5; -s_min = p(5) / 5; - -sx = linspace(s_max, s_min, size_z); -sy = linspace(s_min, s_max, size_z); -sz = 8; - -delta_s = sx(1) - sx(2); -sx = [sx ,s_min - 1*delta_s, s_min - 2*delta_s, s_min - 3*delta_s]; -sy = [sy ,s_max + 1*delta_s, s_max + 2*delta_s, s_max + 3*delta_s]; - -x = linspace(0,pi,size_z); -a = zeros(numel(x), 4); -a(:,1) = sin(x + pi * 0/4) * 0.5 + 0.5; -a(:,2) = sin(x + pi * 1/4) * 0.5 + 0.5; -a(:,3) = sin(x + pi * 2/4) * 0.5 + 0.5; -a(:,4) = sin(x + pi * 3/4) * 0.5 + 0.5; - -f = zeros(size_x, size_y, size_z, n_channels, 'single'); - -for ch = 0 : n_channels-1 - for xi = 0 : size_x-1 - for yi = 0 : size_y-1 - for zi = 0 : size_z-1 - - arg_x = exp(-1/2*((xi-p(2))/sx(zi+1+ch)).^2); - arg_y = exp(-1/2*((yi-p(3))/sy(zi+1+ch)).^2); - arg_z = exp(-1/2*((zi+ch-p(4))/sz ).^2); - - f(xi+1,yi+1,zi+1,ch+1) = a(zi+1,ch+1) * p(1) * arg_x * arg_y * arg_z + p(6); - - end - end - end -end - -end - -function f = spline_values_3d_multichannel(coefficients, n_intervals, n_channels, p) - -coefficients = reshape(coefficients, 64, n_intervals(1), n_intervals(2), n_intervals(3), n_channels); -x = (0:n_intervals(1))-p(2); -y = (0:n_intervals(2))-p(3); -z = -p(4); -f = p(1) * spline_values(coefficients, x, y, z) + p(5); - -end \ No newline at end of file