diff options
5 files changed, 535 insertions, 27 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py index 0836324..4178e6c 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -20,11 +20,14 @@ class BlockFunction(Function): ''' def __init__(self, *functions): - + + super(BlockFunction, self).__init__() self.functions = functions self.length = len(self.functions) - super(BlockFunction, self).__init__() + self.L = self.functions.L + + def __call__(self, x): diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py index 6c18ebf..ace0ba7 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py +++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py @@ -22,6 +22,7 @@ from ccpi.optimisation.functions import Function import numpy +from ccpi.framework import ImageData class IndicatorBox(Function): '''Box constraints indicator function. @@ -39,7 +40,7 @@ class IndicatorBox(Function): def __call__(self,x): - + if (numpy.all(x.array>=self.lower) and numpy.all(x.array <= self.upper) ): val = 0 @@ -51,8 +52,9 @@ class IndicatorBox(Function): return ValueError('Not Differentiable') def convex_conjugate(self,x): - # support function sup <x^*, x> - return 0 + # support function sup <x, z>, z \in [lower, upper] + # ???? + return x.maximum(0).sum() def proximal(self, x, tau, out=None): @@ -78,7 +80,7 @@ class IndicatorBox(Function): if __name__ == '__main__': - from ccpi.framework import ImageGeometry + from ccpi.framework import ImageGeometry, BlockDataContainer N, M = 2,3 ig = ImageGeometry(voxel_num_x = N, voxel_num_y = M) @@ -109,6 +111,23 @@ if __name__ == '__main__': numpy.testing.assert_array_equal(p.as_array(), u.as_array()) + d = f.convex_conjugate(u) + print(d) + + + + # what about n-dimensional Block + #uB = BlockDataContainer(u,u,u) + #lowerB = BlockDataContainer(1,2,3) + #upperB = BlockDataContainer(10,21,30) + + #fB = IndicatorBox(lowerB, upperB) + + #z1B = fB.proximal(uB, tau) + + + + diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index 4dd77cf..0d3c8f5 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -20,33 +20,42 @@ import numpy from ccpi.optimisation.functions import Function from ccpi.optimisation.functions.ScaledFunction import ScaledFunction -from ccpi.framework import ImageData, ImageGeometry import functools +import scipy.special class KullbackLeibler(Function): - ''' Assume that data >= 0 + ''' + + KL_div(x, y + back) = int x * log(x/(y+back)) - x + (y+back) + + Assumption: y>=0 + back>=0 ''' - def __init__(self,data, **kwargs): + def __init__(self, data, **kwargs): super(KullbackLeibler, self).__init__() - self.b = data - self.bnoise = kwargs.get('bnoise', 0) - - + self.b = data + self.bnoise = 0 + + def __call__(self, x): - - res_tmp = numpy.zeros(x.shape) - tmp = x + self.bnoise - ind = x.as_array()>0 - res_tmp[ind] = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind]) + ''' - return res_tmp.sum() + x - y * log( x + bnoise) + y * log(y) - y + bnoise + + + ''' + + # TODO avoid scipy import ???? + tmp = scipy.special.kl_div(self.b.as_array(), x.as_array()) + return numpy.sum(tmp) + def log(self, datacontainer): '''calculates the in-place log of the datacontainer''' @@ -58,7 +67,6 @@ class KullbackLeibler(Function): def gradient(self, x, out=None): - #TODO Division check if out is None: return 1 - self.b/(x + self.bnoise) else: @@ -70,12 +78,10 @@ class KullbackLeibler(Function): def convex_conjugate(self, x): - tmp = self.b/(1-x) - ind = tmp.as_array()>0 - - return (self.b.as_array()[ind] * (numpy.log(tmp.as_array()[ind])-1)).sum() - - + # TODO avoid scipy import ???? + xlogy = scipy.special.xlogy(self.b.as_array(), 1 - x.as_array()) + return numpy.sum(-xlogy) + def proximal(self, x, tau, out=None): if out is None: @@ -140,8 +146,86 @@ if __name__ == '__main__': 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) + b = ig.allocate('random_int') + u.as_array()[1,1]=0 + u.as_array()[2,0]=0 + b.as_array()[1,1]=0 + b.as_array()[2,0]=0 + + f = KullbackLeibler(b) + + +# longest = reduce(lambda x, y: len(x) if len(x) > len(y) else len(y), strings) + + +# tmp = functools.reduce(lambda x, y: \ +# 0 if x==0 and not numpy.isnan(y) else x * numpy.log(y), \ +# zip(b.as_array().ravel(), u.as_array().ravel()),0) + + +# np.multiply.reduce(X, 0) + + +# sf = reduce(lambda x,y: x + y[0]*y[1], +# zip(self.as_array().ravel(), +# other.as_array().ravel()), +# 0) +#cdef inline number_t xlogy(number_t x, number_t y) nogil: +# if x == 0 and not zisnan(y): +# return 0 +# else: +# return x * zlog(y) + +# if npy_isnan(x): +# return x +# elif x > 0: +# return -x * log(x) +# elif x == 0: +# return 0 +# else: +# return -inf + +# cdef inline double kl_div(double x, double y) nogil: +# if npy_isnan(x) or npy_isnan(y): +# return nan +# elif x > 0 and y > 0: +# return x * log(x / y) - x + y +# elif x == 0 and y >= 0: +# return y +# else: +# return inf + + + + +# def xlogy(self, dc1, dc2): + +# return numpy.sum(numpy.where(dc1.as_array() != 0, dc2.as_array() * numpy.log(dc2.as_array() / dc1.as_array()), 0)) + + +# f.xlog(u, b) + + + + +# tmp1 = b.as_array() +# tmp2 = u.as_array() +# +# zz = scipy.special.xlogy(tmp1, tmp2) +# +# print(np.sum(zz)) + + +# ww = f.xlogy(b, u) + +# print(ww) + + +#cdef inline double kl_div(double x, double y) nogil: + + + diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py new file mode 100644 index 0000000..5b1bb16 --- /dev/null +++ b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py @@ -0,0 +1,184 @@ +#======================================================================== +# 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 for Poisson denoising using FISTA algorithm: + +Problem: min_x, x>0 \alpha * ||\nabla x||_{2} + \int x - g * log(x) + + \alpha: Regularization parameter + + \nabla: Gradient operator + + g: Noisy Data with Poisson Noise + + +""" + +from ccpi.framework import ImageData + +import numpy as np +import numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import FISTA + +from ccpi.optimisation.operators import Gradient, BlockOperator, Identity +from ccpi.optimisation.functions import KullbackLeibler, IndicatorBox, BlockFunction, \ + L2NormSquared, IndicatorBox, FunctionOperatorComposition + +from ccpi.framework import TestData +import os, sys +from skimage.util import random_noise + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + +# Load Data +N = 50 +M = 50 +data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) + +ig = data.geometry +ag = ig + +# Create Noisy data. Add Gaussian noise +n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) +noisy_data = ImageData(n1) + +# 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 = 20 + +# Setup and run the FISTA algorithm +op1 = Gradient(ig) +op2 = BlockOperator(Identity(ig), Identity(ig), shape=(2,1)) + +tmp_function = BlockFunction( KullbackLeibler(noisy_data), IndicatorBox(lower=0) ) + +fid = tmp +reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator) + +x_init = ig.allocate() +opt = {'memopt':True} +fista = FISTA(x_init=x_init , f=reg, g=fid, opt=opt) +fista.max_iteration = 2000 +fista.update_objective_interval = 500 +fista.run(2000, verbose=True) + +#%% +# Show results +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(fista.get_output().as_array()) +plt.title('Reconstruction') +plt.colorbar() +plt.show() + +plt.plot(np.linspace(0,N,M), data.as_array()[int(N/2),:], label = 'GTruth') +plt.plot(np.linspace(0,N,M), fista.get_output().as_array()[int(N/2),:], label = 'Reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +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: + + ##Construct problem + u1 = Variable(ig.shape) + q = Variable() + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0)) + fidelity = sum(kl_div(noisy_data.as_array(), u1)) + + constraints = [q>=fidelity, u1>=0] + + solver = SCS + obj = Minimize( regulariser + q) + prob = Problem(obj, constraints) + result = prob.solve(verbose = True, solver = solver) + + diff_cvx = numpy.abs( fista.get_output().as_array() - u1.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(fista.get_output().as_array()) + plt.title('FISTA 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,M), fista.get_output().as_array()[int(N/2),:], label = 'FISTA') + plt.plot(np.linspace(0,N,M), u1.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 (FISTA) {} '.format(fista.objective[-1][0])) + + + + + + + diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py new file mode 100644 index 0000000..830ea00 --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Tomo2D_poisson.py @@ -0,0 +1,218 @@ +#======================================================================== +# 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 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, KullbackLeibler, \ + MixedL21Norm, BlockFunction, IndicatorBox + +from ccpi.astra.ops import AstraProjectorSimple +from ccpi.framework import TestData +import os, sys + +""" + +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 + +""" + +loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi')) + +# Load Data +N = 50 +M = 50 +data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1)) + +ig = data.geometry +ag = ig + +#Create Acquisition Data and apply poisson noise + +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.5 +eta = 0 +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 = IndicatorBox(lower=0) + + +# 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) +pdhg.max_iteration = 3000 +pdhg.update_objective_interval = 500 +pdhg.run(3000, 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() +# +# fidelity = sum(kl_div(tmp, ProjMat * u)) +# +# constraints = [q>=fidelity, u>=0] +# solver = SCS +# obj = Minimize( regulariser + q) +# prob = Problem(obj, constraints) +# result = prob.solve(verbose = True, solver = solver) +# +# diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N))) +# +# 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), np.reshape(u.value, (N, N))[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 |