From d6f9bb8e7d9eae126c90d2cf74110ac53c859d37 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 23 Mar 2018 15:14:34 +0000 Subject: initial revision --- Wrappers/Python/wip/DemoRecIP.py | 110 ++++++++++++++++++ Wrappers/Python/wip/astra_test.py | 87 ++++++++++++++ Wrappers/Python/wip/demo_sophiabeads.py | 65 +++++++++++ Wrappers/Python/wip/desktop.ini | 4 + Wrappers/Python/wip/simple_demo_astra.py | 189 +++++++++++++++++++++++++++++++ Wrappers/Python/wip/simple_mc_demo.py | 142 +++++++++++++++++++++++ 6 files changed, 597 insertions(+) create mode 100755 Wrappers/Python/wip/DemoRecIP.py create mode 100755 Wrappers/Python/wip/astra_test.py create mode 100755 Wrappers/Python/wip/demo_sophiabeads.py create mode 100755 Wrappers/Python/wip/desktop.ini create mode 100755 Wrappers/Python/wip/simple_demo_astra.py create mode 100755 Wrappers/Python/wip/simple_mc_demo.py (limited to 'Wrappers/Python/wip') 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 . +# +# ----------------------------------------------------------------------- + +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 -- cgit v1.2.3