From 011eaa3e313f9d84ce7698ddb57d84028b1d1749 Mon Sep 17 00:00:00 2001 From: Thierry Zamofing Date: Mon, 29 Aug 2022 17:38:40 +0200 Subject: [PATCH] towards least_square_trf --- ModuleFixTarget.py | 50 ++++++++++++++++------------ geometry.py | 81 ++++++++++++++++++++++++++++++++++++++++++++++ swissmx.py | 33 +++++++++++++++++++ 3 files changed, 144 insertions(+), 20 deletions(-) diff --git a/ModuleFixTarget.py b/ModuleFixTarget.py index 63f363f..4d211df 100644 --- a/ModuleFixTarget.py +++ b/ModuleFixTarget.py @@ -175,7 +175,7 @@ class WndFixTarget(QWidget): self._cbType=cb=QComboBox() #cb.addItem("C") #cb.addItem("C++") - cb.addItems(["Fiducial", "FixTarget(12.5x12.5)", "FixTarget(23.0x23.0)", "FixTarget()", "Grid()"]) + #cb.addItems(["Fiducial", "FixTarget(12.5x12.5)", "FixTarget(23.0x23.0)", "FixTarget()", "Grid()"]) #cb.currentIndexChanged.connect(self.selectionchange) @@ -191,6 +191,9 @@ class WndFixTarget(QWidget): self._btnDelAll=btnDelAll = QPushButton("del all") #btnDelAll.clicked.connect(lambda x: _log.warning("TODO: IMPLEMENT") ) bl.addWidget(btnDelAll, 2, 2,1,1) + self._btnFit=btnFit = QPushButton("fit with fiducials") + #btnDelAll.clicked.connect(lambda x: _log.warning("TODO: IMPLEMENT") ) + bl.addWidget(btnFit, 2, 1,1,1) @@ -267,21 +270,25 @@ class WndFixTarget(QWidget): data_folder='' if filename is None: filename, _ = QFileDialog.getOpenFileName(self,"Load data file",None, - "json files (*.json);;yaml files (*.yaml);;pickle files (*.pkl);;text files (*.txt);;all files (*)",) + "json files (*.json);;all files (*)",) + #"json files (*.json);;yaml files (*.yaml);;pickle files (*.pkl);;text files (*.txt);;all files (*)",) if not filename: # cancelled dialog return ext=filename.rsplit('.',1)[1].lower() - if ext=='pkl': - with open(filename, 'rb') as f: - data=pickle.load(f) - if ext=='yaml': - with open(filename, 'r') as f: - data = yaml.load(f) - elif ext=='json': + + if ext=='json': with open(filename, 'r') as f: data=json.load(f,object_hook=MyJsonDecoder) - print(data) + else: + raise(IOError('unsupported file type')) + #elif ext=='yaml': + # with open(filename, 'r') as f: + # data = yaml.load(f) + #elif ext=='pkl': + # with open(filename, 'rb') as f: + # data=pickle.load(f) + self._data=data self._tree.setData(data) try: @@ -354,7 +361,8 @@ class WndFixTarget(QWidget): data_folder='' if filename is None: filename, _ = QFileDialog.getSaveFileName(self,"Save data file",data_folder, - "json files (*.json);;yaml files (*.yaml);;pickle files (*.pkl);;text files (*.txt);;all files (*)",) + "json files (*.json);;all files (*)",) + #"json files (*.json);;yaml files (*.yaml);;pickle files (*.pkl);;text files (*.txt);;all files (*)",) if not filename: return @@ -365,16 +373,18 @@ class WndFixTarget(QWidget): #df.to_csv(filename, float_format="%.6f") #import numpy as np ext=filename.rsplit('.',1)[1].lower() - if ext=='pkl': - with open(filename, 'wb') as f: - pickle.dump(self._data, f) - elif ext=='yaml': - with open(filename, 'w') as f: - yaml.dump(self._data, f) - elif ext=='json': + if ext=='json': with open(filename, 'w') as f: json.dump(self._data, f,cls=MyJsonEncoder, indent=2)#separators=(',', ':') - print(self._data) + else: + raise(IOError('unsupported file type')) + #elif ext=='yaml': + # with open(filename, 'w') as f: + # yaml.dump(self._data, f) + #elif ext=='pkl': + # with open(filename, 'wb') as f: + # pickle.dump(self._data, f) + #print(self._data) def delete_selected(self): row = self._current_row @@ -729,7 +739,7 @@ void itemSelectionChanged() if __name__ == "__main__": import sys - class Monster(yaml.YAMLObject): + class Monster: yaml_tag = u'!Monster' def __init__(self, name, hp, ac, attacks): self.name = name diff --git a/geometry.py b/geometry.py index 498edb9..3720043 100755 --- a/geometry.py +++ b/geometry.py @@ -267,6 +267,56 @@ class geometry: # returns the vector m(x,y) of the optical center to the xray pass + @staticmethod + def least_square_trf(points): + # 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.float) + + + y=points # sort points p00 p01 p10 p11 + s=y[:,1].argsort() + y=y[s,:] + s=y[:2,0].argsort() + y[:2,:]=y[:2,:][s,:] + s=y[2:,0].argsort() + y[2:,:]=y[2:,:][s,:] + 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 + + @staticmethod + def least_square_plane(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 + pass + if __name__=="__main__": import argparse @@ -408,4 +458,35 @@ if __name__=="__main__": 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) diff --git a/swissmx.py b/swissmx.py index 8959b1c..a907f8c 100755 --- a/swissmx.py +++ b/swissmx.py @@ -1541,8 +1541,12 @@ class SwissMxWnd(QMainWindow, Ui_MainWindow): self._moduleFixTarget =mft = ModuleFixTarget.WndFixTarget(self,data) tab = self._tabs_daq_methods.insertTab(0,mft,'Fix Target') + self._tabs_daq_methods.setCurrentWidget(mft) #set this as the active tabs + + mft._cbType.addItems(["Fiducial", "FixTarget(12.5x12.5)", "FixTarget(23.0x23.0)", "FixTarget()", "Grid()"]) mft._btnAdd.clicked.connect(self.module_fix_target_add_obj) mft._btnDelAll.clicked.connect(self.module_fix_target_del_all_obj) + mft._btnFit.clicked.connect(self.module_fix_target_fit_fiducial) #tab.layout().addWidget(mft) mft.prefixSelected.connect(lambda prefix: self._le_prefix.setText(prefix)) @@ -1659,6 +1663,35 @@ class SwissMxWnd(QMainWindow, Ui_MainWindow): objLst.clear() mft._tree.setData(objLst) + def module_fix_target_fit_fiducial(self): + mft=self._moduleFixTarget + vb=self.vb + objLst=self._goTracked['objLst'] + n=0 + for go in objLst: + if type(go)==UsrGO.Fiducial: + n+=1 + ptFitTrf=np.ndarray((4,2)) # 4 (x,y) points to fit a transformation of a parallelogram + ptFitPlane=np.ndarray((n,3)) #n (x,y,z) points to fit a plane transformation z=ax+by+c + i=j=0 + for go in reversed(objLst): + if type(go)==UsrGO.FixTargetFrame: + if j==4: + trf=geometry.geometry.least_square_trf(ptFitTrf) + tr=go.transform() + tr.setMatrix(100, 0, 0, + 0, 100, 0, + 0, 0, 1) + go.setTransform(tr) + go.setPos(trf[:,2]) + go.setSize((1,1)) + j=0 + elif type(go)==UsrGO.Fiducial: + ptFitTrf[j]=go.pos()+go.size()/2 + ptFitPlane[i]=go._xyz + i+=1;j+=1 + plane=geometry.geometry.least_square_plane(ptFitPlane) + # **************** OBSOLETE AND/OR OLD STUFF ****************