Move lprec iterative schemes from wrappers.py to the lprec package (#349

)

* Move iterative schemes from wrappers.py to lprec package, add CG method

* Move iterative schemes from wrappers.py to lprec package, add CG method

* Move iterative schemes from wrappers.py to lprec package, add CG method

* Move iterative schemes from wrappers.py to lprec package, add CG method

* Move iterative schemes from wrappers.py to lprec package, add CG method

doc/source/ipynb/lprec.rst

@@ -5,19 +5,13 @@ LPrec
 Here is an example on how to use the log-polar based method
 (https://github.com/math-vrn/lprec) for reconstruction in Tomopy
 
-.. code:: ipython2
+.. code:: ipython3
 
     %pylab inline
 
-
-.. parsed-literal::
-
-    Populating the interactive namespace from numpy and matplotlib
-
-
 Install lprec from github, then
 
-.. code:: ipython2
+.. code:: ipython3
 
     import tomopy
 
@@ -27,27 +21,27 @@ from all major
 `synchrotron <http://dxchange.readthedocs.io/en/latest/source/demo.html>`__
 facilities are supported.
 
-.. code:: ipython2
+.. code:: ipython3
 
     import dxchange
 
 matplotlib provide plotting of the result in this notebook.
 `Paraview <http://www.paraview.org/>`__ or other tools are available for
 more sophisticated 3D rendering.
 
-.. code:: ipython2
+.. code:: ipython3
 
     import matplotlib.pyplot as plt
 
 Set the path to the micro-CT data to reconstruct.
 
-.. code:: ipython2
+.. code:: ipython3
 
     fname = '../../tomopy/data/tooth.h5'
 
 Select the sinogram range to reconstruct.
 
-.. code:: ipython2
+.. code:: ipython3
 
     start = 0
     end = 2
@@ -57,13 +51,13 @@ beamline `2-BM and 32-ID <https://www1.aps.anl.gov/Imaging>`__
 definition. Other file format readers are available at
 `DXchange <http://dxchange.readthedocs.io/en/latest/source/api/dxchange.exchange.html>`__.
 
-.. code:: ipython2
+.. code:: ipython3
 
-    proj, flat, dark, theta = dxchange.read_aps_32id(fname, sino=(start, end))
+    proj, flat, dark = dxchange.read_aps_32id(fname, sino=(start, end))
 
 Plot the sinogram:
 
-.. code:: ipython2
+.. code:: ipython3
 
     plt.imshow(proj[:, 0, :], cmap='Greys_r')
     plt.show()
@@ -77,45 +71,45 @@ If the angular information is not avaialable from the raw data you need
 to set the data collection angles. In this case theta is set as equally
 spaced between 0-180 degrees.
 
-.. code:: ipython2
+.. code:: ipython3
 
     theta = tomopy.angles(proj.shape[0])
 
 Perform the flat-field correction of raw data:
 
 .. math::  \frac{proj - dark} {flat - dark} 
 
-.. code:: ipython2
+.. code:: ipython3
 
     proj = tomopy.normalize(proj, flat, dark)
 
 Select the rotation center manually
 
-.. code:: ipython2
+.. code:: ipython3
 
     rot_center = 296
 
 Calculate
 
 .. math::  -log(proj) 
 
-.. code:: ipython2
+.. code:: ipython3
 
     proj = tomopy.minus_log(proj)
 
 Reconstruction using FBP method with the log-polar coordinates
 
-.. code:: ipython2
+.. code:: ipython3
 
-    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='lpfbp', filter_name='parzen')
+    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='fbp', filter_name='parzen')
 
 Mask each reconstructed slice with a circle.
 
-.. code:: ipython2
+.. code:: ipython3
 
     recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
 
-.. code:: ipython2
+.. code:: ipython3
 
     plt.imshow(recon[0, :,:], cmap='Greys_r')
     plt.show()
@@ -128,9 +122,9 @@ Mask each reconstructed slice with a circle.
 Reconstruction using the gradient descent method with the log-polar
 coordinates
 
-.. code:: ipython2
+.. code:: ipython3
 
-    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='lpgrad', ncore=1, num_iter=64, reg_par=-1)
+    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='grad', ncore=1, num_iter=64, reg_par=-1)
     recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
     plt.imshow(recon[0, :,:], cmap='Greys_r')
     plt.show()
@@ -140,11 +134,12 @@ coordinates
 .. image:: lprec_files/output_30_0.png
 
 
