diff options
Diffstat (limited to 'Wrappers/Python')
3 files changed, 143 insertions, 8 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 752fd21..484dc61 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -8,7 +8,7 @@ Created on Thu Feb 14 12:36:40 2019 import numpy from numbers import Number import functools -from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer +from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer, DataContainer from ccpi.optimisation.operators import Operator, LinearOperator from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator from ccpi.framework import BlockGeometry @@ -224,7 +224,7 @@ class BlockOperator(Operator): res = [] for row in range(self.shape[0]): - for col in range(self.shape[1]): + for col in range(self.shape[1]): if col == 0: prod = self.get_item(row,col).sum_abs_row() else: @@ -269,8 +269,8 @@ if __name__ == '__main__': B = BlockOperator(G, Id) - print(B.sum_abs_row().as_array()) - + print(B.sum_abs_row()) +# Gx = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') Gy = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') @@ -282,17 +282,17 @@ if __name__ == '__main__': d_res = numpy.reshape(d1 + d2 + d3, ig.shape, 'F') print(d_res) - +# z1 = abs(Gx.matrix()).toarray().sum(axis=1) z2 = abs(Gy.matrix()).toarray().sum(axis=1) z3 = abs(Id.matrix()).toarray().sum(axis=1) - +# z_res = BlockDataContainer(BlockDataContainer(ImageData(numpy.reshape(z2, ig.shape, 'F')),\ ImageData(numpy.reshape(z1, ig.shape, 'F'))),\ ImageData(numpy.reshape(z3, ig.shape, 'F'))) - +# ttt = B.sum_abs_col() - +# numpy.testing.assert_array_almost_equal(z_res[0][0].as_array(), ttt[0][0].as_array(), decimal=4) numpy.testing.assert_array_almost_equal(z_res[0][1].as_array(), ttt[0][1].as_array(), decimal=4) numpy.testing.assert_array_almost_equal(z_res[1].as_array(), ttt[1].as_array(), decimal=4) diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py index 87723f0..cd65ee4 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py @@ -122,6 +122,7 @@ if __name__ == '__main__': d2 = G_neum.sum_abs_col() print(d2) + d1 * d2 diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py new file mode 100644 index 0000000..518ead2 --- /dev/null +++ b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py @@ -0,0 +1,134 @@ +#!/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 ZeroFun, L2NormSquared, \ + MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction + +from skimage.util import random_noise + + + +# ############################################################################ +# 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', seed=10) +noisy_data = ImageData(n1) + + +#%% + +# 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) ) + + #### Create functions +# f = FunctionComposition_new(operator, mixed_L12Norm(alpha), \ +# L2NormSq(0.5, b = noisy_data) ) + + f1 = alpha * MixedL21Norm() + f2 = 0.5 * L2NormSquared(b = noisy_data) + + f = BlockFunction(f1, f2 ) + g = ZeroFun() + +else: + + ########################################################################### + # No Composite # + ########################################################################### + operator = Gradient(ig) + f = alpha * FunctionOperatorComposition(operator, MixedL21Norm()) + g = 0.5 * L2NormSquared(b = noisy_data) + ########################################################################### +#%% + +diag_precon = True + +if diag_precon: + tmp_tau = operator.sum_abs_row() + tmp_sigma = operator.sum_abs_col() + + tmp_sigma[0][0].as_array()[tmp_sigma[0][0].as_array()==0]=1 + tmp_sigma[0][1].as_array()[tmp_sigma[0][1].as_array()==0]=1 + tmp_sigma[1].as_array()[tmp_sigma[1].as_array()==0]=1 + + tau = 1/tmp_tau + sigma = 1/tmp_sigma + +else: + # Compute operator Norm + normK = operator.norm() + print ("normK", normK) + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + +#%% + +#opt = {'niter':2000} + +#res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) + + +#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() +# + +#%% +# |