From 0847a315ce744e52be3dade398fb16c58323084e Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 18 Oct 2017 13:44:57 +0100 Subject: linked demos to TomoPhantom, cone beam demo, some FISTA modificastions --- demos/Demo1.m | 174 ---------------------------------------- demos/Demo_Phantom3D_Cone.m | 66 +++++++++++++++ demos/Demo_Phantom3D_Parallel.m | 49 +++++++++++ main_func/FISTA_REC.m | 155 +++++++++++++++-------------------- 4 files changed, 179 insertions(+), 265 deletions(-) delete mode 100644 demos/Demo1.m create mode 100644 demos/Demo_Phantom3D_Cone.m create mode 100644 demos/Demo_Phantom3D_Parallel.m diff --git a/demos/Demo1.m b/demos/Demo1.m deleted file mode 100644 index 15e2e5b..0000000 --- a/demos/Demo1.m +++ /dev/null @@ -1,174 +0,0 @@ -% Demonstration of tomographic reconstruction from noisy and corrupted by -% artifacts undersampled projection data using Students't penalty -% Optimisation problem is solved using FISTA algorithm (see Beck & Teboulle) - -% see Readme file for instructions -%% -% compile MEX-files ones -% cd .. -% cd main_func -% compile_mex -% cd .. -% cd demos -%% - -close all;clc;clear all; -% adding paths -addpath('../data/'); -addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); -addpath('../supp/'); - -load phantom_bone512.mat % load the phantom -load my_red_yellowMAP.mat % load the colormap -% load sino1.mat; % load noisy sinogram - -N = 512; % the size of the tomographic image NxN -theta = 1:1:180; % acquisition angles (in parallel beam from 0 to Pi) -theta_rad = theta*(pi/180); % conversion to radians -P = 2*ceil(N/sqrt(2))+1; % the size of the detector array -ROI = find(phantom > 0); - -% using ASTRA to set the projection geometry -% potentially parallel geometry can be replaced with a divergent one -Z_slices = 1; -det_row_count = Z_slices; -proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, P, theta_rad); -vol_geom = astra_create_vol_geom(N,N,Z_slices); - -zing_rings_add; % generating data, adding zingers and stripes -%% -fprintf('%s\n', 'Direct reconstruction using FBP...'); -FBP_1 = iradon(sino_zing_rings', theta, N); - -fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:))); - -figure(1); -subplot_tight(1,2,1, [0.05 0.05]); imshow(FBP_1,[0 0.6]); title('FBP reconstruction of noisy and corrupted by artifacts sinogram'); colorbar; -subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - FBP_1).^2,[0 0.1]); title('residual: (ideal phantom - FBP)^2'); colorbar; -colormap(cmapnew); - -%% -fprintf('%s\n', 'Reconstruction using FISTA-PWLS without regularization...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; % sinogram -params.iterFISTA = 45; %max number of outer iterations -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.show = 1; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -params.weights = Dweights; % statistical weighting -tic; [X_FISTA, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:))); -error_FISTA = output.Resid_error; obj_FISTA = output.objective; - -figure(2); clf -%set(gcf, 'Position', get(0,'Screensize')); -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); -figure(3); clf -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; -params.iterFISTA = 45; % max number of outer iterations -params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.show = 1; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_TV, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction is:', min(error_FISTA_TV(:))); -error_FISTA_TV = output.Resid_error; obj_FISTA_TV = output.objective; - -figure(4); clf -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_TV).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); -figure(5); clf -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_TV); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_TV); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; -params.iterFISTA = 50; % max number of outer iterations -params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.Ring_LambdaR_L1 = 0.002; % parameter to sparsify the "rings vector" -params.Ring_Alpha = 20; % to accelerate ring-removal procedure -params.show = 0; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_GH_TV, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction is:', min(error_FISTA_GH_TV(:))); -error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective; - -figure(6); clf -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]);imshow((phantom - X_FISTA_GH_TV).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); - -figure(7); clf -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_GH_TV); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_GH_TV); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; -params.iterFISTA = 55; % max number of outer iterations -params.L_const = 0.1; % Lipshitz constant (can be chosen manually to accelerate convergence) -params.Regul_LambdaTV = 0.00152; % regularization parameter for TV problem -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.fidelity = 'student'; % selecting students t fidelity -params.show = 1; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -params.initilize = 1; % warm start with SIRT -tic; [X_FISTA_student_TV, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction is:', min(error_FISTA_student_TV(:))); -error_FISTA_student_TV = output.Resid_error; obj_FISTA_student_TV = output.objective; - -figure(8); -set(gcf, 'Position', get(0,'Screensize')); -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_student_TV).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); - -figure(9); -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_student_TV); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -% print all RMSE's -fprintf('%s\n', '--------------------------------------------'); -fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction:', min(error_FISTA(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction:', min(error_FISTA_TV(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction:', min(error_FISTA_GH_TV(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction:', min(error_FISTA_student_TV(:))); -% \ No newline at end of file diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m new file mode 100644 index 0000000..6419386 --- /dev/null +++ b/demos/Demo_Phantom3D_Cone.m @@ -0,0 +1,66 @@ +% A demo script to reconstruct 3D synthetic data using FISTA method for +% CONE BEAM geometry +% requirements: ASTRA-toolbox and TomoPhantom toolbox + +close all;clc;clear all; +% adding paths +addpath('../data/'); +addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('../main_func/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/'); +addpath('../supp/'); + + +%% +% build 3D phantom using TomoPhantom +modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom +N = 256; % x-y-z size (cubic image) +angles = 0:1.5:360; % angles vector in degrees +angles_rad = angles*(pi/180); % conversion to radians +det_size = round(sqrt(2)*N); % detector size +% in order to run functions you have to go to the directory: +cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ +TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom +%% +% using ASTRA-toolbox to set the projection geometry (cone beam) +% eg: astra.create_proj_geom('cone', 1.0 (resol), 1.0 (resol), detectorRowCount, detectorColCount, angles, originToSource, originToDetector) +vol_geom = astra_create_vol_geom(N,N,N); +proj_geom = astra_create_proj_geom('cone', 1.0, 1.0, N, det_size, angles_rad, 2000, 2160); +%% +% do forward projection using ASTRA +% inverse crime data generation +[sino_id, SinoCone3D] = astra_create_sino3d_cuda(TomoPhantom, proj_geom, vol_geom); +astra_mex_data3d('delete', sino_id); +%% +fprintf('%s\n', 'Reconstructing with CGLS using ASTRA-toolbox ...'); +vol_id = astra_mex_data3d('create', '-vol', vol_geom, 0); +proj_id = astra_mex_data3d('create', '-proj3d', proj_geom, SinoCone3D); +cfg = astra_struct('CGLS3D_CUDA'); +cfg.ProjectionDataId = proj_id; +cfg.ReconstructionDataId = vol_id; +cfg.option.MinConstraint = 0; +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('iterate', alg_id, 15); +reconASTRA_3D = astra_mex_data3d('get', vol_id); +%% +fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); +clear params +% define parameters +params.proj_geom = proj_geom; % pass geometry to the function +params.vol_geom = vol_geom; +params.sino = single(SinoCone3D); % sinogram +params.iterFISTA = 30; %max number of outer iterations +params.X_ideal = TomoPhantom; % ideal phantom +params.show = 1; % visualize reconstruction on each iteration +params.slice = round(N/2); params.maxvalplot = 1; +tic; [X_FISTA, output] = FISTA_REC(params); toc; + +error_FISTA = output.Resid_error; obj_FISTA = output.objective; +fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction is:', min(error_FISTA(:))); + +Resid3D = (TomoPhantom - X_FISTA).^2; +figure(2); +subplot(1,2,1); imshow(X_FISTA(:,:,params.slice),[0 params.maxvalplot]); title('FISTA-LS reconstruction'); colorbar; +subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]); title('residual'); colorbar; +figure(3); +subplot(1,2,1); plot(error_FISTA); title('RMSE plot'); colorbar; +subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); colorbar; +%% \ No newline at end of file diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m new file mode 100644 index 0000000..fd8096a --- /dev/null +++ b/demos/Demo_Phantom3D_Parallel.m @@ -0,0 +1,49 @@ +% A demo script to reconstruct 3D synthetic data using FISTA method for +% PARALLEL BEAM geometry +% requirements: ASTRA-toolbox and TomoPhantom toolbox + +close all;clc;clear all; +% adding paths +addpath('../data/'); +addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('../main_func/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/'); +addpath('../supp/'); + +%% +% build 3D phantom using TomoPhantom and generate projection data +modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom +N = 256; % x-y-z size (cubic image) +angles = 1:0.5:180; % angles vector in degrees +angles_rad = angles*(pi/180); % conversion to radians +det_size = round(sqrt(2)*N); % detector size +% in order to run functions you have to go to the directory: +cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ +TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom +sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles)); % generate data +%% +% using ASTRA-toolbox to set the projection geometry (parallel beam) +proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad); +vol_geom = astra_create_vol_geom(N,N); +%% +fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); +clear params +% define parameters +params.proj_geom = proj_geom; % pass geometry to the function +params.vol_geom = vol_geom; +params.sino = single(sino_tomophan3D); % sinogram +params.iterFISTA = 5; %max number of outer iterations +params.X_ideal = TomoPhantom; % ideal phantom +params.show = 1; % visualize reconstruction on each iteration +params.slice = round(N/2); params.maxvalplot = 1; +tic; [X_FISTA, output] = FISTA_REC(params); toc; + +error_FISTA = output.Resid_error; obj_FISTA = output.objective; +fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:))); + +Resid3D = (TomoPhantom - X_FISTA).^2; +figure(2); +subplot(1,2,1); imshow(X_FISTA(:,:,params.slice),[0 params.maxvalplot]); title('FISTA-LS reconstruction'); colorbar; +subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]); title('residual'); colorbar; +figure(3); +subplot(1,2,1); plot(error_FISTA); title('RMSE plot'); colorbar; +subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); colorbar; +%% \ No newline at end of file diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index 6987dca..dde0e73 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -15,7 +15,7 @@ function [X, output] = FISTA_REC(params) %----------------General Parameters------------------------ % - .proj_geom (geometry of the projector) [required] % - .vol_geom (geometry of the reconstructed object) [required] -% - .sino (vectorized in 2D or 3D sinogram) [required] +% - .sino (2D or 3D sinogram) [required] % - .iterFISTA (iterations for the main loop, default 40) % - .L_const (Lipschitz constant, default Power method) ) % - .X_ideal (ideal image, if given) @@ -94,45 +94,36 @@ else weights = ones(size(sino)); end if (isfield(params,'fidelity')) - studentt = 0; + studentt = 0; if (strcmp(params.fidelity,'studentt') == 1) - studentt = 1; - lambdaR_L1 = 0; - end + studentt = 1; + end else - studentt = 0; + studentt = 0; end if (isfield(params,'L_const')) L_const = params.L_const; else % using Power method (PM) to establish L constant - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'parallel3d')) - % for parallel geometry we can do just one slice - fprintf('%s \n', 'Calculating Lipshitz constant for parallel beam geometry...'); + fprintf('%s %s %s \n', 'Calculating Lipshitz constant for',proj_geom.type, 'beam geometry...'); + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % for 2D geometry we can do just one selected slice niter = 15; % number of iteration for the PM x1 = rand(N,N,1); sqweight = sqrt(weights(:,:,1)); - proj_geomT = proj_geom; - proj_geomT.DetectorRowCount = 1; - vol_geomT = vol_geom; - vol_geomT.GridSliceCount = 1; - [sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT); - y = sqweight.*y; - astra_mex_data3d('delete', sino_id); - + [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); + y = sqweight.*y'; + astra_mex_data2d('delete', sino_id); for i = 1:niter - [id,x1] = astra_create_backprojection3d_cuda(sqweight.*y, proj_geomT, vol_geomT); + [x1] = astra_create_backprojection_cuda((sqweight.*y)', proj_geom, vol_geom); s = norm(x1(:)); x1 = x1./s; - [sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT); - y = sqweight.*y; - astra_mex_data3d('delete', sino_id); - astra_mex_data3d('delete', id); - end - %clear proj_geomT vol_geomT - else - % divergen beam geometry - fprintf('%s \n', 'Calculating Lipshitz constant for divergen beam geometry... will take some time!'); + [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); + y = sqweight.*y'; + astra_mex_data2d('delete', sino_id); + end + elseif (strcmp(proj_geom.type,'cone') || strcmp(proj_geom.type,'parallel3d') || strcmp(proj_geom.type,'parallel3d_vec') || strcmp(proj_geom.type,'cone_vec')) + % 3D geometry niter = 8; % number of iteration for PM x1 = rand(N,N,SlicesZ); sqweight = sqrt(weights); @@ -150,6 +141,8 @@ else astra_mex_data3d('delete', id); end clear x1 + else + error('%s \n', 'No suitable geometry has been found!'); end L_const = s; end @@ -272,28 +265,10 @@ else slice = 1; end if (isfield(params,'initialize')) - % a 'warm start' with SIRT method - % Create a data object for the reconstruction - rec_id = astra_mex_data3d('create', '-vol', vol_geom); - - sinogram_id = astra_mex_data3d('create', '-proj3d', proj_geom, sino); - - % Set up the parameters for a reconstruction algorithm using the GPU - cfg = astra_struct('SIRT3D_CUDA'); - cfg.ReconstructionDataId = rec_id; - cfg.ProjectionDataId = sinogram_id; - - % Create the algorithm object from the configuration structure - alg_id = astra_mex_algorithm('create', cfg); - astra_mex_algorithm('iterate', alg_id, 35); - % Get the result - X = astra_mex_data3d('get', rec_id); - - % Clean up. Note that GPU memory is tied up in the algorithm object, - % and main RAM in the data objects. - astra_mex_algorithm('delete', alg_id); - astra_mex_data3d('delete', rec_id); - astra_mex_data3d('delete', sinogram_id); + X = params.initialize; + if ((size(X,1) ~= N) || (size(X,2) ~= N) || (size(X,3) ~= SlicesZ)) + error('%s \n', 'The initialized volume has different dimensions!'); + end else X = zeros(N,N,SlicesZ, 'single'); % storage for the solution end @@ -345,16 +320,18 @@ if (subsets == 0) t_old = t; r_old = r; - % if the geometry is parallel use slice-by-slice projection-backprojection routine - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'parallel3d')) + % if the geometry is 2D use slice-by-slice projection-backprojection routine + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) sino_updt = zeros(size(sino),'single'); - for kkk = 1:SlicesZ - [sino_id, sino_updt(:,:,kkk)] = astra_create_sino3d_cuda(X_t(:,:,kkk), proj_geomT, vol_geomT); - astra_mex_data3d('delete', sino_id); + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geom, vol_geom); + sino_updt(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); end else - % for divergent 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom); + astra_mex_data3d('delete', sino_id); end if (lambdaR_L1 > 0) @@ -368,40 +345,37 @@ if (subsets == 0) end r = r_x - (1./L_const).*vec; objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output - else - if (studentt == 1) - % artifacts removal with Students t penalty - residual = weights.*(sino_updt - sino); - for kkk = 1:SlicesZ - res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram + elseif (studentt > 0) + % artifacts removal with Students t penalty + residual = weights.*(sino_updt - sino); + for kkk = 1:SlicesZ + res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram %s = 100; %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); [ff, gr] = studentst(res_vec, 1); residual(:,:,kkk) = reshape(gr, Detectors, anglesNumb); - end - objective(i) = ff; % for the objective function output - else + end + objective(i) = ff; % for the objective function output + else % no ring removal (LS model) residual = weights.*(sino_updt - sino); objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output - end - end + end + - % if the geometry is parallel use slice-by-slice projection-backprojection routine - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'parallel3d')) + % if the geometry is 2D use slice-by-slice projection-backprojection routine + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) x_temp = zeros(size(X),'single'); for kkk = 1:SlicesZ - [id, x_temp(:,:,kkk)] = astra_create_backprojection3d_cuda(squeeze(residual(:,:,kkk)), proj_geomT, vol_geomT); - astra_mex_data3d('delete', id); + [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residual(:,:,kkk))', proj_geom, vol_geom); end else [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom); + astra_mex_data3d('delete', id); end - X = X_t - (1/L_const).*x_temp; - astra_mex_data3d('delete', sino_id); - astra_mex_data3d('delete', id); + X = X_t - (1/L_const).*x_temp; - % regularization + % ----------------Regularization part------------------------ if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -572,28 +546,27 @@ else end r = r_x - (1./L_const).*vec; else - [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); if (studentt == 1) - % artifacts removal with Students t penalty - residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); - - for kkk = 1:SlicesZ - res_vec = reshape(residual(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram - %s = 100; - %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); - [ff, gr] = studentst(res_vec, 1); - residual(:,:,kkk) = reshape(gr, Detectors, numProjSub); - end - objective(i) = ff; % for the objective function output - else - % no ring removal (LS model) - residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); - end + % artifacts removal with Students t penalty + residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); + + for kkk = 1:SlicesZ + res_vec = reshape(residual(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram + %s = 100; + %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); + [ff, gr] = studentst(res_vec, 1); + residual(:,:,kkk) = reshape(gr, Detectors, numProjSub); + end + objective(i) = ff; % for the objective function output + else + % no ring removal (LS model) + residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); + end end [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geomSUB, vol_geom); - X = X_t - (1/L_const).*x_temp; astra_mex_data3d('delete', sino_id); astra_mex_data3d('delete', id); -- cgit v1.2.3 From cb8ef11f00e897b6f5b3049126dd32baa3c50cf9 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 18 Oct 2017 21:46:29 +0100 Subject: ordered subsets fix for GH term --- demos/Demo_Phantom3D_Parallel.m | 5 +- main_func/FISTA_REC.m | 179 +++++++++++++++++++++++----------------- 2 files changed, 108 insertions(+), 76 deletions(-) diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index fd8096a..ac9827c 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -33,6 +33,7 @@ params.sino = single(sino_tomophan3D); % sinogram params.iterFISTA = 5; %max number of outer iterations params.X_ideal = TomoPhantom; % ideal phantom params.show = 1; % visualize reconstruction on each iteration +params.subsets = 12; params.slice = round(N/2); params.maxvalplot = 1; tic; [X_FISTA, output] = FISTA_REC(params); toc; @@ -44,6 +45,6 @@ figure(2); subplot(1,2,1); imshow(X_FISTA(:,:,params.slice),[0 params.maxvalplot]); title('FISTA-LS reconstruction'); colorbar; subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]); title('residual'); colorbar; figure(3); -subplot(1,2,1); plot(error_FISTA); title('RMSE plot'); colorbar; -subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); colorbar; +subplot(1,2,1); plot(error_FISTA); title('RMSE plot'); +subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); %% \ No newline at end of file diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index dde0e73..bea1860 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -106,14 +106,14 @@ if (isfield(params,'L_const')) else % using Power method (PM) to establish L constant fprintf('%s %s %s \n', 'Calculating Lipshitz constant for',proj_geom.type, 'beam geometry...'); - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % for 2D geometry we can do just one selected slice + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % for 2D geometry we can do just one selected slice niter = 15; % number of iteration for the PM x1 = rand(N,N,1); sqweight = sqrt(weights(:,:,1)); [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); y = sqweight.*y'; - astra_mex_data2d('delete', sino_id); + astra_mex_data2d('delete', sino_id); for i = 1:niter [x1] = astra_create_backprojection_cuda((sqweight.*y)', proj_geom, vol_geom); s = norm(x1(:)); @@ -121,9 +121,9 @@ else [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); y = sqweight.*y'; astra_mex_data2d('delete', sino_id); - end + end elseif (strcmp(proj_geom.type,'cone') || strcmp(proj_geom.type,'parallel3d') || strcmp(proj_geom.type,'parallel3d_vec') || strcmp(proj_geom.type,'cone_vec')) - % 3D geometry + % 3D geometry niter = 8; % number of iteration for PM x1 = rand(N,N,SlicesZ); sqweight = sqrt(weights); @@ -268,7 +268,7 @@ if (isfield(params,'initialize')) X = params.initialize; if ((size(X,1) ~= N) || (size(X,2) ~= N) || (size(X,3) ~= SlicesZ)) error('%s \n', 'The initialized volume has different dimensions!'); - end + end else X = zeros(N,N,SlicesZ, 'single'); % storage for the solution end @@ -320,10 +320,11 @@ if (subsets == 0) t_old = t; r_old = r; - % if the geometry is 2D use slice-by-slice projection-backprojection routine - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine sino_updt = zeros(size(sino),'single'); - for kkk = 1:SlicesZ + for kkk = 1:SlicesZ [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geom, vol_geom); sino_updt(:,:,kkk) = sinoT'; astra_mex_data2d('delete', sino_id); @@ -359,12 +360,11 @@ if (subsets == 0) else % no ring removal (LS model) residual = weights.*(sino_updt - sino); - objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output + objective(i) = 0.5*norm(residual(:)); % for the objective function output end - % if the geometry is 2D use slice-by-slice projection-backprojection routine - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) x_temp = zeros(size(X),'single'); for kkk = 1:SlicesZ [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residual(:,:,kkk))', proj_geom, vol_geom); @@ -373,9 +373,9 @@ if (subsets == 0) [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom); astra_mex_data3d('delete', id); end - X = X_t - (1/L_const).*x_temp; + X = X_t - (1/L_const).*x_temp; - % ----------------Regularization part------------------------ + % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -484,94 +484,128 @@ if (subsets == 0) end end else - % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than classical) + % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than the classical version) t = 1; X_t = X; - proj_geomSUB = proj_geom; - + proj_geomSUB = proj_geom; r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors r_x = r; % another ring variable residual2 = zeros(size(sino),'single'); + sino_updt_FULL = zeros(size(sino),'single'); % Outer FISTA iterations loop - for i = 1:iterFISTA + for i = 1:iterFISTA - % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required - % one solution is to work with a full sinogram at times - if ((i >= 3) && (lambdaR_L1 > 0)) - [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom); - astra_mex_data3d('delete', sino_id2); + if ((i > 1) && (lambdaR_L1 > 0)) + % in order to make Group-Huber fidelity work with ordered subsets + % we still need to work with full sinogram + + % the offset variable must be calculated for the whole + % updated sinogram - sino_updt_FULL + for kkk = 1:anglesNumb + residual2(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt_FULL(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x)); + end + + r_old = r; + vec = sum(residual2,2); + if (SlicesZ > 1) + vec = squeeze(vec(:,1,:)); + end + r = r_x - (1./L_const).*vec; % update ring variable end % subsets loop counterInd = 1; for ss = 1:subsets X_old = X; - t_old = t; - r_old = r; + t_old = t; numProjSub = binsDiscr(ss); % the number of projections per subset CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces); + sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single'); if (lambdaR_L1 > 0) - - % the ring removal part (Group-Huber fidelity) - % first 2 iterations do additional work reconstructing whole dataset to ensure - % the stablility - if (i < 3) - [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom); - astra_mex_data3d('delete', sino_id2); + % Group-Huber fidelity (ring removal) + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); + end else - [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - end - - for kkk = 1:anglesNumb - residual2(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt2(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x)); + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); end - residual = zeros(Detectors, numProjSub, SlicesZ,'single'); + residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub indC = CurrSubIndeces(kkk); - if (i < 3) - residual(:,kkk,:) = squeeze(residual2(:,indC,:)); - else - residual(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); + residualSub(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); + sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram + end + + elseif (studentt > 0) + % student t data fidelity + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); end + else + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); end - vec = sum(residual2,2); - if (SlicesZ > 1) - vec = squeeze(vec(:,1,:)); + + % artifacts removal with Students t penalty + residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); + + for kkk = 1:SlicesZ + res_vec = reshape(residualSub(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram + %s = 100; + %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); + [ff, gr] = studentst(res_vec, 1); + residualSub(:,:,kkk) = reshape(gr, Detectors, numProjSub); end - r = r_x - (1./L_const).*vec; + objective(i) = ff; % for the objective function output else - [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - - if (studentt == 1) - % artifacts removal with Students t penalty - residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); - + % PWLS model + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine for kkk = 1:SlicesZ - res_vec = reshape(residual(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram - %s = 100; - %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); - [ff, gr] = studentst(res_vec, 1); - residual(:,:,kkk) = reshape(gr, Detectors, numProjSub); + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); end - objective(i) = ff; % for the objective function output else - % no ring removal (LS model) - residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); end + residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); + objective(i) = 0.5*norm(residualSub(:)); % for the objective function output end - [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geomSUB, vol_geom); - X = X_t - (1/L_const).*x_temp; - astra_mex_data3d('delete', sino_id); - astra_mex_data3d('delete', id); + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + x_temp = zeros(size(X),'single'); + for kkk = 1:SlicesZ + [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residualSub(:,:,kkk))', proj_geomSUB, vol_geom); + end + else + [id, x_temp] = astra_create_backprojection3d_cuda(residualSub, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', id); + end + + X = X_t - (1/L_const).*x_temp; - % regularization + % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -653,20 +687,17 @@ else end end - if (lambdaR_L1 > 0) - r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector - end - t = (1 + sqrt(1 + 4*t^2))/2; % updating t X_t = X + ((t_old-1)/t).*(X - X_old); % updating X - - if (lambdaR_L1 > 0) - r_x = r + ((t_old-1)/t).*(r - r_old); % updating r - end - counterInd = counterInd + numProjSub; end + % working with a 'ring vector' + if (lambdaR_L1 > 0) + r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector + r_x = r + ((t_old-1)/t).*(r - r_old); % updating r + end + if (show == 1) figure(10); imshow(X(:,:,slice), [0 maxvalplot]); if (lambdaR_L1 > 0) -- cgit v1.2.3 From dd6e415991f312bf54cdd69d1dd09fb8bcdebd2a Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 19 Oct 2017 12:38:58 +0100 Subject: moved code out of if-the-else closes issue #5 --- main_func/FISTA_REC.m | 54 ++++++++++++++------------------------------------- 1 file changed, 15 insertions(+), 39 deletions(-) diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index bea1860..d177748 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -517,6 +517,18 @@ else % subsets loop counterInd = 1; + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); + end + else + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); + end for ss = 1:subsets X_old = X; t_old = t; @@ -528,18 +540,7 @@ else if (lambdaR_L1 > 0) % Group-Huber fidelity (ring removal) - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end + residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub @@ -550,18 +551,7 @@ else elseif (studentt > 0) % student t data fidelity - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end + % artifacts removal with Students t penalty residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); @@ -576,21 +566,7 @@ else objective(i) = ff; % for the objective function output else % PWLS model - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end - residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); - objective(i) = 0.5*norm(residualSub(:)); % for the objective function output - end + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) % if geometry is 2D use slice-by-slice projection-backprojection routine -- cgit v1.2.3 From 99cfddc7764a292c438868a02c9d512dabefbdc4 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 19 Oct 2017 12:47:02 +0100 Subject: moved code out of if-then-else fixed #5 --- main_func/FISTA_REC.m | 54 ++++++++++++++------------------------------------- 1 file changed, 15 insertions(+), 39 deletions(-) diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index bea1860..d177748 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -517,6 +517,18 @@ else % subsets loop counterInd = 1; + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); + end + else + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); + end for ss = 1:subsets X_old = X; t_old = t; @@ -528,18 +540,7 @@ else if (lambdaR_L1 > 0) % Group-Huber fidelity (ring removal) - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end + residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub @@ -550,18 +551,7 @@ else elseif (studentt > 0) % student t data fidelity - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end + % artifacts removal with Students t penalty residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); @@ -576,21 +566,7 @@ else objective(i) = ff; % for the objective function output else % PWLS model - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end - residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); - objective(i) = 0.5*norm(residualSub(:)); % for the objective function output - end + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) % if geometry is 2D use slice-by-slice projection-backprojection routine -- cgit v1.2.3 From ff543b21e175d49b95ea6da9c03f21a1789f6833 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 19 Oct 2017 13:08:13 +0100 Subject: OS loop now fixed --- demos/Demo_Phantom3D_Cone.m | 4 ++-- demos/Demo_Phantom3D_Parallel.m | 24 +++++++++++++++++++---- main_func/FISTA_REC.m | 43 +++++++++++++++++++++-------------------- 3 files changed, 44 insertions(+), 27 deletions(-) diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m index 6419386..3a9178b 100644 --- a/demos/Demo_Phantom3D_Cone.m +++ b/demos/Demo_Phantom3D_Cone.m @@ -17,8 +17,8 @@ angles = 0:1.5:360; % angles vector in degrees angles_rad = angles*(pi/180); % conversion to radians det_size = round(sqrt(2)*N); % detector size % in order to run functions you have to go to the directory: -cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ -TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom +pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file +TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % generate 3D phantom %% % using ASTRA-toolbox to set the projection geometry (cone beam) % eg: astra.create_proj_geom('cone', 1.0 (resol), 1.0 (resol), detectorRowCount, detectorColCount, angles, originToSource, originToDetector) diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index ac9827c..6a54450 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -10,19 +10,35 @@ addpath('../supp/'); %% % build 3D phantom using TomoPhantom and generate projection data -modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom +modelNo = 2; % see Phantom3DLibrary.dat file in TomoPhantom N = 256; % x-y-z size (cubic image) angles = 1:0.5:180; % angles vector in degrees angles_rad = angles*(pi/180); % conversion to radians det_size = round(sqrt(2)*N); % detector size % in order to run functions you have to go to the directory: -cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ -TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom -sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles)); % generate data +pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file +TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % generate 3D phantom +sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles),pathTP); % generate ideal data +% Adding noise and distortions if required +sino_artifacts = sino_add_artifacts(sino_tomophan3D,'rings'); +% %% % using ASTRA-toolbox to set the projection geometry (parallel beam) proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad); vol_geom = astra_create_vol_geom(N,N); +%% +fprintf('%s\n', 'Reconstructing with FBP using ASTRA-toolbox ...'); +for i = 1:k +vol_id = astra_mex_data2d('create', '-vol', vol_geom, 0); +proj_id = astra_mex_data2d('create', '-proj3d', proj_geom, sino_artifacts(:,:,k)); +cfg = astra_struct('FBP_CUDA'); +cfg.ProjectionDataId = proj_id; +cfg.ReconstructionDataId = vol_id; +cfg.option.MinConstraint = 0; +alg_id = astra_mex_algorithm('create', cfg); +astra_mex_algorithm('iterate', alg_id, 15); +reconASTRA_3D = astra_mex_data2d('get', vol_id); +end %% fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); clear params diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index d177748..bdaeb18 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -487,7 +487,7 @@ else % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than the classical version) t = 1; X_t = X; - proj_geomSUB = proj_geom; + proj_geomSUB = proj_geom; r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors r_x = r; % another ring variable @@ -495,7 +495,7 @@ else sino_updt_FULL = zeros(size(sino),'single'); % Outer FISTA iterations loop - for i = 1:iterFISTA + for i = 1:iterFISTA if ((i > 1) && (lambdaR_L1 > 0)) % in order to make Group-Huber fidelity work with ordered subsets @@ -517,31 +517,31 @@ else % subsets loop counterInd = 1; - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end for ss = 1:subsets X_old = X; - t_old = t; + t_old = t; numProjSub = binsDiscr(ss); % the number of projections per subset + sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single'); CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces); - sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single'); + + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); + end + else + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); + end if (lambdaR_L1 > 0) % Group-Huber fidelity (ring removal) - residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub indC = CurrSubIndeces(kkk); @@ -551,8 +551,6 @@ else elseif (studentt > 0) % student t data fidelity - - % artifacts removal with Students t penalty residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); @@ -566,8 +564,11 @@ else objective(i) = ff; % for the objective function output else % PWLS model - + residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); + objective(i) = 0.5*norm(residualSub(:)); % for the objective function output + end + % perform backprojection of a subset if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) % if geometry is 2D use slice-by-slice projection-backprojection routine x_temp = zeros(size(X),'single'); @@ -579,7 +580,7 @@ else astra_mex_data3d('delete', id); end - X = X_t - (1/L_const).*x_temp; + X = X_t - (1/L_const).*x_temp; % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) -- cgit v1.2.3 From fe4e60b66c1aeccb373b013bfd8396ce413d9d3e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 19 Oct 2017 17:01:55 +0100 Subject: First commit of CMakeLists.txt attempt to use CMake to create the build environment --- CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..7eb7185 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,4 @@ +# +cmake_minimum_required(3.0.2 FATAL_ERROR) + +project(FISTA) \ No newline at end of file -- cgit v1.2.3 From b6c314b371ef3081828fa007cd3fcaf1dc820477 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 19 Oct 2017 17:07:48 +0100 Subject: First commit of CMakeLists.txt attempting to locate conda python environment --- src/Python/CMakeLists.txt | 79 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 src/Python/CMakeLists.txt diff --git a/src/Python/CMakeLists.txt b/src/Python/CMakeLists.txt new file mode 100644 index 0000000..a75c062 --- /dev/null +++ b/src/Python/CMakeLists.txt @@ -0,0 +1,79 @@ +cmake_minimum_required (VERSION 3.0) + +project(FISTA) + +# The version number. +set (FISTA_VERSION_MAJOR 1) +set (FISTA_VERSION_MINOR 0) + +set (CIL_VERSION_MAJOR 0) +set (CIL_VERSION_MINOR 9) +set (CIL_VERSION_PATCH 1) + +set (CIL_VERSION '${CIL_VERSION_MAJOR}.${CIL_VERSION_MINOR}.${CIL_VERSION_PATCH}') + +message("CIL VERSION " ${CIL_VERSION}) + +# variables we need to run the conda build +#PREFIX=C:\Apps\Miniconda2\envs\cil\Library +#LIBRARY_INC=C:\\Apps\\Miniconda2\\envs\\cil\\Library\\include + +set (NUMPY_VERSION 1.12) +#set (PYTHON_VERSION 3.5) + +#https://groups.google.com/a/continuum.io/forum/#!topic/anaconda/R9gWjl09UFs +set (CONDA_ENVIRONMENT "C:\\Apps\\Miniconda2\\envs\\cil27" CACHE PATH "env dir") + +function (findPythonForAnacondaEnvironment env) + + file(TO_CMAKE_PATH ${env}/python.exe PYTHON_EXECUTABLE) + + message("Found " ${PYTHON_EXECUTABLE}) + execute_process(COMMAND ${PYTHON_EXECUTABLE} pythonversion.py major + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} + OUTPUT_VARIABLE major + ERROR_QUIET) + execute_process(COMMAND ${PYTHON_EXECUTABLE} pythonversion.py minor + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} + OUTPUT_VARIABLE minor + ERROR_QUIET) + execute_process(COMMAND ${PYTHON_EXECUTABLE} pythonversion.py patch + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} + OUTPUT_VARIABLE patch + ERROR_QUIET) + execute_process(COMMAND ${PYTHON_EXECUTABLE} pythonversion.py + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} + OUTPUT_VARIABLE version + ERROR_QUIET) + + + set (PYTHON_EXECUTABLE ${PYTHON_EXECUTABLE} PARENT_SCOPE) + set (PYTHONINTERP_FOUND "ON" PARENT_SCOPE) + set (PYTHON_VERSION_STRING ${version}) + message("My version found " ${PYTHON_VERSION_STRING}) + +endfunction() + +findPythonForAnacondaEnvironment(${CONDA_ENVIRONMENT}) + +set(Python_ADDITIONAL_VERSIONS 3) + +find_package(PythonInterp) +if (PYTHONINTERP_FOUND) + + message("Found interpret " ${PYTHON_EXECUTABLE}) + + foreach(pv ${PYTHON_VERSION_STRING}) + message("Found interpret " ${pv}) + endforeach() +endif() + +find_package(PythonLibs) +if (PYTHONLIB_FOUND) + message("Found PythonLibs PYTHON_LIBRARIES " ${PYTHON_LIBRARIES}) + message("Found PythonLibs PYTHON_INCLUDE_PATH " ${PYTHON_INCLUDE_PATH}) + message("Found PythonLibs PYTHON_INCLUDE_DIRS " ${PYTHON_INCLUDE_DIRS}) + message("Found PythonLibs PYTHONLIBS_VERSION_STRING " ${PYTHONLIBS_VERSION_STRING} ) + +endif() + -- cgit v1.2.3 From 8b427d82acfaeb4671484bc459343c5e2e412736 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 19 Oct 2017 17:01:55 +0100 Subject: Initial revision of build environment made with CMake Initial revision of build environment made with CMake First commit of CMakeLists.txt attempt to use CMake to create the build environment First commit of CMakeLists.txt attempting to locate conda python environment Added a few files for CMake Many changes for the CMake compilation. Tested CMake build Bugfixes --- CMakeLists.txt | 30 ++++++ src/CMakeLists.txt | 14 +++ src/Python/CMakeLists.txt | 58 +++++++++++ src/Python/FindAnacondaEnvironment.cmake | 166 +++++++++++++++++++++++++++++++ src/Python/compile.bat.in | 4 + src/Python/compile.sh.in | 6 ++ src/Python/conda-recipe/bld.bat | 14 +++ src/Python/conda-recipe/build.sh | 14 +++ src/Python/conda-recipe/meta.yaml | 30 ++++++ 9 files changed, 336 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 src/CMakeLists.txt create mode 100644 src/Python/CMakeLists.txt create mode 100644 src/Python/FindAnacondaEnvironment.cmake create mode 100644 src/Python/compile.bat.in create mode 100644 src/Python/compile.sh.in create mode 100644 src/Python/conda-recipe/bld.bat create mode 100644 src/Python/conda-recipe/build.sh create mode 100644 src/Python/conda-recipe/meta.yaml diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..c9e8cb7 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,30 @@ +# Copyright 2017 Edoardo Pasca +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +cmake_minimum_required (VERSION 3.0) + +project(FISTA) +#https://stackoverflow.com/questions/13298504/using-cmake-with-setup-py + +# The version number. +set (FISTA_VERSION_MAJOR 1) +set (FISTA_VERSION_MINOR 0) + +set (CIL_VERSION_MAJOR 0) +set (CIL_VERSION_MINOR 9) +set (CIL_VERSION_PATCH 1) + +set (CIL_VERSION '${CIL_VERSION_MAJOR}.${CIL_VERSION_MINOR}.${CIL_VERSION_PATCH}') + +add_subdirectory(src) \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..cbe2fec --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,14 @@ +# Copyright 2017 Edoardo Pasca +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +add_subdirectory(Python) \ No newline at end of file diff --git a/src/Python/CMakeLists.txt b/src/Python/CMakeLists.txt new file mode 100644 index 0000000..fd377cc --- /dev/null +++ b/src/Python/CMakeLists.txt @@ -0,0 +1,58 @@ +# Copyright 2017 Edoardo Pasca +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +message("CIL VERSION " ${CIL_VERSION}) + +# variables that are set by conda +#PREFIX=C:\Apps\Miniconda2\envs\cil\Library +#LIBRARY_INC=C:\\Apps\\Miniconda2\\envs\\cil\\Library\\include + +set (NUMPY_VERSION 1.12) +#set (PYTHON_VERSION 3.5) + +#https://groups.google.com/a/continuum.io/forum/#!topic/anaconda/R9gWjl09UFs +set (CONDA_ENVIRONMENT "cil") +set (CONDA_ENVIRONMENT_PATH "C:\\Apps\\Miniconda2\\envs\\${CONDA_ENVIRONMENT}" CACHE PATH "env dir") + +message("CIL VERSION " ${CIL_VERSION}) + +# set the Python variables for the Conda environment +include(FindAnacondaEnvironment.cmake) +findPythonForAnacondaEnvironment(${CONDA_ENVIRONMENT_PATH}) +message("Python found " ${PYTHON_VERSION_STRING}) +findPythonPackagesPath() +message("PYTHON_PACKAGES_FOUND " ${PYTHON_PACKAGES_PATH}) + +# copy the Pyhon files of the package +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/ccpi/imaging/) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/__init__.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/imaging/__init__.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi/imaging) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/imaging/Regularizer.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi/imaging) + + +# Copy and configure the relative conda build and recipes +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in ${CMAKE_CURRENT_BINARY_DIR}/setup.py) +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/conda-recipe) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/conda-recipe/meta.yaml DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/conda-recipe) + +if (WIN32) + file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/conda-recipe/bld.bat DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/conda-recipe/) + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/compile.bat.in ${CMAKE_CURRENT_BINARY_DIR}/compile.bat) +elseif(UNIX) + file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/conda-recipe/build.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/conda-recipe/) + # assumes we will use bash + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/compile.sh.in ${CMAKE_CURRENT_BINARY_DIR}/compile.sh) +endif() + + diff --git a/src/Python/FindAnacondaEnvironment.cmake b/src/Python/FindAnacondaEnvironment.cmake new file mode 100644 index 0000000..3abb5d1 --- /dev/null +++ b/src/Python/FindAnacondaEnvironment.cmake @@ -0,0 +1,166 @@ +# Copyright 2017 Edoardo Pasca +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# #.rst: +# FindAnacondaEnvironment +# -------------- +# +# Find Python executable and library for a specific Anaconda environment +# +# This module finds the Python interpreter for a specific Anaconda enviroment, +# if installed and determines where the include files and libraries are. +# This code sets the following variables: +# +# :: +# PYTHONINTERP_FOUND - if the Python interpret has been found +# PYTHON_EXECUTABLE - the Python interpret found +# PYTHON_LIBRARY - path to the python library +# PYTHON_INCLUDE_PATH - path to where Python.h is found (deprecated) +# PYTHON_INCLUDE_DIRS - path to where Python.h is found +# PYTHONLIBS_VERSION_STRING - version of the Python libs found (since CMake 2.8.8) +# PYTHON_VERSION_MAJOR - major Python version +# PYTHON_VERSION_MINOR - minor Python version +# PYTHON_VERSION_PATCH - patch Python version + + + +function (findPythonForAnacondaEnvironment env) + + file(TO_CMAKE_PATH ${env}/python.exe PYTHON_EXECUTABLE) + + message("Found " ${PYTHON_EXECUTABLE}) + ####### FROM FindPythonInterpr ######## + # determine python version string + if(PYTHON_EXECUTABLE) + execute_process(COMMAND "${PYTHON_EXECUTABLE}" -c + "import sys; sys.stdout.write(';'.join([str(x) for x in sys.version_info[:3]]))" + OUTPUT_VARIABLE _VERSION + RESULT_VARIABLE _PYTHON_VERSION_RESULT + ERROR_QUIET) + if(NOT _PYTHON_VERSION_RESULT) + string(REPLACE ";" "." _PYTHON_VERSION_STRING "${_VERSION}") + list(GET _VERSION 0 _PYTHON_VERSION_MAJOR) + list(GET _VERSION 1 _PYTHON_VERSION_MINOR) + list(GET _VERSION 2 _PYTHON_VERSION_PATCH) + if(PYTHON_VERSION_PATCH EQUAL 0) + # it's called "Python 2.7", not "2.7.0" + string(REGEX REPLACE "\\.0$" "" _PYTHON_VERSION_STRING "${PYTHON_VERSION_STRING}") + endif() + else() + # sys.version predates sys.version_info, so use that + execute_process(COMMAND "${PYTHON_EXECUTABLE}" -c "import sys; sys.stdout.write(sys.version)" + OUTPUT_VARIABLE _VERSION + RESULT_VARIABLE _PYTHON_VERSION_RESULT + ERROR_QUIET) + if(NOT _PYTHON_VERSION_RESULT) + string(REGEX REPLACE " .*" "" _PYTHON_VERSION_STRING "${_VERSION}") + string(REGEX REPLACE "^([0-9]+)\\.[0-9]+.*" "\\1" _PYTHON_VERSION_MAJOR "${PYTHON_VERSION_STRING}") + string(REGEX REPLACE "^[0-9]+\\.([0-9])+.*" "\\1" _PYTHON_VERSION_MINOR "${PYTHON_VERSION_STRING}") + if(PYTHON_VERSION_STRING MATCHES "^[0-9]+\\.[0-9]+\\.([0-9]+)") + set(PYTHON_VERSION_PATCH "${CMAKE_MATCH_1}") + else() + set(PYTHON_VERSION_PATCH "0") + endif() + else() + # sys.version was first documented for Python 1.5, so assume + # this is older. + set(PYTHON_VERSION_STRING "1.4" PARENT_SCOPE) + set(PYTHON_VERSION_MAJOR "1" PARENT_SCOPE) + set(PYTHON_VERSION_MINOR "4" PARENT_SCOPE) + set(PYTHON_VERSION_PATCH "0" PARENT_SCOPE) + endif() + endif() + unset(_PYTHON_VERSION_RESULT) + unset(_VERSION) + endif() + ############################################### + + set (PYTHON_EXECUTABLE ${PYTHON_EXECUTABLE} PARENT_SCOPE) + set (PYTHONINTERP_FOUND "ON" PARENT_SCOPE) + set (PYTHON_VERSION_STRING ${_PYTHON_VERSION_STRING} PARENT_SCOPE) + set (PYTHON_VERSION_MAJOR ${_PYTHON_VERSION_MAJOR} PARENT_SCOPE) + set (PYTHON_VERSION_MINOR ${_PYTHON_VERSION_MINOR} PARENT_SCOPE) + set (PYTHON_VERSION_PATCH ${_PYTHON_VERSION_PATCH} PARENT_SCOPE) + message("My version found " ${PYTHON_VERSION_STRING}) + +endfunction() + + + +set(Python_ADDITIONAL_VERSIONS 3.5) + +find_package(PythonInterp) +if (PYTHONINTERP_FOUND) + + message("Found interpret " ${PYTHON_EXECUTABLE}) + message("Python Library " ${PYTHON_LIBRARY}) + message("Python Include Dir " ${PYTHON_INCLUDE_DIR}) + message("Python Include Path " ${PYTHON_INCLUDE_PATH}) + + foreach(pv ${PYTHON_VERSION_STRING}) + message("Found interpret " ${pv}) + endforeach() +endif() + + + +find_package(PythonLibs) +if (PYTHONLIB_FOUND) + message("Found PythonLibs PYTHON_LIBRARIES " ${PYTHON_LIBRARIES}) + message("Found PythonLibs PYTHON_INCLUDE_PATH " ${PYTHON_INCLUDE_PATH}) + message("Found PythonLibs PYTHON_INCLUDE_DIRS " ${PYTHON_INCLUDE_DIRS}) + message("Found PythonLibs PYTHONLIBS_VERSION_STRING " ${PYTHONLIBS_VERSION_STRING} ) +else() + message("No PythonLibs Found") +endif() + + + + +function(findPythonPackagesPath) +### https://openlab.ncl.ac.uk/gitlab/john.shearer/clappertracker/raw/549885e5decd37f7b23e9c1fd39e86f207156795/src/3rdparty/opencv/cmake/OpenCVDetectPython.cmake +### +if(CMAKE_HOST_UNIX) + execute_process(COMMAND ${PYTHON_EXECUTABLE} -c "from distutils.sysconfig import *; print get_python_lib()" + RESULT_VARIABLE PYTHON_CVPY_PROCESS + OUTPUT_VARIABLE PYTHON_STD_PACKAGES_PATH + OUTPUT_STRIP_TRAILING_WHITESPACE) + if("${PYTHON_STD_PACKAGES_PATH}" MATCHES "site-packages") + set(_PYTHON_PACKAGES_PATH "python${PYTHON_VERSION_MAJOR_MINOR}/site-packages") + else() #debian based assumed, install to the dist-packages. + set(_PYTHON_PACKAGES_PATH "python${PYTHON_VERSION_MAJOR_MINOR}/dist-packages") + endif() + if(EXISTS "${CMAKE_INSTALL_PREFIX}/lib${LIB_SUFFIX}/${PYTHON_PACKAGES_PATH}") + set(_PYTHON_PACKAGES_PATH "lib${LIB_SUFFIX}/${_PYTHON_PACKAGES_PATH}") + else() + set(_PYTHON_PACKAGES_PATH "lib/${_PYTHON_PACKAGES_PATH}") + endif() + elseif(CMAKE_HOST_WIN32) + get_filename_component(PYTHON_PATH "${PYTHON_EXECUTABLE}" PATH) + file(TO_CMAKE_PATH "${PYTHON_PATH}" PYTHON_PATH) + if(NOT EXISTS "${PYTHON_PATH}/Lib/site-packages") + unset(PYTHON_PATH) + get_filename_component(PYTHON_PATH "[HKEY_LOCAL_MACHINE\\SOFTWARE\\Python\\PythonCore\\${PYTHON_VERSION_MAJOR_MINOR}\\InstallPath]" ABSOLUTE) + if(NOT PYTHON_PATH) + get_filename_component(PYTHON_PATH "[HKEY_CURRENT_USER\\SOFTWARE\\Python\\PythonCore\\${PYTHON_VERSION_MAJOR_MINOR}\\InstallPath]" ABSOLUTE) + endif() + file(TO_CMAKE_PATH "${PYTHON_PATH}" PYTHON_PATH) + endif() + set(_PYTHON_PACKAGES_PATH "${PYTHON_PATH}/Lib/site-packages") + endif() + SET(PYTHON_PACKAGES_PATH "${_PYTHON_PACKAGES_PATH}" PARENT_SCOPE) + +endfunction() + + diff --git a/src/Python/compile.bat.in b/src/Python/compile.bat.in new file mode 100644 index 0000000..d4ddc92 --- /dev/null +++ b/src/Python/compile.bat.in @@ -0,0 +1,4 @@ +set CIL_VERSION=@CIL_VERSION@ + +activate @CONDA_ENVIRONMENT@ +conda build conda-recipe --python=@PYTHON_VERSION_MAJOR@.@PYTHON_VERSION_MINOR@ --numpy=@NUMPY_VERSION@ -c ccpi \ No newline at end of file diff --git a/src/Python/compile.sh.in b/src/Python/compile.sh.in new file mode 100644 index 0000000..dd29973 --- /dev/null +++ b/src/Python/compile.sh.in @@ -0,0 +1,6 @@ +#!/bin/sh + +export CIL_VERSION=@CIL_VERSION@ +module load python/anaconda +source activate @CONDA_ENVIRONMENT@ +conda build conda-recipe --python=@PYTHON_VERSION_MAJOR@.@PYTHON_VERSION_MINOR@ --numpy=@NUMPY_VERSION@ -c ccpi \ No newline at end of file diff --git a/src/Python/conda-recipe/bld.bat b/src/Python/conda-recipe/bld.bat new file mode 100644 index 0000000..69491de --- /dev/null +++ b/src/Python/conda-recipe/bld.bat @@ -0,0 +1,14 @@ +IF NOT DEFINED CIL_VERSION ( +ECHO CIL_VERSION Not Defined. +exit 1 +) + +mkdir "%SRC_DIR%\ccpi" +xcopy /e "%RECIPE_DIR%\..\.." "%SRC_DIR%\ccpi" + +cd %SRC_DIR%\ccpi\Python + +%PYTHON% setup.py build_ext +if errorlevel 1 exit 1 +%PYTHON% setup.py install +if errorlevel 1 exit 1 diff --git a/src/Python/conda-recipe/build.sh b/src/Python/conda-recipe/build.sh new file mode 100644 index 0000000..855047f --- /dev/null +++ b/src/Python/conda-recipe/build.sh @@ -0,0 +1,14 @@ + +if [ -z "$CIL_VERSION" ]; then + echo "Need to set CIL_VERSION" + exit 1 +fi +mkdir "$SRC_DIR/ccpi" +cp -r "$RECIPE_DIR/../.." "$SRC_DIR/ccpi" + +cd $SRC_DIR/ccpi/Python + +$PYTHON setup.py build_ext +$PYTHON setup.py install + + diff --git a/src/Python/conda-recipe/meta.yaml b/src/Python/conda-recipe/meta.yaml new file mode 100644 index 0000000..c5b7a89 --- /dev/null +++ b/src/Python/conda-recipe/meta.yaml @@ -0,0 +1,30 @@ +package: + name: ccpi-fista + version: {{ environ['CIL_VERSION'] }} + + +build: + preserve_egg_dir: False + script_env: + - CIL_VERSION +# number: 0 + +requirements: + build: + - python + - numpy + - setuptools + - boost ==1.64 + - boost-cpp ==1.64 + - cython + + run: + - python + - numpy + - boost ==1.64 + + +about: + home: http://www.ccpi.ac.uk + license: BSD license + summary: 'CCPi Core Imaging Library Quantification Toolbox' -- cgit v1.2.3