diff options
| -rw-r--r-- | Wrappers/Python/wip/demo_sophiabeads.py | 65 | ||||
| -rw-r--r-- | Wrappers/Python/wip/simple_demo.py | 189 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/simple_demo_ccpi.py | 270 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/simple_mc_demo.py | 142 |
4 files changed, 0 insertions, 666 deletions
diff --git a/Wrappers/Python/wip/demo_sophiabeads.py b/Wrappers/Python/wip/demo_sophiabeads.py deleted file mode 100644 index e3c7f3a..0000000 --- a/Wrappers/Python/wip/demo_sophiabeads.py +++ /dev/null @@ -1,65 +0,0 @@ - -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/simple_demo.py b/Wrappers/Python/wip/simple_demo.py deleted file mode 100644 index d5d6826..0000000 --- a/Wrappers/Python/wip/simple_demo.py +++ /dev/null @@ -1,189 +0,0 @@ -#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_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py deleted file mode 100755 index d2cd6f5..0000000 --- a/Wrappers/Python/wip/simple_demo_ccpi.py +++ /dev/null @@ -1,270 +0,0 @@ -#import sys -#sys.path.append("..") - -from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry -from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D -#from ccpi.astra.astra_ops import AstraProjectorSimple -from ccpi.reconstruction.ops import CCPiProjectorSimple -from ccpi.reconstruction.parallelbeam import alg as pbalg -from ccpi.processors import CCPiForwardProjector, CCPiBackwardProjector - -import numpy as np -import matplotlib.pyplot as plt - -test_case = 1 # 1=parallel2D, 2=cone2D, 3=parallel3D - -# Set up phantom -N = 128 -vert = 1 -# Set up measurement geometry -angles_num = 20; # angles number -det_w = 1.0 -det_num = N -SourceOrig = 200 -OrigDetec = 0 - -if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi - #nangles = angles_num - #angles = np.linspace(0,360, nangles, dtype=np.float32) - -elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) -elif test_case == 3: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) -else: - NotImplemented - - -def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter): - vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) - Phantom_ccpi = ImageData(geometry=vg, - dimension_labels=['horizontal_x','horizontal_y','vertical']) - #.subset(['horizontal_x','horizontal_y','vertical']) - # ask the ccpi code what dimensions it would like - - voxel_per_pixel = 1 - geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), - angles, - voxel_per_pixel ) - - pg = AcquisitionGeometry('parallel', - '3D', - angles, - geoms['n_h'],det_w, - geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick - ) - - center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 - ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) - geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), - angles, - center_of_rotation, - voxel_per_pixel ) - - #print (counter) - counter+=1 - #print (geoms , geoms_i) - if counter < 4: - if (not ( geoms_i == geoms )): - print ("not equal and {0}".format(counter)) - X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) - Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) - Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) - return setupCCPiGeometries(X,Y,Z,angles, counter) - else: - print ("return geoms {0}".format(geoms)) - return geoms - else: - print ("return geoms_i {0}".format(geoms_i)) - return geoms_i - -geoms = setupCCPiGeometries(N,N,vert,angles,0) -#%% -#geoms = {'n_v': 12, 'output_volume_y': 128, 'n_angles': 20, -# 'output_volume_x': 128, 'output_volume_z': 12, 'n_h': 128} -vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], - voxel_num_y=geoms['output_volume_y'], - voxel_num_z=geoms['output_volume_z']) -Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1 - - -#x = Phantom.as_array() -i = 0 -while i < geoms['n_v']: - #x = Phantom.subset(vertical=i, dimensions=['horizontal_x','horizontal_y']).array - x = Phantom.subset(vertical=i).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 - Phantom.fill(x, vertical=i) - i += 1 - -plt.imshow(Phantom.subset(vertical=0).as_array()) -plt.show() - - - -# Parallelbeam geometry test -if test_case==1: - #Phantom_ccpi = Phantom.subset(['horizontal_x','horizontal_y','vertical']) - #Phantom_ccpi.geometry = vg.clone() - center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2 - - pg = AcquisitionGeometry('parallel', - '3D', - angles, - geoms['n_h'],det_w, - geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick - ) -elif test_case==2: - raise NotImplemented('cone beam projector not yet available') - pg = AcquisitionGeometry('cone', - '2D', - angles, - det_num, - det_w, - vert, det_w, #2D in 3D is a slice 1 pixel thick - dist_source_center=SourceOrig, - dist_center_detector=OrigDetec) - -# ASTRA operator using volume and sinogram geometries -#Aop = AstraProjectorSimple(vg, pg, 'cpu') -Cop = CCPiProjectorSimple(vg, pg) - -# Try forward and backprojection -b = Cop.direct(Phantom) -out2 = Cop.adjoint(b) - -#%% -for i in range(b.get_dimension_size('vertical')): - plt.imshow(b.subset(vertical=i).array) - #plt.imshow(Phantom.subset( vertical=i).array) - #plt.imshow(b.array[:,i,:]) - plt.show() -#%% - -plt.imshow(out2.subset( vertical=0).array) -plt.show() - -# Create least squares object instance with projector and data. -f = Norm2sq(Cop,b,c=0.5) - -# Initial guess -x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical']) -#invL = 0.5 -#g = f.grad(x_init) -#print (g) -#u = x_init - invL*f.grad(x_init) - -#%% -# Run FISTA for least squares without regularization -opt = {'tol': 1e-4, 'iter': 100} -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) - -plt.imshow(x_fista0.subset(vertical=0).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,opt=opt) - -plt.imshow(x_fista0.subset(vertical=0).array) -plt.title('FISTA1') -plt.show() - -plt.semilogy(criter1) -plt.show() - -# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm -x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) - -plt.imshow(x_fbpd1.subset(vertical=0).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.subset(vertical=0).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, Cop, b, opt=opt) - -plt.imshow(x_CGLS.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).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.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) - -plt.show() -#%% -#current = current + 1 -#a=fig.add_subplot(rows,cols,current) -#a.set_title('FBPD TV') -#imgplot = plt.imshow(x_fbpdtv.subset(vertical=0).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 deleted file mode 100755 index f77a678..0000000 --- a/Wrappers/Python/wip/simple_mc_demo.py +++ /dev/null @@ -1,142 +0,0 @@ -#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 |