diff options
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py | 26 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py (renamed from Wrappers/Python/demos/PDHG_TGV_Denoising_SaltPepper.py) | 67 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py (renamed from Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian.py) | 35 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py (renamed from Wrappers/Python/demos/PDHG_TV_Denoising_Poisson.py) | 98 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py (renamed from Wrappers/Python/demos/PDHG_TV_Denoising_SaltPepper.py) | 69 | ||||
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py (renamed from Wrappers/Python/demos/PDHG_Tikhonov_Denoising.py) | 74 |
6 files changed, 239 insertions, 130 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index cf1a244..3096191 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -106,17 +106,27 @@ class KullbackLeibler(Function): z = x + tau * self.bnoise return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) else: - - z_m = x + tau * self.bnoise -1 - self.b.multiply(4*tau, out=out) - z_m.multiply(z_m, out=z_m) - out += z_m - + + tmp = x + tau * self.bnoise + self.b.multiply(4*tau, out=out) + out.add((tmp-1)**2, out=out) out.sqrt(out=out) - out *= -1 - out += tmp2 + out.add(tmp+1, out=out) out *= 0.5 + + + +# z_m = x + tau * self.bnoise -1 +# self.b.multiply(4*tau, out=out) +# z_m.multiply(z_m, out=z_m) +# out += z_m +# +# out.sqrt(out=out) +# +# out *= -1 +# out += tmp2 +# out *= 0.5 diff --git a/Wrappers/Python/demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py index 7b65c31..57f6fcd 100644 --- a/Wrappers/Python/demos/PDHG_TGV_Denoising_SaltPepper.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TGV_Denoising_SaltPepper.py @@ -1,9 +1,48 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Feb 22 14:53:03 2019 - -@author: evangelos +#======================================================================== +# 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 Generalised Variation (TGV) Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x - w||_{2,1} + + \beta * || E w ||_{2,1} + + \frac{1}{2} * || x - g ||_{2}^{2} + + \alpha: Regularization parameter + \alpha: Regularization parameter + + \nabla: Gradient operator + E: Symmetrized Gradient operator + + g: Noisy Data with Salt & Pepper Noise + + Method = 0 ( PDHG - split ) : K = [ \nabla, - Identity + ZeroOperator, E + Identity, ZeroOperator] + + + Method = 1 (PDHG - explicit ): K = [ \nabla, - Identity + ZeroOperator, E ] + """ from ccpi.framework import ImageData, ImageGeometry @@ -43,6 +82,18 @@ ag = ig n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) noisy_data = ImageData(n1) +# Show Ground Truth and Noisy Data +plt.figure(figsize=(15,15)) +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 Parameters alpha = 0.8 beta = numpy.sqrt(2)* alpha @@ -94,10 +145,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.run(2000) +pdhg.run(2000, verbose = False) #%% plt.figure(figsize=(15,15)) diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py index c830025..afdb6a2 100644 --- a/Wrappers/Python/demos/PDHG_TV_Denoising_Gaussian.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Gaussian.py @@ -18,14 +18,10 @@ # limitations under the License. # #========================================================================= - - """ Total Variation Denoising using PDHG algorithm: - min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) - Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2} @@ -35,12 +31,12 @@ Problem: min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2} g: Noisy Data with Gaussian Noise - Method = 0: K = [ \nabla, - Identity] + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + - Method = 1: K = \nabla - - + Method = 1 (PDHG - explicit ): K = \nabla + """ from ccpi.framework import ImageData, ImageGeometry @@ -54,20 +50,17 @@ from ccpi.optimisation.algorithms import PDHG from ccpi.optimisation.operators import BlockOperator, Identity, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction - - -from ccpi.data import camera - - -# Load Data -data = camera(size=(256,256)) - -N, M = data.shape - -# Image and Acquitition Geometries + +# Load Data +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 +data = ImageData(data) ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) ag = ig - + # Create Noisy data. Add Gaussian noise np.random.seed(10) noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.1, size=ig.shape) ) diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py index 70f6b9b..4d53635 100644 --- a/Wrappers/Python/demos/PDHG_TV_Denoising_Poisson.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_Poisson.py @@ -1,41 +1,43 @@ -# -*- 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 STFC, 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 - -# 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 Denoising using PDHG algorithm: - min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) +Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + \int x - g * log(x) - -Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + \int x - g * log(x) - - \nabla: Gradient operator - g: Noisy Data with Poisson Noise \alpha: Regularization parameter - Method = 0: K = [ \nabla, - Identity] - - Method = 1: K = \nabla + \nabla: Gradient operator + + g: Noisy Data with Poisson Noise + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + + + Method = 1 (PDHG - explicit ): K = \nabla + """ from ccpi.framework import ImageData, ImageGeometry @@ -66,10 +68,25 @@ ag = ig n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) noisy_data = ImageData(n1) + +# Show Ground Truth and Noisy Data +plt.figure(figsize=(15,15)) +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 -method = '1' +method = '0' if method == '0': @@ -99,26 +116,15 @@ else: # Compute operator Norm normK = operator.norm() -# Primal & dual stepsizes +# Primal & Dual stepsizes sigma = 1 tau = 1/(sigma*normK**2) -opt = {'niter':2000, 'memopt': True} - -# 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 - -def pdgap_objectives(niter, objective, solution): - - - print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\ - format(niter, pdhg.max_iteration,'', \ - objective[0],'',\ - objective[1],'',\ - objective[2])) -pdhg.run(2000, callback = pdgap_objectives) +# Setup and Run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) +pdhg.max_iteration = 3000 +pdhg.update_objective_interval = 200 +pdhg.run(3000, verbose=False) plt.figure(figsize=(15,15)) diff --git a/Wrappers/Python/demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py index f5d4ce4..c5709c3 100644 --- a/Wrappers/Python/demos/PDHG_TV_Denoising_SaltPepper.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Denoising_SaltPepper.py @@ -1,39 +1,43 @@ -# -*- 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 Denoising using PDHG algorithm: - min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) - -Problem: min_x, x>0 \alpha * ||\nabla x||_{1} + ||x-g||_{1} +Problem: min_x, x>0 \alpha * ||\nabla x||_{2,1} + ||x-g||_{1} + \alpha: Regularization parameter + \nabla: Gradient operator + g: Noisy Data with Salt & Pepper Noise - \alpha: Regularization parameter - Method = 0: K = [ \nabla, - Identity] + + Method = 0 ( PDHG - split ) : K = [ \nabla, + Identity] + - Method = 1: K = \nabla + Method = 1 (PDHG - explicit ): K = \nabla """ @@ -66,6 +70,18 @@ ag = ig n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) noisy_data = ImageData(n1) +# Show Ground Truth and Noisy Data +plt.figure(figsize=(15,15)) +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 @@ -80,8 +96,7 @@ if method == '0': # Create BlockOperator operator = BlockOperator(op1, op2, shape=(2,1) ) - # Create functions - + # Create functions f1 = alpha * MixedL21Norm() f2 = L1Norm(b = noisy_data) f = BlockFunction(f1, f2) diff --git a/Wrappers/Python/demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py index 041d4ee..7b73c1a 100644 --- a/Wrappers/Python/demos/PDHG_Tikhonov_Denoising.py +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_Tikhonov_Denoising.py @@ -1,21 +1,43 @@ -# -*- 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. +# +#========================================================================= +""" + +Tikhonov Denoising using PDHG algorithm: + + +Problem: min_{x} \alpha * ||\nabla x||_{2}^{2} + \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 @@ -41,13 +63,25 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) ag = ig # Create noisy data. Apply Salt & Pepper noise -n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10) +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) noisy_data = ImageData(n1) +# Show Ground Truth and Noisy Data +plt.figure(figsize=(15,15)) +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 = 4 -method = '0' +method = '1' if method == '0': |