00001
00002
00003 """:py:class:`GeometryAccess` - holds and access hierarchical geometry for generic pixel detector
00004
00005 Usage::
00006
00007 from PSCalib.GeometryAccess import GeometryAccess, img_from_pixel_arrays
00008
00009 fname_geometry = '/reg/d/psdm/CXI/cxitut13/calib/CsPad::CalibV1/CxiDs1.0:Cspad.0/geometry/0-end.data'
00010 geometry = GeometryAccess(fname_geometry, 0377)
00011
00012 # load constants from geometry file
00013 geometry.load_pars_from_file(path=None)
00014
00015 # load constants from '\n'-separated str / text object
00016 geometry.load_pars_from_str(s)
00017
00018 # check if geometry info is available, returns bool
00019 status = geometry.is_valid()
00020
00021 # get pixel coordinate [um] arrays
00022 X, Y, Z = geometry.get_pixel_coords(oname=None, oindex=0, do_tilt=True)
00023
00024 # get pixel x,y coordinate [um] arrays projected toward origin on specified zplane, if zplane=None then zplane=Z.mean()
00025 X, Y = geometry.get_pixel_xy_at_z(zplane=None, oname=None, oindex=0, do_tilt=True)
00026
00027 # print a portion of pixel X, Y, and Z coordinate arrays
00028 geometry.print_pixel_coords(oname=None, oindex=0)
00029
00030 # get pixel area array; A=1 for regular pixels, =2.5 for wide.
00031 area = geometry.get_pixel_areas(oname=None, oindex=0)
00032
00033 # returns (smallest) pixel size [um]
00034 pixel_size = geometry.get_pixel_scale_size(oname=None, oindex=0)
00035
00036 # returns dictionary of comments associated with geometry (file)
00037 dict_of_comments = geometry.get_dict_of_comments()
00038
00039 # print comments associated with geometry (file)
00040 geometry.print_comments_from_dict()
00041
00042 # print list of geometry objects
00043 geometry.print_list_of_geos()
00044
00045 # print list of geometry object children
00046 geometry.print_list_of_geos_children()
00047
00048 # get pixel mask array;
00049 # mbits = +1-mask edges, +2-wide pixels, +4-non-bonded pixels, +8/+16 - four/eight neighbours of non-bonded
00050 mask = geometry.get_pixel_mask(oname=None, oindex=0, mbits=0377)
00051
00052 # get index arrays for entire detector
00053 iX, iY = geometry.get_pixel_coord_indexes(do_tilt=True)
00054
00055 # get index arrays for specified quad with offset
00056 iX, iY = geometry.get_pixel_coord_indexes('QUAD:V1', 1, pix_scale_size_um=None, xy0_off_pix=(1000,1000), do_tilt=True)
00057
00058 # get index arrays for pixel coordinates projected toward origin on specified zplane
00059 iX, iY = geometry.get_pixel_xy_inds_at_z(zplane=None, oname=None, oindex=0, pix_scale_size_um=None, xy0_off_pix=None, do_tilt=True)
00060
00061 # get ix and iy indexes for specified point in [um]. By default p_um=(0,0) - detector origin coordinates (center).
00062 ix, iy = geometry.point_coord_indexes(p_um=(0,0))
00063 # all other parameters should be the same as in get_pixel_coord_indexes method
00064 ix, iy = geometry.point_coord_indexes(p_um=(0,0), 'QUAD:V1', 1, pix_scale_size_um=None, xy0_off_pix=(1000,1000), do_tilt=True)
00065
00066 # get 2-d image from index arrays
00067 img = img_from_pixel_arrays(iX,iY,W=arr)
00068
00069 # Get specified object of the class GeometryObject, all objects are kept in the list self.list_of_geos
00070 geo = geometry.get_geo('QUAD:V1', 1)
00071 # Get top GeometryObject - the object which includes all other geometry objects
00072 geo = geometry.get_top_geo()
00073
00074 # modify currect geometry objects' parameters
00075 geometry.set_geo_pars('QUAD:V1', 1, x0, y0, z0, rot_z, rot_y, rot_x, tilt_z, tilt_y, tilt_x)
00076 geometry.move_geo('QUAD:V1', 1, 10, 20, 0)
00077 geometry.tilt_geo('QUAD:V1', 1, 0.01, 0, 0)
00078
00079 # save current geometry parameters in file
00080 geometry.save_pars_in_file(fname_geometry_new)
00081
00082 # change verbosity bit-control word; to print everythisg use pbits = 0xffff
00083 geometry.set_print_bits(pbits=0)
00084
00085 # return geometry parameters in "psf" format as a tuple psf[32][3][3]
00086 psf = geometry.get_psf()
00087 geometry.print_psf()
00088
00089 @see :py:class:`PSCalib.GeometryObject`, :py:class:`PSCalib.SegGeometry`, :py:class:`PSCalib.SegGeometryCspad2x1V1`, :py:class:`PSCalib.SegGeometryStore`
00090
00091 For more detail see `Detector Geometry <https://confluence.slac.stanford.edu/display/PSDM/Detector+Geometry>`_.
00092
00093 This software was developed for the SIT project. If you use all or
00094 part of it, please give an appropriate acknowledgment.
00095
00096 $Revision: 12898 $
00097
00098 @version $Id: GeometryAccess.py 12898 2016-11-22 19:27:56Z dubrovin@SLAC.STANFORD.EDU $
00099
00100 @author Mikhail S. Dubrovin
00101 """
00102
00103 __version__ = "$Revision: 12898 $"
00104
00105
00106 import os
00107 import sys
00108 from PSCalib.GeometryObject import GeometryObject
00109
00110 import numpy as np
00111 from math import floor, fabs
00112
00113
00114
00115 def divide_protected(num, den, vsub_zero=0) :
00116 """Returns result of devision of numpy arrays num/den with substitution of value vsub_zero for zero den elements.
00117 """
00118 pro_num = np.select((den!=0,), (num,), default=vsub_zero)
00119 pro_den = np.select((den!=0,), (den,), default=1)
00120 return pro_num / pro_den
00121
00122
00123
00124 class GeometryAccess :
00125 """ :py:class:`GeometryAccess`
00126 """
00127
00128 def __init__(self, path=None, pbits=0) :
00129 """Constructor of the class :py:class:`GeometryAccess`
00130 """
00131 self.path = path
00132 self.pbits = pbits
00133 self.valid = False
00134
00135 if path is None or not os.path.exists(path) :
00136 if pbits : print '%s: geometry file "%s" does not exist' % (self.__class__.__name__, path)
00137 return
00138
00139 self.load_pars_from_file()
00140
00141 if self.pbits & 1 : self.print_list_of_geos()
00142 if self.pbits & 2 : self.print_list_of_geos_children()
00143 if self.pbits & 4 : self.print_comments_from_dict()
00144
00145
00146
00147 def reset_cash(self) :
00148
00149 self.geo_old = None
00150 self.oname_old = None
00151 self.oindex_old = None
00152 self.tilt_old = None
00153 self.X_old = None
00154 self.Y_old = None
00155 self.Z_old = None
00156 self.iX_old = None
00157 self.iY_old = None
00158 self.ipx_old = None
00159 self.ipy_old = None
00160 self.p_um_old = None
00161
00162
00163
00164 def is_valid(self) :
00165 """returns True if geometry is loaded and presumably valid, otherwise False
00166 """
00167 return self.valid
00168
00169
00170
00171 def load_pars_from_file(self, path=None) :
00172 """Reads input "geometry" file, discards empty lines and comments, fills the list of geometry objects for data lines
00173 """
00174 self.valid = False
00175 if path is not None : self.path = path
00176
00177 self.reset_cash()
00178 self.dict_of_comments = {}
00179 self.list_of_geos = []
00180
00181 if self.pbits & 32 : print 'Load file: %s' % self.path
00182
00183 f=open(self.path,'r')
00184 for linef in f :
00185 line = linef.strip('\n')
00186 if self.pbits & 128 : print line
00187 if not line : continue
00188 if line[0] == '#' :
00189 self._add_comment_to_dict(line)
00190 continue
00191
00192 self.list_of_geos.append(self._parse_line(line))
00193
00194 f.close()
00195
00196 self._set_relations()
00197 self.valid = True
00198
00199
00200
00201 def load_pars_from_str(self, s) :
00202 """Reads input geometry from str, discards empty lines and comments, fills the list of geometry objects for data lines
00203 """
00204 self.valid = False
00205 if not isinstance(s, str) :
00206 if pbits : print '%s.load_pars_from_str input parameter is not a str, s: %s' % (self.__class__.__name__, str(s))
00207 return
00208
00209 self.reset_cash()
00210 self.dict_of_comments = {}
00211 self.list_of_geos = []
00212
00213 if self.pbits & 32 : print 'Load text: %s' % s
00214
00215 for linef in s.split('\n') :
00216 line = linef.strip('\n')
00217 if self.pbits & 128 : print line
00218 if not line : continue
00219 if line[0] == '#' :
00220 self._add_comment_to_dict(line)
00221 continue
00222
00223 self.list_of_geos.append(self._parse_line(line))
00224
00225 self._set_relations()
00226 self.valid = True
00227
00228
00229
00230 def save_pars_in_file(self, path) :
00231 """Save geometry file with current content
00232 """
00233 if not self.valid : return
00234
00235 if self.pbits & 32 : print 'Save file: %s' % path
00236
00237 txt = ''
00238
00239 for k in sorted(self.dict_of_comments) :
00240
00241 txt += '# %s\n' % (self.dict_of_comments[k])
00242
00243 txt += '\n'
00244
00245
00246 for geo in self.list_of_geos :
00247 if geo.get_parent_name() is None : continue
00248 txt += '%s\n' % (geo.str_data())
00249
00250 f=open(path,'w')
00251 f.write(txt)
00252 f.close()
00253
00254 if self.pbits & 64 : print txt
00255
00256
00257
00258 def _add_comment_to_dict(self, line) :
00259 """Splits the line of comments for keyward and value and store it in the dictionary
00260 """
00261
00262 cmt = line.lstrip('#').lstrip(' ')
00263 if len(cmt)<1 : return
00264 ind = len(self.dict_of_comments)
00265 if len(cmt)==1 :
00266
00267 self.dict_of_comments[ind] = ''
00268 return
00269
00270
00271
00272
00273
00274
00275 self.dict_of_comments[ind] = cmt.strip()
00276
00277
00278
00279 def _parse_line(self, line) :
00280 """Gets the string line with data from input file,
00281 creates and returns the geometry object for this string.
00282 """
00283 keys = ['pname','pindex','oname','oindex','x0','y0','z0','rot_z','rot_y','rot_x','tilt_z','tilt_y','tilt_x']
00284 f = line.split()
00285 if len(f) != len(keys) :
00286 print 'The list length for fields from file: %d is not equal to expected: %d' % (len(f), len(keys))
00287 return
00288
00289 vals = [str (f[0]),
00290 int (f[1]),
00291 str (f[2]),
00292 int (f[3]),
00293 float(f[4]),
00294 float(f[5]),
00295 float(f[6]),
00296 float(f[7]),
00297 float(f[8]),
00298 float(f[9]),
00299 float(f[10]),
00300 float(f[11]),
00301 float(f[12])
00302 ]
00303
00304
00305
00306
00307 d = dict(zip(keys, vals))
00308
00309
00310 return GeometryObject(**d)
00311
00312
00313
00314 def _find_parent(self, geobj) :
00315 """Finds and returns parent for geobj geometry object
00316 """
00317 for geo in self.list_of_geos :
00318 if geo == geobj : continue
00319 if geo.oindex == geobj.pindex \
00320 and geo.oname == geobj.pname :
00321 return geo
00322
00323
00324
00325 if geobj.pname is not None :
00326 top_parent = GeometryObject(pname=None, pindex=0, oname=geobj.pname, oindex=geobj.pindex)
00327 self.list_of_geos.append(top_parent)
00328 return top_parent
00329
00330 return None
00331
00332
00333
00334 def _set_relations(self) :
00335 """Set relations between geometry objects in the list_of_geos
00336 """
00337 for geo in self.list_of_geos :
00338
00339 parent = self._find_parent(geo)
00340
00341 if parent is None : continue
00342
00343 geo.set_parent(parent)
00344 parent.add_child(geo)
00345
00346 if self.pbits & 16 : print 'geo:%s:%d has parent:%s:%d' % (geo.oname, geo.oindex, parent.oname, parent.oindex)
00347
00348
00349
00350 def get_geo(self, oname, oindex) :
00351 """Returns specified geometry object
00352 """
00353 if not self.valid : return None
00354
00355 if oindex == self.oindex_old and oname == self.oname_old : return self.geo_old
00356
00357 for geo in self.list_of_geos :
00358 if geo.oindex == oindex \
00359 and geo.oname == oname :
00360 self.oindex_old = oindex
00361 self.oname_old = oname
00362 self.geo_old = geo
00363 return geo
00364 return None
00365
00366
00367
00368 def get_top_geo(self) :
00369 """Returns top geometry object
00370 """
00371 if not self.valid : return None
00372 return self.list_of_geos[-1]
00373
00374
00375
00376 def get_pixel_coords(self, oname=None, oindex=0, do_tilt=True) :
00377 """Returns three pixel X,Y,Z coordinate arrays for top or specified geometry object
00378 """
00379 if not self.valid : return None
00380
00381 if oindex == self.oindex_old\
00382 and oname == self.oname_old\
00383 and do_tilt == self.tilt_old :
00384 return self.X_old, self.Y_old, self.Z_old
00385
00386 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00387 if self.pbits & 8 :
00388 print 'get_pixel_coords(...) for geo:',
00389 geo.print_geo_children();
00390
00391 self.X_old, self.Y_old, self.Z_old = geo.get_pixel_coords(do_tilt)
00392 self.tilt = do_tilt
00393
00394 return self.X_old, self.Y_old, self.Z_old
00395
00396
00397
00398 def get_pixel_xy_at_z(self, zplane=None, oname=None, oindex=0, do_tilt=True) :
00399 """Returns pixel coordinate arrays XatZ, YatZ, for specified zplane and geometry object.
00400
00401 This method projects pixel X, Y coordinates in 3-D
00402 on the specified Z plane along direction to origin.
00403 """
00404 if not self.valid : return None, None
00405
00406 X, Y, Z = self.get_pixel_coords(oname, oindex, do_tilt)
00407 if X is None : return None, None
00408 Z0 = Z.mean() if zplane is None else zplane
00409 if fabs(Z0) < 1000 : return X, Y
00410
00411 XatZ = Z0 * divide_protected(X,Z)
00412 YatZ = Z0 * divide_protected(Y,Z)
00413 return XatZ, YatZ
00414
00415
00416
00417 def get_pixel_areas(self, oname=None, oindex=0) :
00418 """Returns pixel areas array for top or specified geometry object
00419 """
00420 if not self.valid : return None
00421 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00422 return geo.get_pixel_areas()
00423
00424
00425
00426 def get_pixel_mask(self, oname=None, oindex=0, mbits=0377) :
00427 """Returns pixel mask array for top or specified geometry object.
00428 mbits =+1 - mask edges
00429 +2 - two wide-pixel central columns
00430 +4 - non-bonded pixels
00431 +8 - four nearest neighbours of non-bonded pixels
00432 +16- eight neighbours of non-bonded pixels
00433 """
00434 if not self.valid : return None
00435 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00436 return geo.get_pixel_mask(mbits)
00437
00438
00439
00440 def get_pixel_scale_size(self, oname=None, oindex=0) :
00441 """Returns pixel scale size for top or specified geometry object
00442 """
00443 if not self.valid : return None
00444 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00445 return geo.get_pixel_scale_size()
00446
00447
00448
00449 def get_dict_of_comments(self) :
00450 """Returns dictionary of comments
00451 """
00452 if not self.valid : return None
00453 return self.dict_of_comments
00454
00455
00456
00457 def set_geo_pars(self, oname=None, oindex=0, x0=0, y0=0, z0=0, rot_z=0, rot_y=0, rot_x=0, tilt_z=0, tilt_y=0, tilt_x=0) :
00458 """Sets geometry parameters for specified or top geometry object
00459 """
00460 if not self.valid : return None
00461 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00462 return geo.set_geo_pars(x0, y0, z0, rot_z, rot_y, rot_x, tilt_z, tilt_y, tilt_x)
00463
00464
00465
00466 def move_geo(self, oname=None, oindex=0, dx=0, dy=0, dz=0) :
00467 """Moves specified or top geometry object by dx, dy, dz
00468 """
00469 if not self.valid : return None
00470 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00471 return geo.move_geo(dx, dy, dz)
00472
00473
00474
00475 def tilt_geo(self, oname=None, oindex=0, dt_x=0, dt_y=0, dt_z=0) :
00476 """Tilts specified or top geometry object by dt_x, dt_y, dt_z
00477 """
00478 if not self.valid : return None
00479 geo = self.get_top_geo() if oname is None else self.get_geo(oname, oindex)
00480 return geo.tilt_geo(dt_x, dt_y, dt_z)
00481
00482
00483
00484 def print_list_of_geos(self) :
00485 ss = '\nprint_list_of_geos():'
00486 if len(self.list_of_geos) == 0 : print '%s List_of_geos is empty...' % ss
00487 if not self.valid : return
00488 for geo in self.list_of_geos : geo.print_geo()
00489
00490
00491
00492 def print_list_of_geos_children(self) :
00493 ss = '\nprint_list_of_geos_children():'
00494 if len(self.list_of_geos) == 0 : print '%s List_of_geos is empty...' % ss
00495 if not self.valid : return
00496 for geo in self.list_of_geos : geo.print_geo_children()
00497
00498
00499
00500 def print_comments_from_dict(self) :
00501 print '\nprint_comments_from_dict():'
00502 if not self.valid : return
00503
00504 for k in sorted(self.dict_of_comments):
00505 print 'key: %3d val: %s' % (k, self.dict_of_comments[k])
00506
00507
00508
00509 def print_pixel_coords(self, oname=None, oindex=0) :
00510 """Partial print of pixel coordinate X,Y,Z arrays for selected or top(by default) geo
00511 """
00512 if not self.valid : return
00513 X, Y, Z = self.get_pixel_coords(oname, oindex, do_tilt=True)
00514
00515 print 'size=', X.size
00516 print 'X: %s...'% ', '.join(['%10.1f'%v for v in X.flatten()[0:9]])
00517 print 'Y: %s...'% ', '.join(['%10.1f'%v for v in Y.flatten()[0:9]])
00518 print 'Z: %s...'% ', '.join(['%10.1f'%v for v in Z.flatten()[0:9]])
00519
00520
00521
00522 def get_pixel_coord_indexes(self, oname=None, oindex=0, pix_scale_size_um=None, xy0_off_pix=None, do_tilt=True) :
00523 """Returns two pixel X,Y coordinate index arrays for top or specified geometry object
00524 """
00525 if not self.valid : return None, None
00526
00527 if oindex == self.oindex_old\
00528 and oname == self.oname_old\
00529 and do_tilt == self.tilt_old\
00530 and pix_scale_size_um is None\
00531 and xy0_off_pix is None\
00532 and self.iX_old is not None:
00533 return self.iX_old, self.iY_old
00534
00535 X, Y, Z = self.get_pixel_coords(oname, oindex, do_tilt)
00536
00537 pix_size = self.get_pixel_scale_size() if pix_scale_size_um is None else pix_scale_size_um
00538 pix_half = 0.5*pix_size
00539
00540 xmin, ymin = X.min()-pix_half, Y.min()-pix_half
00541
00542 if xy0_off_pix is not None :
00543
00544 x_off_um = xy0_off_pix[0] * pix_size
00545 y_off_um = xy0_off_pix[1] * pix_size
00546
00547 xmin += x_off_um
00548 ymin += y_off_um
00549 x_off_um = x_off_um + pix_half if xmin>0 else x_off_um - xmin
00550 y_off_um = y_off_um + pix_half if ymin>0 else y_off_um - ymin
00551 self.iX_old, self.iY_old = np.array((X+x_off_um)/pix_size, dtype=np.uint), np.array((Y+y_off_um)/pix_size, dtype=np.uint)
00552
00553 else :
00554 self.iX_old, self.iY_old = np.array((X-xmin)/pix_size, dtype=np.uint), np.array((Y-ymin)/pix_size, dtype=np.uint)
00555
00556 return self.iX_old, self.iY_old
00557
00558
00559
00560 def get_pixel_xy_inds_at_z(self, zplane=None, oname=None, oindex=0, pix_scale_size_um=None, xy0_off_pix=None, do_tilt=True) :
00561 """Returns pixel coordinate index arrays iX, iY of size for specified zplane and geometry object
00562 """
00563 if not self.valid : return None, None
00564
00565 X, Y = self.get_pixel_xy_at_z(zplane, oname, oindex, do_tilt)
00566
00567 if X is None : return None, None
00568
00569 pix_size = self.get_pixel_scale_size() if pix_scale_size_um is None else pix_scale_size_um
00570 pix_half = 0.5*pix_size
00571
00572 xmin, ymin = X.min()-pix_half, Y.min()-pix_half
00573
00574 if xy0_off_pix is not None :
00575
00576 x_off_um = xy0_off_pix[0] * pix_size
00577 y_off_um = xy0_off_pix[1] * pix_size
00578
00579 xmin += x_off_um
00580 ymin += y_off_um
00581 x_off_um = x_off_um + pix_half if xmin>0 else x_off_um - xmin
00582 y_off_um = y_off_um + pix_half if ymin>0 else y_off_um - ymin
00583 self.iX_old, self.iY_old = np.array((X+x_off_um)/pix_size, dtype=np.uint), np.array((Y+y_off_um)/pix_size, dtype=np.uint)
00584
00585 else :
00586 self.iX_old, self.iY_old = np.array((X-xmin)/pix_size, dtype=np.uint), np.array((Y-ymin)/pix_size, dtype=np.uint)
00587
00588 return self.iX_old, self.iY_old
00589
00590
00591
00592 def point_coord_indexes(self, p_um=(0,0), oname=None, oindex=0, pix_scale_size_um=None, xy0_off_pix=None, do_tilt=True) :
00593 """Converts point (x_um, y_um) corrdinates [um] to pixel (ix, iy) indexes.
00594 All other parameters are the same as in get_pixel_coord_indexes.
00595 WARNING: indexes are not required to be inside the image. They are integer, may be negative or exceed pixel maximal index.
00596 """
00597 if not self.valid : return None, None
00598
00599 if oindex == self.oindex_old\
00600 and oname == self.oname_old\
00601 and do_tilt == self.tilt_old\
00602 and p_um == self.p_um_old\
00603 and pix_scale_size_um is None\
00604 and xy0_off_pix is None\
00605 and self.ipx_old is not None:
00606 return self.ipx_old, self.ipy_old
00607
00608 X, Y, Z = self.get_pixel_coords(oname, oindex, do_tilt)
00609
00610 pix_size = self.get_pixel_scale_size() if pix_scale_size_um is None else pix_scale_size_um
00611 pix_half = 0.5*pix_size
00612
00613 xmin, ymin = X.min()-pix_half, Y.min()-pix_half
00614 x_um, y_um = self.p_um_old = p_um
00615
00616 if xy0_off_pix is not None :
00617
00618 x_off_um = xy0_off_pix[0] * pix_size
00619 y_off_um = xy0_off_pix[1] * pix_size
00620
00621 xmin += x_off_um
00622 ymin += y_off_um
00623 x_off_um = x_off_um + pix_half if xmin>0 else x_off_um - xmin
00624 y_off_um = y_off_um + pix_half if ymin>0 else y_off_um - ymin
00625
00626 self.ipx_old, self.ipy_old = int(floor((x_um+x_off_um)/pix_size)), int(floor((y_um+y_off_um)/pix_size))
00627
00628 else :
00629 self.ipx_old, self.ipy_old = int(floor((x_um-xmin)/pix_size)), int(floor((y_um-ymin)/pix_size))
00630
00631 return self.ipx_old, self.ipy_old
00632
00633
00634
00635 def set_print_bits(self, pbits=0) :
00636 """ Sets printout control bitword
00637 """
00638 self.pbits = pbits
00639
00640
00641
00642 def get_psf(self) :
00643 """Returns array of vectors in TJ format (psf stands for position-slow-fast vectors)
00644 """
00645 if not self.valid : return None
00646 X, Y, Z = self.get_pixel_coords()
00647 if X.size != 32*185*388 : return None
00648
00649 shape_cspad = (32,185,388)
00650 X.shape, Y.shape, Z.shape, = shape_cspad, shape_cspad, shape_cspad
00651
00652 psf = []
00653
00654 for s in range(32) :
00655 vp = (X[s,0,0], Y[s,0,0], Z[s,0,0])
00656
00657 vs = (X[s,1,0]-X[s,0,0], \
00658 Y[s,1,0]-Y[s,0,0], \
00659 Z[s,1,0]-Z[s,0,0])
00660
00661 vf = (X[s,0,1]-X[s,0,0], \
00662 Y[s,0,1]-Y[s,0,0], \
00663 Z[s,0,1]-Z[s,0,0])
00664
00665 psf.append((vp,vs,vf))
00666
00667 return psf
00668
00669
00670
00671
00672 def print_psf(self) :
00673 """ Gets and prints psf array for test purpose
00674 """
00675 if not self.valid : return None
00676 psf = np.array(self.get_psf())
00677 print 'print_psf(): psf.shape: %s \npsf vectors:' % (str(psf.shape))
00678 for (px,py,pz), (sx,xy,xz), (fx,fy,fz) in psf:
00679 print ' p=(%12.2f, %12.2f, %12.2f), s=(%8.2f, %8.2f, %8.2f) f=(%8.2f, %8.2f, %8.2f)' \
00680 % (px,py,pz, sx,xy,xz, fx,fy,fz)
00681
00682
00683
00684
00685
00686 def img_default(shape=(10,10), dtype = np.float32) :
00687 """Returns default image
00688 """
00689 arr = np.arange(shape[0]*shape[1], dtype=dtype)
00690 arr.shape = shape
00691 return arr
00692
00693
00694
00695 def img_from_pixel_arrays(iX, iY, W=None, dtype=np.float32, vbase=0) :
00696 """Returns image from iX, iY coordinate index arrays and associated weights W.
00697 """
00698 if iX.size != iY.size \
00699 or (W is not None and iX.size != W.size) :
00700 msg = 'img_from_pixel_arrays(): WARNING input array sizes are different;' \
00701 + ' iX.size=%d, iY.size=%d, W.size=%d' % (iX.size, iY.size, W.size)
00702 print msg
00703 return img_default()
00704
00705 iXfl = iX.flatten()
00706 iYfl = iY.flatten()
00707
00708 xsize = iXfl.max()+1
00709 ysize = iYfl.max()+1
00710
00711 weight = W.flatten() if W is not None else np.ones_like(iXfl)
00712 img = vbase*np.ones((xsize,ysize), dtype=dtype)
00713 img[iXfl,iYfl] = weight
00714 return img
00715
00716
00717
00718
00719
00720
00721
00722 if __name__ == "__main__" :
00723 from time import time
00724
00725
00726 import pyimgalgos.GlobalGraphics as gg
00727 import pyimgalgos.TestImageGenerator as tig
00728
00729
00730
00731 def test_access(geometry) :
00732 """ Tests geometry acess methods of the class GeometryAccess
00733 """
00734 geometry.print_list_of_geos()
00735 geometry.print_list_of_geos_children()
00736
00737 print '\nTOP GEO:'
00738 top_geo = geometry.get_top_geo()
00739 top_geo.print_geo_children()
00740
00741 print '\nINTERMEDIATE GEO (QUAD):'
00742 geo = geometry.get_geo('QUAD:V1', 0)
00743
00744 geo.print_geo_children()
00745
00746 t0_sec = time()
00747 X,Y,Z = geo.get_pixel_coords(do_tilt=True)
00748
00749 print 'X:\n', X
00750 print 'Consumed time to get 3d pixel coordinates = %7.3f sec' % (time()-t0_sec)
00751 print 'Geometry object: %s:%d X.shape:%s' % (geo.oname, geo.oindex, str(X.shape))
00752
00753 print '\nTest of print_pixel_coords() for quad:'
00754 geometry.print_pixel_coords('QUAD:V1', 1)
00755 print '\nTest of print_pixel_coords() for CSPAD:'
00756 geometry.print_pixel_coords()
00757
00758 print '\nTest of get_pixel_areas() for QUAD:'
00759 A = geo.get_pixel_areas()
00760 print 'Geometry object: %s:%d A.shape:%s' % (geo.oname, geo.oindex, str(A.shape))
00761 print 'A[0,0:5,190:198]:\n', A[0,0:5,190:198]
00762
00763 print '\nTest of get_pixel_areas() for CSPAD:'
00764 A = top_geo.get_pixel_areas()
00765 print 'Geometry object: %s:%d A.shape:%s' % (geo.oname, geo.oindex, str(A.shape))
00766 print 'A[0,0,0:5,190:198]:\n', A[0,0,0:5,190:198]
00767
00768 print '\nTest of get_size_geo_array()'
00769 print 'for QUAD : %d' % geo.get_size_geo_array()
00770 print 'for CSPAD: %d' % top_geo.get_size_geo_array()
00771
00772 print '\nTest of get_pixel_scale_size()'
00773 print 'for QUAD : %8.2f' % geo.get_pixel_scale_size()
00774 print 'for CSPAD : %8.2f' % top_geo.get_pixel_scale_size()
00775 print 'for geometry: %8.2f' % geometry.get_pixel_scale_size()
00776
00777 print '\nTest of get_dict_of_comments():'
00778 d = geometry.get_dict_of_comments()
00779 print "d[0] = %s" % d[0]
00780
00781
00782
00783 def test_plot_quad(geometry) :
00784 """ Tests geometry acess methods of the class GeometryAccess object for CSPAD quad
00785 """
00786
00787 iX, iY = geometry.get_pixel_coord_indexes('QUAD:V1', 1, pix_scale_size_um=None, xy0_off_pix=None, do_tilt=True)
00788
00789
00790 arr = tig.cspad_nparr(n2x1=iX.shape[0])
00791 arr.shape = (8,185,388)
00792 amp_range = (0,185+388)
00793
00794 print 'iX, iY, W shape:', iX.shape, iY.shape, arr.shape
00795 img = img_from_pixel_arrays(iX,iY,W=arr)
00796
00797 gg.plotImageLarge(img,amp_range=amp_range)
00798 gg.move(500,10)
00799 gg.show()
00800
00801
00802
00803 def test_mask_quad(geometry, mbits) :
00804 """ Tests geometry acess methods of the class GeometryAccess object for CSPAD quad
00805 """
00806
00807 iX, iY = geometry.get_pixel_coord_indexes('QUAD:V1', 1, pix_scale_size_um=None, xy0_off_pix=None, do_tilt=True)
00808
00809
00810 arr = geometry.get_pixel_mask('QUAD:V1', 1, mbits)
00811 arr.shape = (8,185,388)
00812 amp_range = (-1,2)
00813
00814 print 'iX, iY, W shape:', iX.shape, iY.shape, arr.shape
00815 img = img_from_pixel_arrays(iX, iY, W=arr, vbase=0.5)
00816
00817 gg.plotImageLarge(img,amp_range=amp_range)
00818 gg.move(500,10)
00819 gg.show()
00820
00821
00822
00823 def test_plot_cspad(geometry, fname_data, amp_range=(0,0.5)) :
00824 """ The same test as previous, but use get_pixel_coord_indexes(...) method
00825 """
00826
00827
00828 rad1 = 655
00829 rad2 = 670
00830
00831
00832 xc, yc = 1000, 1000
00833 xyc = xc, yc
00834
00835
00836 iX, iY = geometry.get_pixel_coord_indexes(xy0_off_pix=xyc, do_tilt=True)
00837
00838 ixo, iyo = geometry.point_coord_indexes(xy0_off_pix=xyc, do_tilt=True)
00839 print 'Detector origin indexes ixo, iyo:', ixo, iyo
00840
00841 root, ext = os.path.splitext(fname_data)
00842 arr = np.load(fname_data) if ext == '.npy' else np.loadtxt(fname_data, dtype=np.float)
00843 arr.shape= (4,8,185,388)
00844
00845 print 'iX, iY, W shape:', iX.shape, iY.shape, arr.shape
00846
00847 arr.shape = iX.shape
00848 img = img_from_pixel_arrays(iX, iY, W=arr)
00849
00850 xyc_ring = (yc, xc)
00851 axim = gg.plotImageLarge(img,amp_range=amp_range)
00852 gg.drawCircle(axim, xyc_ring, rad1, linewidth=1, color='w', fill=False)
00853 gg.drawCircle(axim, xyc_ring, rad2, linewidth=1, color='w', fill=False)
00854 gg.drawCenter(axim, xyc_ring, rad1, linewidth=1, color='w')
00855 gg.move(500,10)
00856 gg.show()
00857
00858
00859
00860 def test_img_default() :
00861 """ Test default image
00862 """
00863 axim = gg.plotImageLarge( img_default() )
00864 gg.move(500,10)
00865 gg.show()
00866
00867
00868
00869 def test_save_pars_in_file(geometry) :
00870 """ Test default image
00871 """
00872
00873 x0, y0, z0, rot_z, rot_y, rot_x, tilt_z, tilt_y, tilt_x = -3500, 5800, 0, 0.123, 0.123, 0.123, 1, 2, 3
00874 geometry.set_geo_pars('QUAD:V1', 1, x0, y0, z0, rot_z, rot_y, rot_x, tilt_z, tilt_y, tilt_x)
00875
00876 geometry.set_print_bits(32)
00877 geometry.save_pars_in_file('./test.txt')
00878
00879
00880
00881 def test_load_pars_from_file(geometry) :
00882 """ Test default image
00883 """
00884 geometry.set_print_bits(32+64)
00885 geometry.load_pars_from_file('./test.txt')
00886 geometry.print_list_of_geos()
00887
00888
00889
00890 def test_cspad2x2() :
00891 """ Test cspad2x2 geometry table
00892 """
00893
00894 basedir = '/reg/neh/home1/dubrovin/LCLS/CSPad2x2Alignment/calib-cspad2x2-01-2013-02-13/'
00895 fname_geometry = basedir + 'calib/CsPad2x2::CalibV1/MecTargetChamber.0:Cspad2x2.1/geometry/0-end.data'
00896 fname_data = basedir + 'cspad2x2.1-ndarr-ave-meca6113-r0028.dat'
00897
00898
00899
00900
00901
00902
00903 geometry = GeometryAccess(fname_geometry, 0377)
00904 amp_range = (0,15000)
00905
00906
00907
00908
00909
00910 iX, iY = geometry.get_pixel_coord_indexes(do_tilt=True)
00911
00912 root, ext = os.path.splitext(fname_data)
00913 arr = np.load(fname_data) if ext == '.npy' else np.loadtxt(fname_data, dtype=np.float)
00914 arr.shape= (185,388,2)
00915
00916 print 'iX, iY, W shape:', iX.shape, iY.shape, arr.shape
00917 img = img_from_pixel_arrays(iX,iY,W=arr)
00918
00919 axim = gg.plotImageLarge(img,amp_range=amp_range)
00920 gg.move(500,10)
00921 gg.show()
00922
00923
00924
00925 def test_epix100a() :
00926 """ Test test_epix100a geometry table
00927 """
00928 basedir = '/reg/g/psdm/detector/alignment/cspad/calib-cxi-ds1-2014-05-15/'
00929 fname_geometry = basedir + 'calib/CsPad::CalibV1/CxiDs1.0:Cspad.0/geometry/2-end.data'
00930 fname_data = basedir + 'cspad-arr-cxid2714-r0023-lysozyme-rings.txt'
00931
00932
00933
00934
00935
00936 geometry = GeometryAccess(fname_geometry, 0177777)
00937 amp_range = (-4,10)
00938
00939 iX, iY = geometry.get_pixel_coord_indexes()
00940
00941 root, ext = os.path.splitext(fname_data)
00942 arr = np.load(fname_data) if ext == '.npy' else np.loadtxt(fname_data, dtype=np.float)
00943
00944 print 'iX, iY, W shape:', iX.shape, iY.shape, arr.shape
00945 img = img_from_pixel_arrays(iX,iY,W=arr)
00946
00947 axim = gg.plotImageLarge(img,amp_range=amp_range)
00948 gg.move(500,10)
00949 gg.show()
00950
00951
00952
00953 def test_cspad_xy_at_z() :
00954 """ Test cspad geometry table
00955 """
00956
00957 basedir = '/reg/g/psdm/detector/alignment/cspad/calib-cxi-camera1-2014-09-24/'
00958 fname_geometry = basedir + '2016-06-03-geometry-cxi06216-r25-camera1-z175mm.txt'
00959 fname_data = basedir + '2016-06-03-chun-cxi06216-0025-DscCsPad-max.txt'
00960
00961 geometry = GeometryAccess(fname_geometry, 0377)
00962
00963
00964 xyc = xc, yc = 1000, 1000
00965
00966
00967
00968 iX, iY = geometry.get_pixel_xy_inds_at_z(zplane=150000)
00969
00970 root, ext = os.path.splitext(fname_data)
00971 arr = np.load(fname_data) if ext == '.npy' else np.loadtxt(fname_data, dtype=np.float)
00972
00973
00974 arr.shape= (32,185,388)
00975
00976
00977
00978 amp_range = (0, 1000)
00979 print 'amp_range', amp_range
00980
00981 print 'iX, iY, W shape:', iX.shape, iY.shape, arr.shape
00982 img = img_from_pixel_arrays(iX,iY,W=arr)
00983
00984 axim = gg.plotImageLarge(img,amp_range=amp_range)
00985 gg.move(500,10)
00986 gg.show()
00987
00988
00989
00990
00991
00992
00993 if __name__ == "__main__" :
00994
00995
00996
00997
00998
00999
01000
01001
01002 basedir = '/reg/g/psdm/detector/alignment/cspad/calib-cxi-ds1-2014-03-19/'
01003 fname_data = basedir + 'cspad-ndarr-ave-cxii0114-r0227.dat'
01004 fname_geometry = basedir + 'calib/CsPad::CalibV1/CxiDs1.0:Cspad.0/geometry/0-end.data'
01005
01006 amp_range = (0,500)
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027 print '%s\nfname_geometry: %s\nfname_data: %s' %(120*'_', fname_geometry, fname_geometry)
01028
01029 geometry = GeometryAccess(fname_geometry, 0)
01030
01031 msg = 'Use command: sys.argv[0] <num>, wher num=1,2,3,...,10'
01032
01033 if len(sys.argv)==1 : print 'App needs in input parameter.' + msg
01034 elif sys.argv[1]=='1' : test_access(geometry)
01035 elif sys.argv[1]=='2' : test_plot_quad(geometry)
01036 elif sys.argv[1]=='3' : test_plot_cspad(geometry, fname_data, amp_range)
01037 elif sys.argv[1]=='4' : test_img_default()
01038 elif sys.argv[1]=='5' :
01039 print 'Init GeometryAccess is silent? (see below)'
01040 ga0 = GeometryAccess(fname_geometry, 0)
01041 elif sys.argv[1]=='6' : ga0377 = GeometryAccess(fname_geometry, 0377)
01042 elif sys.argv[1]=='7' : test_save_pars_in_file(geometry)
01043 elif sys.argv[1]=='8' : test_load_pars_from_file(geometry)
01044 elif sys.argv[1]=='9' : test_mask_quad(geometry, 1+2+8)
01045 elif sys.argv[1]=='10': geometry.print_psf()
01046 elif sys.argv[1]=='11': test_cspad2x2()
01047 elif sys.argv[1]=='12': test_epix100a()
01048 elif sys.argv[1]=='13': geometry.print_comments_from_dict()
01049 elif sys.argv[1]=='14': test_cspad_xy_at_z()
01050 else : print 'Wrong input parameter.' + msg
01051
01052 sys.exit ('End of %s' % sys.argv[0])
01053
01054
01055
01056