diff options
22 files changed, 973 insertions, 2324 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py index df8dc89..6c18ebf 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py +++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py @@ -1,21 +1,25 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= -# Copyright 2018-2019 Jakob Jorgensen, Daniil Kazantsev and 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. from ccpi.optimisation.functions import Function import numpy @@ -43,23 +47,69 @@ class IndicatorBox(Function): val = numpy.inf return val - def prox(self,x,tau=None): - return (x.maximum(self.lower)).minimum(self.upper) + def gradient(self,x): + return ValueError('Not Differentiable') + + def convex_conjugate(self,x): + # support function sup <x^*, x> + return 0 def proximal(self, x, tau, out=None): + if out is None: - return self.prox(x, tau) + return (x.maximum(self.lower)).minimum(self.upper) + else: + x.maximum(self.lower, out=out) + out.minimum(self.upper, out=out) + + def proximal_conjugate(self, x, tau, out=None): + + if out is None: + + return x - tau * self.proximal(x/tau, tau) + else: - if not x.shape == out.shape: - raise ValueError('Norm1 Incompatible size:', - x.shape, out.shape) - #(x.abs() - tau*self.gamma).maximum(0) * x.sign() - x.abs(out = out) - out.__isub__(tau*self.gamma) - out.maximum(0, out=out) - if self.sign_x is None or not x.shape == self.sign_x.shape: - self.sign_x = x.sign() - else: - x.sign(out=self.sign_x) - - out.__imul__( self.sign_x ) + + self.proximal(x/tau, tau, out=out) + out *= -1*tau + out += x + + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + + N, M = 2,3 + ig = ImageGeometry(voxel_num_x = N, voxel_num_y = M) + + u = ig.allocate('random_int') + tau = 2 + + f = IndicatorBox(2, 3) + + lower = 10 + upper = 30 + + z1 = f.proximal(u, tau) + + z2 = f.proximal_conjugate(u/tau, 1/tau) + + z = z1 + tau * z2 + + numpy.testing.assert_array_equal(z.as_array(), u.as_array()) + + out1 = ig.allocate() + out2 = ig.allocate() + + f.proximal(u, tau, out=out1) + f.proximal_conjugate(u/tau, 1/tau, out = out2) + + p = out1 + tau * out2 + + numpy.testing.assert_array_equal(p.as_array(), u.as_array()) + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 2de11ce..4dd77cf 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -25,7 +25,7 @@ import functools class KullbackLeibler(Function): - ''' Assume that data > 0 + ''' Assume that data >= 0 ''' @@ -129,7 +129,18 @@ class KullbackLeibler(Function): ''' - return ScaledFunction(self, scalar) + return ScaledFunction(self, scalar) + + +if __name__ == '__main__': + + from ccpi.framework import ImageGeometry + import numpy + + M, N = 2,3 + ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N) + u = ig.allocate('random_int') + b = np.random.normal(0, 0.1, size=ig.shape) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py index 4e53f2c..79040a0 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py @@ -1,21 +1,23 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and 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. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= from ccpi.optimisation.functions import Function diff --git a/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py index 39f0907..0875c20 100644 --- a/Wrappers/Python/wip/Compare_Algs/LeastSq_CGLS_FISTA_PDHG.py +++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py @@ -1,26 +1,28 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and 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. - +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= """ -Compare Least Squares minimization problem using FISTA, PDHG, CGLS classes -and Astra Built-in CGLS +Compare solutions of FISTA & PDHG + & CGLS & Astra Built-in algorithms for Least Squares + Problem: min_x || A x - g ||_{2}^{2} @@ -29,6 +31,7 @@ Problem: min_x || A x - g ||_{2}^{2} """ + from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry import numpy as np @@ -54,8 +57,14 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) detectors = N angles = np.linspace(0, np.pi, N, dtype=np.float32) +device = input('Available device: GPU==1 / CPU==0 ') ag = AcquisitionGeometry('parallel','2D', angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') +if device=='1': + dev = 'gpu' +else: + dev = 'cpu' + +Aop = AstraProjectorSimple(ig, ag, dev) sin = Aop.direct(data) noisy_data = sin @@ -63,7 +72,6 @@ noisy_data = sin # Setup and run Astra CGLS algorithm vol_geom = astra.create_vol_geom(N, N) proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) - proj_id = astra.create_projector('strip', proj_geom, vol_geom) # Create a sinogram from a phantom @@ -85,23 +93,18 @@ astra.algorithm.run(alg_id, 1000) recon_cgls_astra = astra.data2d.get(rec_id) # Setup and run the CGLS algorithm - -x_init = ig.allocate() - +x_init = ig.allocate() cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) cgls.max_iteration = 1000 cgls.update_objective_interval = 200 -cgls.run(1000) +cgls.run(1000, verbose=False) # Setup and run the PDHG algorithm - operator = Aop f = L2NormSquared(b = noisy_data) g = ZeroFunction() - ## Compute operator Norm normK = operator.norm() - ## Primal & dual stepsizes sigma = 1 tau = 1/(sigma*normK**2) @@ -109,10 +112,9 @@ tau = 1/(sigma*normK**2) pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) pdhg.max_iteration = 1000 pdhg.update_objective_interval = 200 -pdhg.run(1000) +pdhg.run(1000, verbose=False) # Setup and run the FISTA algorithm - fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop) regularizer = ZeroFunction() @@ -120,12 +122,12 @@ opt = {'memopt':True} fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt) fista.max_iteration = 1000 fista.update_objective_interval = 200 -fista.run(1000, verbose=True) +fista.run(1000, verbose=False) #%% Show results -plt.figure(figsize=(15,15)) -plt.suptitle('Reconstructions ') +plt.figure(figsize=(10,10)) +plt.suptitle('Reconstructions ', fontsize=16) plt.subplot(2,2,1) plt.imshow(cgls.get_output().as_array()) @@ -143,15 +145,10 @@ plt.subplot(2,2,4) plt.imshow(recon_cgls_astra) plt.title('CGLS astra') -#%% diff1 = pdhg.get_output() - cgls.get_output() diff2 = fista.get_output() - cgls.get_output() diff3 = ImageData(recon_cgls_astra) - cgls.get_output() -print( diff1.squared_norm()) -print( diff2.squared_norm()) -print( diff3.squared_norm()) - plt.figure(figsize=(15,15)) plt.subplot(3,1,1) diff --git a/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py index 984fca4..942d328 100644 --- a/Wrappers/Python/wip/Compare_Algs/Tikhonov_CGLS_PDHG.py +++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_PDHG_Tikhonov.py @@ -1,28 +1,31 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and 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. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= """ -Compare Tikhonov with PDHG, CGLS classes +Compare solutions of PDHG & "Block CGLS" algorithms for Problem: min_x alpha * ||\grad x ||^{2}_{2} + || A x - g ||_{2}^{2} + A: Projection operator g: Sinogram @@ -39,11 +42,9 @@ import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG, CGLS from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import ZeroFunction, BlockFunction, L2NormSquared -from ccpi.astra.ops import AstraProjectorSimple +from ccpi.optimisation.functions import ZeroFunction, BlockFunction, L2NormSquared # Create Ground truth phantom and Sinogram - N = 128 x = np.zeros((N,N)) x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 @@ -56,19 +57,22 @@ detectors = N angles = np.linspace(0, np.pi, N, dtype=np.float32) ag = AcquisitionGeometry('parallel','2D', angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') +device = input('Available device: GPU==1 / CPU==0 ') +ag = AcquisitionGeometry('parallel','2D', angles, detectors) +if device=='1': + dev = 'gpu' +else: + dev = 'cpu' + sin = Aop.direct(data) noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,3,ig.shape)) # Setup and run the CGLS algorithm - alpha = 50 Grad = Gradient(ig) # Form Tikhonov as a Block CGLS structure - -#%% op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1)) block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate()) @@ -77,16 +81,15 @@ x_init = ig.allocate() cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data) cgls.max_iteration = 1000 cgls.update_objective_interval = 200 -cgls.run(1000) +cgls.run(1000,verbose=False) + -#%% #Setup and run the PDHG algorithm # Create BlockOperator op_PDHG = BlockOperator(Grad, Aop, shape=(2,1) ) -# Create functions - +# Create functions f1 = 0.5 * alpha**2 * L2NormSquared() f2 = 0.5 * L2NormSquared(b = noisy_data) f = BlockFunction(f1, f2) @@ -102,18 +105,18 @@ tau = 1/(sigma*normK**2) pdhg = PDHG(f=f,g=g,operator=op_PDHG, tau=tau, sigma=sigma, memopt=True) pdhg.max_iteration = 1000 pdhg.update_objective_interval = 200 -pdhg.run(1000) +pdhg.run(1000, verbose=False) + #%% # Show results +plt.figure(figsize=(10,10)) -plt.figure(figsize=(15,15)) - -plt.subplot(1,2,1) +plt.subplot(2,1,1) plt.imshow(cgls.get_output().as_array()) plt.title('CGLS reconstruction') -plt.subplot(1,2,2) +plt.subplot(2,1,2) plt.imshow(pdhg.get_output().as_array()) plt.title('PDHG reconstruction') @@ -126,8 +129,6 @@ plt.title('Diff PDHG vs CGLS') plt.colorbar() plt.show() -#%% - plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') plt.plot(np.linspace(0,N,N), cgls.get_output().as_array()[int(N/2),:], label = 'CGLS') plt.legend() diff --git a/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py b/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py index eb62761..c24ebac 100644 --- a/Wrappers/Python/wip/Compare_Algs/FISTA_vs_PDHG.py +++ b/Wrappers/Python/demos/CompareAlgorithms/FISTA_vs_PDHG.py @@ -1,31 +1,36 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and 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. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= """ -Compare FISTA & PDHG classes +Compare solutions of FISTA & PDHG algorithms Problem: min_x alpha * ||\grad x ||^{2}_{2} + || x - g ||_{1} - A: Projection operator - g: Sinogram + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Salt & Pepper Noise """ @@ -44,8 +49,7 @@ from ccpi.optimisation.functions import L2NormSquared, L1Norm, \ from skimage.util import random_noise -# Create Ground truth and noisy data - +# Create Ground truth and Noisy data N = 100 data = np.zeros((N,N)) @@ -61,7 +65,7 @@ noisy_data = ImageData(n1) # Regularisation Parameter alpha = 5 - +############################################################################### # Setup and run the FISTA algorithm operator = Gradient(ig) fidelity = L1Norm(b=noisy_data) @@ -73,7 +77,10 @@ fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt) fista.max_iteration = 2000 fista.update_objective_interval = 50 fista.run(2000, verbose=False) +############################################################################### + +############################################################################### # Setup and run the PDHG algorithm op1 = Gradient(ig) op2 = Identity(ig, ag) @@ -89,19 +96,19 @@ tau = 1/(sigma*normK**2) pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 +pdhg.update_objective_interval = 200 pdhg.run(2000, verbose=False) +############################################################################### -#%% # Show results -plt.figure(figsize=(15,15)) +plt.figure(figsize=(10,10)) -plt.subplot(1,2,1) +plt.subplot(2,1,1) plt.imshow(pdhg.get_output().as_array()) plt.title('PDHG reconstruction') -plt.subplot(1,2,2) +plt.subplot(2,1,2) plt.imshow(fista.get_output().as_array()) plt.title('FISTA reconstruction') @@ -110,7 +117,7 @@ plt.show() diff1 = pdhg.get_output() - fista.get_output() plt.imshow(diff1.abs().as_array()) -plt.title('Diff PDHG vs CGLS') +plt.title('Diff PDHG vs FISTA') plt.colorbar() plt.show() diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py new file mode 100644 index 0000000..72d0670 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_TV_Tomo2D_poisson.py @@ -0,0 +1,258 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction, IndicatorBox + +from ccpi.astra.ops import AstraProjectorSimple + +""" + +Total Variation Denoising using PDHG algorithm: + + +Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + int A x -g log(Ax + \eta) + + \nabla: Gradient operator + + A: Projection Matrix + g: Noisy sinogram corrupted with Poisson Noise + + \eta: Background Noise + \alpha: Regularization parameter + +""" + +# Create phantom for TV 2D tomography +N = 50 +x = np.zeros((N,N)) +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +detectors = N +angles = np.linspace(0, np.pi, N) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data. Apply Poisson noise +scale = 0.25 +eta = 0 #np.random.randint(0, sin.as_array().max()/2, ag.shape) +n1 = scale * np.random.poisson(eta + sin.as_array()/scale) + +noisy_data = AcquisitionData(n1, ag) + + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,10)) +plt.subplot(2,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(2,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() + +#%% + +# Regularisation Parameter +alpha = 2 + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = KullbackLeibler(noisy_data) +f = BlockFunction(f1, f2) + +#g = ZeroFunction() +g = IndicatorBox(lower=0) + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 2 +tau = 1/(sigma*normK**2) +#sigma = 1/normK +#tau = 1/normK + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 500 +pdhg.run(2000, verbose = True) + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(pdhg.get_output().as_array()) +plt.title('TV Reconstruction') +plt.colorbar() +plt.show() +## +plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + +##%% Check with CVX solution +# +#from ccpi.optimisation.operators import SparseFiniteDiff +#import astra +#import numpy +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +# +#if cvx_not_installable: +# +# +# ##Construct problem +# u = Variable(N*N) +# #q = Variable() +# +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) +# +# # create matrix representation for Astra operator +# +# vol_geom = astra.create_vol_geom(N, N) +# proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) +# +# proj_id = astra.create_projector('strip', proj_geom, vol_geom) +# +# matrix_id = astra.projector.matrix(proj_id) +# +# ProjMat = astra.matrix.get(matrix_id) +# +# tmp = noisy_data.as_array().ravel('F') +# +## fidelity = sum( ProjMat * u - tmp * log(ProjMat * u + 1e-6)) +# #constraints = [q>= fidelity, u>=0] +# constraints = [] +# +# fidelity = sum(kl_div(tmp, ProjMat * u + 1e-6)) +## fidelity = kl_div(cp.multiply(alpha, W), +## cp.multiply(alpha, W + cp.multiply(beta, P))) - \ +## cp.multiply(alpha, cp.multiply(beta, P)) +# +# +# +# solver = SCS +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj, constraints) +# result = prob.solve(verbose = True, solver = solver) +# +# +####%% Check with CVX solution +## +##from ccpi.optimisation.operators import SparseFiniteDiff +## +##try: +## from cvxpy import * +## cvx_not_installable = True +##except ImportError: +## cvx_not_installable = False +## +## +##if cvx_not_installable: +## +## ##Construct problem +## u = Variable(ig.shape) +## +## DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +## DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +## +## # Define Total Variation as a regulariser +## regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) +## fidelity = pnorm( u - noisy_data.as_array(),1) +## +## # choose solver +## if 'MOSEK' in installed_solvers(): +## solver = MOSEK +## else: +## solver = SCS +## +## obj = Minimize( regulariser + fidelity) +## prob = Problem(obj) +## result = prob.solve(verbose = True, solver = solver) +## +## +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(pdhg.get_output().as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# plt.subplot(3,1,2) +# plt.imshow(np.reshape(u.value, (N, N))) +# plt.title('CVX solution') +# plt.colorbar() +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.legend() +# plt.title('Middle Line Profiles') +# plt.show() +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py new file mode 100644 index 0000000..045458a --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_2D_time_denoising.py @@ -0,0 +1,169 @@ +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018-2019 Evangelos Papoutsellis and 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. + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorMC + +import os +import tomophantom +from tomophantom import TomoP2D + +# Create phantom for TV 2D dynamic tomography + +model = 102 # note that the selected model is temporal (2D + time) +N = 50 # set dimension of the phantom +# one can specify an exact path to the parameters file +# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' +path = os.path.dirname(tomophantom.__file__) +path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +#This will generate a N_size x N_size x Time frames phantom (2D + time) +phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D) + +plt.close('all') +plt.figure(1) +plt.rcParams.update({'font.size': 21}) +plt.title('{}''{}'.format('2D+t phantom using model no.',model)) +for sl in range(0,np.shape(phantom_2Dt)[0]): + im = phantom_2Dt[sl,:,:] + plt.imshow(im, vmin=0, vmax=1) + plt.pause(.1) + plt.draw + + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0]) +data = ImageData(phantom_2Dt, geometry=ig) + +detectors = N +angles = np.linspace(0,np.pi,N) + +ag = AcquisitionGeometry('parallel','2D', angles, detectors, channels = np.shape(phantom_2Dt)[0]) +Aop = AstraProjectorMC(ig, ag, 'gpu') +sin = Aop.direct(data) + +scale = 2 +n1 = scale * np.random.poisson(sin.as_array()/scale) +noisy_data = AcquisitionData(n1, ag) + +tindex = [3, 6, 10] + +fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10)) +plt.subplot(1,3,1) +plt.imshow(noisy_data.as_array()[tindex[0],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) +plt.subplot(1,3,2) +plt.imshow(noisy_data.as_array()[tindex[1],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) +plt.subplot(1,3,3) +plt.imshow(noisy_data.as_array()[tindex[2],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + +fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8, + wspace=0.02, hspace=0.02) + +plt.show() + +#%% +# Regularisation Parameter +alpha = 5 + +# Create operators +#op1 = Gradient(ig) +op1 = Gradient(ig, correlation='SpaceChannels') +op2 = Aop + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = KullbackLeibler(noisy_data) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 200 +pdhg.run(2000) + + +#%% +fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) + +plt.subplot(2,3,1) +plt.imshow(phantom_2Dt[tindex[0],:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) + +plt.subplot(2,3,2) +plt.imshow(phantom_2Dt[tindex[1],:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) + +plt.subplot(2,3,3) +plt.imshow(phantom_2Dt[tindex[2],:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + + +plt.subplot(2,3,4) +plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:]) +plt.axis('off') +plt.subplot(2,3,5) +plt.imshow(pdhg.get_output().as_array()[tindex[1],:,:]) +plt.axis('off') +plt.subplot(2,3,6) +plt.imshow(pdhg.get_output().as_array()[tindex[2],:,:]) +plt.axis('off') +im = plt.imshow(pdhg.get_output().as_array()[tindex[0],:,:]) + + +fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8, + wspace=0.02, hspace=0.02) + +cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8]) +cbar = fig.colorbar(im, cax=cb_ax) + + +plt.show() + diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_2D_time.py new file mode 100644 index 0000000..14608db --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_2D_time.py @@ -0,0 +1,192 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= + +from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG + +from ccpi.optimisation.operators import BlockOperator, Gradient, Identity +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorMC + +import os +import tomophantom +from tomophantom import TomoP2D + +# Create phantom for TV 2D dynamic tomography + +model = 102 # note that the selected model is temporal (2D + time) +N = 128 # set dimension of the phantom +# one can specify an exact path to the parameters file +# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' +path = os.path.dirname(tomophantom.__file__) +path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +#This will generate a N_size x N_size x Time frames phantom (2D + time) +phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D) + +plt.close('all') +plt.figure(1) +plt.rcParams.update({'font.size': 21}) +plt.title('{}''{}'.format('2D+t phantom using model no.',model)) +for sl in range(0,np.shape(phantom_2Dt)[0]): + im = phantom_2Dt[sl,:,:] + plt.imshow(im, vmin=0, vmax=1) +# plt.pause(.1) +# plt.draw + + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0]) +data = ImageData(phantom_2Dt, geometry=ig) +ag = ig + +# Create Noisy data. Add Gaussian noise +np.random.seed(10) +noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) ) + +tindex = [3, 6, 10] + +fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10)) +plt.subplot(1,3,1) +plt.imshow(noisy_data.as_array()[tindex[0],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) +plt.subplot(1,3,2) +plt.imshow(noisy_data.as_array()[tindex[1],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) +plt.subplot(1,3,3) +plt.imshow(noisy_data.as_array()[tindex[2],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + +fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8, + wspace=0.02, hspace=0.02) + +plt.show() + +#%% +# Regularisation Parameter +alpha = 0.3 + +# Create operators +#op1 = Gradient(ig) +op1 = Gradient(ig, correlation='Space') +op2 = Gradient(ig, correlation='SpaceChannels') + +op3 = Identity(ig, ag) + +# Create BlockOperator +operator1 = BlockOperator(op1, op3, shape=(2,1) ) +operator2 = BlockOperator(op2, op3, shape=(2,1) ) + +# Create functions + +f1 = alpha * MixedL21Norm() +f2 = 0.5 * L2NormSquared(b = noisy_data) +f = BlockFunction(f1, f2) + +g = ZeroFunction() + +# Compute operator Norm +normK1 = operator1.norm() +normK2 = operator2.norm() + +#%% +# Primal & dual stepsizes +sigma1 = 1 +tau1 = 1/(sigma1*normK1**2) + +sigma2 = 1 +tau2 = 1/(sigma2*normK2**2) + +# Setup and run the PDHG algorithm +pdhg1 = PDHG(f=f,g=g,operator=operator1, tau=tau1, sigma=sigma1) +pdhg1.max_iteration = 2000 +pdhg1.update_objective_interval = 200 +pdhg1.run(2000) + +# Setup and run the PDHG algorithm +pdhg2 = PDHG(f=f,g=g,operator=operator2, tau=tau2, sigma=sigma2) +pdhg2.max_iteration = 2000 +pdhg2.update_objective_interval = 200 +pdhg2.run(2000) + + +#%% + +tindex = [3, 6, 10] +fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) + +plt.subplot(3,3,1) +plt.imshow(phantom_2Dt[tindex[0],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[0])) + +plt.subplot(3,3,2) +plt.imshow(phantom_2Dt[tindex[1],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[1])) + +plt.subplot(3,3,3) +plt.imshow(phantom_2Dt[tindex[2],:,:]) +plt.axis('off') +plt.title('Time {}'.format(tindex[2])) + +plt.subplot(3,3,4) +plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:]) +plt.axis('off') +plt.subplot(3,3,5) +plt.imshow(pdhg1.get_output().as_array()[tindex[1],:,:]) +plt.axis('off') +plt.subplot(3,3,6) +plt.imshow(pdhg1.get_output().as_array()[tindex[2],:,:]) +plt.axis('off') + + +plt.subplot(3,3,7) +plt.imshow(pdhg2.get_output().as_array()[tindex[0],:,:]) +plt.axis('off') +plt.subplot(3,3,8) +plt.imshow(pdhg2.get_output().as_array()[tindex[1],:,:]) +plt.axis('off') +plt.subplot(3,3,9) +plt.imshow(pdhg2.get_output().as_array()[tindex[2],:,:]) +plt.axis('off') + +#%% +im = plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:]) + + +fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8, + wspace=0.02, hspace=0.02) + +cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8]) +cbar = fig.colorbar(im, cax=cb_ax) + + +plt.show() + diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian_3D.py index c86ddc9..03dc2ef 100644 --- a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian_3D.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian_3D.py @@ -1,26 +1,46 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and 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. +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= +""" + +Total Variation (3D) Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Gaussian Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + +""" from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy + import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG @@ -47,6 +67,8 @@ path_library3D = os.path.join(path, "Phantom3DLibrary.dat") #This will generate a N x N x N phantom (3D) phantom_tm = TomoP3D.Model(model, N, path_library3D) +#%% + # Create noisy data. Add Gaussian noise ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N) ag = ig @@ -69,6 +91,7 @@ plt.title('Sagittal View') plt.colorbar() plt.show() +#%% # Regularisation Parameter alpha = 0.05 @@ -107,14 +130,17 @@ normK = operator.norm() sigma = 1 tau = 1/(sigma*normK**2) - # Setup and run the PDHG algorithm pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) pdhg.max_iteration = 2000 pdhg.update_objective_interval = 200 -pdhg.run(2000) +pdhg.run(2000, verbose = True) + +#%% fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8)) +fig.suptitle('TV Reconstruction',fontsize=20) + plt.subplot(2,3,1) plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1) diff --git a/Wrappers/Python/demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_gaussian.py index 87d5328..6acbfcc 100644 --- a/Wrappers/Python/demos/PDHG_TV_Tomo2D.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_gaussian.py @@ -1,21 +1,39 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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. +# +#========================================================================= -# Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca +""" + +Total Variation 2D Tomography Reconstruction using PDHG algorithm: -# 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 +Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + \frac{1}{2}||Ax - g||^{2} -# 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. + \nabla: Gradient operator + + A: Projection Matrix + g: Noisy sinogram corrupted with Gaussian Noise + + \alpha: Regularization parameter + +""" from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData @@ -25,33 +43,16 @@ import matplotlib.pyplot as plt from ccpi.optimisation.algorithms import PDHG -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction from ccpi.astra.ops import AstraProjectorSimple -""" - -Total Variation Denoising using PDHG algorithm: - - min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) - -Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + int A x -g log(Ax + \eta) - - \nabla: Gradient operator - - A: Projection Matrix - g: Noisy sinogram corrupted with Poisson Noise - - \eta: Background Noise - \alpha: Regularization parameter - -""" # Create phantom for TV 2D tomography -N = 75 +N = 100 x = np.zeros((N,N)) x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 @@ -60,19 +61,38 @@ data = ImageData(x) ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) detectors = N -angles = np.linspace(0, np.pi, N, dtype=np.float32) +angles = np.linspace(0, np.pi, N) ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') + +device = input('Available device: GPU==1 / CPU==0 ') + +if device=='1': + dev = 'gpu' +else: + dev = 'cpu' + +Aop = AstraProjectorSimple(ig, ag, dev) sin = Aop.direct(data) # Create noisy data. Apply Poisson noise -scale = 2 -n1 = scale * np.random.poisson(sin.as_array()/scale) -noisy_data = AcquisitionData(n1, ag) +n1 = np.random.normal(0, 3, size=ig.shape) +noisy_data = AcquisitionData(n1 + sin.as_array(), ag) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(10,10)) +plt.subplot(2,1,1) +plt.imshow(data.as_array()) +plt.title('Ground Truth') +plt.colorbar() +plt.subplot(2,1,2) +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Data') +plt.colorbar() +plt.show() # Regularisation Parameter -alpha = 5 +alpha = 50 # Create operators op1 = Gradient(ig) @@ -84,40 +104,24 @@ operator = BlockOperator(op1, op2, shape=(2,1) ) # Create functions f1 = alpha * MixedL21Norm() -f2 = KullbackLeibler(noisy_data) +f2 = 0.5 * L2NormSquared(b=noisy_data) f = BlockFunction(f1, f2) g = ZeroFunction() - -diag_precon = True -if diag_precon: +# Compute operator Norm +normK = operator.norm() - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - # Compute operator Norm - normK = operator.norm() - # Primal & dual stepsizes - sigma = 10 - tau = 1/(sigma*normK**2) - +# Primal & dual stepsizes +sigma = 10 +tau = 1/(sigma*normK**2) # Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 +pdhg.update_objective_interval = 200 pdhg.run(2000) - -#%% plt.figure(figsize=(15,15)) plt.subplot(3,1,1) plt.imshow(data.as_array()) @@ -151,14 +155,11 @@ try: cvx_not_installable = True except ImportError: cvx_not_installable = False - - + if cvx_not_installable: - ##Construct problem u = Variable(N*N) - #q = Variable() DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') @@ -166,7 +167,6 @@ if cvx_not_installable: regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) # create matrix representation for Astra operator - vol_geom = astra.create_vol_geom(N, N) proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) @@ -176,49 +176,16 @@ if cvx_not_installable: ProjMat = astra.matrix.get(matrix_id) - fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) - #constraints = [q>= fidelity, u>=0] - constraints = [u>=0] - - solver = SCS - obj = Minimize( regulariser + fidelity) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - -##%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) + tmp = noisy_data.as_array().ravel() - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = pnorm( u - noisy_data.as_array(),1) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - + fidelity = 0.5 * sum_squares(ProjMat * u - tmp) + + solver = MOSEK obj = Minimize( regulariser + fidelity) prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - np.reshape(u.value, (N,N) )) plt.figure(figsize=(15,15)) plt.subplot(3,1,1) @@ -236,7 +203,7 @@ if cvx_not_installable: plt.show() plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') + plt.plot(np.linspace(0,N,N), np.reshape(u.value, (N,N) )[int(N/2),:], label = 'CVX') plt.legend() plt.title('Middle Line Profiles') plt.show() diff --git a/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py b/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py deleted file mode 100644 index 3155654..0000000 --- a/Wrappers/Python/wip/Compare_Algs/PDHG_vs_CGLS.py +++ /dev/null @@ -1,127 +0,0 @@ -# -*- coding: utf-8 -*- -# This work is part of the Core Imaging Library developed by -# Visual Analytics and Imaging System Group of the Science Technology -# Facilities Council, STFC - -# Copyright 2018-2019 Evangelos Papoutsellis and 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. - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, CGLS - -from ccpi.optimisation.operators import Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition -from skimage.util import random_noise -from ccpi.astra.ops import AstraProjectorSimple - -#%% - -N = 128 -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -detectors = N -angles = np.linspace(0, np.pi, N, dtype=np.float32) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -noisy_data = sin - -x_init = ig.allocate() - -## Setup and run the CGLS algorithm -cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data) -cgls.max_iteration = 500 -cgls.update_objective_interval = 50 -cgls.run(500, verbose=True) - -# Create BlockOperator -operator = Aop -f = 0.5 * L2NormSquared(b = noisy_data) -g = ZeroFunction() - -## Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes -sigma = 0.1 -tau = 1/(sigma*normK**2) -# -# -## Setup and run the PDHG algorithm -pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -pdhg.max_iteration = 2000 -pdhg.update_objective_interval = 50 -pdhg.run(2000) - -#%% - -diff = pdhg.get_output() - cgls.get_output() -print( diff.norm()) -# -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(pdhg.get_output().as_array()) -plt.title('PDHG reconstruction') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(cgls.get_output().as_array()) -plt.title('CGLS reconstruction') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow(diff.abs().as_array()) -plt.title('Difference reconstruction') -plt.colorbar() -plt.show() - - - - - - - - - - - - - - - - - - - - - - -# -# -# -# -# -# -# -# diff --git a/Wrappers/Python/wip/demo_box_constraints_FISTA.py b/Wrappers/Python/wip/demo_box_constraints_FISTA.py index 2f9e0c6..b15dd45 100644 --- a/Wrappers/Python/wip/demo_box_constraints_FISTA.py +++ b/Wrappers/Python/wip/demo_box_constraints_FISTA.py @@ -72,7 +72,7 @@ else: # Set up Operator object combining the ImageGeometry and AcquisitionGeometry # wrapping calls to ASTRA as well as specifying whether to use CPU or GPU. -Aop = AstraProjectorSimple(ig, ag, 'gpu') +Aop = AstraProjectorSimple(ig, ag, 'cpu') Aop = Identity(ig,ig) diff --git a/Wrappers/Python/wip/fista_TV_denoising.py b/Wrappers/Python/wip/fista_TV_denoising.py deleted file mode 100644 index a9712fc..0000000 --- a/Wrappers/Python/wip/fista_TV_denoising.py +++ /dev/null @@ -1,72 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old, FISTA - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction - -from ccpi.optimisation.algs import FISTA - -from skimage.util import random_noise - -from timeit import default_timer as timer -def dt(steps): - return steps[-1] - steps[-2] - -# Create phantom for TV denoising - -N = 100 - -data = np.zeros((N,N)) -data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -# Regularisation Parameter -alpha = 2 - -operator = Gradient(ig) -g = alpha * MixedL21Norm() -f = 0.5 * L2NormSquared(b = noisy_data) - -x_init = ig.allocate() -opt = {'niter':2000} - - -x = FISTA(x_init, f, g, opt) - -#fista = FISTA() -#fista.set_up(x_init, f, g, opt ) -#fista.max_iteration = 10 -# -#fista.run(2000) -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(fista.get_output().as_array()) -#plt.title('no memopt class') - - - diff --git a/Wrappers/Python/wip/pdhg_TGV_denoising.py b/Wrappers/Python/wip/pdhg_TGV_denoising.py deleted file mode 100644 index 1c570cb..0000000 --- a/Wrappers/Python/wip/pdhg_TGV_denoising.py +++ /dev/null @@ -1,230 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TGV Gaussian denoising - -N = 100 - -data = np.zeros((N,N)) - -x1 = np.linspace(0, int(N/2), N) -x2 = np.linspace(int(N/2), 0., N) -xv, yv = np.meshgrid(x1, x2) - -xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T - -data = xv -data = data/data.max() - -plt.imshow(data) -plt.show() - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.005, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -alpha, beta = 0.1, 0.5 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '1' - -if method == '0': - - # Create operators - op11 = Gradient(ig) - op12 = Identity(op11.range_geometry()) - - op22 = SymmetrizedGradient(op11.domain_geometry()) - op21 = ZeroOperator(ig, op22.range_geometry()) - - op31 = Identity(ig, ag) - op32 = ZeroOperator(op22.domain_geometry(), ag) - - operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) - - - f1 = alpha * MixedL21Norm() - f2 = beta * MixedL21Norm() - f3 = 0.5 * L2NormSquared(b = noisy_data) - f = BlockFunction(f1, f2, f3) - g = ZeroFunction() - - -else: - - # Create operators - op11 = Gradient(ig) - op12 = Identity(op11.range_geometry()) - op22 = SymmetrizedGradient(op11.domain_geometry()) - op21 = ZeroOperator(ig, op22.range_geometry()) - - operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) ) - - f1 = alpha * MixedL21Norm() - f2 = beta * MixedL21Norm() - - f = BlockFunction(f1, f2) - g = BlockFunction(0.5 * L2NormSquared(b = noisy_data), ZeroFunction()) - -## Compute operator Norm -normK = operator.norm() -# -## Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -## -opt = {'niter':500} -opt1 = {'niter':500, 'memopt': True} -# -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() -# -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() -# -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res[0].as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1[0].as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1[0] - res[0]).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() - -print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) - - -###### - -#%% - -#def update_plot(it_update, objective, x): -# -## sol = pdhg.get_output() -# plt.figure() -# plt.imshow(x[0].as_array()) -# plt.show() -# -# -##def stop_criterion(x,) -# -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 100 -# -#pdhg.run(4000, verbose=False, callback=update_plot) - - -#%% - - - - - - - - -#%% Check with CVX solution - -#from ccpi.optimisation.operators import SparseFiniteDiff -# -#try: -# from cvxpy import * -# cvx_not_installable = True -#except ImportError: -# cvx_not_installable = False -# -#if cvx_not_installable: -# -# u = Variable(ig.shape) -# w1 = Variable((N, N)) -# w2 = Variable((N, N)) -# -# # create TGV regulariser -# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -# -# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ -# DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ -# beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ -# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ -# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) -# -# constraints = [] -# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) -# solver = MOSEK -# -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj) -# result = prob.solve(verbose = True, solver = solver) -# -# diff_cvx = numpy.abs( res[0].as_array() - u.value ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res[0].as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# -# plt.subplot(3,1,2) -# plt.imshow(u.value) -# plt.title('CVX solution') -# plt.colorbar() -# -# plt.subplot(3,1,3) -# plt.imshow(diff_cvx) -# plt.title('Difference') -# plt.colorbar() -# plt.show() -# -# plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.legend() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1])) -# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) -# - - diff --git a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py b/Wrappers/Python/wip/pdhg_TGV_tomography2D.py deleted file mode 100644 index ee3b089..0000000 --- a/Wrappers/Python/wip/pdhg_TGV_tomography2D.py +++ /dev/null @@ -1,200 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, \ - Gradient, SymmetrizedGradient, ZeroOperator -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -from ccpi.astra.ops import AstraProjectorSimple - -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TGV Gaussian denoising - -N = 100 - -data = np.zeros((N,N)) - -x1 = np.linspace(0, int(N/2), N) -x2 = np.linspace(int(N/2), 0., N) -xv, yv = np.meshgrid(x1, x2) - -xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T - -data = xv -data = ImageData(data/data.max()) - -plt.imshow(data.as_array()) -plt.show() - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -detectors = N -angles = np.linspace(0,np.pi,N) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -plt.imshow(sin.as_array()) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - -#%% - -alpha, beta = 20, 50 - - -# Create operators -op11 = Gradient(ig) -op12 = Identity(op11.range_geometry()) - -op22 = SymmetrizedGradient(op11.domain_geometry()) -op21 = ZeroOperator(ig, op22.range_geometry()) - -op31 = Aop -op32 = ZeroOperator(op22.domain_geometry(), ag) - -operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) - - -f1 = alpha * MixedL21Norm() -f2 = beta * MixedL21Norm() -f3 = 0.5 * L2NormSquared(b = noisy_data) -f = BlockFunction(f1, f2, f3) -g = ZeroFunction() - -## Compute operator Norm -normK = operator.norm() -# -## Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -## -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} -# -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() -# -plt.imshow(res[0].as_array()) -plt.show() - - -#t3 = timer() -#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -#t4 = timer() -# -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(res[0].as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1[0].as_array()) -#plt.title('memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((res1[0] - res[0]).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() -# -#print("NoMemopt/Memopt is {}/{}".format(t2-t1, t4-t3)) - - -#%% Check with CVX solution - -#from ccpi.optimisation.operators import SparseFiniteDiff -# -#try: -# from cvxpy import * -# cvx_not_installable = True -#except ImportError: -# cvx_not_installable = False -# -#if cvx_not_installable: -# -# u = Variable(ig.shape) -# w1 = Variable((N, N)) -# w2 = Variable((N, N)) -# -# # create TGV regulariser -# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -# -# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u) - vec(w1), \ -# DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \ -# beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \ -# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \ -# 0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0 ) ) -# -# constraints = [] -# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) -# solver = MOSEK -# -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj) -# result = prob.solve(verbose = True, solver = solver) -# -# diff_cvx = numpy.abs( res[0].as_array() - u.value ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res[0].as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# -# plt.subplot(3,1,2) -# plt.imshow(u.value) -# plt.title('CVX solution') -# plt.colorbar() -# -# plt.subplot(3,1,3) -# plt.imshow(diff_cvx) -# plt.title('Difference') -# plt.colorbar() -# plt.show() -# -# plt.plot(np.linspace(0,N,N), res[0].as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.legend() -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1])) -# print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max())) - - - diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py deleted file mode 100755 index b16e8f9..0000000 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ /dev/null @@ -1,204 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -#def dt(steps): -# return steps[-1] - steps[-2] - -# Create phantom for TV Gaussian denoising - -N = 100 - -data = np.zeros((N,N)) -data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10) -noisy_data = ImageData(n1) - - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy data') -plt.show() - -# Regularisation Parameter -alpha = 0.5 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '1' - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Form Composite Operator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - #### Create functions - - f1 = alpha * MixedL21Norm() - f2 = 0.5 * L2NormSquared(b = noisy_data) - f = BlockFunction(f1, f2) - - g = ZeroFunction() - -else: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = 0.5 * L2NormSquared(b = noisy_data) - - -diag_precon = False - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( res.as_array() - u.value ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.as_array()) - plt.title('PDHG solution') - plt.colorbar() - - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - - - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - - - - diff --git a/Wrappers/Python/wip/pdhg_TV_denoising3D.py b/Wrappers/Python/wip/pdhg_TV_denoising3D.py deleted file mode 100644 index 06ecfa2..0000000 --- a/Wrappers/Python/wip/pdhg_TV_denoising3D.py +++ /dev/null @@ -1,360 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer -def dt(steps): - return steps[-1] - steps[-2] - -#%% - -# Create phantom for TV denoising - -import timeit -import os -from tomophantom import TomoP3D -import tomophantom - -print ("Building 3D phantom using TomoPhantom software") -tic=timeit.default_timer() -model = 13 # select a model number from the library -N_size = 64 # Define phantom dimensions using a scalar value (cubic phantom) -path = os.path.dirname(tomophantom.__file__) -path_library3D = os.path.join(path, "Phantom3DLibrary.dat") -#This will generate a N_size x N_size x N_size phantom (3D) -phantom_tm = TomoP3D.Model(model, N_size, path_library3D) -#toc=timeit.default_timer() -#Run_time = toc - tic -#print("Phantom has been built in {} seconds".format(Run_time)) -# -#sliceSel = int(0.5*N_size) -##plt.gray() -#plt.figure() -#plt.subplot(131) -#plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) -#plt.title('3D Phantom, axial view') -# -#plt.subplot(132) -#plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1) -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(133) -#plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1) -#plt.title('3D Phantom, sagittal view') -#plt.show() - -#%% - -N = N_size -ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N) - -n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10) -noisy_data = ImageData(n1) -#plt.imshow(noisy_data.as_array()[:,:,32]) - -#%% - -# Regularisation Parameter -alpha = 0.02 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '0' - -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig) - - # Form Composite Operator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - #### Create functions - - f1 = alpha * MixedL21Norm() - f2 = 0.5 * L2NormSquared(b = noisy_data) - f = BlockFunction(f1, f2) - - g = ZeroFunction() - -else: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * FunctionOperatorComposition(operator, MixedL21Norm()) - g = L2NormSquared(b = noisy_data) - - ########################################################################### -#%% - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} - -#t1 = timer() -#res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -#t2 = timer() - - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -#import cProfile -#cProfile.run('res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) ') -### -print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) -# -## -##%% -# -#plt.figure(figsize=(10,10)) -#plt.subplot(311) -#plt.imshow(res1.as_array()[sliceSel,:,:]) -#plt.colorbar() -#plt.title('3D Phantom, axial view') -# -#plt.subplot(312) -#plt.imshow(res1.as_array()[:,sliceSel,:]) -#plt.colorbar() -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(313) -#plt.imshow(res1.as_array()[:,:,sliceSel]) -#plt.colorbar() -#plt.title('3D Phantom, sagittal view') -#plt.show() -# -#plt.figure(figsize=(10,10)) -#plt.subplot(311) -#plt.imshow(res.as_array()[sliceSel,:,:]) -#plt.colorbar() -#plt.title('3D Phantom, axial view') -# -#plt.subplot(312) -#plt.imshow(res.as_array()[:,sliceSel,:]) -#plt.colorbar() -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(313) -#plt.imshow(res.as_array()[:,:,sliceSel]) -#plt.colorbar() -#plt.title('3D Phantom, sagittal view') -#plt.show() -# -#diff = (res1 - res).abs() -# -#plt.figure(figsize=(10,10)) -#plt.subplot(311) -#plt.imshow(diff.as_array()[sliceSel,:,:]) -#plt.colorbar() -#plt.title('3D Phantom, axial view') -# -#plt.subplot(312) -#plt.imshow(diff.as_array()[:,sliceSel,:]) -#plt.colorbar() -#plt.title('3D Phantom, coronal view') -# -#plt.subplot(313) -#plt.imshow(diff.as_array()[:,:,sliceSel]) -#plt.colorbar() -#plt.title('3D Phantom, sagittal view') -#plt.show() -# -# -# -# -##%% -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 100 -#### -#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -#pdhgo.max_iteration = 2000 -#pdhgo.update_objective_interval = 100 -#### -#steps = [timer()] -#pdhgo.run(2000) -#steps.append(timer()) -#t1 = dt(steps) -## -#pdhg.run(2000) -#steps.append(timer()) -#t2 = dt(steps) -# -#print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1)) -#res = pdhg.get_output() -#res1 = pdhgo.get_output() - -#%% -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(res.as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1.as_array()) -#plt.title('memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((res1 - res).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() - - -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(pdhg.get_output().as_array()) -#plt.title('no memopt class') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res.as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((pdhg.get_output() - res).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() -# -# -# -#plt.figure(figsize=(15,15)) -#plt.subplot(3,1,1) -#plt.imshow(pdhgo.get_output().as_array()) -#plt.title('memopt class') -#plt.colorbar() -#plt.subplot(3,1,2) -#plt.imshow(res1.as_array()) -#plt.title('no memopt') -#plt.colorbar() -#plt.subplot(3,1,3) -#plt.imshow((pdhgo.get_output() - res1).abs().as_array()) -#plt.title('diff') -#plt.colorbar() -#plt.show() - - - - - -# print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1)) -# res = pdhg.get_output() -# res1 = pdhgo.get_output() -# -# diff = (res-res1) -# print ("diff norm {} max {}".format(diff.norm(), diff.abs().as_array().max())) -# print ("Sum ( abs(diff) ) {}".format(diff.abs().sum())) -# -# -# plt.figure(figsize=(5,5)) -# plt.subplot(1,3,1) -# plt.imshow(res.as_array()) -# plt.colorbar() -# #plt.show() -# -# #plt.figure(figsize=(5,5)) -# plt.subplot(1,3,2) -# plt.imshow(res1.as_array()) -# plt.colorbar() - -#plt.show() - - - -#======= -## opt = {'niter':2000, 'memopt': True} -# -## res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -# -#>>>>>>> origin/pdhg_fix -# -# -## opt = {'niter':2000, 'memopt': False} -## res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -# -## plt.figure(figsize=(5,5)) -## plt.subplot(1,3,1) -## plt.imshow(res.as_array()) -## plt.title('memopt') -## plt.colorbar() -## plt.subplot(1,3,2) -## plt.imshow(res1.as_array()) -## plt.title('no memopt') -## plt.colorbar() -## plt.subplot(1,3,3) -## plt.imshow((res1 - res).abs().as_array()) -## plt.title('diff') -## plt.colorbar() -## plt.show() -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 100 -# -# -#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True) -#pdhgo.max_iteration = 2000 -#pdhgo.update_objective_interval = 100 -# -#steps = [timer()] -#pdhgo.run(200) -#steps.append(timer()) -#t1 = dt(steps) -# -#pdhg.run(200) -#steps.append(timer()) -#t2 = dt(steps) -# -#print ("Time difference {} {} {}".format(t1,t2,t2-t1)) -#sol = pdhg.get_output().as_array() -##sol = result.as_array() -## -#fig = plt.figure() -#plt.subplot(1,3,1) -#plt.imshow(noisy_data.as_array()) -#plt.colorbar() -#plt.subplot(1,3,2) -#plt.imshow(sol) -#plt.colorbar() -#plt.subplot(1,3,3) -#plt.imshow(pdhgo.get_output().as_array()) -#plt.colorbar() -# -#plt.show() -### -## -#### -##plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth') -##plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') -##plt.legend() -##plt.show() -# -# -##%% -## diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py deleted file mode 100644 index bd330fc..0000000 --- a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py +++ /dev/null @@ -1,268 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction - -from skimage.util import random_noise - -from timeit import default_timer as timer - - - -# ############################################################################ -# Create phantom for TV denoising - -N = 100 -data = np.zeros((N,N)) -data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) -noisy_data = ImageData(n1) - -plt.imshow(noisy_data.as_array()) -plt.colorbar() -plt.show() - -#%% - -# Regularisation Parameter -alpha = 2 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '0' -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - operator = BlockOperator(op1, op2, shape=(2,1) ) - - f1 = alpha * MixedL21Norm() - f2 = L1Norm(b = noisy_data) - - f = BlockFunction(f1, f2 ) - g = ZeroFunction() - -else: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = L1Norm(b = noisy_data) - ########################################################################### -#%% - -diag_precon = False - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 10 -# -#pdhg.run(2000) - - - -#sol = pdhg.get_output().as_array() -##sol = result.as_array() -## -#fig = plt.figure() -#plt.subplot(1,2,1) -#plt.imshow(noisy_data.as_array()) -##plt.colorbar() -#plt.subplot(1,2,2) -#plt.imshow(sol) -##plt.colorbar() -#plt.show() -## - -## -#plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth') -#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') -#plt.legend() -#plt.show() - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - fidelity = pnorm( u - noisy_data.as_array(),1) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( res.as_array() - u.value ) - -# Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.as_array()) - plt.title('PDHG solution') - plt.colorbar() - - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - - - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - - - -#try_cvx = input("Do you want CVX comparison (0/1)") -# -#if try_cvx=='0': -# -# from cvxpy import * -# import sys -# sys.path.insert(0,'/Users/evangelos/Desktop/Projects/CCPi/CCPi-Framework/Wrappers/Python/ccpi/optimisation/cvx_scripts') -# from cvx_functions import TV_cvx -# -# u = Variable((N, N)) -# fidelity = pnorm( u - noisy_data.as_array(),1) -# regulariser = alpha * TV_cvx(u) -# solver = MOSEK -# obj = Minimize( regulariser + fidelity) -# constraints = [] -# prob = Problem(obj, constraints) -# -# # Choose solver (SCS is fast but less accurate than MOSEK) -# result = prob.solve(verbose = True, solver = solver) -# -# print('Objective value is {} '.format(obj.value)) -# -# diff_pdhg_cvx = np.abs(u.value - res.as_array()) -# plt.imshow(diff_pdhg_cvx) -# plt.colorbar() -# plt.title('|CVX-PDHG|') -# plt.show() -# -# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') -# plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG') -# plt.legend() -# plt.show() -# -#else: -# print('No CVX solution available') - - - - diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py deleted file mode 100644 index e123739..0000000 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ /dev/null @@ -1,231 +0,0 @@ -# -*- coding: utf-8 -*- - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction - -from ccpi.astra.operators import AstraProjectorSimple -from timeit import default_timer as timer - - -#%%############################################################################### -# Create phantom for TV tomography - -#import os -#import tomophantom -#from tomophantom import TomoP2D -#from tomophantom.supp.qualitymetrics import QualityTools - -#model = 1 # select a model number from the library -#N = 150 # set dimension of the phantom -## one can specify an exact path to the parameters file -## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' -#path = os.path.dirname(tomophantom.__file__) -#path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -##This will generate a N_size x N_size phantom (2D) -#phantom_2D = TomoP2D.Model(model, N, path_library2D) -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -#data = ImageData(phantom_2D, geometry=ig) - -N = 75 -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -data = ImageData(x) -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - - -detectors = N -angles = np.linspace(0,np.pi,N) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'cpu') -sin = Aop.direct(data) - -plt.imshow(sin.as_array()) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Form Composite Operator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -alpha = 10 -f = BlockFunction( alpha * MixedL21Norm(), \ - 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFunction() - -# Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes -diag_precon = True - -if diag_precon: - - def tau_sigma_precond(operator): - - tau = 1/operator.sum_abs_row() - sigma = 1/ operator.sum_abs_col() - - return tau, sigma - - tau, sigma = tau_sigma_precond(operator) - -else: - normK = operator.norm() - sigma = 1 - tau = 1/(sigma*normK**2) - -# Compute operator Norm - - -# Primal & dual stepsizes - -opt = {'niter':200} -opt1 = {'niter':200, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() -# -# -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() -# -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() -# -print(" Max of abs difference is {}".format(diff)) - - -#%% Check with CVX solution - -#from ccpi.optimisation.operators import SparseFiniteDiff -#import astra -#import numpy -# -#try: -# from cvxpy import * -# cvx_not_installable = True -#except ImportError: -# cvx_not_installable = False -# -# -#if cvx_not_installable: -# -# -# u = Variable( N*N) -# -# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') -# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') -# -# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) -# -# # create matrix representation for Astra operator -# -# vol_geom = astra.create_vol_geom(N, N) -# proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) -# -# proj_id = astra.create_projector('strip', proj_geom, vol_geom) -# -# matrix_id = astra.projector.matrix(proj_id) -# -# ProjMat = astra.matrix.get(matrix_id) -# -# fidelity = 0.5 * sum_squares(ProjMat * u - noisy_data.as_array().ravel()) -# -# # choose solver -# if 'MOSEK' in installed_solvers(): -# solver = MOSEK -# else: -# solver = SCS -# -# obj = Minimize( regulariser + fidelity) -# prob = Problem(obj) -# result = prob.solve(verbose = True, solver = solver) -# -##%% -# -# u_rs = np.reshape(u.value, (N,N)) -# -# diff_cvx = numpy.abs( res.as_array() - u_rs ) -# -# # Show result -# plt.figure(figsize=(15,15)) -# plt.subplot(3,1,1) -# plt.imshow(res.as_array()) -# plt.title('PDHG solution') -# plt.colorbar() -# -# plt.subplot(3,1,2) -# plt.imshow(u_rs) -# plt.title('CVX solution') -# plt.colorbar() -# -# plt.subplot(3,1,3) -# plt.imshow(diff_cvx) -# plt.title('Difference') -# plt.colorbar() -# plt.show() -# -# plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') -# plt.plot(np.linspace(0,N,N), u_rs[int(N/2),:], label = 'CVX') -# plt.legend() -# -# -# print('Primal Objective (CVX) {} '.format(obj.value)) -# print('Primal Objective (PDHG) {} '.format(primal[-1]))
\ No newline at end of file diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py deleted file mode 100644 index 5423b22..0000000 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py +++ /dev/null @@ -1,152 +0,0 @@ -# -*- coding: utf-8 -*- - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData - -import numpy as np -import matplotlib.pyplot as plt - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction, ScaledFunction - -from ccpi.astra.ops import AstraProjectorSimple, AstraProjectorMC -from skimage.util import random_noise - - -#%%############################################################################### -# Create phantom for TV tomography - -import numpy as np -import matplotlib.pyplot as plt -import os -import tomophantom -from tomophantom import TomoP2D - -model = 102 # note that the selected model is temporal (2D + time) -N = 150 # set dimension of the phantom -# one can specify an exact path to the parameters file -# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' -path = os.path.dirname(tomophantom.__file__) -path_library2D = os.path.join(path, "Phantom2DLibrary.dat") -#This will generate a N_size x N_size x Time frames phantom (2D + time) -phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D) - -plt.close('all') -plt.figure(1) -plt.rcParams.update({'font.size': 21}) -plt.title('{}''{}'.format('2D+t phantom using model no.',model)) -for sl in range(0,np.shape(phantom_2Dt)[0]): - im = phantom_2Dt[sl,:,:] - plt.imshow(im, vmin=0, vmax=1) - plt.pause(.1) - plt.draw - -#N = 150 -#x = np.zeros((N,N)) -#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -#%% -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0]) -data = ImageData(phantom_2Dt, geometry=ig) - - - -detectors = 150 -angles = np.linspace(0,np.pi,100) - -ag = AcquisitionGeometry('parallel','2D',angles, detectors, channels = np.shape(phantom_2Dt)[0]) -Aop = AstraProjectorMC(ig, ag, 'gpu') -sin = Aop.direct(data) - -plt.imshow(sin.as_array()[10]) -plt.title('Sinogram') -plt.colorbar() -plt.show() - -# Add Gaussian noise to the sinogram data -np.random.seed(10) -n1 = np.random.random(sin.shape) - -noisy_data = sin + ImageData(5*n1) - -plt.imshow(noisy_data.as_array()[10]) -plt.title('Noisy Sinogram') -plt.colorbar() -plt.show() - - -#%% Works only with Composite Operator Structure of PDHG - -#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) - -# Create operators -op1 = Gradient(ig) -op2 = Aop - -# Form Composite Operator -operator = BlockOperator(op1, op2, shape=(2,1) ) - -alpha = 50 -f = BlockFunction( alpha * MixedL21Norm(), \ - 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFunction() - -# Compute operator Norm -normK = operator.norm() - -## Primal & dual stepsizes - -sigma = 1 -tau = 1/(sigma*normK**2) - -#sigma = 1/normK -#tau = 1/normK - -opt = {'niter':2000} - -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) - -plt.figure(figsize=(5,5)) -plt.imshow(res.as_array()) -plt.colorbar() -plt.show() - -#sigma = 10 -#tau = 1/(sigma*normK**2) -# -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 5000 -#pdhg.update_objective_interval = 20 -# -#pdhg.run(5000) -# -##%% -#sol = pdhg.get_output().as_array() -#fig = plt.figure() -#plt.subplot(1,2,1) -#plt.imshow(noisy_data.as_array()) -##plt.colorbar() -#plt.subplot(1,2,2) -#plt.imshow(sol) -##plt.colorbar() -#plt.show() - - -#%% -plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') -plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') -plt.legend() -plt.show() - - diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py deleted file mode 100644 index cb37598..0000000 --- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py +++ /dev/null @@ -1,187 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos -""" - -import numpy as np -import numpy -import matplotlib.pyplot as plt - -from ccpi.framework import ImageData, ImageGeometry - -from ccpi.optimisation.algorithms import PDHG, PDHG_old - -from ccpi.optimisation.operators import BlockOperator, Identity, Gradient -from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ - MixedL21Norm, BlockFunction - - -from skimage.util import random_noise -from timeit import default_timer as timer - - - -# ############################################################################ -# Create phantom for TV Poisson denoising - -N = 200 -data = np.zeros((N,N)) -data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1 - -ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -ag = ig - -# Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'poisson', seed = 10) -noisy_data = ImageData(n1) - -plt.imshow(noisy_data.as_array()) -plt.colorbar() -plt.show() - -# Regularisation Parameter -alpha = 2 - -#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") - -method = '0' -if method == '0': - - # Create operators - op1 = Gradient(ig) - op2 = Identity(ig, ag) - - # Form Composite Operator - operator = BlockOperator(op1, op2, shape=(2,1) ) - - f1 = alpha * MixedL21Norm() - f2 = KullbackLeibler(noisy_data) - - f = BlockFunction(f1, f2 ) - g = ZeroFunction() - -else: - - ########################################################################### - # No Composite # - ########################################################################### - operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = KullbackLeibler(noisy_data) - ########################################################################### - -# Compute operator Norm -normK = operator.norm() - -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) - -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} - -t1 = timer() -res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) -t2 = timer() - -t3 = timer() -res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) -t4 = timer() - -print(pdgap[-1]) - - -plt.figure(figsize=(15,15)) -plt.subplot(3,1,1) -plt.imshow(res.as_array()) -plt.title('no memopt') -plt.colorbar() -plt.subplot(3,1,2) -plt.imshow(res1.as_array()) -plt.title('memopt') -plt.colorbar() -plt.subplot(3,1,3) -plt.imshow((res1 - res).abs().as_array()) -plt.title('diff') -plt.colorbar() -plt.show() -# -plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') -plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') -plt.legend() -plt.show() - -print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) -diff = (res1 - res).abs().as_array().max() - -print(" Max of abs difference is {}".format(diff)) - - -#%% Check with CVX solution - -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u1 = Variable(ig.shape) - q = Variable() - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) - - fidelity = sum( u1 - multiply(noisy_data.as_array(), log(u1)) ) - constraints = [q>= fidelity, u1>=0] - - solver = ECOS - obj = Minimize( regulariser + q) - prob = Problem(obj, constraints) - result = prob.solve(verbose = True, solver = solver) - - - diff_cvx = numpy.abs( res.as_array() - u1.value ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.as_array()) - plt.title('PDHG solution') - plt.colorbar() - - plt.subplot(3,1,2) - plt.imshow(u1.value) - plt.title('CVX solution') - plt.colorbar() - - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX') - plt.legend() - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) - - - - - |