Commit 85de1a5c authored by Nick Sauerwein's avatar Nick Sauerwein
Browse files

change tomography to tomopy

parent 115e259c
......@@ -9,9 +9,10 @@ import matplotlib.pyplot as plt
from scipy import misc
from scipy import optimize
from scipy import ndimage
import os
import scipy as sp
import math
def sampleavg(path, Basler, g, N=500.):
#return avg. Image of N Samples, g is an int indicating the summation:
......@@ -206,6 +207,10 @@ def phaseunwrap(Z1, Z2):
def Ga(x):
return np.exp(-1./x**2*(6.*(x[1]-x[0]))**2)
def Ga_gpu(x):
import arrayfire as af
xx = af.algorithm.max(x[1] - x[0])
return af.exp(-1./x**2*(6.*xx)**2)
def ft(D):
......@@ -213,11 +218,26 @@ def ft(D):
kx = np.r_[-lx/2.:lx/2.]/lx
F0=np.fft.fft(D,axis=0)
F0=np.fft.fftshift(F0, axes=0)
F0[int(lx/2):,:]=0.
F0[lx//2:,:]=0.
F0=np.transpose(np.transpose(F0)*Ga(kx))
IF0=np.fft.ifft(F0,axis=0)
return IF0
def ft_gpu(D):
import arrayfire as af
lx = D.dims()[0]
ly = D.dims()[1]
kx = np.r_[-lx/2.:lx/2.]/lx
kxx = af.np_to_af_array(kx)
F0 = af.signal.fft(D, dim0 = 0)
F0 = af.data.shift(F0, lx / 2)
F0[lx//2:,:] = 0.
Ga_arr = af.data.tile(af.array.transpose(Ga_gpu(kxx)), ly)
F0=af.array.transpose(af.array.transpose(F0)*Ga_arr)
IF0=af.signal.ifft(F0)
return IF0
def unwrap(D0,D1,fringe_dir='ver'):
if fringe_dir=="hor":
F0=ft(D0)
......@@ -229,6 +249,19 @@ def unwrap(D0,D1,fringe_dir='ver'):
if fringe_dir=="hor": return np.angle(F1/F0)
if fringe_dir=="ver": return np.transpose(np.angle(F1/F0))
def unwrap_gpu(D0, D1, fringe_dir='ver'):
import arrayfire as af
if fringe_dir=="hor":
F0=ft_gpu(D0)
F1=ft_gpu(D1)
if fringe_dir=="ver":
F0=ft_gpu(af.transpose(D0))
F1=ft_gpu(af.transpose(D1))
if fringe_dir=="hor": return af.arith.arg(F1/F0)
if fringe_dir=="ver": return af.transpose(af.arith.arg(F1/F0))
#if fringe_dir=="ver": return F1/F0
def correct_phase_jumps(P):
Pc=np.copy(P)
lx=len(P[:,0])
......@@ -241,37 +274,87 @@ def correct_phase_jumps(P):
g=np.round(diff/np.pi)
Pc[j:,i+1]+=g*np.pi
return Pc
return Pc
def phase(A,As,shift=217):
def phase(A,As,shift=217, y_calib = 200):
lx=len(A[:,0])
ly=len(A[0,:])
Ph=np.copy(A)
Phs=np.copy(As)
for i in range(lx):
Ph[i,:]-=np.average(Ph[i,0:50])
Phs[i,:]-=np.average(Ph[i,0:50])
Ph[i,:]-=np.average(Ph[i,:y_calib])
for j in range(shift,ly):
Phs[i,j]=Phs[i,j]+Phs[i,j-shift]
Ph[i,j]=Ph[i,j]+Phs[i,j-shift]
return Ph
Ph[i,j]=Ph[i,j]+Ph[i,j-shift]
return Ph
def MLEM(X,theta, meas):
from scipy import ndimage
#input: X: current guess of distribution, theta: projection angle, meas: corresponding measurement data
#ouptu: new guess for X
#rotate
Xp=ndimage.rotate(X, -theta, reshape=False)
if (theta != 0):
Xp=ndimage.rotate(X, -theta, reshape=False)
else:
Xp = X
Lxt=len(X[0,:])
X1=np.zeros((Lxt,Lxt))
ysim=np.sum(Xp,axis=0)
ZZ=(len(ysim)-len(meas))/2
meas=np.append(ysim[:ZZ],meas)
ysim=np.sum(Xp,axis=0) # forward projected data from current guess
ZZ=int((len(ysim)-len(meas))/2) # field of view correction
meas=np.append(ysim[:ZZ],meas)# append simulation where out of view
meas=np.append(meas,ysim[-ZZ:])
T=np.abs(meas/ysim)
T=np.abs(meas/ysim) #R, Ratio, correction factor
Lxt4 = int(Lxt/4)
X1[:,Lxt4:-Lxt4]=Xp[:,Lxt4:-Lxt4]*T[Lxt4:-Lxt4]# new guess, eq. 4.14
if (theta != 0):
return ndimage.rotate(X1,theta, reshape=False)
else:
return X1
def MLEM_gpu(X,theta, meas):
from arrayfire.image import rotate as afrotate
from arrayfire import data as afdata
from arrayfire import array as afarray
from arrayfire import arith as afarith
from arrayfire import library as aflib
from arrayfire import algorithm as afalg
#rotate X by theta
if (theta != 0):
thetar = theta * math.pi / 180
Xp = afrotate(X, -thetar, is_crop=True, method=aflib.INTERP.BICUBIC_SPLINE)
else:
Xp = X
#get X lenght and create a temp array wiht zeros
Lxt=X.dims()[1]
X1 = afdata.constant(0.0, Lxt, Lxt, 1, 1, aflib.Dtype.f32)
ysim = afalg.sum(Xp, dim=0) #forward projected data from current guess
ZZ=(ysim.dims()[1] - meas.dims()[0]) / 2
ysim = afarray.transpose(ysim)
meas = afarray.transpose(meas)
meas = afdata.flat(meas)
X1[:,Lxt/4:-Lxt/4]=Xp[:,Lxt/4:-Lxt/4]*T[Lxt/4:-Lxt/4]
return ndimage.rotate(X1,theta, reshape=False)
meas = afdata.join(0, ysim[:ZZ, :], meas)
meas = afdata.join(0, meas, ysim[-ZZ:, :])
T = afarith.abs(meas / ysim)
T = afarray.transpose(T)
T = afdata.tile(T, Lxt)
Lxt4 = int(Lxt/4)
X1[:,Lxt4:-Lxt4]=Xp[:,Lxt4:-Lxt4]*T[:, Lxt4:-Lxt4]# new guess, eq. 4.14
if (theta != 0):
return afrotate(X1, thetar, is_crop=True, method=aflib.INTERP.BICUBIC_SPLINE)
else:
return X1
def iabel(dfdx, x):
nx = len(x)
......@@ -334,4 +417,4 @@ def gauss_fit(x, y):
mean = n/2
sigma = n/4
popt,pcov = sp.optimize.curve_fit(gauss,x,y,p0=[1,mean,sigma])
return popt, pcov
\ No newline at end of file
return popt, pcov
import tomopy
import numpy as np
import matplotlib.pyplot as plt
import time
def calculate_phase(images_angles, y1, y2, Lx,center = 1110, y_calib=200, slices=None, plot_images = False, plot_phases = False):
from routines import phase, unwrap, correct_phase_jumps
if plot_images:
plt.figure()
plt.imshow(images_angles[-1][0])
plt.show()
if slices==None:
slices=range(len(images_angles[0][0]))
phase_angles = []
print ("Calculate phase:")
center = center - int(y1) #if calibration is available. 0-180
for idx, images_angle in enumerate(images_angles):
#loop trough the needed heights
print (idx)
alpha = unwrap(images_angle[0][slices,:], images_angle[1][slices,:])[:,int(y1):int(y2)]
alpha = correct_phase_jumps(alpha) #correct for phases larger than 2pi
sh=217
phase_l = phase(alpha, alpha, shift=sh)
phase_r = -phase(alpha[:,::-1], alpha[:,::-1], shift=sh)[:,::-1]
phase_l[:,:-sh] = (phase_l[:,:-sh] + phase_r[:,sh:])/2.
phase_l = phase_l[:,center-Lx//2:center+Lx//2] - np.average(phase_l[:,0:y_calib])
if plot_phases:
plt.figure()
plt.plot(phase_l[0])
plt.show()
phase_angles.append(phase_l)
return phase_angles
def create_images(x1, x2):
from skimage import io
print ("Create images")
Ang=np.array([0,15,30,45,60,75,90])
Data = []
for R in Ang:
path1="data/1R" + str(R) + ".bmp"
path2="data/2R" + str(R) + ".bmp"
D = []
im1 = io.imread(path1, as_grey=True)
im2 = io.imread(path2, as_grey=True)
im1 = im1[x1:x2,:]
im2 = im2[x1:x2,:]
D.append(im1)
D.append(im2)
Data.append(D)
np.save("data/images_angles.npy",Data)
def reconstruct(images_angles, angles, x1, x2, y1, y2, Ly,slices = [1100]):
#create_images(x1, x2)
phase_angles = calculate_phase(images_angles ,y1 ,y2 , Ly, center = 1311, slices = slices, plot_images = False)
phase_angles = np.array(list(phase_angles) + [np.fliplr(ph) for ph in phase_angles[::-1]])
proj = tomopy.normalize(phase_angles, np.ones_like(phase_angles), np.zeros_like(phase_angles))
start = time.time()
rot_center = tomopy.find_center(proj, angles, init = int(proj.shape[2]/2), ind=0, tol=0.5)
print ('Center of rotation: ', rot_center)
start = time.time()
recon = tomopy.recon(proj, angles, center=rot_center, algorithm='sirt', ncore = 8)
recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
print ('Time for reconstruction: ',time.time() - start)
save_file(recon, plot = True)
......@@ -52,7 +52,7 @@ class ParabolicMirror:
for motor in self.motors:
motor.home()
def optimize(self, focus, laser, Cz, wanted_focus_um, n_av = 10, sleep_time = 0.2, output = True, maxiter = 100):
def optimize(self, focus, laser, Cz, wanted_focus_um, n_av = 1, sleep_time = 0.2, output = True, maxiter = 100):
#this variable is necessary to compensate the attenuation of the laser
global intmult
intmult = 1.
......@@ -77,11 +77,12 @@ class ParabolicMirror:
del measurement
if output:
print ('FWHMxs: ',FWHMxs)
print ('FWHMys: ',FWHMys)
print ('maxints: ',maxints)
print ('mean(maxits): ',np.mean(maxints))
res = np.sqrt(np.array(FWHMxs)**2 + np.array(FWHMys)**2)/np.sqrt(2)*4096/np.array(maxints)*intmult
mean, error = np.mean(res), np.std(res, ddof = -1)
res = (np.array(FWHMxs)**2 + np.array(FWHMys)**2)/2#*4096/np.array(maxints)*intmult
mean, error = np.mean(res), np.std(res, ddof = -1)/np.sqrt(len(res))
global reduce_int
if np.mean(maxints) > 4e3:
reduce_int = True
......@@ -116,7 +117,7 @@ class ParabolicMirror:
def opt_f(x):
y = x[0]
z = x[1]
zc = x[2]
zc = x[2]*30.
if output:
print ('-------------------------------------------')
change_parameter(y,z, zc)
......@@ -134,7 +135,7 @@ class ParabolicMirror:
focus.set_config(config0_focus)
x0 = [config0['pos_Pay'],
config0['pos_Paz'],
Cz.position()]
Cz.position()/30.]
#config_per0['deltaz_Cz']]
def after_it(res):
......@@ -149,21 +150,23 @@ class ParabolicMirror:
reduce_int = False
print ('laser attenuated')
# directly Nedler-Mead
import scipy.optimize
res = scipy.optimize.minimize(opt_f, x0, method='Nelder-Mead',callback = after_it, options={'disp': output,
'initial_simplex': None,
'maxiter': maxiter,
'xatol': 0.00005,
'return_all': False,
'fatol': wanted_focus_um})
'fatol': wanted_focus_um**2})
#import noisyopt
#res = noisyopt.minimizeSPSA(opt_f, x0,bounds = [[0.,6.],[0.,6.],[0.,2.5]], c = 0.001,a = 0.001, niter=maxiter,paired = False, disp=False, callback = after_it)
if output:
print (res)
return res
......@@ -32,7 +32,8 @@ im1: first image
im2: second image
'''
mode = {False: 'Off', True: 'On'}
modei = {'Off': False, 'On': True}
class PlasmaCamHorizontal:
#mendatory functions
def __init__(self, config):
......@@ -50,14 +51,14 @@ class PlasmaCamHorizontal:
def set_config(self,config):
self.cam.properties['ExposureTimeAbs']=config['ExposureTimeAbs']
self.cam.properties['TriggerMode']=config['TriggerMode']
self.cam.properties['TriggerMode']=mode[config['TriggerMode']]
self.y0 = config['y0']
self.z0 = config['z0']
def get_config(self):
config = {}
config['ExposureTimeAbs'] = self.cam.properties['ExposureTimeAbs']
config['TriggerMode'] = self.cam.properties['TriggerMode']
config['TriggerMode'] = modei[self.cam.properties['TriggerMode']]
config['y0'] = self.y0
config['z0'] = self.z0
return config
......
No preview for this file type
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment