#!/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 modes: 0x01: update_pix2pos 0x02: update_optical_center ''' import logging _log=logging.getLogger(__name__) import numpy as np import PIL.Image try: from scipy import ndimage, signal except ImportError as e: _log.warning(e) try: import cv2 as cv except ImportError as e: _log.warning(e) 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 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 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.float64) K=np.array(tuple(meas.keys()), np.float64) 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.float64) 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.float64) 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.float64)) # A[:, 1]=1 # AA=np.ndarray((numPoints, 2), np.float64) # 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.float64) # 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.float64) # k=0 # D=np.ndarray((2, 2), np.float64) # 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.float64)) 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(f'{zoom} -> {pix2pos.ravel()}') 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] [px] [mx] # [ ]*[ ]=[ ] # [c d] [py] [my] # # [a*px1 b*py1] [mx1] # [ ... ... ] [...] # [a*pxn b*pyn] =[mxn] # [c*px1 d*py1] [my1] # [ ... ... ] [...] # [c*pxn d*pyn] [myn] # # [px1 py1 0 0 ] [a] [mx1] # [... .. 0 0 ] [b] [...] # [pxn pyn 0 0 ]*[c]=[mxn] # [ 0 0 px1 py1] [d] [my1] # [ 0 0 ... .. ] [...] # [ 0 0 pxn pyn] [myn] # # 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) _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.float64) K=np.array(tuple(meas.keys()), np.float64) 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.float64) A[:d.shape[0], :2]=A[d.shape[0]:, 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) if np.diff(K).min()<0: #zoom need to be from lowest to highest ordered _log.debug('resorting') s=K.argsort() K=K[s] AA=AA[s] #resort A if needed self._lut_pix2pos=(K, AA) _log.debug('least square data:\nK:{}\nAA:{}'.format(K, AA)) @staticmethod def find_fiducial(cam,sz=(210,210),brd=(20,20)): if type(cam)==str: img=PIL.Image.open(cam) img=np.asarray(img) else: img=cam.get_image() # get latest image if img.dtype==np.uint16: img=np.array(img,np.uint8) fid=np.ones((sz[1]+2*brd[1],sz[0]+2*brd[0]),dtype=np.uint8)*255 fid[brd[1]:sz[1]+brd[1],brd[0]:sz[0]+brd[0]]=0 mask=np.ones((sz[1]+2*brd[1],sz[0]+2*brd[0]),dtype=np.uint8)*255 mask[2*brd[1]:sz[1],2*brd[0]:sz[0]]=0 #https://docs.opencv.org/4.5.2/d4/dc6/tutorial_py_template_matching.html #res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED ) res = cv.matchTemplate(img,fid,cv.TM_CCORR_NORMED,mask=mask) h,w=img.shape fh2,fw2=fid.shape fw2//=2 fh2//=2 mtr=np.ndarray((5,2),np.uint16) corr=np.ndarray((5,),np.float32) for i in range(5): p=np.unravel_index(res.argmax(), res.shape) corr[i]=res[p] mtr[i,:]=p y0=max(p[0]-fh2,0) y1=min(p[0]+fh2,h) x0=max(p[1]-fw2,0) x1=min(p[1]+fw2,w) res[y0:y1,x0:x1]*=.5 pos=mtr.mean(0)[::-1]+(fw2,fh2) crm=corr.mean() _log.debug(f'position: {pos} correlation:{crm}') return (pos,crm) 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.asmatrix(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.asmatrix(p).T).ravel() def optctr2xray(self): # returns the vector m(x,y) of the optical center to the xray pass @staticmethod def least_square_trf(points,fid=None,sort=True): # inputs are 4 points of a parallelogram # this function returns the least square transformation #[ px] [a b c ] [ q ] #[ py] =[d e f ]*[ r ] # [0 0 1 ] [ 1 ] # with (q,r) in [0,0],[0,1],[1,0],[1,1] # a b c d e f # [q00 r00 1 0 0 0 ] [a] [px00] # [0 0 0 q00 r00 1 ] [b] [py00] # [q01 r01 1 0 0 0 ]*[c]=[px01] # [0 0 0 q01 r01 1 ] [d] [py01] # [q10 r10 1 0 0 0 ] [e] [px10] # [0 0 0 q10 r10 1 ] [f] [py10] # [q11 r11 1 0 0 0 ] [px11] # [0 0 0 q11 r11 1 ] [py11] # # A *aa = y # A=np.array(( # (0,0,1,0,0,0), # (0,0,0,0,0,1), # (0,1,1,0,0,0), # (0,0,0,0,1,1), # (1,0,1,0,0,0), # (0,0,0,1,0,1), # (1,1,1,0,0,0), # (0,0,0,1,1,1)), np.float64) if fid is None: fid=np.array(((0,0),(0,1),(1,0),(1,1))) elif type(fid)!=np.ndarray: fid=np.array(fid) y=points # sort points p00 p01 p10 p11 if sort: # p01 ----------- p11 # | | # | | # p00-------------p10 def sort(a): s=a[:,0].argsort() a=a[s,:] s=a[:2,1].argsort() a[:2,:]=a[:2,:][s,:] s=a[2:,1].argsort() a[2:,:]=a[2:,:][s,:] return a y=sort(y) fid=sort(fid) A=np.ndarray((len(fid)*2,6)) i=0 for q,r in fid: A[i] =(q,r,1,0,0,0) A[i+1]=(0,0,0,q,r,1) i+=2 y=np.asmatrix(y.ravel()).T A=np.asmatrix(A) aa=(A.T*A).I*A.T*y aa=aa.reshape((2,3)) return aa def least_square_plane(self,points): #inputs are multiple (x,y,z) points # this function returns the parameters of least square fitted plane x=x,y=y,z=ax+bx+c # a b c # [px1 py1 1 ] [a] [pz1] # [px2 py2 1 ] [b] [pz2] # [... ... 1 ]*[c]=[...] # [pxn pyn 1 ] [pzn] A=np.ndarray((points.shape[0],3)) y=points[:,2] A[:,2]=1 A[:,0:2]=points[:,:2] y=np.asmatrix(y.ravel()).T A=np.asmatrix(A) try: aa=(A.T*A).I*A.T*y except np.linalg.LinAlgError as e: _log.warning(e) return aa=aa.A.ravel() if _log.level==logging.DEBUG: s=f'plane={aa[0]}X+{aa[1]}Y+{aa[2]}\n' for p in points: s+=f' {p}->{aa[0]*p[0]+aa[1]*p[1]+aa[2]}\n' _log.debug(s) self._fitPlane=aa class autofocus: def __init__(self,cam,mot=None): # cam camera object or image list # mot motor object (e.g. SimMotorTweak) or None # rng region (min max relative to current position) to seek self._cam=cam self._mot=mot def save_img(self,rng=(-1,1),n=30,progressDlg=None): mot=self._mot cam=self._cam try: p0=mot.get_rbv() except AttributeError: _log.warning('needs a real motor') return for i,p in enumerate(np.linspace(p0+rng[0],p0+rng[1],n)): mot.move_abs(p,wait=True) try: pic=cam._pic except AttributeError: pic=cam.get_image() mx=pic.max() if pic.max()>255: scl=2**int(round(np.log2(mx)-8)) pic=np.array(pic/scl, dtype=np.uint8) elif pic.dtype!=np.uint8: pic=np.array(pic, dtype=np.uint8) img=PIL.Image.fromarray(pic) fn=f'/tmp/image{i:03d}.png' img.save(fn) _log.debug(f'{fn} {pic.dtype} {pic.min()} {pic.max()}') mot.move_abs(p0) def cb_autofocus(self,*args,**kwargs): cam=self._cam mot=self._mot mtr=self._mtr msk=self._msk idx=self._idx p=mot.get_rbv() # the position is only polled at slow rate. # ŧherefore the received positions are sometimes the same # could try direct deltatau communications img16=np.array(cam._pic, np.int16) sb1=signal.convolve2d(img16, msk, mode='same', boundary='fill', fillvalue=0) sb2=signal.convolve2d(img16, msk.T, mode='same', boundary='fill', fillvalue=0) sb=np.abs(sb1)+np.abs(sb2) m=sb.sum() mtr[idx,:]=(p,m) _log.debug(f'{args} {kwargs} p:{p} mtr={m}') #print(p,cam._timestamp,m) if abs(p-self._pDst)<0.01 or idx>mtr.shape[0]: _log.debug(f'DONE->{self._idx}') del cam.process print(mtr[:idx,:]) mtr[:idx,0]=np.linspace(mtr[0,0],mtr[idx-1,0],idx) # smoothing the positions with equally spaced distances print(mtr[:idx,:]) mx=mtr[:,1].argmax() self._af=self.interpolate(mtr[:idx,0], mtr[:idx,1]) return self._idx+=1 #if p==mot.move_abs(p0+rng[1], wait=False) # move to sharp location # and unregister callback def run_continous(self,rng=(-1,1),velo=0.5): #do not stop and go, use ballback for processing mot=self._mot cam=self._cam p0=mot.get_rbv() v=mot._motor.get('VELO') p0=mot.get_rbv() mot._motor.put('VELO',2) mot.move_abs(p0+rng[0], wait=True) mot._motor.put('VELO',velo) n=int((rng[1]-rng[0])/velo)*20 # assume max 20fps self._mtr=mtr=np.ndarray(shape=(n,2)) self._idx=0 self._pDst=p0+rng[1] self._msk=np.array(((1, 0, -1), (2, 0, -2), (1, 0, -1)), np.int16) cam.process=self.cb_autofocus mot.move_abs(p0+rng[1], wait=False) import pyqtgraph as pg import time with pg.ProgressDialog('Progress', 0, n) as dlg: while True: try: p=self._af except AttributeError: time.sleep(.1) dlg.setValue(self._idx) continue break _log.debug(f' Final AF:{p}') mot._motor.put('VELO',2) mot.move_abs(p, wait=False) @staticmethod def interpolate(posLst,mtr): mx=mtr[:].argmax() _log.debug(f'best focus at idx:{mx}= pos:{posLst[mx]} = metric:{mtr[mx]:.6g}') if mx>1 and mx+2<=mtr.shape[0]: #fit parabola and interpolate: # y=ax2+bx+c, at positions x=-1, 0, 1, y'= 2a+b == 0 (top of parabola) # calc a,b,c: # y(-1)=a-b+c # y( 0)= +c # y( 1)=a+b+c # c=y(0) # b=(y(1)-y(-1))/2 # a=(y(1)+y(-1))/2-y(0) # x=-b/2a=(y(-1)-y(1))/2(y(-1)+y(1)-2y(0)) u,v,w=mtr[mx-1:mx+2] d=posLst[1]-posLst[0] p=posLst[mx]+d*.5*(u-w)/(u+w-2*v) else: p=posLst[mx] return p def run(self,rng=(-1,1),n=30,progressDlg=None,saveImg=False): # rng region (min max relative to current position) to seek # n number of images to take in region # progressDlg a pg.ProgressDialog to allow to abort the function call mot=self._mot cam=self._cam if mot is not None: p0=mot.get_rbv() else: p0=0 if type(cam) == list: imgLst=cam n=len(imgLst) mtr=np.ndarray(shape=(n,)) posLst=np.linspace(p0+rng[0], p0+rng[1], n) for i,p in enumerate(posLst): if type(cam)==list: img=PIL.Image.open(imgLst[i]) img=np.asarray(img) else: mot.move_abs(p, wait=True) img=cam.get_image() # get latest image img16=np.array(img, np.int16) msk=np.array(((1, 0, -1), (2, 0, -2), (1, 0, -1)), np.int16) sb1=signal.convolve2d(img16, msk, mode='same', boundary='fill', fillvalue=0) sb2=signal.convolve2d(img16, msk.T, mode='same', boundary='fill', fillvalue=0) sb=np.abs(sb1)+np.abs(sb2) mtr[i]=sb.sum() try: progressDlg+=1 if progressDlg.wasCanceled(): return except TypeError as e: pass _log.debug(f'{i}/{p:.4g} -> {mtr[i]:.4g}') p=autofocus.interpolate(posLst,mtr) if mot is not None: mot.move_abs(p) return p if __name__=="__main__": import argparse logging.basicConfig(level=logging.DEBUG, format='%(levelname)s:%(module)s:%(lineno)d:%(funcName)s:%(message)s ') #(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__)) logging.getLogger('matplotlib').setLevel(logging.INFO) import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import json class MyJsonEncoder(json.JSONEncoder): """ Special json encoder for numpy types """ def default(self, obj): if isinstance(obj, np.integer): return int(obj) elif isinstance(obj, np.floating): return float(obj) elif isinstance(obj, np.ndarray): return obj.tolist() elif type(obj) not in (dict, list, str, int): _log.error('dont know how to json') return repr(obj) return json.JSONEncoder.default(self, obj) np.set_printoptions(suppress=True) np.set_printoptions(linewidth=196) obj=geometry() 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)]} measure=\ {1000: [(-0.7310599999999997, 0.643119999999999, 1068.294887300626, 37.80427622442167), (-0.8105999999999997, 0.6437299218749977, 325.69398272322496, 66.36584947739863), (-0.8073999999999999, 0.548, 394.71778475125257, 954.1547514240993), (-0.7257699218749997, 0.5511799999999984, 1155.1696726117643, 900.6018015747674)], 800: [(-0.7087599999999997, 0.6409199999999983, 1056.3942317785522, 65.17578392519127), (-0.8374999999999998, 0.6427399999999998, 354.2555559762018, 79.45657055167973), (-0.8393599999999997, 0.4836299999999992, 378.0568670203493, 955.3448169763066), (-0.6977599999999998, 0.487619999999999, 1159.9299348205936, 906.5521293358044)], 600: [(-0.6560399999999997, 0.6448999999999997, 1112.3273127322987, 61.60558726856914), (-0.8914799999999998, 0.6495299218750006, 324.5039171710175, 67.55591502960601), (-0.8760799999999997, 0.3893999999999978, 407.80850582553353, 930.3534403799518), (-0.6621199999999997, 0.39789999999999964, 1128.988230463202, 873.2302938739978)], 400: [(-0.6044999999999998, 0.6441200000000008, 1074.2452150616627, 63.98571837298388), (-0.9739299218749998, 0.6475599999999995, 310.2231305445289, 73.50624279064286), (-0.9546199999999997, 0.22334999999999855, 398.28798140787467, 926.7832437233297), (-0.5859199999999999, 0.23333999999999833, 1144.4590826418978, 880.3706871872422)], 200: [(-0.4556199999999997, 0.6560399999999991, 1119.4677060455429, 44.9446695376659), (-1.0542399999999998, 0.6175399999999991, 381.6270636769714, 108.01814380465672), (-1.0568399999999998, -0.05702000000000044, 413.7588335865705, 932.7335714843665), (-0.47519999999999984, -0.03876000000000022, 1136.1286237764461, 881.5607527394494)], 1: [(-0.3154599999999996, 0.6150199999999987, 1112.3273127322987, 78.26650499947237), (-1.1583899218749996, 0.5868999999999996, 425.65948910864427, 129.43932374438944), (-1.2067399218749997, -0.3514000000000015, 400.6681125122893, 905.3620637835969), (-0.3347399999999998, -0.36506000000000133, 1136.1286237764461, 880.3706871872422)]} obj.update_pix2pos(measure) print(json.dumps(obj._lut_pix2pos, cls=MyJsonEncoder)) obj.interp_zoom(1) print(obj._pix2pos) print(obj._pix2pos.I) 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)]} opt_ctr_meas=\ {1000: [(943.5814614822486, 271.79551650350834), (310.28543090843203, 317.87448013362246), (357.2157641345181, 945.3504169712918), (982.8856744171983, 907.5691713113356)], # 800: [(804.3925790602868, 377.45059798190084), # (428.94991325159833, 401.78484483987137), # (456.760481089279, 776.3584304036324), # (829.595906163185, 752.0241835456618)], # 600: [(725.3062767718826, 435.6789743920446), # (497.52754628726467, 448.43934430925407), # (514.6798917766387, 674.3986678435945), # (740.7483129345815, 659.5001702939031)], 400: [(676.5616510826158, 467.6400729767663), (539.6440684653923, 477.20915075936176), (550.068542225442, 612.9457099265683), (686.4200153677604, 603.8457508612842)], # 200: [(647.9882475277134, 488.31261219978296), # (564.236743827639, 494.4770028644908), # (571.9870755133423, 577.6833939808159), # (654.4623032630916, 571.5416217015793)], # 1: [(633.8759316189822, 496.84436992197425), # (578.3851177374523, 500.9472972316709), # (582.6205681144841, 559.3272486023482), # (638.3177769017573, 554.8006203263507)] } obj.update_optical_center(opt_ctr_meas, True) if args.mode&0x04: pts=np.array([[ -633.75367631, -561.87640769], [ -636.87852469, -258.90639846], [-1098.22395446, -351.56498398], [-1096.62487933, -694.59151602]]) pts=np.array([[-699.81571798, -700.30880259], [-700.33262571, -500.3524493], [-1000.63771992, -501.08112667], [-999.69062995, -701.32290775]]) pts=np.array([[10,15], [40, 35], [10, 35], [40, 15]]) #pts=-pts for p in pts: print(p) trf=geometry.least_square_trf(pts) for p in ((0,0),(0,1),(1,0),(1,1)): fit=(trf*np.matrix((p[0],p[1],1)).T).A.ravel() print(p,fit) if args.mode&0x08: pts=np.array([[3.95620203, 3.17424935, 3.1415], [3.92067666, 2.43961212, 3.1415], [2.80894834, 2.71536886, 3.1415], [2.84446241, 3.54734879, 3.1415]]) plane=geometry.least_square_plane(pts) if args.mode&0x08: import glob imgLst=sorted(glob.glob("scratch/autofocus2/image*.png")) af=autofocus(imgLst,None) af.run() if args.mode&0x10: geometry.find_fiducial("scratch/fiducial/image001.png") # pix2pos="[[1.0, 200.0, 400.0, 600.0, 800.0], [[[0.001182928623952055, -2.6941995127711305e-05], [-4.043716694634124e-05, -0.0011894314084263825]], [[0.0007955995220142541, -3.175003727901119e-05], [-2.0896601103372113e-05, -0.0008100805094631365]], [[0.00048302539335378367, -1.1661121407652543e-05], [-2.0673746995751222e-05, -0.0004950857899461772]], [[0.00028775903460576395, -1.3762555219494508e-05], [-9.319936861519268e-06, -0.0002889214488565999]], [[0.0001788819256630411, -6.470841493681516e-06], [-2.0336605088889967e-06, -0.0001831131753499113]]]]" #lut_pix2pos=(np.array([1., 200., 400., 600., 800., 1000.]), # np.array([[[ 2.42827273e-03, -9.22117396e-05],[-1.10489804e-04, -2.42592492e-03]], # [[ 1.64346103e-03, -7.52341417e-05],[-6.60711165e-05, -1.64190224e-03]], # [[ 9.91307639e-04, -4.02008751e-05],[-4.23878232e-05, -9.91563507e-04]], # [[ 5.98443038e-04, -2.54046255e-05],[-2.76831563e-05, -6.02738142e-04]], # [[ 3.64418977e-04, -1.41389267e-05],[-1.55708176e-05, -3.66233567e-04]], # [[ 2.16526433e-04, -8.23070130e-06],[-9.29894004e-06, -2.16842976e-04]]])) #np.array([[[ 0.0010, 0.0001],[-0.0001, -0.0010]], # [[ 0.00163, 0.00000],[-0.00000, -0.00163]], # [[ 0.00099, 0.00000],[-0.00000, -0.00099]], # [[ 0.00060, 0.00000],[-0.00000, -0.00060]], # [[ 0.00036, 0.00000],[-0.00000, -0.00036]], # [[ 0.00022, 0.00000],[-0.00000, -0.00022]]]))