-Reconstruction using the TV method with the log-polar coordinates
+Reconstruction using the conjugate gradient method with the log-polar
+coordinates
 
-.. code:: ipython2
+.. code:: ipython3
 
-    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='lptv', ncore=1, num_iter=256, reg_par=2e-4)
+    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='cg', ncore=1, num_iter=16, reg_par=-1)
     recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
     plt.imshow(recon[0, :,:], cmap='Greys_r')
     plt.show()
@@ -154,11 +149,11 @@ Reconstruction using the TV method with the log-polar coordinates
 .. image:: lprec_files/output_32_0.png
 
 
-Reconstruction using the MLEM method with the log-polar coordinates
+Reconstruction using the TV method with the log-polar coordinates
 
-.. code:: ipython2
+.. code:: ipython3
 
-    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='lpem', ncore=1, num_iter=64, reg_par=0.05)
+    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='tv', ncore=1, num_iter=256, reg_par=1e-3)
     recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
     plt.imshow(recon[0, :,:], cmap='Greys_r')
     plt.show()
@@ -167,3 +162,17 @@ Reconstruction using the MLEM method with the log-polar coordinates
 
 .. image:: lprec_files/output_34_0.png
 
+
+Reconstruction using the MLEM method with the log-polar coordinates
+
+.. code:: ipython3
+
+    recon = tomopy.recon(proj, theta, center=rot_center, algorithm=tomopy.lprec, lpmethod='em', ncore=1, num_iter=64, reg_par=0.05)
+    recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
+    plt.imshow(recon[0, :,:], cmap='Greys_r')
+    plt.show()
+
+
+
+.. image:: lprec_files/output_36_0.png
+

tomopy/recon/wrappers.py

