sf-mx/SwissMX
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

towards new pix2pos calibration part 2

Browse Source
This commit is contained in:
zamofing_t
2022-07-20 10:50:01 +02:00
parent 7209243f2f
commit 4a6db36805
5 changed files with 251 additions and 159 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
PrelocatedCoordinatesModel.py
View File

@@ -6,7 +6,7 @@ from typing import List
import numpy as np
import pandas as pd
import transformations as tfs
#import transformations as tfs
from PyQt5.QtCore import Qt, QFileInfo, pyqtSignal, pyqtSlot
from PyQt5.QtWidgets import (

4
Readme.md
View File

@@ -155,3 +155,7 @@ python swissmx.py
caput SAR-EXPMX:MOT_FX.RBV -2.700
caput SAR-EXPMX:MOT_FY.RBV -7.450
caput SAR-EXPMX:MOT_CZ.RBV 1.509
#ZAC: orig. code
git dt 98297263 swissmx.py
git dt 7445a5aa CustomROI.py app_config.py app_utils.py epics_widgets/MotorTweak.py epics_widgets/SmaractMotorTweak.py

10
app_config.py
View File

@@ -1,3 +1,11 @@
#TODO:
# currently 2 settings/configs are used
# QSettings -> cat ~/.config/Paul\ Scherrer\ Institut/SwissMX.conf
# yaml -> swissmx.yaml
# QSettings are changed by program
# #yaml is fixed and not altened by program
import yaml
from PyQt5.QtCore import QSettings
from pathlib import Path
@@ -6,7 +14,7 @@ import logging
logger = logging.getLogger(__name__)
settings = QSettings("Paul Scherrer Institut", "SwissMX")
settings = QSettings("PSI", "SwissMX")
inst_folder = Path(__file__).absolute().parent

341
geometry.py Normal file → Executable file
View File

@@ -11,16 +11,22 @@ coordinate systems, optical center, xray axis, pixel sizes etc.
ppm= pixel per mm
update_pix2pos
modes:
0x01: update_pix2pos
0x02: update_optical_center
'''
import logging
import numpy as np
_log = logging.getLogger(__name__)
_log=logging.getLogger(__name__)
class geometry:
def __init__(self):
pass
def find_optical_center(self,p):
def find_optical_center(self, p):
# p is an array of
# at zoom out: (p1x,p1y),(p2x,p2y),(p3x,p3y),...
# at zoom in : (p1x,p1y),(p2x,p2y),(p3x,p3y),...
@@ -31,27 +37,144 @@ class geometry:
# this coordinate represents also the origin of other coordinates
pass
def interp_zoom(self,zoom):
# this calculates the interpoated pix2pos coordinate transformation matrix
def update_optical_center(self, meas, debug=False):
if debug:
import sys
import matplotlib.pyplot as plt
if sys.gettrace(): # in debuggingmode
plt.interactive(True) #uncomment to debug
print('DEBUG MODE')
fig, ax=plt.subplots()
ax.invert_yaxis()
ax.grid()
ax.axis('equal')
for k, v in opt_ctr_meas.items():
v=np.array(v)
# ax.plot(v[:, 0], v[:, 1])
ax.fill(v[:, 0], v[:, 1], fill=False)
numPoints=len(next(iter(meas.values())))
numZoom=len(meas)
M=np.ndarray((numPoints, numZoom, 2), np.float)
K=np.array(tuple(meas.keys()), np.float)
for i, (k, v) in enumerate(meas.items()):
M[:, i, :]=np.array(v)
if debug:
for i in range(numPoints):
ax.plot(M[i, :, 0], M[i, :, 1])
# parametrizex lines:
# (USE ONLY 2 zoom levels (first and last!!!)
# (x,y)=(p0x+a*q, p0y+b*q, q=0 -> (p0x,p0y), q=1 (p1x,p1y)
# -> p1x=p0x+a*1 -> a=p1x-p0x
# -> p1y=p0y+b*1 -> b=p1y-p0y
AB=M[:,-1,:]-M[:,0,:] # parameters a1..an b1..bn
# least square fitting for x,y,q
# x=p0x+a0*q -> x-a0*q=p0x
# y=p0y+b0*q
# x=pnx+an*q
# y=pny+bn*q
# [1 0 -a0] [x] [p0x]
# [.. .. ..] * [y] = [..]
# [1 0 -an] [q] [pnx]
# [0 1 -b0] [p0y]
# [.. .. ..] [..]
# [0 1 -bn] [pny]
# A * aa = y
A=np.ndarray((numPoints*2, 3), np.float)
A[:numPoints,:2]=(1,0)
A[numPoints:,:2]=(0,1)
A[:,2]=-AB.T.ravel()
A=np.asmatrix(A)
y=np.asmatrix(M[:, 0, :].T.ravel()).T
aa=(A.T*A).I*A.T*y
if debug:
# plot least square line fitting
ax.autoscale(False)
P=np.ndarray((2, 2), np.float)
for i in range(numPoints):
q=aa[2]
P[0, :]=M[i, 0, :]
P[1, :]=M[i, 0, :]+q*AB[i,:]
ax.plot(P[:, 0], P[:, 1])
ax.plot(aa[0], aa[1], marker='*',color='r')
plt.show()
self._opt_ctr=np.asarray(aa[:2]).ravel()
_log.debug('optical center:{}'.format(tuple(self._opt_ctr)))
# # line fitting y= ax+b
# # calculate line fitting (least square)
# # y= a*x+b
# # [x1 1] [a] [y1]
# # [.. 1] * [b] = [..]
# # [xn 1] [yn]
# # A * aa = y
# # calculate least square line fitting for each point
# A=np.asmatrix(np.ndarray((numZoom, 2), np.float))
# A[:, 1]=1
# AA=np.ndarray((numPoints, 2), np.float) # fitting results a1..an b1..bn
# for i in range(numPoints):
# A[:, 0]=M[i, :, 0].reshape(-1, 1)
# y=np.asmatrix(M[i, :, 1]).T
# aa=(A.T*A).I*A.T*y
# AA[i]=aa.reshape(2)
# if debug:
# # plot least square line fitting
# ax.autoscale(False)
# rngx=(M[:, :, 0].min(), M[:, :, 0].max())
# Y=np.stack((AA[:, 0]*rngx[0]+AA[:, 1], AA[:, 0]*rngx[1]+AA[:, 1]))
# P=np.ndarray((2, 2), np.float)
# P[:, 0]=rngx
# for i in range(numPoints):
# P[:, 1]=Y[:, i]
# ax.plot(P[:, 0], P[:, 1])
# plt.show()
# # TODO: calculate all cross points of least square line fitting
# # a1*x + b1 = y
# # a2*x + b2 = y
# # [a1 -1] * [x] [-b1]
# # [a2 -1] * [y] = [-b2]
# P=np.ndarray(((numPoints-1)*numPoints//2, 2), np.float)
# k=0
# D=np.ndarray((2, 2), np.float)
# D[:, 1]=-1
# # for i in range(numPoints):
# # for j in range(i+1,numPoints):
# # P[k,:]=
# # k+=1
def interp_zoom(self, zoom):
# this calculates the interpolated pix2pos coordinate transformation matrix
# the returned value is a 2x2 matrix:
# [a b] [mx] [px]
# [ ]*[ ]=[ ]
# [c d] [my] [py]
K,AA=self._lut_pix2pos
pix2pos=self._pix2pos=np.asmatrix(np.ndarray((2,2),np.float))
pix2pos[0,0]=np.interp(zoom,K,AA[:,0,0])
pix2pos[0,1]=np.interp(zoom,K,AA[:,0,1])
pix2pos[1,0]=np.interp(zoom,K,AA[:,1,0])
pix2pos[1,1]=np.interp(zoom,K,AA[:,1,1])
K, AA=self._lut_pix2pos
pix2pos=self._pix2pos=np.asmatrix(np.ndarray((2, 2), np.float))
pix2pos[0, 0]=np.interp(zoom, K, AA[:, 0, 0])
pix2pos[0, 1]=np.interp(zoom, K, AA[:, 0, 1])
pix2pos[1, 0]=np.interp(zoom, K, AA[:, 1, 0])
pix2pos[1, 1]=np.interp(zoom, K, AA[:, 1, 1])
_log.debug('interpolation array:{}'.format(tuple(pix2pos.ravel())))
def update_pix2pos(self,meas):
#calculates _lut_z2p out of measurements
# the _lut_z2p is dictionaty a lookuptable
# zoom {1,200,400,600,800,1000}
#[pxx pxy]
#[pyx pyy]
def update_pix2pos(self, meas):
# calculates _lut_pix2pos out of measurements
# meas is a structure as in the sample code
# [pxx pxy]
# [pyx pyy]
#
# [a b] [mx] [px]
# [ ]*[ ]=[ ]
@@ -90,25 +213,32 @@ class geometry:
# a,b,c,d=np.asarray(aa).ravel()
#
# aa=aa.reshape(2,-1)
AA=np.ndarray((len(meas),2,2),np.float)
K=np.array(tuple(meas.keys()),np.float)
self._lut_pix2pos=(K,AA)
for i,(k,v) in enumerate(meas.items()):
m=np.array(v) # measurements
d=m[1:]-m[0] # distances
_log.debug('raw data\nmeas:{}'.format(meas))
if len(meas)<2:
_log.error('need at least 2 zoom levels:\nmeas:{}'.format(meas))
return
AA=np.ndarray((len(meas), 2, 2), np.float)
K=np.array(tuple(meas.keys()), np.float)
self._lut_pix2pos=(K, AA)
for i, (k, v) in enumerate(meas.items()):
m=np.array(v) # measurements
if m.shape[0]<3:
_log.error('zoom:{}: need at least 3 points per zoom levels:\nmeas:{}'.format(k, v))
return
d=m[1:]-m[0] # distances
A=np.zeros((d.shape[0]*2,d.shape[1]),np.float)
A[:d.shape[1]-1,:2]=A[d.shape[1]-1:,2:]=d[:,:2]
A=np.zeros((d.shape[0]*2, d.shape[1]), np.float)
A[:d.shape[0], :2]=A[d.shape[0]:, 2:]=d[:, :2]
y=d[:,2:].T.ravel()
y=d[:, 2:].T.ravel()
y=np.asmatrix(y).T
A=np.asmatrix(A)
aa=(A.T*A).I*A.T*y
AA[i]=aa.reshape(2,-1)
#bx_coefs = np.polyfit(nbx[0], nbx[1], 3)
#np.poly1d(bx_coefs)
_log.debug('least square data:\nK:{}\nAA:{}'.format(K,AA))
AA[i]=aa.reshape(2, -1)
# bx_coefs = np.polyfit(nbx[0], nbx[1], 3)
# np.poly1d(bx_coefs)
_log.debug('least square data:\nK:{}\nAA:{}'.format(K, AA))
def autofocus(self):
# cam camera object
@@ -119,14 +249,14 @@ class geometry:
# mode mode to calculate sharpness (sum/max-min/hist? of edge detection in roi)
pass
def pix2pos(self,p,zoom=None):
def pix2pos(self, p, zoom=None):
# returns the position m(x,y) in meter relative to the optical center at a given zoom level of the pixel p(x,y)
# if zoom=None, the last zoom value is used
if zoom is not None:
self.interp_zoom(zoom)
return np.asarray(self._pix2pos*np.mat(p).T).ravel()
def pos2pix(self,p,zoom=None):
def pos2pix(self, p, zoom=None):
# returns the pixel p(x,y) of the position m(x,y) in meter relative to the optical center at a given zoom level
# if zoom=None, the last zoom value is used
if zoom is not None:
@@ -138,96 +268,89 @@ class geometry:
pass
if __name__ == "__main__":
if __name__=="__main__":
import argparse
logging.basicConfig(level=logging.DEBUG,format='%(levelname)s:%(module)s:%(lineno)d:%(funcName)s:%(message)s ')
parser = argparse.ArgumentParser()
parser.add_argument('-m', '--mode', help='mode')
parser.add_argument("-t", "--test", help="test sequence", action="store_true")
logging.basicConfig(level=logging.DEBUG, format='%(levelname)s:%(module)s:%(lineno)d:%(funcName)s:%(message)s ')
args = parser.parse_args()
#(h, t)=os.path.split(sys.argv[0]);cmd='\n '+(t if len(h)>20 else sys.argv[0])+' '
#exampleCmd=('', '-m0xf -v0')
epilog=__doc__ #+'\nExamples:'+''.join(map(lambda s:cmd+s, exampleCmd))+'\n'
parser=argparse.ArgumentParser(epilog=epilog, formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("-m", "--mode", type=lambda x: int(x,0), help="mode (see bitmasks) default=0x%(default)x", default=0xff)
args=parser.parse_args()
_log.info('Arguments:{}'.format(args.__dict__))
# recorded data:
# x y x y
#Zoom 1 x+1.2 97 543 -> 1100 507
# y+1.0 607 941 -> 575 106
#Zoom 200 x+0.6 93 504 -> 853 510
# y+0.4 475 897 -> 472 157
#Zoom 400 x+0.5 88 615 -> 1094 579
# y+0.4 705 991 -> 673 190
#Zoom 600 x+0.3 32 460 -> 103 416
# y+0.25 551 937 -> 520 106
#Zoom 800 x+0.18 65 524 -> 1050 484
# y+0.14 632 946 -> 602 168
#Zoom 1000 x+0.1 121 632 -> 1044 592
# y+0.08 593 883 -> 570 145
measure={
1:{'x': (+1.2 , ( 97, 543),(1100, 507)),
'y': (+1.0 , ( 607, 941),( 575, 106))},
200:{'x': (+0.6 , ( 93, 504),( 853, 510)),
'y': (+0.4 , ( 475, 897),( 472, 157))},
400:{'x': (+0.5 , ( 88, 615),(1094, 579)),
'y': (+0.4 , ( 705, 991),( 673, 190))},
600:{'x': (+0.3 , ( 32, 460),( 103, 416)),
'y': (+0.25 , ( 551, 937),( 520, 106))},
800:{'x': (+0.18 , ( 65, 524),(1050, 484)),
'y': (+0.14 , ( 632, 946),( 602, 168))},
1000:{'x': (+0.1 , ( 121, 632),(1044, 592)),
'y': (+0.08 , ( 593, 883),( 570, 145))},
}
measure={
1: [(-3.0 , -7.6 , 116.99110193046, 632.5463827525),
(-2.0 , -7.6 , 934.07501517856, 600.7926167715),
(-2.0 , -7.1 , 916.54131238102, 191.0366615002),
(-3.0 , -7.1 , 103.74668003329, 226.2150231456)],
200: [(-3.1 , -7.3 , 113.66321353086, 824.9041423107),
(-2.3 , -7.3 , 1065.97386092697, 792.2851118419),
(-2.3 , -6.7 , 1033.68452410347, 74.0336610693),
(-3.1 , -6.7 , 84.62681572700, 116.6832971512)],
400: [(-3.4 , -6.7 , 155.00053674203, 601.3838942136),
(-3.0 , -6.7 , 957.95919656052, 573.0827012272),
(-3.2 , -6.5 , 541.08684037200, 187.9171307943),
(-3.2 , -6.8 , 564.32152887203, 789.1146957326)],
600: [(-3.3 , -6.8 , 328.27244399903, 509.5061192017),
(-3.1 , -6.8 , 992.78996735279, 488.0323963092),
(-3.2 , -6.9 , 672.03111567934, 832.4122409755),
(-3.2 , -6.7 , 645.70960116180, 164.2534779331)],
800: [(-3.2 , -6.7 , 639.52253576449, 53.4455632943),
(-3.2 , -6.85, 671.47023245203, 882.6335091391),
(-3.3 , -6.75, 105.12470026379, 361.3051859197),
(-3.1 , -6.75, 1195.96864609255, 313.1068618673)],
1000:[(-3.25, -6.75, 195.05641095116, 353.3492286375),
(-3.15, -6.75, 1117.27204644084, 314.9636405871),
(-3.2 , -6.8 , 675.10991143017, 790.3040145281),
(-3.2 , -6.72, 638.98580653116, 59.3803912957)]}
logging.getLogger('matplotlib').setLevel(logging.INFO)
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
np.set_printoptions(suppress=True)
np.set_printoptions(linewidth=196)
obj=geometry()
obj.update_pix2pos(measure)
obj.interp_zoom(1)
obj.interp_zoom(200)
obj.interp_zoom(100)
print(np.asarray(obj._pix2pos*np.mat((1,0)).T).ravel())
print(obj._pix2pos*np.mat((0,1)).T)
print(obj._pix2pos*np.mat((2,.1)).T)
if args.mode&0x01:
measure={
1:[(-3.0, -7.6, 116.99110193046, 632.5463827525),
(-2.0, -7.6, 934.07501517856, 600.7926167715),
(-2.0, -7.1, 916.54131238102, 191.0366615002),
(-3.0, -7.1, 103.74668003329, 226.2150231456)],
200:[(-3.1, -7.3, 113.66321353086, 824.9041423107),
(-2.3, -7.3, 1065.97386092697, 792.2851118419),
(-2.3, -6.7, 1033.68452410347, 74.0336610693),
(-3.1, -6.7, 84.62681572700, 116.6832971512)],
400:[(-3.4, -6.7, 155.00053674203, 601.3838942136),
(-3.0, -6.7, 957.95919656052, 573.0827012272),
(-3.2, -6.5, 541.08684037200, 187.9171307943),
(-3.2, -6.8, 564.32152887203, 789.1146957326)],
600:[(-3.3, -6.8, 328.27244399903, 509.5061192017),
(-3.1, -6.8, 992.78996735279, 488.0323963092),
(-3.2, -6.9, 672.03111567934, 832.4122409755),
(-3.2, -6.7, 645.70960116180, 164.2534779331)],
800:[(-3.2, -6.7, 639.52253576449, 53.4455632943),
(-3.2, -6.85, 671.47023245203, 882.6335091391),
(-3.3, -6.75, 105.12470026379, 361.3051859197),
(-3.1, -6.75, 1195.96864609255, 313.1068618673)],
1000:[(-3.25, -6.75, 195.05641095116, 353.3492286375),
(-3.15, -6.75, 1117.27204644084, 314.9636405871),
(-3.2, -6.8, 675.10991143017, 790.3040145281),
(-3.2, -6.72, 638.98580653116, 59.3803912957)]}
obj.update_pix2pos(measure)
obj.interp_zoom(1)
obj.interp_zoom(200)
obj.interp_zoom(100)
if args.mode&0x02:
opt_ctr_meas={ # x,y = 0.02, -4.89
1000:[(1057.4251530483375, 116.10122290395591),
(117.84916300310408, 190.27827474963223),
(184.2181041281829, 963.2812360887852),
(1092.5616512910262, 899.514998537239)],
800:[(888.2494207687248, 203.2917926172947),
(329.96950424600305, 248.83910515411347),
(372.9141132092893, 708.2162858826),
(906.4683457834523, 675.6824912134438)],
600:[(781.5385742538922, 251.44180872764602),
(447.09116505496564, 264.4553265953085),
(471.81684900352445, 554.6567750441825),
(798.4561474818535, 535.1364982426887)],
400:[(722.9777438494109, 286.5783069703348),
(525.1722722609408, 295.68776947769857),
(535.5830865550707, 462.26079818377866),
(729.4845027832422, 450.5486321028824)],
200:[(689.1425973934884, 308.70128734536104),
(565.5141776506945, 307.39993555859473),
(574.6236401580583, 406.30267135282986),
(693.0466527537872, 399.79591241899857)],
1:[(673.5263759522934, 307.39993555859473),
(591.5412133860195, 308.70128734536104),
(595.4452687463182, 376.3715802572061),
(672.2250241655271, 373.7688766836736)]}
obj.update_optical_center(opt_ctr_meas, True)

53
swissmx.py
View File

@@ -75,8 +75,7 @@ TASK_PRELOCATED = "prelocated"
TASK_HELICAL = "helical"
TASK_EMBL = "embl"
import PrelocatedCoordinatesModel
import PrelocatedCoordinatesModel # ZAC: orig. code
from EmblModule import EmblWidget #ZAC: orig. code
from HelicalTable import HelicalTableWidget #ZAC: orig. code
#from Wigis import Spinner, Checkbox #ZAC: orig. code
@@ -798,7 +797,6 @@ class Main(QMainWindow, Ui_MainWindow):
grp.layout().addWidget(self._ppm_feature_size_spinbox, 0, 1)
self._ppm_calibration = but = QPushButton("Start calibration")
but.setCheckable(True)
#but.clicked.connect(lambda x: self.update_ppm_fitters())
but.clicked.connect(self.update_ppm_fitters)
grp.layout().addWidget(but, 1, 0, 1, 2)
help = QTextBrowser()
@@ -973,46 +971,6 @@ class Main(QMainWindow, Ui_MainWindow):
del app._ppmRaw
return
if calib is not None:
_log.info("received new calibration for PPM")
self._zoom_to_ppm = calib
# self._zoom_to_ppm = { # FIXME: eventually automate
# 1: 830,
# 100: 940,
# 200: 1220,
# 400: 1950,
# 600: 3460,
# 800: 5400,
# 900: 7150,
# 1000: 9200,
# }
num_points = len(self._zoom_to_ppm)
if num_points < 2:
return
elif num_points < 4:
order = 2
elif num_points < 6:
order = 3
else:
order = 4
_log.debug("polynomial fitting using {} data points of order {}".format(num_points, order))
bx = [(z, ppm) for z, ppm in self._zoom_to_ppm.items()]
nbx = np.asarray(bx).T
bx_coefs = np.polyfit(nbx[0], nbx[1], order)
_log.debug(".... ppm fitting coeficients {}".format(bx_coefs))
self.ppm_fitter = np.poly1d(bx_coefs)
def getFastX(self):
return self.tweakers["fast_x"].get_position()
def getFastY(self):
return self.tweakers["fast_y"].get_position()
def zoom_changed_cb(self, value):
self.zoomChanged.emit(value)
try:
@@ -1125,8 +1083,6 @@ class Main(QMainWindow, Ui_MainWindow):
raw=app._ppmRaw[zoom]
except (AttributeError,KeyError) as e:
if type(e) is AttributeError:
_log.info("")
msg="SwissMX *POST-SAMPLE-TUBE* configuration is incomplete!!!"
QMessageBox.warning(self, "pix2pos calibration", "Press 'Start calibration' button first")
return
app._ppmRaw=dict()
@@ -1967,7 +1923,8 @@ class Main(QMainWindow, Ui_MainWindow):
def daq_grid_add_grid(self, gx=None, gy=None):
grid_index = len(self._grids)
if gx in (False, None):
gx, gy = self.getFastX(), self.getFastY()
gx=self.tweakers["fast_x"].get_position()
gy=self.tweakers["fast_y"].get_position()
xstep = self._sb_grid_x_step.value()
ystep = self._sb_grid_y_step.value()
xoffset = self._sb_grid_x_offset.value()
@@ -3174,8 +3131,8 @@ class Main(QMainWindow, Ui_MainWindow):
if SKIP_ESCAPE_TRANSITIONS_IF_SAFE not in keys:
settings.setValue(SKIP_ESCAPE_TRANSITIONS_IF_SAFE, False)
if EXPERIMENT_PGROUP not in keys:
self.update_user_and_storage()
#if EXPERIMENT_PGROUP not in keys:
# self.update_user_and_storage()
if "hits/marker_size" not in keys:
settings.setValue("hits/marker_size", 10)