6 files changed, 597 insertions, 0 deletions
diff --git a/Wrappers/Python/wip/DemoRecIP.py b/Wrappers/Python/wip/DemoRecIP.py
new file mode 100755
index 0000000..442e40e
--- /dev/null
+++ b/Wrappers/Python/wip/DemoRecIP.py
@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Reading multi-channel data and reconstruction using FISTA modular
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+#import sys
+#sys.path.append('../../../data/')
+from read_IPdata import read_IPdata
+
+from ccpi.astra.astra_ops import AstraProjectorSimple, AstraProjectorMC
+from ccpi.reconstruction.funcs import Norm2sq, Norm1, BaseFunction
+from ccpi.reconstruction.algs import FISTA
+#from ccpi.reconstruction.funcs import BaseFunction
+
+from ccpi.framework import ImageData, AcquisitionData, AcquisitionGeometry, ImageGeometry
+
+# read IP paper data into a dictionary
+dataDICT = read_IPdata('..\..\..\data\IP_data70channels.mat')
+
+# Set ASTRA Projection-backprojection class (fan-beam geometry)
+DetWidth = dataDICT.get('im_size')[0] * dataDICT.get('det_width')[0] / \
+            dataDICT.get('detectors_numb')[0]
+SourceOrig = dataDICT.get('im_size')[0] * dataDICT.get('src_to_rotc')[0] / \
+            dataDICT.get('dom_width')[0]
+OrigDetec = dataDICT.get('im_size')[0] * \
+            (dataDICT.get('src_to_det')[0] - dataDICT.get('src_to_rotc')[0]) /\
+            dataDICT.get('dom_width')[0]
+
+N = dataDICT.get('im_size')[0]
+
+vg = ImageGeometry(voxel_num_x=dataDICT.get('im_size')[0],
+                   voxel_num_y=dataDICT.get('im_size')[0],
+                   channels=1)
+
+pg = AcquisitionGeometry('cone',
+                         '2D',
+                         angles=(np.pi/180)*dataDICT.get('theta')[0],
+                         pixel_num_h=dataDICT.get('detectors_numb')[0],
+                         pixel_size_h=DetWidth,
+                         dist_source_center=SourceOrig, 
+                         dist_center_detector=OrigDetec,
+                         channels=1)
+
+
+sino = dataDICT.get('data_norm')[0][:,:,34] # select mid-channel 
+b = AcquisitionData(sino,geometry=pg)
+
+# Initial guess
+x_init = ImageData(np.zeros((N, N)),geometry=vg)
+
+
+
+
+
+Aop = AstraProjectorSimple(vg,pg,'gpu')
+f = Norm2sq(Aop,b,c=0.5)
+
+# Run FISTA for least squares without regularization
+opt = {'tol': 1e-4, 'iter': 10}
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+
+plt.imshow(x_fista0.array)
+plt.show()
+
+# Now least squares plus 1-norm regularization
+g1 = Norm1(10)
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt)
+
+plt.imshow(x_fista1.array)
+plt.show()
+
+# Multiple channels
+sino_mc = dataDICT.get('data_norm')[0][:,:,32:37] # select mid-channel 
+
+vg_mc = ImageGeometry(voxel_num_x=dataDICT.get('im_size')[0],
+                   voxel_num_y=dataDICT.get('im_size')[0],
+                   channels=5)
+
+pg_mc = AcquisitionGeometry('cone',
+                         '2D',
+                         angles=(np.pi/180)*dataDICT.get('theta')[0],
+                         pixel_num_h=dataDICT.get('detectors_numb')[0],
+                         pixel_size_h=DetWidth,
+                         dist_source_center=SourceOrig, 
+                         dist_center_detector=OrigDetec,
+                         channels=5)
+
+b_mc = AcquisitionData(np.transpose(sino_mc,(2,0,1)),
+                       geometry=pg_mc, 
+                       dimension_labels=("channel","angle","horizontal"))
+
+# ASTRA operator using volume and sinogram geometries
+Aop_mc = AstraProjectorMC(vg_mc, pg_mc, 'gpu')
+
+f_mc = Norm2sq(Aop_mc,b_mc,c=0.5)
+
+# Initial guess
+x_init_mc = ImageData(np.zeros((5, N, N)),geometry=vg_mc)
+
+
+x_fista0_mc, it0_mc, timing0_mc, criter0_mc = FISTA(x_init_mc, f_mc, None, opt)
+
+plt.imshow(x_fista0_mc.as_array()[4,:,:])
+plt.show()
+\ No newline at end of file
diff --git a/Wrappers/Python/wip/astra_test.py b/Wrappers/Python/wip/astra_test.py
new file mode 100755
index 0000000..c0a7359
--- /dev/null
+++ b/Wrappers/Python/wip/astra_test.py
@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar  7 15:07:16 2018
+
+@author: ofn77899
+"""
+
+# -----------------------------------------------------------------------
+# Copyright: 2010-2018, imec Vision Lab, University of Antwerp
+#            2013-2018, CWI, Amsterdam
+#
+# Contact: astra@astra-toolbox.com
+# Website: http://www.astra-toolbox.com/
+#
+# This file is part of the ASTRA Toolbox.
+#
+#
+# The ASTRA Toolbox is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The ASTRA Toolbox is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+#
+# -----------------------------------------------------------------------
+
+import astra
+import numpy as np
+
+vol_geom = astra.create_vol_geom(256, 256)
+proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False))
+
+# For CPU-based algorithms, a "projector" object specifies the projection
+# model used. In this case, we use the "strip" model.
+proj_id = astra.create_projector('strip', proj_geom, vol_geom)
+
+# Create a sinogram from a phantom
+import scipy.io
+P = scipy.io.loadmat('phantom.mat')['phantom256']
+sinogram_id, sinogram = astra.create_sino(P, proj_id)
+
+import pylab
+pylab.gray()
+pylab.figure(1)
+pylab.imshow(P)
+pylab.figure(2)
+pylab.imshow(sinogram)
+
+# Create a data object for the reconstruction
+rec_id = astra.data2d.create('-vol', vol_geom)
+
+# Set up the parameters for a reconstruction algorithm using the CPU
+# The main difference with the configuration of a GPU algorithm is the
+# extra ProjectorId setting.
+cfg = astra.astra_dict('SIRT')
+cfg['ReconstructionDataId'] = rec_id
+cfg['ProjectionDataId'] = sinogram_id
+cfg['ProjectorId'] = proj_id
+
+# Available algorithms:
+# ART, SART, SIRT, CGLS, FBP
+
+
+# Create the algorithm object from the configuration structure
+alg_id = astra.algorithm.create(cfg)
+
+# Run 20 iterations of the algorithm
+# This will have a runtime in the order of 10 seconds.
+astra.algorithm.run(alg_id, 20)
+
+# Get the result
+rec = astra.data2d.get(rec_id)
+pylab.figure(3)
+pylab.imshow(rec)
+pylab.show()
+
+# Clean up.
+astra.algorithm.delete(alg_id)
+astra.data2d.delete(rec_id)
+astra.data2d.delete(sinogram_id)
+astra.projector.delete(proj_id)
diff --git a/Wrappers/Python/wip/demo_sophiabeads.py b/Wrappers/Python/wip/demo_sophiabeads.py
new file mode 100755
index 0000000..e3c7f3a
--- /dev/null
+++ b/Wrappers/Python/wip/demo_sophiabeads.py
@@ -0,0 +1,65 @@
+
+from ccpi.io.reader import XTEKReader
+import numpy as np
+import matplotlib.pyplot as plt
+from ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData, ImageData
+from ccpi.astra.astra_ops import AstraProjectorSimple
+from ccpi.reconstruction.algs import CGLS
+
+# Set up reader object and read the data
+datareader = XTEKReader("C:/Users/mbbssjj2/Documents/SophiaBeads_256_averaged/SophiaBeads_256_averaged.xtekct")
+data = datareader.getAcquisitionData()
+
+# Extract central slice, scale and negative-log transform
+sino = -np.log(data.as_array()[:,:,1000]/60000.0)
+
+# Apply centering correction by zero padding, amount found manually
+cor_pad = 30
+sino_pad = np.zeros((sino.shape[0],sino.shape[1]+cor_pad))
+sino_pad[:,cor_pad:] = sino
+
+# Simple beam hardening correction as done in SophiaBeads coda
+sino_pad = sino_pad**2
+
+# Extract AcquisitionGeometry for central slice for 2D fanbeam reconstruction
+ag2d = AcquisitionGeometry('cone',
+                          '2D',
+                          angles=-np.pi/180*data.geometry.angles,
+                          pixel_num_h=data.geometry.pixel_num_h + cor_pad,
+                          pixel_size_h=data.geometry.pixel_size_h,
+                          dist_source_center=data.geometry.dist_source_center, 
+                          dist_center_detector=data.geometry.dist_center_detector)
+
+# Set up AcquisitionData object for central slice 2D fanbeam
+data2d = AcquisitionData(sino_pad,geometry=ag2d)
+
+# Choose the number of voxels to reconstruct onto as number of detector pixels
+N = data.geometry.pixel_num_h
+
+# Geometric magnification
+mag = (np.abs(data.geometry.dist_center_detector) + \
+      np.abs(data.geometry.dist_source_center)) / \
+      np.abs(data.geometry.dist_source_center)
+
+# Voxel size is detector pixel size divided by mag
+voxel_size_h = data.geometry.pixel_size_h / mag
+
+# Construct the appropriate ImageGeometry
+ig2d = ImageGeometry(voxel_num_x=N,
+                   voxel_num_y=N,
+                   voxel_size_x=voxel_size_h, 
+                   voxel_size_y=voxel_size_h)
+
+# Set up the Projector (AcquisitionModel) using ASTRA on GPU
+Aop = AstraProjectorSimple(ig2d, ag2d,"gpu")
+
+# Set initial guess for CGLS reconstruction
+x_init = ImageData(np.zeros((N,N)),geometry=ig2d)
+
+# Run CGLS reconstruction
+num_iter = 50
+x, it, timing, criter = CGLS(Aop,data2d,num_iter,x_init)
+
+# Display reconstruction
+plt.figure()
+plt.imshow(x.as_array())
+\ No newline at end of file
diff --git a/Wrappers/Python/wip/desktop.ini b/Wrappers/Python/wip/desktop.ini
new file mode 100755
index 0000000..bb9f3d6
--- /dev/null
+++ b/Wrappers/Python/wip/desktop.ini
@@ -0,0 +1,4 @@
+[ViewState]
+Mode=
+Vid=
+FolderType=Generic
diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py
new file mode 100755
index 0000000..3365504
--- /dev/null
+++ b/Wrappers/Python/wip/simple_demo_astra.py
@@ -0,0 +1,189 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry
+from ccpi.reconstruction.algs import FISTA, FBPD, CGLS
+from ccpi.reconstruction.funcs import Norm2sq, Norm1 , TV2D
+from ccpi.astra.astra_ops import AstraProjectorSimple
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+
+
+vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = ImageData(geometry=vg)
+
+x = Phantom.as_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)] = 1
+
+plt.imshow(x)
+plt.show()
+
+# Set up measurement geometry
+angles_num = 20; # angles number
+
+if test_case==1:
+    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
+elif test_case==2:
+    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+else:
+    NotImplemented
+
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+# Parallelbeam geometry test
+if test_case==1:
+    pg = AcquisitionGeometry('parallel',
+                          '2D',
+                          angles,
+                          det_num,det_w)
+elif test_case==2:
+    pg = AcquisitionGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec)
+
+# ASTRA operator using volume and sinogram geometries
+Aop = AstraProjectorSimple(vg, pg, 'cpu')
+
+# Unused old astra projector without geometry
+# Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
+#                      OrigDetec, angles, 
+#                      N,'fanbeam','gpu') 
+
+# Try forward and backprojection
+b = Aop.direct(Phantom)
+out2 = Aop.adjoint(b)
+
+plt.imshow(b.array)
+plt.show()
+
+plt.imshow(out2.array)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(np.zeros(x.shape),geometry=vg)
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None)
+
+plt.imshow(x_fista0.array)
+plt.title('FISTA0')
+plt.show()
+
+# Now least squares plus 1-norm regularization
+lam = 0.1
+g0 = Norm1(lam)
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0)
+
+plt.imshow(x_fista1.array)
+plt.title('FISTA1')
+plt.show()
+
+plt.semilogy(criter1)
+plt.show()
+
+# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
+opt = {'tol': 1e-4, 'iter': 100}
+x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
+
+plt.imshow(x_fbpd1.array)
+plt.title('FBPD1')
+plt.show()
+
+plt.semilogy(criter_fbpd1)
+plt.show()
+
+# Now FBPD for least squares plus TV
+#lamtv = 1
+#gtv = TV2D(lamtv)
+
+#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
+
+#plt.imshow(x_fbpdtv.array)
+#plt.show()
+
+#plt.semilogy(criter_fbpdtv)
+#plt.show()
+
+
+# Run CGLS, which should agree with the FISTA0
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt )
+
+plt.imshow(x_CGLS.array)
+plt.title('CGLS')
+#plt.title('CGLS recon, compare FISTA0')
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+plt.show()
+
+
+#%%
+
+clims = (0,1)
+cols = 3
+rows = 2
+current = 1
+fig = plt.figure()
+# projections row
+a=fig.add_subplot(rows,cols,current)
+a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
+
+imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA0')
+imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA1')
+imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FBPD1')
+imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('CGLS')
+imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
+
+#current = current + 1
+#a=fig.add_subplot(rows,cols,current)
+#a.set_title('FBPD TV')
+#imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1])
+
+fig = plt.figure()
+# projections row
+b=fig.add_subplot(1,1,1)
+b.set_title('criteria')
+imgplot = plt.loglog(criter0 , label='FISTA0')
+imgplot = plt.loglog(criter1 , label='FISTA1')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
+imgplot = plt.loglog(criter_CGLS, label='CGLS')
+#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
+b.legend(loc='right')
+plt.show()
+#%%
+\ No newline at end of file
diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py
new file mode 100755
index 0000000..f77a678
--- /dev/null
+++ b/Wrappers/Python/wip/simple_mc_demo.py
@@ -0,0 +1,142 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry
+from ccpi.reconstruction.algs import FISTA
+from ccpi.reconstruction.funcs import Norm2sq, Norm1
+from ccpi.astra.astra_ops import AstraProjectorMC
+
+import numpy
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+M = 100
+numchannels = 3
+
+vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels)
+Phantom = ImageData(geometry=vg)
+
+x = Phantom.as_array()
+x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.0
+x[0 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.0
+
+x[1 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 0.7
+x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2
+
+x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.5
+x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2
+
+#x = numpy.zeros((N,N,1,numchannels))
+
+#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0
+#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0
+
+#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7
+#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2
+
+#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5
+#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2
+
+f, axarr = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarr[k].imshow(x[k],vmin=0,vmax=2.5)
+plt.show()
+
+#vg = ImageGeometry(N,N,None, 1,1,None,channels=numchannels)
+
+
+#Phantom = ImageData(x,geometry=vg)
+
+# Set up measurement geometry
+angles_num = 20; # angles number
+
+if test_case==1:
+    angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False)
+elif test_case==2:
+    angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False)
+else:
+    NotImplemented
+
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+# Parallelbeam geometry test
+if test_case==1:
+    pg = AcquisitionGeometry('parallel',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          channels=numchannels)
+elif test_case==2:
+    pg = AcquisitionGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec,
+                          channels=numchannels)
+
+# ASTRA operator using volume and sinogram geometries
+Aop = AstraProjectorMC(vg, pg, 'cpu')
+
+
+# Try forward and backprojection
+b = Aop.direct(Phantom)
+
+fb, axarrb = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrb[k].imshow(b.as_array()[k],vmin=0,vmax=250)
+plt.show()
+
+out2 = Aop.adjoint(b)
+
+fo, axarro = plt.subplots(1,numchannels)
+for k in range(3):
+    axarro[k].imshow(out2.as_array()[k],vmin=0,vmax=3500)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(numpy.zeros(x.shape),geometry=vg)
+
+# FISTA options
+opt = {'tol': 1e-4, 'iter': 200}
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+
+
+ff0, axarrf0 = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5)
+plt.show()
+
+plt.semilogy(criter0)
+plt.title('Criterion vs iterations, least squares')
+plt.show()
+
+# Now least squares plus 1-norm regularization
+lam = 0.1
+g0 = Norm1(lam)
+
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
+
+ff1, axarrf1 = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5)
+plt.show()
+
+plt.semilogy(criter1)
+plt.title('Criterion vs iterations, least squares plus 1-norm regu')
+plt.show()
+\ No newline at end of file