#!/usr/bin/env python # *-----------------------------------------------------------------------* # | | # | Copyright (c) 2022 by Paul Scherrer Institute (http://www.psi.ch) | # | Based on Zac great first implementation | # | Author Thierry Zamofing (thierry.zamofing@psi.ch) | # *-----------------------------------------------------------------------* ''' coordinate systems, optical center, xray axis, pixel sizes etc. ppm= pixel per mm update_pix2pos ''' import logging import numpy as np _log = logging.getLogger(__name__) class geometry: def __init__(self): pass 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),... # # the pixel positions are given in chip pixel coordinates (0/0= top/right) # and ignore roi and binning # the returned (cx,cy) is the pixel koordinate that does not change with zoom # this coordinate represents also the origin of other coordinates pass def interp_zoom(self,zoom): # this calculates the interpoated 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]) _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] # # [a b] [mx] [px] # [ ]*[ ]=[ ] # [c d] [my] [py] # # [a*mx1 b*my1] [px1] # [ ... ... ] [...] # [a*mxn b*myn] =[pxn] # [c*mx1 d*my1] [py1] # [ ... ... ] [...] # [c*mxn d*myn] [pyn] # # [mx1 my1 0 0 ] [a] [px1] # [... .. 0 0 ] [b] [...] # [mxn myn 0 0 ]*[c]=[pxn] # [ 0 0 mx1 my1] [d] [py1] # [ 0 0 ... .. ] [...] # [ 0 0 mxn myn] [pyn] # # A *aa = y # https://de.wikipedia.org/wiki/Methode_der_kleinsten_Quadrate # r= A*aa-y # mit aa=(A.T*A)^-1*A.T*y # A=np.zeros((d.shape[0]*2,d.shape[1])) # A[:d.shape[1]-1,:2]=A[d.shape[1]-1:,2:]=d[:,:2] # # y=d[:,2:].T.ravel() # y=np.asmatrix(y).T # # A=np.asmatrix(A) # aa=(A.T*A).I*A.T*y # a,b,c,d=aa # 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 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] 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)) def autofocus(self): # cam camera object # mot motor object # rng region (min max relative to current position) to seek # n number of images to take in region # roi region of interrest to calculate sharpness # mode mode to calculate sharpness (sum/max-min/hist? of edge detection in roi) pass 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): # 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: self.interp_zoom(zoom) return np.asarray(self._pix2pos.I*np.mat(p).T).ravel() def optctr2xray(self): # returns the vector m(x,y) of the optical center to the xray pass 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") 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)]} import numpy as np 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)