diff options
| -rw-r--r-- | CMakeLists.txt | 30 | ||||
| -rw-r--r-- | demos/Demo1.m | 174 | ||||
| -rw-r--r-- | demos/Demo_Phantom3D_Cone.m | 66 | ||||
| -rw-r--r-- | demos/Demo_Phantom3D_Parallel.m | 66 | ||||
| -rw-r--r-- | main_func/FISTA_REC.m | 267 | ||||
| -rw-r--r-- | src/CMakeLists.txt | 14 | ||||
| -rw-r--r-- | src/Python/CMakeLists.txt | 59 | ||||
| -rw-r--r-- | src/Python/FindAnacondaEnvironment.cmake | 166 | ||||
| -rw-r--r-- | src/Python/compile.bat.in | 4 | ||||
| -rw-r--r-- | src/Python/compile.sh.in | 6 | ||||
| -rw-r--r-- | src/Python/conda-recipe/bld.bat | 14 | ||||
| -rw-r--r-- | src/Python/conda-recipe/build.sh | 14 | ||||
| -rw-r--r-- | src/Python/conda-recipe/meta.yaml | 30 |
13 files changed, 601 insertions, 309 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..d05cdd9 --- /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) 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..3a9178b --- /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: +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) +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..6a54450 --- /dev/null +++ b/demos/Demo_Phantom3D_Parallel.m @@ -0,0 +1,66 @@ +% 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 = 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:
+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
+% 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.subsets = 12;
+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');
+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 6987dca..1e4228d 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); + [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); + y = sqweight.*y'; + astra_mex_data2d('delete', sino_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!'); + 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,19 @@ 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 (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 - [sino_id, sino_updt(:,:,kkk)] = astra_create_sino3d_cuda(X_t(:,:,kkk), proj_geomT, vol_geomT); - astra_mex_data3d('delete', sino_id); + [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 +346,36 @@ 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 + objective(i) = 0.5*norm(residual(:)); % for the objective function output + 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); - % regularization + % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -510,95 +484,121 @@ 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; - 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 - % 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; + 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; - r_old = r; 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); - 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); - 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)); + 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) + - 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)); - end + 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 - vec = sum(residual2,2); - if (SlicesZ > 1) - vec = squeeze(vec(:,1,:)); + + elseif (studentt > 0) + % student t data fidelity + + % 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,:))); - - 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 + % PWLS model + + 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); + % 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'); + 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; - astra_mex_data3d('delete', sino_id); - astra_mex_data3d('delete', id); - % regularization + % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -680,20 +680,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) 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..3eb4158 --- /dev/null +++ b/src/Python/CMakeLists.txt @@ -0,0 +1,59 @@ +# 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' |