diff options
41 files changed, 3599 insertions, 934 deletions
diff --git a/Wrappers/Python/.DS_Store b/Wrappers/Python/.DS_Store Binary files differnew file mode 100644 index 0000000..5008ddf --- /dev/null +++ b/Wrappers/Python/.DS_Store diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index 386cd87..166014b 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -270,6 +270,7 @@ class BlockDataContainer(object): def squared_norm(self): y = numpy.asarray([el.squared_norm() for el in self.containers]) return y.sum() + def norm(self): return numpy.sqrt(self.squared_norm()) @@ -442,9 +443,12 @@ class BlockDataContainer(object): '''Inline truedivision''' return self.__idiv__(other) + def dot(self, other): +# + tmp = [ self.containers[i].dot(other.containers[i]) for i in range(self.shape[0])] + return sum(tmp) - - + if __name__ == '__main__': @@ -456,6 +460,7 @@ if __name__ == '__main__': BG = BlockGeometry(ig, ig) U = BG.allocate('random_int') + V = BG.allocate('random_int') print ("test sum BDC " ) @@ -468,10 +473,10 @@ if __name__ == '__main__': z1 = sum(U**2).sqrt().as_array() numpy.testing.assert_array_equal(z, z1) - - z2 = U.pnorm(2) + zzz = U.dot(V) + diff --git a/Wrappers/Python/ccpi/framework/BlockGeometry.py b/Wrappers/Python/ccpi/framework/BlockGeometry.py index 5dd6750..ed44d99 100755 --- a/Wrappers/Python/ccpi/framework/BlockGeometry.py +++ b/Wrappers/Python/ccpi/framework/BlockGeometry.py @@ -31,7 +31,50 @@ class BlockGeometry(object): '''returns the Geometry in the BlockGeometry located at position index'''
return self.geometries[index]
- def allocate(self, value=0, dimension_labels=None):
+ def allocate(self, value=0, dimension_labels=None, **kwargs):
+
+ symmetry = kwargs.get('symmetry',False)
containers = [geom.allocate(value) for geom in self.geometries]
+
+ if symmetry == True:
+
+ # for 2x2
+ # [ ig11, ig12\
+ # ig21, ig22]
+
+ # Row-wise Order
+
+ if len(containers)==4:
+ containers[1]=containers[2]
+
+ # for 3x3
+ # [ ig11, ig12, ig13\
+ # ig21, ig22, ig23\
+ # ig31, ig32, ig33]
+
+ elif len(containers)==9:
+ containers[1]=containers[3]
+ containers[2]=containers[6]
+ containers[5]=containers[7]
+
+ # for 4x4
+ # [ ig11, ig12, ig13, ig14\
+ # ig21, ig22, ig23, ig24\
+ # ig31, ig32, ig33, ig34
+ # ig41, ig42, ig43, ig44]
+
+ elif len(containers) == 16:
+ containers[1]=containers[4]
+ containers[2]=containers[8]
+ containers[3]=containers[12]
+ containers[6]=containers[9]
+ containers[7]=containers[10]
+ containers[11]=containers[15]
+
+
+
+
return BlockDataContainer(*containers)
+
+
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py index ed95c3f..c923a30 100755 --- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py @@ -34,7 +34,7 @@ class Algorithm(object): method will stop when the stopping cryterion is met. ''' - def __init__(self): + def __init__(self, **kwargs): '''Constructor Set the minimal number of parameters: @@ -48,7 +48,7 @@ class Algorithm(object): when evaluating the objective is computationally expensive. ''' self.iteration = 0 - self.__max_iteration = 0 + self.__max_iteration = kwargs.get('max_iteration', 0) self.__loss = [] self.memopt = False self.timing = [] @@ -91,9 +91,11 @@ class Algorithm(object): if self.iteration % self.update_objective_interval == 0: self.update_objective() self.iteration += 1 + def get_output(self): '''Returns the solution found''' return self.x + def get_last_loss(self): '''Returns the last stored value of the loss function @@ -145,13 +147,14 @@ class Algorithm(object): if self.should_stop(): print ("Stop cryterion has been reached.") i = 0 + for _ in self: - if verbose and self.iteration % self.update_objective_interval == 0: - print ("Iteration {}/{}, objective {}".format(self.iteration, + if (self.iteration -1) % self.update_objective_interval == 0: + if verbose: + print ("Iteration {}/{}, = {}".format(self.iteration-1, self.max_iteration, self.get_last_objective()) ) - else: if callback is not None: - callback(self.iteration, self.get_last_objective()) + callback(self.iteration -1, self.get_last_objective(), self.x) i += 1 if i == iterations: break diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index d1b5351..39b092b 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -18,93 +18,76 @@ class PDHG(Algorithm): '''Primal Dual Hybrid Gradient''' def __init__(self, **kwargs): - super(PDHG, self).__init__() + super(PDHG, self).__init__(max_iteration=kwargs.get('max_iteration',0)) self.f = kwargs.get('f', None) self.operator = kwargs.get('operator', None) self.g = kwargs.get('g', None) self.tau = kwargs.get('tau', None) self.sigma = kwargs.get('sigma', None) - self.memopt = kwargs.get('memopt', False) - + if self.f is not None and self.operator is not None and \ self.g is not None: print ("Calling from creator") self.set_up(self.f, + self.g, self.operator, - self.g, self.tau, self.sigma) def set_up(self, f, g, operator, tau = None, sigma = None, opt = None, **kwargs): # algorithmic parameters - + self.operator = operator + self.f = f + self.g = g + self.tau = tau + self.sigma = sigma + self.opt = opt if sigma is None and tau is None: raise ValueError('Need sigma*tau||K||^2<1') - - + self.x_old = self.operator.domain_geometry().allocate() - self.y_old = self.operator.range_geometry().allocate() - - self.xbar = self.x_old.copy() - + self.x_tmp = self.x_old.copy() self.x = self.x_old.copy() + + self.y_old = self.operator.range_geometry().allocate() + self.y_tmp = self.y_old.copy() self.y = self.y_old.copy() - if self.memopt: - self.y_tmp = self.y_old.copy() - self.x_tmp = self.x_old.copy() - #y = y_tmp - + + self.xbar = self.x_old.copy() + # relaxation parameter self.theta = 1 def update(self): - if self.memopt: - # Gradient descent, Dual problem solution - # self.y_old += self.sigma * self.operator.direct(self.xbar) - self.operator.direct(self.xbar, out=self.y_tmp) - self.y_tmp *= self.sigma - self.y_old += self.y_tmp - - #self.y = self.f.proximal_conjugate(self.y_old, self.sigma) - self.f.proximal_conjugate(self.y_old, self.sigma, out=self.y) - - # Gradient ascent, Primal problem solution - self.operator.adjoint(self.y, out=self.x_tmp) - self.x_tmp *= self.tau - self.x_old -= self.x_tmp - - self.g.proximal(self.x_old, self.tau, out=self.x) - - #Update - self.x.subtract(self.x_old, out=self.xbar) - self.xbar *= self.theta - self.xbar += self.x + + # Gradient descent, Dual problem solution + self.operator.direct(self.xbar, out=self.y_tmp) + self.y_tmp *= self.sigma + self.y_tmp += self.y_old - self.x_old.fill(self.x) - self.y_old.fill(self.y) + #self.y = self.f.proximal_conjugate(self.y_old, self.sigma) + self.f.proximal_conjugate(self.y_tmp, self.sigma, out=self.y) + + # Gradient ascent, Primal problem solution + self.operator.adjoint(self.y, out=self.x_tmp) + self.x_tmp *= -1*self.tau + self.x_tmp += self.x_old - else: - # Gradient descent, Dual problem solution - self.y_old += self.sigma * self.operator.direct(self.xbar) - self.y = self.f.proximal_conjugate(self.y_old, self.sigma) - # Gradient ascent, Primal problem solution - self.x_old -= self.tau * self.operator.adjoint(self.y) - self.x = self.g.proximal(self.x_old, self.tau) + self.g.proximal(self.x_tmp, self.tau, out=self.x) - #Update - #xbar = x + theta * (x - x_old) - self.xbar.fill(self.x) - self.xbar -= self.x_old - self.xbar *= self.theta - self.xbar += self.x + #Update + self.x.subtract(self.x_old, out=self.xbar) + self.xbar *= self.theta + self.xbar += self.x - self.x_old = self.x - self.y_old = self.y + self.x_old.fill(self.x) + self.y_old.fill(self.y) def update_objective(self): + p1 = self.f(self.operator.direct(self.x)) + self.g(self.x) - d1 = -(self.f.convex_conjugate(self.y) + self.g(-1*self.operator.adjoint(self.y))) + d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y))) self.loss.append([p1,d1,p1-d1]) @@ -148,63 +131,44 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): for i in range(niter): + - if not memopt: - - y_old += sigma * operator.direct(xbar) - y = f.proximal_conjugate(y_old, sigma) - - x_old -= tau*operator.adjoint(y) - x = g.proximal(x_old, tau) - - x.subtract(x_old, out=xbar) - xbar *= theta - xbar += x - - x_old.fill(x) - y_old.fill(y) - - -# if i%100==0: -# -# p1 = f(operator.direct(x)) + g(x) -# d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) ) -# primal.append(p1) -# dual.append(d1) -# pdgap.append(p1-d1) -# print(p1, d1, p1-d1) - - - else: - + if memopt: operator.direct(xbar, out = y_tmp) y_tmp *= sigma - y_tmp += y_old - f.proximal_conjugate(y_tmp, sigma, out=y) - + y_tmp += y_old + else: + y_tmp = y_old + sigma * operator.direct(xbar) + + f.proximal_conjugate(y_tmp, sigma, out=y) + + if memopt: operator.adjoint(y, out = x_tmp) - x_tmp *= -tau + x_tmp *= -1*tau x_tmp += x_old - - g.proximal(x_tmp, tau, out = x) + else: + x_tmp = x_old - tau*operator.adjoint(y) - x.subtract(x_old, out=xbar) - xbar *= theta - xbar += x - - x_old.fill(x) - y_old.fill(y) + g.proximal(x_tmp, tau, out=x) + + x.subtract(x_old, out=xbar) + xbar *= theta + xbar += x + + x_old.fill(x) + y_old.fill(y) + + if i%10==0: -# if i%100==0: -# -# p1 = f(operator.direct(x)) + g(x) -# d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) ) -# primal.append(p1) -# dual.append(d1) -# pdgap.append(p1-d1) -# print(p1, d1, p1-d1) + p1 = f(operator.direct(x)) + g(x) + d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) ) + primal.append(p1) + dual.append(d1) + pdgap.append(p1-d1) + print(p1, d1, p1-d1) + t_end = time.time() diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py index bf627a5..0836324 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py @@ -93,16 +93,28 @@ class BlockFunction(Function): ''' + + if out is None: - out = [None]*self.length - if isinstance(tau, Number): - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau) + out = [None]*self.length + if isinstance(tau, Number): + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau) + else: + for i in range(self.length): + out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i)) + + return BlockDataContainer(*out) + else: - for i in range(self.length): - out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i)) - - return BlockDataContainer(*out) + if isinstance(tau, Number): + for i in range(self.length): + self.functions[i].proximal(x.get_item(i), tau, out[i]) + else: + for i in range(self.length): + self.functions[i].proximal(x.get_item(i), tau.get_item(i), out[i]) + + def gradient(self,x, out=None): diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py index d4370a8..36ba146 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py @@ -38,22 +38,15 @@ class KullbackLeibler(Function): def __call__(self, x): + + res_tmp = numpy.zeros(x.shape) - # TODO check + tmp = x + self.bnoise + ind = x.as_array()>0 - self.sum_value = x + self.bnoise - if (self.sum_value.as_array()<0).any(): - self.sum_value = numpy.inf + res_tmp[ind] = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind]) - if self.sum_value==numpy.inf: - return numpy.inf - else: - tmp = self.sum_value.copy() - #tmp.fill( numpy.log(tmp.as_array()) ) - self.log(tmp) - return (x - self.b * tmp ).sum() - -# return numpy.sum( x.as_array() - self.b.as_array() * numpy.log(self.sum_value.as_array())) + return res_tmp.sum() def log(self, datacontainer): '''calculates the in-place log of the datacontainer''' @@ -61,6 +54,7 @@ class KullbackLeibler(Function): datacontainer.as_array().ravel(), True): raise ValueError('KullbackLeibler. Cannot calculate log of negative number') datacontainer.fill( numpy.log(datacontainer.as_array()) ) + def gradient(self, x, out=None): @@ -68,6 +62,7 @@ class KullbackLeibler(Function): if out is None: return 1 - self.b/(x + self.bnoise) else: + x.add(self.bnoise, out=out) self.b.divide(out, out=out) out.subtract(1, out=out) @@ -75,16 +70,18 @@ class KullbackLeibler(Function): def convex_conjugate(self, x): - tmp = self.b/( 1 - x ) - self.log(tmp) - return (self.b * ( tmp - 1 ) - self.bnoise * (x - 1)).sum() -# return self.b * ( ImageData(numpy.log(self.b/(1-x)) - 1 )) - self.bnoise * (x - 1) + tmp = self.b/(1-x) + ind = tmp.as_array()>0 + + return (self.b.as_array()[ind] * (numpy.log(tmp.as_array()[ind])-1)).sum() + def proximal(self, x, tau, out=None): if out is None: return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() ) else: + tmp = 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b ) .sqrt() ) @@ -101,26 +98,26 @@ class KullbackLeibler(Function): out += tmp out *= 0.5 - - - + def proximal_conjugate(self, x, tau, out=None): if out is None: z = x + tau * self.bnoise - return (z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt() + return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) else: - z_m = x + tau * self.bnoise - 1 + + 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) - # z = z_m + 2 - z_m.sqrt(out=z_m) + z_m.sqrt(out=z_m) z_m += 2 out *= -1 out += z_m + out *= 0.5 + def __rmul__(self, scalar): @@ -137,19 +134,20 @@ class KullbackLeibler(Function): if __name__ == '__main__': + + + from ccpi.framework import ImageGeometry + import numpy N, M = 2,3 ig = ImageGeometry(N, M) - data = ImageData(numpy.random.randint(-10, 100, size=(M, N))) - x = ImageData(numpy.random.randint(-10, 100, size=(M, N))) + data = ImageData(numpy.random.randint(-10, 10, size=(M, N))) + x = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - bnoise = ImageData(numpy.random.randint(-100, 100, size=(M, N))) + bnoise = ImageData(numpy.random.randint(-10, 10, size=(M, N))) - f = KullbackLeibler(data, bnoise=bnoise) - print(f.sum_value) - - print(f(x)) - + f = KullbackLeibler(data) + print(f(x)) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 2bb4ca7..8740376 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -71,7 +71,7 @@ class L2NormSquared(Function): def convex_conjugate(self, x): ''' Evaluate convex conjugate of L2NormSquared at x: f^{*}(x)''' - + tmp = 0 if self.b is not None: @@ -94,28 +94,19 @@ class L2NormSquared(Function): return x/(1+2*tau) else: tmp = x.subtract(self.b) - #tmp -= self.b tmp /= (1+2*tau) tmp += self.b return tmp # return (x-self.b)/(1+2*tau) + self.b - -# if self.b is not None: -# out=x -# if self.b is not None: -# out -= self.b -# out /= (1+2*tau) -# if self.b is not None: -# out += self.b -# return out + + else: - out.fill(x) - if self.b is not None: - out -= self.b - out /= (1+2*tau) if self.b is not None: - out += self.b - + x.subtract(self.b, out=out) + out /= (1+2*tau) + out += self.b + else: + x.divide((1+2*tau), out=out) def proximal_conjugate(self, x, tau, out=None): @@ -145,7 +136,13 @@ class L2NormSquared(Function): ''' - return ScaledFunction(self, scalar) + return ScaledFunction(self, scalar) + + + def operator_composition(self, operator): + + return FunctionOperatorComposition + if __name__ == '__main__': @@ -287,17 +284,3 @@ if __name__ == '__main__': numpy.testing.assert_array_almost_equal(res1.as_array(), \ res2.as_array(), decimal=4) - - - - - - - - - - - - - - diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index 7e6b6e7..e8f6da4 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py @@ -43,19 +43,9 @@ class MixedL21Norm(Function): ''' if not isinstance(x, BlockDataContainer): raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) - - if self.SymTensor: - - #TODO fix this case - param = [1]*x.shape[0] - param[-1] = 2 - tmp = [param[i]*(x[i] ** 2) for i in range(x.shape[0])] - res = sum(tmp).sqrt().sum() - - else: - - tmp = [ el**2 for el in x.containers ] - res = sum(tmp).sqrt().sum() + + tmp = [ el**2 for el in x.containers ] + res = sum(tmp).sqrt().sum() return res @@ -67,7 +57,12 @@ class MixedL21Norm(Function): ''' This is the Indicator function of ||\cdot||_{2, \infty} which is either 0 if ||x||_{2, \infty} or \infty ''' + return 0.0 + + #tmp = [ el**2 for el in x.containers ] + #print(sum(tmp).sqrt().as_array().max()) + #return sum(tmp).sqrt().as_array().max() def proximal(self, x, tau, out=None): @@ -80,34 +75,24 @@ class MixedL21Norm(Function): def proximal_conjugate(self, x, tau, out=None): - if self.SymTensor: - - param = [1]*x.shape[0] - param[-1] = 2 - tmp = [param[i]*(x[i] ** 2) for i in range(x.shape[0])] - frac = [x[i]/(sum(tmp).sqrt()).maximum(1.0) for i in range(x.shape[0])] - res = BlockDataContainer(*frac) - - return res - - else: - if out is None: - tmp = [ el*el for el in x.containers] - res = sum(tmp).sqrt().maximum(1.0) - frac = [el/res for el in x.containers] - return BlockDataContainer(*frac) - - #TODO this is slow, why??? + if out is None: + tmp = [ el*el for el in x.containers] + res = sum(tmp).sqrt().maximum(1.0) + frac = [el/res for el in x.containers] + return BlockDataContainer(*frac) + + + #TODO this is slow, why??? # return x.divide(x.pnorm().maximum(1.0)) - else: - - res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ) - res = res1.sqrt().maximum(1.0) - x.divide(res, out=out) - + else: + + res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ) + res = res1.sqrt().maximum(1.0) + x.divide(res, out=out) + # x.divide(sum([el*el for el in x.containers]).sqrt().maximum(1.0), out=out) - #TODO this is slow, why ??? + #TODO this is slow, why ??? # x.divide(x.pnorm().maximum(1.0), out=out) diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py index cce519a..a019815 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py +++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py @@ -39,14 +39,8 @@ class ZeroFunction(Function): indicator function for the set = {0} So 0 if x=0, or inf if x neq 0 ''' + return x.maximum(0).sum() - if x.shape[0]==1: - return x.maximum(0).sum() - else: - if isinstance(x, BlockDataContainer): - return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum() - else: - return x.maximum(0).sum() + x.maximum(0).sum() def proximal(self, x, tau, out=None): if out is None: diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 1d77510..c8bd546 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -266,15 +266,30 @@ class BlockOperator(Operator): # column BlockOperator return self.get_item(0,0).domain_geometry() else: - shape = (self.shape[0], 1) - return BlockGeometry(*[el.domain_geometry() for el in self.operators], - shape=shape) + # get the geometries column wise + # we need only the geometries from the first row + # since it is compatible from __init__ + tmp = [] + for i in range(self.shape[1]): + tmp.append(self.get_item(0,i).domain_geometry()) + return BlockGeometry(*tmp) + + #shape = (self.shape[0], 1) + #return BlockGeometry(*[el.domain_geometry() for el in self.operators], + # shape=self.shape) def range_geometry(self): '''returns the range of the BlockOperator''' - shape = (self.shape[1], 1) - return BlockGeometry(*[el.range_geometry() for el in self.operators], - shape=shape) + + tmp = [] + for i in range(self.shape[0]): + tmp.append(self.get_item(i,0).range_geometry()) + return BlockGeometry(*tmp) + + + #shape = (self.shape[1], 1) + #return BlockGeometry(*[el.range_geometry() for el in self.operators], + # shape=shape) def sum_abs_row(self): @@ -312,7 +327,8 @@ class BlockOperator(Operator): if __name__ == '__main__': from ccpi.framework import ImageGeometry - from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff + from ccpi.optimisation.operators import Gradient, Identity, \ + SparseFiniteDiff, SymmetrizedGradient, ZeroOperator M, N = 4, 3 @@ -363,4 +379,39 @@ if __name__ == '__main__': + ########################################################################### + # Block Operator for TGV reconstruction + + M, N = 2,3 + ig = ImageGeometry(M, N) + ag = ig + + 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) ) + + z1 = operator.domain_geometry() + z2 = operator.range_geometry() + + print(z1.shape) + print(z2.shape) + + + + + + + + + +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py index 2c42532..f459334 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py @@ -41,7 +41,7 @@ class FiniteDiff(LinearOperator): #self.voxel_size = kwargs.get('voxel_size',1) # this wrongly assumes a homogeneous voxel size - self.voxel_size = self.gm_domain.voxel_size_x +# self.voxel_size = self.gm_domain.voxel_size_x def direct(self, x, out=None): diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py index 3fc43b8..243809b 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py @@ -14,6 +14,12 @@ from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff class SymmetrizedGradient(Gradient): + ''' Symmetrized Gradient, denoted by E: V --> W + where V is the Range of the Gradient Operator + and W is the Range of the Symmetrized Gradient. + ''' + + def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs): super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) @@ -21,165 +27,218 @@ class SymmetrizedGradient(Gradient): ''' Domain of SymGrad is the Range of Gradient ''' + self.gm_domain = self.gm_range self.bnd_cond = bnd_cond self.channels = self.gm_range.get_item(0).channels - if self.correlation=='Space': - if self.channels>1: - pass - else: -# # 2D image ---> Dx v1, Dyv2, Dx - tmp = self.gm_domain.geometries + (self.gm_domain.get_item(0),) - self.gm_range = BlockGeometry(*tmp ) - self.ind1 = range(self.gm_domain.get_item(0).length) - self.ind2 = range(self.gm_domain.get_item(0).length-1, -1, -1) -# self.order = myorder = [0,1,2 3] - - elif self.correlation=='SpaceChannels': - if self.channels>1: - pass - else: - raise ValueError('No channels to correlate') - - - def direct(self, x, out=None): + tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries -# tmp = numpy.zeros(self.gm_range) -# tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) -# tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) -# tmp[2] = 0.5 * (FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) + -# FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) ) -# -# return type(x)(tmp) - - tmp = [[None]*2]*2 - for i in range(2): - for j in range(2): - tmp[i][j]=FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).adjoint(x.get_item(j)) - tmp = numpy.array(tmp) - z = 0.5 * (tmp.T + tmp) - z = z.to + self.gm_range = BlockGeometry(*tmp_gm) - return BlockDataContainer(*z.tolist()) - - - def adjoint(self, x, out=None): - pass + self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond) - res = [] - for i in range(2): - tmp = ImageData(np.zeros(x.get_item(0))) - for j in range(2): - tmp += FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).direct(x.get_item(j)) - res.append(tmp) - return res - + if self.gm_domain.shape[0]==2: + self.order_ind = [0,2,1,3] + else: + self.order_ind = [0,3,6,1,4,7,2,5,8] -# for + def direct(self, x, out=None): + + if out is None: + + tmp = [] + for i in range(self.gm_domain.shape[0]): + for j in range(x.shape[0]): + self.FD.direction = i + tmp.append(self.FD.adjoint(x.get_item(j))) + + tmp1 = [tmp[i] for i in self.order_ind] + + res = [0.5 * sum(x) for x in zip(tmp, tmp1)] + + return BlockDataContainer(*res) + + else: + + ind = 0 + for i in range(self.gm_domain.shape[0]): + for j in range(x.shape[0]): + self.FD.direction = i + self.FD.adjoint(x.get_item(j), out=out[ind]) + ind+=1 + out1 = BlockDataContainer(*[out[i] for i in self.order_ind]) + out.fill( 0.5 * (out + out1) ) + + + def adjoint(self, x, out=None): -# tmp = numpy.zeros(self.gm_domain) -# -# tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) + \ -# FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) -# -# tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) + \ -# FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1]) -# -# return type(x)(tmp) + if out is None: - def alloc_domain_dim(self): - return ImageData(numpy.zeros(self.gm_domain)) - - def alloc_range_dim(self): - return ImageData(numpy.zeros(self.range_dim)) - - def domain_dim(self): + tmp = [None]*self.gm_domain.shape[0] + i = 0 + + for k in range(self.gm_domain.shape[0]): + tmp1 = 0 + for j in range(self.gm_domain.shape[0]): + self.FD.direction = j + tmp1 += self.FD.direct(x[i]) + i+=1 + tmp[k] = tmp1 + return BlockDataContainer(*tmp) + + + else: + + tmp = self.gm_domain.allocate() + i = 0 + for k in range(self.gm_domain.shape[0]): + tmp1 = 0 + for j in range(self.gm_domain.shape[0]): + self.FD.direction = j + self.FD.direct(x[i], out=tmp[j]) + i+=1 + tmp1+=tmp[j] + out[k].fill(tmp1) +# tmp = self.adjoint(x) +# out.fill(tmp) + + + def domain_geometry(self): return self.gm_domain - def range_dim(self): + def range_geometry(self): return self.gm_range def norm(self): -# return np.sqrt(4*len(self.domainDim())) - #TODO this takes time for big ImageData - # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12) - x0 = ImageData(np.random.random_sample(self.domain_dim())) - self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0) - return self.s1 + +# TODO need dot method for BlockDataContainer +# return numpy.sqrt(4*self.gm_domain.shape[0]) + +# x0 = self.gm_domain.allocate('random') + self.s1, sall, svec = LinearOperator.PowerMethod(self, 50) + return self.s1 if __name__ == '__main__': ########################################################################### - ## Symmetrized Gradient + ## Symmetrized Gradient Tests from ccpi.framework import DataContainer from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff import numpy as np N, M = 2, 3 K = 2 + C = 2 ig1 = ImageGeometry(N, M) - ig2 = ImageGeometry(N, M, channels=K) + ig2 = ImageGeometry(N, M, channels=C) E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann') - E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Periodic') - print(E1.domain_geometry().shape) - print(E2.domain_geometry().shape) + try: + E1 = SymmetrizedGradient(ig1, correlation = 'SpaceChannels', bnd_cond='Neumann') + except: + print("No Channels to correlate") + + E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Neumann') + + print(E1.domain_geometry().shape, E1.range_geometry().shape) + print(E2.domain_geometry().shape, E2.range_geometry().shape) + + #check Linear operator property + + u1 = E1.domain_geometry().allocate('random_int') + u2 = E2.domain_geometry().allocate('random_int') + + # Need to allocate random_int at the Range of SymGradient + + #a1 = ig1.allocate('random_int') + #a2 = ig1.allocate('random_int') + #a3 = ig1.allocate('random_int') + + #a4 = ig1.allocate('random_int') + #a5 = ig1.allocate('random_int') + #a6 = ig1.allocate('random_int') + + # TODO allocate has to create this symmetry by default!!!!! + #w1 = BlockDataContainer(*[a1, a2, \ + # a2, a3]) + w1 = E1.range_geometry().allocate('random_int',symmetry=True) + + LHS = (E1.direct(u1) * w1).sum() + RHS = (u1 * E1.adjoint(w1)).sum() + + numpy.testing.assert_equal(LHS, RHS) - u1 = E1.gm_domain.allocate('random_int') u2 = E2.gm_domain.allocate('random_int') - - res = E1.direct(u1) + #aa1 = ig2.allocate('random_int') + #aa2 = ig2.allocate('random_int') + #aa3 = ig2.allocate('random_int') + #aa4 = ig2.allocate('random_int') + #aa5 = ig2.allocate('random_int') + #aa6 = ig2.allocate('random_int') - res1 = E1.adjoint(res) + #w2 = BlockDataContainer(*[aa1, aa2, aa3, \ + # aa2, aa4, aa5, \ + # aa3, aa5, aa6]) + w2 = E2.range_geometry().allocate('random_int',symmetry=True) + -# Dx = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann') -# Dy = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann') -# -# B = BlockOperator(Dy, Dx) -# V = BlockDataContainer(u1,u2) -# -# res = B.adjoint(V) + LHS1 = (E2.direct(u2) * w2).sum() + RHS1 = (u2 * E2.adjoint(w2)).sum() -# ig = (N,M) -# ig2 = (2,) + ig -# ig3 = (3,) + ig -# u1 = ig.allocate('random_int') -# w1 = E.gm_range.allocate('random_int') -# DataContainer(np.random.randint(10, size=ig3)) + numpy.testing.assert_equal(LHS1, RHS1) + out = E1.range_geometry().allocate() + E1.direct(u1, out=out) + a1 = E1.direct(u1) + numpy.testing.assert_array_equal(a1[0].as_array(), out[0].as_array()) + numpy.testing.assert_array_equal(a1[1].as_array(), out[1].as_array()) + numpy.testing.assert_array_equal(a1[2].as_array(), out[2].as_array()) + numpy.testing.assert_array_equal(a1[3].as_array(), out[3].as_array()) -# d1 = E.direct(u1) -# d2 = E.adjoint(w1) + out1 = E1.domain_geometry().allocate() + E1.adjoint(w1, out=out1) + b1 = E1.adjoint(w1) -# LHS = (d1.as_array()[0]*w1.as_array()[0] + \ -# d1.as_array()[1]*w1.as_array()[1] + \ -# 2*d1.as_array()[2]*w1.as_array()[2]).sum() -# -# RHS = (u1.as_array()[0]*d2.as_array()[0] + \ -# u1.as_array()[1]*d2.as_array()[1]).sum() -# -# -# print(LHS, RHS, E.norm()) -# + LHS_out = (out * w1).sum() + RHS_out = (u1 * out1).sum() + print(LHS_out, RHS_out) -# + out2 = E2.range_geometry().allocate() + E2.direct(u2, out=out2) + a2 = E2.direct(u2) + + out21 = E2.domain_geometry().allocate() + E2.adjoint(w2, out=out21) + b2 = E2.adjoint(w2) + LHS_out = (out2 * w2).sum() + RHS_out = (u2 * out21).sum() + print(LHS_out, RHS_out) + out = E1.range_geometry().allocate() + E1.direct(u1, out=out) + E1.adjoint(out, out=out1) + print(E1.norm()) + print(E2.norm()) -
\ No newline at end of file + +# +# +# +#
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py index a7c5f09..8168f0b 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py @@ -8,32 +8,37 @@ Created on Wed Mar 6 19:25:53 2019 import numpy as np from ccpi.framework import ImageData -from ccpi.optimisation.operators import Operator +from ccpi.optimisation.operators import LinearOperator -class ZeroOp(Operator): +class ZeroOperator(LinearOperator): - def __init__(self, gm_domain, gm_range): + def __init__(self, gm_domain, gm_range=None): + + super(ZeroOperator, self).__init__() + self.gm_domain = gm_domain - self.gm_range = gm_range - super(ZeroOp, self).__init__() + self.gm_range = gm_range + if self.gm_range is None: + self.gm_range = self.gm_domain + def direct(self,x,out=None): if out is None: - return ImageData(np.zeros(self.gm_range)) + return self.gm_range.allocate() else: - return ImageData(np.zeros(self.gm_range)) + out.fill(self.gm_range.allocate()) def adjoint(self,x, out=None): if out is None: - return ImageData(np.zeros(self.gm_domain)) + return self.gm_domain.allocate() else: - return ImageData(np.zeros(self.gm_domain)) + out.fill(self.gm_domain.allocate()) def norm(self): return 0 - def domain_dim(self): + def domain_geometry(self): return self.gm_domain - def range_dim(self): + def range_geometry(self): return self.gm_range
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py index 811adf6..23222d4 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py @@ -18,6 +18,6 @@ from .FiniteDifferenceOperator import FiniteDiff from .GradientOperator import Gradient
from .SymmetrizedGradientOperator import SymmetrizedGradient
from .IdentityOperator import Identity
-from .ZeroOperator import ZeroOp
+from .ZeroOperator import ZeroOperator
from .LinearOperatorMatrix import LinearOperatorMatrix
diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py index ccef410..936dc05 100755 --- a/Wrappers/Python/ccpi/processors.py +++ b/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py @@ -19,115 +19,9 @@ from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\ AcquisitionGeometry, ImageGeometry, ImageData -from ccpi.reconstruction.parallelbeam import alg as pbalg import numpy from scipy import ndimage -import matplotlib.pyplot as plt - - -class Normalizer(DataProcessor): - '''Normalization based on flat and dark - - This processor read in a AcquisitionData and normalises it based on - the instrument reading with and without incident photons or neutrons. - - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSetn - ''' - - def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5): - kwargs = { - 'flat_field' : flat_field, - 'dark_field' : dark_field, - # very small number. Used when there is a division by zero - 'tolerance' : tolerance - } - - #DataProcessor.__init__(self, **kwargs) - super(Normalizer, self).__init__(**kwargs) - if not flat_field is None: - self.set_flat_field(flat_field) - if not dark_field is None: - self.set_dark_field(dark_field) - - def check_input(self, dataset): - if dataset.number_of_dimensions == 3 or\ - dataset.number_of_dimensions == 2: - return True - else: - raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - .format(dataset.number_of_dimensions)) - - def set_dark_field(self, df): - if type(df) is numpy.ndarray: - if len(numpy.shape(df)) == 3: - raise ValueError('Dark Field should be 2D') - elif len(numpy.shape(df)) == 2: - self.dark_field = df - elif issubclass(type(df), DataContainer): - self.dark_field = self.set_dark_field(df.as_array()) - - def set_flat_field(self, df): - if type(df) is numpy.ndarray: - if len(numpy.shape(df)) == 3: - raise ValueError('Flat Field should be 2D') - elif len(numpy.shape(df)) == 2: - self.flat_field = df - elif issubclass(type(df), DataContainer): - self.flat_field = self.set_flat_field(df.as_array()) - - @staticmethod - def normalize_projection(projection, flat, dark, tolerance): - a = (projection - dark) - b = (flat-dark) - with numpy.errstate(divide='ignore', invalid='ignore'): - c = numpy.true_divide( a, b ) - c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0 - return c - - @staticmethod - def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark): - '''returns the estimated relative error of the normalised projection - - n = (projection - dark) / (flat - dark) - Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d) - ''' - a = (projection - dark) - b = (flat-dark) - df = delta_flat / flat - dd = delta_dark / dark - rel_norm_error = (b + a) / (b * a) * (df + dd) - return rel_norm_error - - def process(self, out=None): - - projections = self.get_input() - dark = self.dark_field - flat = self.flat_field - - if projections.number_of_dimensions == 3: - if not (projections.shape[1:] == dark.shape and \ - projections.shape[1:] == flat.shape): - raise ValueError('Flats/Dark and projections size do not match.') - - - a = numpy.asarray( - [ Normalizer.normalize_projection( - projection, flat, dark, self.tolerance) \ - for projection in projections.as_array() ] - ) - elif projections.number_of_dimensions == 2: - a = Normalizer.normalize_projection(projections.as_array(), - flat, dark, self.tolerance) - y = type(projections)( a , True, - dimension_labels=projections.dimension_labels, - geometry=projections.geometry) - return y - - class CenterOfRotationFinder(DataProcessor): '''Processor to find the center of rotation in a parallel beam experiment diff --git a/Wrappers/Python/ccpi/processors/Normalizer.py b/Wrappers/Python/ccpi/processors/Normalizer.py new file mode 100755 index 0000000..da65e5c --- /dev/null +++ b/Wrappers/Python/ccpi/processors/Normalizer.py @@ -0,0 +1,124 @@ +# -*- 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 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 DataProcessor, DataContainer, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData +import numpy + +class Normalizer(DataProcessor): + '''Normalization based on flat and dark + + This processor read in a AcquisitionData and normalises it based on + the instrument reading with and without incident photons or neutrons. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSetn + ''' + + def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5): + kwargs = { + 'flat_field' : flat_field, + 'dark_field' : dark_field, + # very small number. Used when there is a division by zero + 'tolerance' : tolerance + } + + #DataProcessor.__init__(self, **kwargs) + super(Normalizer, self).__init__(**kwargs) + if not flat_field is None: + self.set_flat_field(flat_field) + if not dark_field is None: + self.set_dark_field(dark_field) + + def check_input(self, dataset): + if dataset.number_of_dimensions == 3 or\ + dataset.number_of_dimensions == 2: + return True + else: + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) + + def set_dark_field(self, df): + if type(df) is numpy.ndarray: + if len(numpy.shape(df)) == 3: + raise ValueError('Dark Field should be 2D') + elif len(numpy.shape(df)) == 2: + self.dark_field = df + elif issubclass(type(df), DataContainer): + self.dark_field = self.set_dark_field(df.as_array()) + + def set_flat_field(self, df): + if type(df) is numpy.ndarray: + if len(numpy.shape(df)) == 3: + raise ValueError('Flat Field should be 2D') + elif len(numpy.shape(df)) == 2: + self.flat_field = df + elif issubclass(type(df), DataContainer): + self.flat_field = self.set_flat_field(df.as_array()) + + @staticmethod + def normalize_projection(projection, flat, dark, tolerance): + a = (projection - dark) + b = (flat-dark) + with numpy.errstate(divide='ignore', invalid='ignore'): + c = numpy.true_divide( a, b ) + c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0 + return c + + @staticmethod + def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark): + '''returns the estimated relative error of the normalised projection + + n = (projection - dark) / (flat - dark) + Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d) + ''' + a = (projection - dark) + b = (flat-dark) + df = delta_flat / flat + dd = delta_dark / dark + rel_norm_error = (b + a) / (b * a) * (df + dd) + return rel_norm_error + + def process(self, out=None): + + projections = self.get_input() + dark = self.dark_field + flat = self.flat_field + + if projections.number_of_dimensions == 3: + if not (projections.shape[1:] == dark.shape and \ + projections.shape[1:] == flat.shape): + raise ValueError('Flats/Dark and projections size do not match.') + + + a = numpy.asarray( + [ Normalizer.normalize_projection( + projection, flat, dark, self.tolerance) \ + for projection in projections.as_array() ] + ) + elif projections.number_of_dimensions == 2: + a = Normalizer.normalize_projection(projections.as_array(), + flat, dark, self.tolerance) + y = type(projections)( a , True, + dimension_labels=projections.dimension_labels, + geometry=projections.geometry) + return y +
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/processors/__init__.py b/Wrappers/Python/ccpi/processors/__init__.py new file mode 100755 index 0000000..f8d050e --- /dev/null +++ b/Wrappers/Python/ccpi/processors/__init__.py @@ -0,0 +1,9 @@ +# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr 30 13:51:09 2019
+
+@author: ofn77899
+"""
+
+from .CenterOfRotationFinder import CenterOfRotationFinder
+from .Normalizer import Normalizer
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index dd3238e..ac5f763 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -26,7 +26,6 @@ requirements: build: - {{ pin_compatible('numpy', max_pin='x.x') }} - python - - numpy - setuptools run: @@ -34,7 +33,6 @@ requirements: - python - numpy - scipy - #- matplotlib - h5py about: diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py new file mode 100644 index 0000000..854f645 --- /dev/null +++ b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py @@ -0,0 +1,238 @@ +# -*- 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.plugins.operators import CCPiProjectorSimple +from timeit import default_timer as timer +from ccpi.reconstruction.parallelbeam import alg as pbalg +import os + +try: + import tomophantom + from tomophantom import TomoP3D + no_tomophantom = False +except ImportError as ie: + no_tomophantom = True + +#%% + +#%%############################################################################### +# 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)) + +vert = 4 +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) + +angles_num = 100 +det_w = 1.0 +det_num = N + +angles = np.linspace(-90.,90.,N, dtype=np.float32) +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +#no_tomophantom = True +if no_tomophantom: + data = ig.allocate() + Phantom = data + # Populate image data by looping over and filling slices + i = 0 + while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + 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)] = 0.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + + Aop = CCPiProjectorSimple(ig, ag, 'cpu') + sin = Aop.direct(data) +else: + + model = 13 # select a model number from the library + N_size = N # 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) + + #%% + Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal) + Vert_det = N_size # detector row count (vertical) (no reason for it to be > N) + #angles_num = int(0.5*np.pi*N_size); # angles number + #angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees + + print ("Building 3D analytical projection data with TomoPhantom") + projData3D_analyt = TomoP3D.ModelSino(model, + N_size, + Horiz_det, + Vert_det, + angles, + path_library3D) + + # tomophantom outputs in [vert,angles,horiz] + # we want [angle,vert,horiz] + data = np.transpose(projData3D_analyt, [1,0,2]) + ag.pixel_num_h = Horiz_det + ag.pixel_num_v = Vert_det + sin = ag.allocate() + sin.fill(data) + ig.voxel_num_y = Vert_det + + Aop = CCPiProjectorSimple(ig, ag, 'cpu') + + +plt.imshow(sin.subset(vertical=0).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.subset(vertical=0).as_array()) +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() + +normK = Aop.norm() + +# Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +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: + sigma = 1 + tau = 1/(sigma*normK**2) + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) +niter = 50 +opt = {'niter':niter} +opt1 = {'niter':niter, 'memopt': True} + + + +pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, max_iteration=niter) +#pdhg1.max_iteration = 2000 +pdhg1.update_objective_interval = 100 + +t1_old = timer() +resold, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) +t2_old = timer() + +pdhg1.run(niter) +print (sum(pdhg1.timing)) +res = pdhg1.get_output().subset(vertical=0) + +#%% +plt.figure() +plt.subplot(1,4,1) +plt.imshow(res.as_array()) +plt.title('Algorithm') +plt.colorbar() +plt.subplot(1,4,2) +plt.imshow(resold.subset(vertical=0).as_array()) +plt.title('function') +plt.colorbar() +plt.subplot(1,4,3) +plt.imshow((res - resold.subset(vertical=0)).abs().as_array()) +plt.title('diff') +plt.colorbar() +plt.subplot(1,4,4) +plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'Algorithm') +plt.plot(np.linspace(0,N,N), resold.subset(vertical=0).as_array()[int(N/2),:], label = 'function') +plt.legend() +plt.show() +# +print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(sum(pdhg1.timing), t2_old -t1_old)) +diff = (res - resold.subset(vertical=0)).abs().as_array().max() +# +print(" Max of abs difference is {}".format(diff)) + diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index a3fde59..95c0dea 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -36,6 +36,7 @@ setup( 'ccpi.optimisation.operators', 'ccpi.optimisation.algorithms', 'ccpi.optimisation.functions', + 'ccpi.processors', 'ccpi.contrib','ccpi.contrib.optimisation', 'ccpi.contrib.optimisation.algorithms'], diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 1c1de3a..3e6a83e 100755 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -11,8 +11,32 @@ from timeit import default_timer as timer from ccpi.framework import AX, CastDataContainer, PixelByPixelDataProcessor
+from ccpi.io.reader import NexusReader
+from ccpi.processors import CenterOfRotationFinder
+import wget
+import os
+
class TestDataProcessor(unittest.TestCase):
+ def setUp(self):
+ wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs')
+ self.filename = '24737_fd.nxs'
+
+ def tearDown(self):
+ os.remove(self.filename)
+ def test_CenterOfRotation(self):
+ reader = NexusReader(self.filename)
+ ad = reader.get_acquisition_data_whole()
+ print (ad.geometry)
+ cf = CenterOfRotationFinder()
+ cf.set_input(ad)
+ print ("Center of rotation", cf.get_output())
+ self.assertAlmostEqual(86.25, cf.get_output())
+ def test_Normalizer(self):
+ pass
+
+
+
def test_DataProcessorChaining(self):
shape = (2,3,4,5)
size = shape[0]
diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py new file mode 100644 index 0000000..7b65c31 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Denoising_SaltPepper.py @@ -0,0 +1,194 @@ +#!/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 + +from ccpi.optimisation.operators import BlockOperator, Identity, \ + Gradient, SymmetrizedGradient, ZeroOperator +from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TGV SaltPepper 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()) + +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Add Gaussian noise +n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameters +alpha = 0.8 +beta = numpy.sqrt(2)* alpha + +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 = L1Norm(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(L1Norm(b=noisy_data), 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 = 50 +pdhg.run(2000) + +#%% +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()[0].as_array()) +plt.title('TGV 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()[0].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 + +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 = pnorm(u - noisy_data.as_array(),1) + solver = MOSEK + + # 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( pdhg.get_output()[0].as_array() - u.value ) + + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(pdhg.get_output()[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), pdhg.get_output()[0].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])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py new file mode 100644 index 0000000..26578bb --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py @@ -0,0 +1,124 @@ +# -*- 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, Identity, \ + SymmetrizedGradient, ZeroOperator +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +N = 75 + +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()) + +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) + +# Create noisy data. Apply Poisson noise +scale = 0.1 +np.random.seed(5) +n1 = scale * np.random.poisson(sin.as_array()/scale) +noisy_data = AcquisitionData(n1, ag) + + +plt.imshow(noisy_data.as_array()) +plt.show() +#%% +# Regularisation Parameters +alpha = 0.7 +beta = 2 + +# 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 = KullbackLeibler(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) + + +# 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) + +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()[0].as_array()) +plt.title('TGV 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()[0].as_array()[int(N/2),:], label = 'TGV reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py new file mode 100644 index 0000000..1a3e0df --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py @@ -0,0 +1,177 @@ +# -*- 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 + +import numpy as np +import numpy +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, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Gaussian denoising +N = 300 + +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 +n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 0.5 + +method = '0' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + 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: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) + + +# 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 = 50 +pdhg.run(2000) + + +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 + +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( pdhg.get_output().as_array() - u.value ) + + 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(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), 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])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py new file mode 100644 index 0000000..c86ddc9 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py @@ -0,0 +1,155 @@ +# -*- 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 + +import numpy as np +import numpy +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, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Gaussian 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 = 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 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 +n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10) +noisy_data = ImageData(n1) + +sliceSel = int(0.5*N) +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1) +plt.title('Axial View') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1) +plt.title('Coronal View') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1) +plt.title('Sagittal View') +plt.colorbar() +plt.show() + + +# Regularisation Parameter +alpha = 0.05 + +method = '0' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + 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: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) + + +# 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(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Axial View') + +plt.subplot(2,3,2) +plt.imshow(noisy_data.as_array()[:,sliceSel,:],vmin=0, vmax=1) +plt.axis('off') +plt.title('Coronal View') + +plt.subplot(2,3,3) +plt.imshow(noisy_data.as_array()[:,:,sliceSel],vmin=0, vmax=1) +plt.axis('off') +plt.title('Sagittal View') + + +plt.subplot(2,3,4) +plt.imshow(pdhg.get_output().as_array()[sliceSel,:,:],vmin=0, vmax=1) +plt.axis('off') +plt.subplot(2,3,5) +plt.imshow(pdhg.get_output().as_array()[:,sliceSel,:],vmin=0, vmax=1) +plt.axis('off') +plt.subplot(2,3,6) +plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1) +plt.axis('off') +im = plt.imshow(pdhg.get_output().as_array()[:,:,sliceSel],vmin=0, vmax=1) + + +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/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py new file mode 100644 index 0000000..ccdabb2 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py @@ -0,0 +1,186 @@ +# -*- 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 + +import numpy as np +import numpy +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, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Poisson 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 +data = ImageData(data) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +ag = ig + +# Create noisy data. Apply Poisson noise +n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 2 + +method = '1' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = KullbackLeibler(noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + 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, '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) + + +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 +# +#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( pdhg.get_output().as_array() - u1.value ) +# +# 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(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), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), 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 (PDHG) {} '.format(pdhg.objective[-1][0])) +# +# +# +# +# diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py new file mode 100644 index 0000000..484fe42 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py @@ -0,0 +1,177 @@ +# -*- 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 + +import numpy as np +import numpy +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, L1Norm, \ + MixedL21Norm, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Salt & Pepper 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 +data = ImageData(data) +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 = 's&p', salt_vs_pepper = 0.9, amount=0.2) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 2 + +method = '1' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * MixedL21Norm() + f2 = L1Norm(b = noisy_data) + f = BlockFunction(f1, f2) + + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * MixedL21Norm() + g = L1Norm(b = noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# 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 +pdhg.run(2000) + + +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 + +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) + + diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value ) + + 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(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), 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])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py new file mode 100644 index 0000000..0711e91 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py @@ -0,0 +1,211 @@ +# -*- 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, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +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, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +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) + +# Regularisation Parameter +alpha = 5 + +# 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() + +# 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 = 50 +pdhg.run(2000) + +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) + + 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) + + 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/wip/Demos/PDHG_TV_Tomo2D_time.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.py new file mode 100644 index 0000000..045458a --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D_time.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/wip/Demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py new file mode 100644 index 0000000..3f275e2 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py @@ -0,0 +1,176 @@ +# -*- 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 + +import numpy as np +import numpy +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, L2NormSquared, BlockFunction + +from skimage.util import random_noise + +# Create phantom for TV Salt & Pepper 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 +data = ImageData(data) +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) +noisy_data = ImageData(n1) + +# Regularisation Parameter +alpha = 4 + +method = '1' + +if method == '0': + + # Create operators + op1 = Gradient(ig) + op2 = Identity(ig, ag) + + # Create BlockOperator + operator = BlockOperator(op1, op2, shape=(2,1) ) + + # Create functions + + f1 = alpha * L2NormSquared() + f2 = 0.5 * L2NormSquared(b = noisy_data) + f = BlockFunction(f1, f2) + g = ZeroFunction() + +else: + + # Without the "Block Framework" + operator = Gradient(ig) + f = alpha * L2NormSquared() + g = 0.5 * L2NormSquared(b = noisy_data) + + +# Compute operator Norm +normK = operator.norm() + +# 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 +pdhg.run(2000) + + +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('Tikhonov 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 = 'Tikhonov 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 + 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_squares(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( pdhg.get_output().as_array() - u.value ) + + 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(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), 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])) + + + + + diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py new file mode 100644 index 0000000..5c03362 --- /dev/null +++ b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Tomo2D.py @@ -0,0 +1,108 @@ +# -*- 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, L2NormSquared, BlockFunction +from skimage.util import random_noise +from ccpi.astra.ops import AstraProjectorSimple + +# Create phantom for TV 2D tomography +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, dtype=np.float32) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +# Create noisy data. Apply Gaussian noise + +np.random.seed(10) +noisy_data = sin + AcquisitionData(np.random.normal(0, 3, sin.shape)) + +# Regularisation Parameter +alpha = 500 + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions + +f1 = alpha * L2NormSquared() +f2 = 0.5 * L2NormSquared(b=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 = 5000 +pdhg.update_objective_interval = 50 +pdhg.run(2000) + +#%% +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('Tikhonov 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 = 'Tikhonov reconstruction') +plt.legend() +plt.title('Middle Line Profiles') +plt.show() + + diff --git a/Wrappers/Python/wip/pdhg_TGV_denoising.py b/Wrappers/Python/wip/pdhg_TGV_denoising.py new file mode 100644 index 0000000..1c570cb --- /dev/null +++ b/Wrappers/Python/wip/pdhg_TGV_denoising.py @@ -0,0 +1,230 @@ +#!/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 new file mode 100644 index 0000000..ee3b089 --- /dev/null +++ b/Wrappers/Python/wip/pdhg_TGV_tomography2D.py @@ -0,0 +1,200 @@ +#!/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 index a5cd1bf..b16e8f9 100755 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising.py @@ -6,24 +6,25 @@ Created on Fri Feb 22 14:53:03 2019 @author: evangelos """ -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer +from ccpi.framework import ImageData, ImageGeometry -import numpy as np +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, FunctionOperatorComposition, BlockFunction, ScaledFunction + MixedL21Norm, BlockFunction from skimage.util import random_noise from timeit import default_timer as timer -def dt(steps): - return steps[-1] - steps[-2] +#def dt(steps): +# return steps[-1] - steps[-2] -# Create phantom for TV denoising +# Create phantom for TV Gaussian denoising N = 100 @@ -44,7 +45,7 @@ plt.title('Noisy data') plt.show() # Regularisation Parameter -alpha = 2 +alpha = 0.5 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") method = '1' @@ -72,30 +73,43 @@ else: # No Composite # ########################################################################### operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = 0.5 * L2NormSquared(b = noisy_data) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) -# Compute operator Norm -normK = operator.norm() + +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 -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) + 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':2000} -opt1 = {'niter':2000, 'memopt': True} + +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() -print(" Run memopt") - 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()) @@ -115,10 +129,76 @@ 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_denoising_salt_pepper.py b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py index cec9770..bd330fc 100644 --- a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py @@ -8,18 +8,20 @@ Created on Fri Feb 22 14:53:03 2019 from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer -import numpy as np +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 ZeroFun, L1Norm, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction - +from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ + MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction from skimage.util import random_noise +from timeit import default_timer as timer + # ############################################################################ @@ -54,18 +56,13 @@ if method == '0': 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 = L1Norm(b = noisy_data) f = BlockFunction(f1, f2 ) - g = ZeroFun() + g = ZeroFunction() else: @@ -73,12 +70,12 @@ else: # No Composite # ########################################################################### operator = Gradient(ig) - f = alpha * FunctionOperatorComposition(operator, MixedL21Norm()) + f = alpha * MixedL21Norm() g = L1Norm(b = noisy_data) ########################################################################### #%% -diag_precon = True +diag_precon = False if diag_precon: @@ -94,20 +91,48 @@ if diag_precon: else: # Compute operator Norm normK = operator.norm() - print ("normK", normK) + # Primal & dual stepsizes - sigma = 1/normK - tau = 1/normK -# tau = 1/(sigma*normK**2) + sigma = 1 + tau = 1/(sigma*normK**2) -opt = {'niter':2000} +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) - -plt.figure(figsize=(5,5)) +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 @@ -137,43 +162,106 @@ plt.show() #plt.show() -#%% Compare with cvx +#%% Check with CVX solution -try_cvx = input("Do you want CVX comparison (0/1)") - -if try_cvx=='0': +from ccpi.optimisation.operators import SparseFiniteDiff +try: 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 + cvx_not_installable = True +except ImportError: + cvx_not_installable = False + - u = Variable((N, N)) +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) - regulariser = alpha * TV_cvx(u) - solver = MOSEK + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + obj = Minimize( regulariser + fidelity) - constraints = [] - prob = Problem(obj, constraints) - - # Choose solver (SCS is fast but less accurate than MOSEK) + prob = Problem(obj) 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) + + 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.title('|CVX-PDHG|') 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.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG') plt.legend() - plt.show() + -else: - print('No CVX solution available') + + + 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 index 3fec34e..cd91409 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -8,16 +8,16 @@ Created on Fri Feb 22 14:53:03 2019 @author: evangelos """ -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData +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, Identity, Gradient +from ccpi.optimisation.operators import BlockOperator, Gradient from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ - MixedL21Norm, BlockFunction, ScaledFunction + MixedL21Norm, BlockFunction from ccpi.astra.ops import AstraProjectorSimple from skimage.util import random_noise @@ -43,7 +43,7 @@ from timeit import default_timer as timer #ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) #data = ImageData(phantom_2D, geometry=ig) -N = 150 +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 @@ -52,11 +52,11 @@ data = ImageData(x) ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -detectors = 150 -angles = np.linspace(0,np.pi,100) +detectors = N +angles = np.linspace(0,np.pi,N) ag = AcquisitionGeometry('parallel','2D',angles, detectors) -Aop = AstraProjectorSimple(ig, ag, 'gpu') +Aop = AstraProjectorSimple(ig, ag, 'cpu') sin = Aop.direct(data) plt.imshow(sin.as_array()) @@ -75,11 +75,6 @@ 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 @@ -87,7 +82,7 @@ op2 = Aop # Form Composite Operator operator = BlockOperator(op1, op2, shape=(2,1) ) -alpha = 50 +alpha = 10 f = BlockFunction( alpha * MixedL21Norm(), \ 0.5 * L2NormSquared(b = noisy_data) ) g = ZeroFunction() @@ -109,19 +104,18 @@ if diag_precon: tau, sigma = tau_sigma_precond(operator) -else: +else: + normK = operator.norm() sigma = 1 tau = 1/(sigma*normK**2) # Compute operator Norm -normK = operator.norm() + # Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} +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) @@ -132,25 +126,107 @@ t3 = timer() res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) t4 = timer() # -print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) - - -#%% -#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.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() # -##%% -#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() +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 index dea8e5c..5423b22 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D_time.py @@ -16,7 +16,7 @@ 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, \ +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ MixedL21Norm, BlockFunction, ScaledFunction from ccpi.astra.ops import AstraProjectorSimple, AstraProjectorMC @@ -100,7 +100,7 @@ operator = BlockOperator(op1, op2, shape=(2,1) ) alpha = 50 f = BlockFunction( alpha * MixedL21Norm(), \ 0.5 * L2NormSquared(b = noisy_data) ) -g = ZeroFun() +g = ZeroFunction() # Compute operator Norm normK = operator.norm() diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py index 9fad6f8..cb37598 100644 --- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py +++ b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py @@ -6,26 +6,28 @@ Created on Fri Feb 22 14:53:03 2019 @author: evangelos """ -from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer - -import numpy as np +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 ZeroFun, KullbackLeibler, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction +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 denoising +# Create phantom for TV Poisson denoising -N = 100 +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 @@ -34,20 +36,19 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) ag = ig # Create noisy data. Add Gaussian noise -n1 = random_noise(data, mode = 'poisson') +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 = 10 +alpha = 2 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '1' + +method = '0' if method == '0': # Create operators @@ -57,15 +58,11 @@ if method == '0': # 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 = KullbackLeibler(b = noisy_data) + f2 = KullbackLeibler(noisy_data) f = BlockFunction(f1, f2 ) - g = ZeroFun() + g = ZeroFunction() else: @@ -73,96 +70,118 @@ else: # No Composite # ########################################################################### operator = Gradient(ig) - f = alpha * FunctionOperatorComposition(operator, MixedL21Norm()) + f = alpha * MixedL21Norm() g = KullbackLeibler(noisy_data) ########################################################################### -#%% # Compute operator Norm normK = operator.norm() -print ("normK", normK) -# Primal & dual stepsizes -#sigma = 1 -#tau = 1/(sigma*normK**2) -sigma = 1/normK -tau = 1/normK +# 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) - -plt.figure(figsize=(5,5)) +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() -#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) -#pdhg.max_iteration = 2000 -#pdhg.update_objective_interval = 10 -# -#pdhg.run(2000) +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)) - -#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() - - -#%% Compare with cvx - -#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') +#%% 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])) + + diff --git a/Wrappers/Python/wip/test_pdhg_profile/profile_pdhg_TV_denoising.py b/Wrappers/Python/wip/test_pdhg_profile/profile_pdhg_TV_denoising.py deleted file mode 100644 index e142d94..0000000 --- a/Wrappers/Python/wip/test_pdhg_profile/profile_pdhg_TV_denoising.py +++ /dev/null @@ -1,273 +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, ScaledFunction - -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 = 500 - -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.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) ) - - #### 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() -# -print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) - -# -#%% -#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/test_symGrad_method1.py b/Wrappers/Python/wip/test_symGrad_method1.py new file mode 100644 index 0000000..36adee1 --- /dev/null +++ b/Wrappers/Python/wip/test_symGrad_method1.py @@ -0,0 +1,178 @@ +#!/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, 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 = 3 + +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.2, 1 + +#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") +method = '1' + + +# 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) + +operator1 = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) + + +f1 = alpha * MixedL21Norm() +f2 = beta * MixedL21Norm() +f3 = ZeroFunction() +f_B3 = BlockFunction(f1, f2, f3) +g_B3 = ZeroFunction() + + + +# Create operators +op11 = Gradient(ig) +op12 = Identity(op11.range_geometry()) + +op22 = SymmetrizedGradient(op11.domain_geometry()) + +op21 = ZeroOperator(ig, op22.range_geometry()) + +operator2 = BlockOperator(op11, -1*op12, \ + op21, op22, \ + shape=(2,2) ) + +#f1 = alpha * MixedL21Norm() +#f2 = beta * MixedL21Norm() +f_B2 = BlockFunction(f1, f2) +g_B2 = 0.5 * L2NormSquared(b = noisy_data) + + +#%% + +x_old1 = operator1.domain_geometry().allocate('random_int') +y_old1 = operator1.range_geometry().allocate() + +xbar1 = x_old1.copy() +x_tmp1 = x_old1.copy() +x1 = x_old1.copy() + +y_tmp1 = y_old1.copy() +y1 = y_tmp1.copy() + +x_old2 = x_old1.copy() +y_old2 = operator2.range_geometry().allocate() + +xbar2 = x_old2.copy() +x_tmp2 = x_old2.copy() +x2 = x_old2.copy() + +y_tmp2 = y_old2.copy() +y2 = y_tmp2.copy() + +sigma = 0.4 +tau = 0.4 + +y_tmp1 = y_old1 + sigma * operator1.direct(xbar1) +y_tmp2 = y_old2 + sigma * operator2.direct(xbar2) + +numpy.testing.assert_array_equal(y_tmp1[0][0].as_array(), y_tmp2[0][0].as_array()) +numpy.testing.assert_array_equal(y_tmp1[0][1].as_array(), y_tmp2[0][1].as_array()) +numpy.testing.assert_array_equal(y_tmp1[1][0].as_array(), y_tmp2[1][0].as_array()) +numpy.testing.assert_array_equal(y_tmp1[1][1].as_array(), y_tmp2[1][1].as_array()) + + +y1 = f_B3.proximal_conjugate(y_tmp1, sigma) +y2 = f_B2.proximal_conjugate(y_tmp2, sigma) + +numpy.testing.assert_array_equal(y1[0][0].as_array(), y2[0][0].as_array()) +numpy.testing.assert_array_equal(y1[0][1].as_array(), y2[0][1].as_array()) +numpy.testing.assert_array_equal(y1[1][0].as_array(), y2[1][0].as_array()) +numpy.testing.assert_array_equal(y1[1][1].as_array(), y2[1][1].as_array()) + + +x_tmp1 = x_old1 - tau * operator1.adjoint(y1) +x_tmp2 = x_old2 - tau * operator2.adjoint(y2) + +numpy.testing.assert_array_equal(x_tmp1[0].as_array(), x_tmp2[0].as_array()) + + + + + + + + + + + +############################################################################## +#x_1 = operator1.domain_geometry().allocate('random_int') +# +#x_2 = BlockDataContainer(x_1[0], x_1[1]) +# +#res1 = operator1.direct(x_1) +#res2 = operator2.direct(x_2) +# +#print(res1[0][0].as_array(), res2[0][0].as_array()) +#print(res1[0][1].as_array(), res2[0][1].as_array()) +# +#print(res1[1][0].as_array(), res2[1][0].as_array()) +#print(res1[1][1].as_array(), res2[1][1].as_array()) +# +##res1 = op11.direct(x1[0]) - op12.direct(x1[1]) +##res2 = op21.direct(x1[0]) - op22.direct(x1[1]) |