项目作者: lucasplagwitz

项目描述 :
Primal-Dual Solver for Inverse Problems
高级语言: Python
项目地址: git://github.com/lucasplagwitz/recon.git
创建时间: 2020-02-15T09:15:01Z
项目社区:https://github.com/lucasplagwitz/recon

开源协议:GNU General Public License v3.0

下载


Recon

A python-based toolbox for solving regularized Inverse Problems using Primal-Dual algorithms.
The project provides an overview of solving regularization problems and is the result of a master’s thesis.
Built as proof of concept.

Overview

  • Reconstruction, Smoothing
  • Class-Based Segmentation
  • Spatially-Adapted Regularization
  • Bregman Iteration
  • Denoising, Deconvolution, Computerized Tomography

Reconstruction

In terms of Inverse Problems one is interested in the reason

of measurment data

with regard to a forward map
.
Due to the fact of measurement inaccuracies, regularization terms

are added and the optimization problem is maintained as



  1. import numpy as np
  2. from recon.operator.ct_radon import CtRt
  3. from recon.interfaces import Recon
  4. from matplotlib.image import imread
  5. gt = imread("../data/phantom.png")
  6. gt = gt/np.max(gt)
  7. theta = np.linspace(0, 180, 180, endpoint=False)
  8. sigma = 0.01
  9. R = CtRt(gt.shape, center=[gt.shape[0]//2, gt.shape[1]//2], theta=theta)
  10. y = R*gt.ravel()
  11. rec = Recon(operator=R, domain_shape=gt.shape, reg_mode='tv', alpha=1, lam=15, extend_pdhgm=True)
  12. x_tv = rec.solve(data=y.ravel(), max_iter=1000, tol=1e-4)

Imaging result for another inverse problem where

is a convolution operator:



Denoising

Image denoising is a special case of regularized reconstruction.



  1. import numpy as np
  2. from scipy import misc
  3. from recon.interfaces import Smoothing
  4. img = misc.ascent()
  5. gt = img/np.max(img)
  6. sigma = 0.2
  7. n = sigma*np.max(gt.ravel()*np.random.uniform(-1, 1, gt.shape))
  8. noise_img = gt + n
  9. tv_smoothing = Smoothing(domain_shape=gt.shape, reg_mode='tv', lam=10, tau='calc')
  10. u0 = tv_smoothing.solve(data=noise_img, maxiter=1500, tol=10**(-4))










Segmentation

Some segmentation methods are implemented as part of regularization approaches and performance measurements.
Through a piecewise constant TV-solution, one quickly obtains a suitable segmentation.

  1. import skimage.data as skd
  2. import numpy as np
  3. from recon.interfaces import Segmentation
  4. gt = rgb2gray(skd.coffee())[:,80:481]
  5. gt = gt/np.max(gt)
  6. gt = gt/np.max(gt)
  7. classes = [0, 50/255, 120/255, 190/255, 220/255]
  8. segmentation = Segmentation(gt.shape, classes=classes, lam=5, tau='calc')
  9. result, _ = segmentation.solve(gt, max_iter=4000)



References

  1. The Repo based on Enhancing joint reconstruction and segmentation with non-convex Bregman iteration - Veronica Corona et al 2019 Inverse Problems 35, and their code on GitHub.
  2. To outsource operator handling PyLops package is used.
  3. Efficient Radon operator Astra Toolbox.