Merge branch 'pythonize' of https://github.com/vais-ral/CCPi-FISTA_Reconstruction into pythonize

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'