@@ -75,7 +75,7 @@
         'gpu_list': None,
     },
     'lprec': {
-        'lpmethod': 'lpfbp',
+        'lpmethod': 'fbp',
         'interp_type': 'cubic',
         'filter_name': 'None',
         'num_iter': '1',
@@ -366,10 +366,11 @@ def lprec(tomo, center, recon, theta, **kwargs):
     ----------
     lpmethod : str
         LP reconsruction method to use
-            - 'lpfbp'
-            - 'lpgrad'
-            - 'lptv'
-            - 'lpem'
+            - 'fbp'
+            - 'grad'
+            - 'cg'
+            - 'tv'
+            - 'em'
     filter_type:
         Filter for backprojection
             - 'ramp'
@@ -399,11 +400,6 @@ def lprec(tomo, center, recon, theta, **kwargs):
     >>> pylab.imshow(rec[64], cmap='gray')
     >>> pylab.show()
     """
-    lpmethods = {'lpfbp': lpfbp,
-                 'lpgrad': lpgrad,
-                 'lptv': lptv,
-                 'lpem': lpem}
-
     from lprec import lpTransform
 
     # set default options
@@ -425,140 +421,27 @@ def lprec(tomo, center, recon, theta, **kwargs):
     lphandle = lpTransform.lpTransform(
         N, Nproj, Nslices0, filter_name, int(center[0]+0.5), interp_type)
 
-    if(lpmethod == 'lpfbp'):
+    if(lpmethod == 'fbp'):
         # precompute only for the adj transform
         lphandle.precompute(0)
         lphandle.initcmem(0)
         for k in range(0, int(np.ceil(Nslices/float(Nslices0)))):
             ids = range(k*Nslices0, min(Nslices, (k+1)*Nslices0))
-            recon[ids] = lpmethods[lpmethod](lphandle, tomo[ids])
+            recon[ids] = lphandle.adj(tomo[ids])
     else:
+        # iterative schemes
         # precompute for both fwd and adj transforms
         lphandle.precompute(1)
         lphandle.initcmem(1)
+
+        from lprec import iterative
+        lpitermethods = {'grad': iterative.grad,
+                         'cg': iterative.cg,
+                         'tv': iterative.tv,
+                         'em': iterative.em
+                         }
         # run
         for k in range(0, int(np.ceil(Nslices/float(Nslices0)))):
             ids = range(k*Nslices0, min(Nslices, (k+1)*Nslices0))
-            recon[ids] = lpmethods[lpmethod](
+            recon[ids] = lpitermethods[lpmethod](
                 lphandle, recon[ids], tomo[ids], num_iter, reg_par)
-
-
-def lpfbp(lp, tomo):
-    """
-    Reconstruction with Filtered backprojection
-    """
-
-    return lp.adj(tomo)
-
-
-def lpgrad(lp, recon, tomo, num_iter, reg_par):
-    """
-    Reconstruction with the gradient descent method
-    with the regularization parameter reg_par,
-    if reg_par<0, then reg_par is computed on each iteration
-    Solving the problem 1/2||R(recon)-tomo||^2_2 -> min
-    """
-
-    recon0 = recon
-    grad = recon*0
-    grad0 = recon*0
-
-    for i in range(0, num_iter):
-        grad = 2*lp.adj(lp.fwd(recon)-tomo)
-        if(reg_par < 0):
-            if(i == 0):
-                lam = np.float32(1e-3*np.ones(tomo.shape[0]))
-            else:
-                lam = np.sum(np.sum((recon-recon0)*(grad-grad0), 1), 1) / \
-                    np.sum(np.sum((grad-grad0)*(grad-grad0), 1), 1)
-        else:
-            lam = np.float32(reg_par*np.ones(tomo.shape[0]))
-        recon0 = recon
-        grad0 = grad
-        recon = recon - np.reshape(lam, [tomo.shape[0], 1, 1])*grad
-
-    return recon
-
-
-def lptv(lp, recon, tomo, num_iter, reg_par):
-    """
-    Reconstruction with the total variation method
-    with the regularization parameter reg_par.
-    Solving the problem 1/2||R(recon)-tomo||^2_2 + reg_par*TV(recon) -> min
-    """
-
-    lam = reg_par
-    c = 0.35  # 1/power_method(lp,tomo,num_iter)
-
-    recon0 = recon
-    prox0x = recon*0
-    prox0y = recon*0
-    div0 = recon*0
-    prox1 = tomo*0
-
-    for i in range(0, num_iter):
-        # forward step
-        # compute proximal prox0
-        prox0x[:, :, :-1] += c*(recon[:, :, 1:]-recon[:, :, :-1])
-        prox0y[:, :-1, :] += c*(recon[:, 1:, :]-recon[:, :-1, :])
-        nprox = np.maximum(1, np.sqrt(prox0x*prox0x+prox0y*prox0y)/lam)
-        prox0x = prox0x/nprox
-        prox0y = prox0y/nprox
-        # compute proximal prox1
-        prox1 = (prox1+c*lp.fwd(recon)-c*tomo)/(1+c)
-
-        # backward step
-        recon = recon0
-        div0[:, :, 1:] = (prox0x[:, :, 1:]-prox0x[:, :, :-1])
-        div0[:, :, 0] = prox0x[:, :, 0]
-        div0[:, 1:, :] += (prox0y[:, 1:, :]-prox0y[:, :-1, :])
-        div0[:, 0, :] += prox0y[:, 0, :]
-        recon0 = recon0-c*lp.adj(prox1)+c*div0
-
-        # update recon
-        recon = 2*recon0 - recon
-
-    return recon
-
-
-def lpem(lp, recon, tomo, num_iter, reg_par):
-    """
-    Reconstruction with the Expectation Maximization algorithm for denoising
-    with parameter reg_par manually chosen for avoiding division by 0.
-    Maximization of the likelihood function L(tomo,rho) 
-    """
-
-    eps = reg_par
-    # R^*(ones)
-    xi = lp.adj(tomo*0+1)
-    # modification for avoiing division by 0
-    xi = xi+eps*np.max(xi)
-    e = np.max(tomo)*eps
-
-    for i in range(0, num_iter):
-        g = lp.fwd(recon)
-        upd = lp.adj(tomo/(g+e))
-        recon = recon*(upd/xi)
-    return recon
-
-
-def power_method(lp, tomo, num_iter):
-    """
-    Power_method for estimating constant c in the lptv method
-    """
-
-    x = lp.adj(tomo)
-    for k in range(0, num_iter):
-        prox0x = x*0
-        prox0y = x*0
-        div0 = x*0
-        prox0x[:, :, :-1] = x[:, :, 1:]-x[:, :, :-1]
-        prox0y[:, :-1, :] = x[:, 1:, :]-x[:, :-1, :]
-        div0[:, :, 1:] = (prox0x[:, :, 1:]-prox0x[:, :, :-1])
-        div0[:, :, 0] = prox0x[:, :, 0]
-        div0[:, 1:, :] += (prox0y[:, 1:, :]-prox0y[:, :-1, :])
-        div0[:, 0, :] += prox0y[:, 0, :]
-        x = lp.adj(lp.fwd(x)) - div0
-        s = np.linalg.norm(x)
-        x = x/s
-    return np.sqrt(s)