Commit 56423c8e authored by Nick Sauerwein's avatar Nick Sauerwein
Browse files

used software successfull in experiment and changed settings for smaract stages

parent 4015efcc
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
File mode changed from 100644 to 100755
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
File mode changed from 100644 to 100755
......@@ -8,15 +8,39 @@ Created on Mon Dec 12 11:10:06 2016
import numpy as np
import matplotlib.pyplot as plt
R0=344.25
theta0 = 0./180 * np.pi
R0 = 1.06
theta1 = 40./180 * np.pi
R1 = 0.94
theta2 = 80./180 * np.pi
R2 = 0.81
z0 = R0 * np.exp(1j * theta0)
z1 = R1 * np.exp(1j * theta1)
z2 = R2 * np.exp(1j * theta2)
z = np.array([z0, z1, z2])
def fit_circle(x, y):
norms = np.array(x)**2 + np.array(y)**2
A = np.array([np.diff(x),np.diff(y)]).T
b = (np.diff(norms))/2.
if len(x) == 3:
z = np.linalg.solve(A,b)
return z
return np.linalg.lstsq(A,b)[0]
e = fit_circle(z.real, z.imag)
def cor(R):
theta=(R0-R)*2*np.pi/360.
ps2=1./268.92857 #mm
e_p=np.array([ 847.30650222, 502.90810925])
c_p=np.array([792,528])
e=e_p-c_p
return np.round((np.cos(theta)*e[0]-np.sin(theta)*e[1])*ps2,6)
theta= R*2*np.pi/360.
return np.round((np.cos(theta)*e[0]-np.sin(theta)*e[1]),6)
#th=R0-np.linspace(0, 360, 100)
......
......@@ -121,13 +121,13 @@ def movetoL(pos, direction):
def init():
cmd("SCLF0,500")#set close-loop frequency!!
cmd("SCLF1,200")
#cmd("SCLS0,100000")
cmd("SCLF1,10000")
cmd("SCLS0,100000000")
cmd("SCLS1,100000000")
print("find rot. ref.")
moverelR(-1)
time.sleep(5)
time.sleep(1)
frmR("fwd")
while cmd("GS0",output=True)[-2]!=str(3):
......@@ -148,6 +148,9 @@ def init():
while cmd("GS1",output=True)[-2]!=str(3):
#print(cmd("GS0",output=True)[-2])
time.sleep(1)
if cmd("GS1",output=True)[-2] == str(0):
print ('lin. stage got stuck... continue')
break
print("L",getposL())
print("lin. ref. done")
......
......@@ -11,6 +11,8 @@ import matplotlib.pyplot as plt
import Devices.InterferometerTools.routines as r
from Devices.CameraTools import binArray
import time as t
# camera setup
......@@ -32,6 +34,7 @@ zmin: minimal z in mm
z0: offset of z to fit blade edge
x0: offset of x to fit blade edge
Ly: width to do tomography in mm
y_calib: length to take for calibration left and right of image in mm
ExposureTimeAbs: exposure time in us
BitsOfImage: 8 or 12 (8 if you want to do fast images, for tomography use 12!)
TriggerMode: On/Off
......@@ -93,6 +96,7 @@ class Interferometer:
self.x_tomo = config['x_tomo']
self.shift = config['shift']
self.live_preview = config['live_preview']
self.y_calib = config['y_calib']
def get_config(self):
......@@ -112,6 +116,7 @@ class Interferometer:
config['x_tomo'] = self.x_tomo
config['shift'] = self.shift
config['live_preview'] = self.live_preview
config['y_calib'] = self.y_calib
return config
......@@ -139,18 +144,21 @@ class Interferometer:
self.plot_server_tomo = Plot2DServer({'title':'Tomography live', 'xlabel':'z [mm]', 'ylabel': 'x [mm]','image': True})
blade_config0 = self.blade.get_config()
global t_begin
t_begin = t.time()
angles_deg =np.round(90 *(np.cos(np.linspace(0, np.pi/2, self.num_dir)) - 1),3)# 90 * np.linspace(0, 1, self.num_dir)#(np.cos(np.linspace(0, np.pi/2, self.num_dir)) - 1)
angles_deg = 90 * np.linspace(0, 1, self.num_dir) # np.round(90 *(np.cos(np.linspace(0, np.pi/2, self.num_dir)) - 1),3)#(np.cos(np.linspace(0, np.pi/2, self.num_dir)) - 1)
if self.take_data:
images_angles = self.tomography_take_data(blade_config0, angles_deg, save = True)
else:
dict_images_angles = np.load("Data/Interferometer/images_angles.npy").item()
images_angles = dict_images_angles['images_angles']
angles_deg = dict_images_angles['angles']
del dict_images_angles
self.save_phases(images_angles, angles_deg)
res = update_tomography(self, images_angles,-angles_deg/360*2*np.pi)
......@@ -224,24 +232,33 @@ class Interferometer:
if self.plot_server_im == None:
self.plot_server_im = Plot2DServer({'title':'Interferometer Image', 'xlabel':'z [mm]', 'ylabel': 'x [mm]','image': True})
if self.plot_server_ph == None:
self.plot_server_ph = Plot2DServer({'title':'Interferometer Phase', 'xlabel':'z [mm]', 'ylabel': 'x [mm]','image': False})
self.plot_server_ph = Plot2DServer({'title':'Interferometer Phase', 'xlabel':'z [mm]', 'ylabel': 'x [mm]','image': True})
im1 = measurement['im1']
im2 = measurement['im2']
shape = measurement['im1'].shape
im1 = binArray(im1, 0, self.n_binning, self.n_binning, np.sum)
im1 = binArray(im1, 1, self.n_binning, self.n_binning, np.sum)
im1bx1 = binArray(im1, 0, self.n_binning, self.n_binning, np.sum)
im1bxy = binArray(im1, 1, self.n_binning, self.n_binning, np.sum)
self.plot_server_im.update(np.flipud(im1bxy),extent=self.psize*np.array([0.,shape[1], 0., shape[0]])-np.array([z0, z0, x0,x0]))
self.plot_server_im.update(np.flipud(im1),extent=self.psize*np.array([0.,shape[1], 0., shape[0]])-np.array([z0, z0, x0,x0]))
alpha = self.Alpha(measurement['im1'],measurement['im2'])
im1bx = binArray(im1[-self.x2p:-self.x1p,:], 0, self.n_binning, self.n_binning, np.sum)
im2bx = binArray(im2[-self.x2p:-self.x1p,:], 0, self.n_binning, self.n_binning, np.sum)
alpha = binArray(alpha, 0, self.n_binning, self.n_binning, np.sum)
alpha = binArray(alpha, 1, self.n_binning, self.n_binning, np.sum)
from Devices.InterferometerTools.tomography import calculate_phase
self.plot_server_ph.update(np.flipud(alpha), extent=[zmin, zmax, xmin, xmax])
phase = calculate_phase([[im1bx, im2bx]], int(self.Ly/self.psize),center = int(self.z0/self.psize),y_calib = int(self.y_calib/self.psize), plot_images = True, plot_phases = False,shift=self.shift)[0]
print (phase)
phase = binArray(phase, 1, self.n_binning, self.n_binning, np.sum)
self.plot_server_ph.update(np.flipud(-phase), extent=[-self.Ly/2,self.Ly/2, xmin, xmax])
def turn_off(self):
self.cam.close()
......@@ -301,10 +318,10 @@ class Interferometer:
def tomography_take_data(self, blade_config0, angles, save = True):
images_angles = []
global t_begin
blade_config1 = blade_config0.copy()
import time
from Devices.InterferometerTools.tomography import reconstruct
from multiprocessing import Process
......@@ -317,10 +334,8 @@ class Interferometer:
print ('Angle ',angle)
blade_config1['pos_R'] = blade_config0['pos_R'] + angle
print ('start movement')
self.blade.set_config(blade_config1)
print ('angle set')
print ('Time [Drive to angle ',angle,'] = ', t.time() - t_begin)
if p:
p.join()
......@@ -328,11 +343,11 @@ class Interferometer:
for j in range(self.num_av):
time.sleep(0.2)
t.sleep(0.2)
print ('Nummer ',j )
self.delay.measure()
D = []
for img in self.cam.grab_images(2):
D.append(np.fliplr(img))
Ds.append(D)
......@@ -344,6 +359,7 @@ class Interferometer:
if self.live_preview and i != len(angles)-1:
p = Process(target = update_tomography, args = (self,images_angles,-angles[:i + 1]/360*2*np.pi))
p.start()
print ('Time [Take images at angle ',angle,'] = ', t.time() - t_begin)
if save:
......@@ -352,14 +368,33 @@ class Interferometer:
self.blade.set_config(blade_config0)
print ('Time [Drive to angle 0] = ', t.time() - t_begin)
return images_angles
def save_phases(self, images_angles, angles):
from Devices.InterferometerTools.tomography import calculate_phase
phase_angles = calculate_phase(images_angles, int(self.Ly/self.psize), center = int(self.z0/self.psize) , plot_images = True, shift = self.shift,y_calib = int(self.y_calib/self.psize)) #1520
del images_angles
if len(angles) >= 2 and angles[-1] != np.pi/2:
phase_angles = np.array(list(phase_angles) + [np.fliplr(ph) for ph in phase_angles[::-1][:]])
angles = list(angles) + list(np.pi - angles[::-1][:])
elif len(angles) >= 2:
phase_angles = np.array(list(phase_angles) + [np.fliplr(ph) for ph in phase_angles[::-1][1:]])
angles = list(angles) + list(np.pi - angles[::-1][1:])
#save
np.save("Data/Interferometer/phase_angles.npy", {'phase_angles': phase_angles, 'angles': angles})
def update_tomography(self, images_angles, angles_rad):
from Devices.InterferometerTools.tomography import reconstruct
res = reconstruct(images_angles,angles_rad, int(self.Ly/self.psize), slices = [images_angles[0][0].shape[0] - int((self.x_tomo + self.x0)/self.psize)], shift = self.shift, output = self.live_preview)
res = reconstruct(images_angles,angles_rad, int(self.Ly/self.psize),center =int(self.z0/self.psize),y_calib = int(self.y_calib/self.psize) , slices = [images_angles[0][0].shape[0] - int((self.x_tomo + self.x0)/self.psize)], shift = self.shift, output = self.live_preview)
A_Ar=4.203e-6 # molar refractifity of Argon, used to go from phase to density
NA= 6e23 #Avogadro const 1/mol
......
......@@ -3,7 +3,7 @@ import tomopy
import numpy as np
import matplotlib.pyplot as plt
import time
import time as t
def calculate_phase(images_angles, Ly,center = 230, y_calib=200, slices=None, plot_images = False, plot_phases = False,shift=217):
......@@ -14,10 +14,10 @@ def calculate_phase(images_angles, Ly,center = 230, y_calib=200, slices=None, pl
if plot_images:
plt.figure()
plt.imshow(images_angles[-1][0])
plt.imshow(images_angles[-1][0][:,int(y1):int(y2)])
plt.show()
if slices==None:
if not slices:
slices=range(len(images_angles[0][0]))
phase_angles = []
......@@ -26,31 +26,19 @@ def calculate_phase(images_angles, Ly,center = 230, y_calib=200, slices=None, pl
center = center - int(y1) #if calibration is available. 0-180
for idx, images_angle in enumerate(images_angles):
#loop trough the needed heights
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
#phase_e = lambda shiftt: np.mean(phase(alpha, alpha, shift=int(shiftt))[:,-y_calib:])**2
from scipy import ndimage
alpha = ndimage.filters.gaussian_filter(alpha,[0,15], mode="nearest", truncate=2)
#shiftarr = np.linspace(shift - 100, shift + 100, 201)
#arg = np.argmin([phase_e(sh) for sh in shiftarr])
#print (shiftarr[arg])
#from scipy import ndimage
#alpha = ndimage.filters.gaussian_filter(alpha,[0,15], mode="nearest", truncate=2)
phase_l = phase(alpha, alpha, shift=shift)
#plt.figure()
#plt.plot(phase_l[0,:])
phase_r = -phase(alpha[:,::-1], alpha[:,::-1], shift=shift)[:,::-1]
#plt.figure()
#plt.plot(phase_r[0,shift:])
#plt.plot(phase_l[0,:-shift])
phase_l = (phase_l[:,:-shift] + phase_r[:,shift:])/2.
x = np.arange(phase_l.shape[1])
......@@ -66,24 +54,23 @@ def calculate_phase(images_angles, Ly,center = 230, y_calib=200, slices=None, pl
phase_l = phase_l[:,y_calib:-y_calib]
if idx == 0 and plot_images:
if idx == len(images_angles) -1 and plot_images:
plt.figure()
plt.title('ta')
plt.plot(phase_l[0,:])
plt.plot(phase_l[0,::-1])
plt.title('phase at half hight x = (xmax - xmin)/2 and last angle ')
plt.plot(phase_l[phase_l.shape[0]//2,:])
plt.plot(phase_l[phase_l.shape[0]//2,::-1])
plt.vlines([len(phase_l[0])//2], -1 , 1)
phase_angles.append(phase_l)
plt.show()
if plot_images or plot_phases:
plt.show()
print (np.sum(phase_angles, axis = 2))
return phase_angles
def reconstruct(images_angles, angles, Ly,slices = [1100], shift = 217, output = False):
phase_angles = calculate_phase(images_angles, Ly, center = 1460 , slices = slices, plot_images = True, shift = shift) #1520
def reconstruct(images_angles, angles, Ly,center, y_calib = 200, slices = [1100], shift = 217, output = False):
phase_angles = calculate_phase(images_angles, Ly, center = center , slices = slices, plot_images = True, plot_phases = True, shift = shift,y_calib = y_calib)
......@@ -116,12 +103,12 @@ def reconstruct(images_angles, angles, Ly,slices = [1100], shift = 217, output =
#np.save("Data/Interferometer/data.npy", {'projections': proj, 'angles': angles})
#print ('saved')
start = time.time()
start = t.time()
recon = tomopy.recon(proj, angles, center=int(proj.shape[2]/2), algorithm='mlem', ncore = 8)
recon = tomopy.circ_mask(recon, axis=0, ratio=0.95)
print ('Time for reconstruction: ',time.time() - start)
print ('Time for reconstruction: ',t.time() - start)
return recon
......@@ -89,18 +89,18 @@ class PlasmaCamHorizontal:
if self.plot_server_im == None:
self.plot_server_im = Plot2DServer({'title':'PlasmaCamHorizontal', 'xlabel':'z [mm]', 'ylabel': 'y [mm]', 'image': True})
im1 = measurement['im']
im1 = measurement['im'].astype('float32')
try:
im_c = measurement['im_c']
im_c = measurement['im_c'].astype('float32')
im1 = im1 - im_c
except:
print ('not calibrated')
shape = im1.shape
im1 = binArray(im1, 0, self.n_binning, self.n_binning, np.sum)
im1 = binArray(im1, 1, self.n_binning, self.n_binning, np.sum)
im1 = binArray(im1, 0, self.n_binning, self.n_binning, np.mean)
im1 = binArray(im1, 1, self.n_binning, self.n_binning, np.mean)
self.plot_server_im.update(np.flipud(im1),extent=self.psize*np.array([0.,shape[1], 0., shape[0]]) - np.array([self.z0, self.z0, self.y0, self.y0]))
......@@ -119,7 +119,7 @@ class PlasmaCamHorizontal:
print(cam)
for availca in available_cameras:
if str(availca)[-3] == '3':
if str(availca)[-3] == '3' and str(availca)[-4] == '3':
self.cam_to_connect = availca
cam = pypylon.factory.create_device(availca)
print ('camera found =) Be happy')
......
......@@ -95,21 +95,21 @@ class ScreenCamera:
if self.plot_server_im == None:
self.plot_server_im = Plot2DServer({'title':'ScreenCamera', 'xlabel':'z [mm]', 'ylabel': 'y [mm]', 'image': False})
im1 = measurement['im']
im1 = measurement['im'].astype('float32')
try:
im_c = measurement['im_c']
im_c = measurement['im_c'].astype('float32')
im1 = im1 - im_c
except:
print ('not calibrated')
im1 = binArray(im1, 0, self.n_binning, self.n_binning, np.sum)
im1 = binArray(im1, 1, self.n_binning, self.n_binning, np.sum)
im1 = binArray(im1, 0, self.n_binning, self.n_binning, np.mean)
im1 = binArray(im1, 1, self.n_binning, self.n_binning, np.mean)
shape = im1.shape
self.plot_server_im.update(np.flipud(im1.astype('int8')),extent=np.array([0.,shape[1], 0., shape[0]]) - np.array([self.z0, self.z0, self.y0, self.y0]))
self.plot_server_im.update(np.flipud(im1),extent=np.array([0.,shape[1], 0., shape[0]]) - np.array([self.z0, self.z0, self.y0, self.y0]))
def turn_off(self):
self.cam.close()
......@@ -128,7 +128,7 @@ class ScreenCamera:
print(cam)
for availca in available_cameras:
if str(availca)[-3] == '2':
if str(availca)[-3] == '3' and str(availca)[-4] == '8':
self.cam_to_connect = availca
cam = pypylon.factory.create_device(availca)
print ('camera found =) Be happy')
......@@ -140,5 +140,5 @@ class ScreenCamera:
#cam.properties['GevSCPSPacketSize']=9014
cam.properties['PixelFormat'] = 'Mono16'
cam.properties['PixelFormat'] = 'Mono12'
self.cam = cam
\ No newline at end of file
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