Commit 48ae962e authored by Thomas Vincent's avatar Thomas Vincent

move qspace conversion with xrayutilities in a separate function

parent 4138114f
......@@ -54,6 +54,73 @@ logger = logging.getLogger(__name__)
disp_times = False
def qspace_conversion(img_size, center_chan, chan_per_deg,
beam_energy, phi, eta, nu, delta):
"""Returns array of Q vectors corresponding to each pixel
This function is based on xrayutilities
:param List[int] img_size:
Size of the detector images in pixels (dim0, dim1)
:param List[float] center_chan:
Calibration center channel in image coordinates (dim0, dim1)
:param List[float] chan_per_deg:
Channel per degree calibration in image coordinates (dim0, dim1)
:param float beam_energy: Energy in eV
:param numpy.ndarray phi: 1D array of phi angle for each image
:param numpy.ndarray eta: 1D array of eta angle for each image
:param numpy.ndarray nu: 1D array of nu angle for each image
:param numpy.ndarray delta: 1D array of delta angle for each image
:return: Array Q vectors of shape: (nb images, image size, 3 (qx, qy, qz))
:rtype: numpy.ndarray
"""
phi = np.array(phi, copy=False, dtype=np.float64)
eta = np.array(eta, copy=False, dtype=np.float64)
nu = np.array(nu, copy=False, dtype=np.float64)
delta = np.array(delta, copy=False, dtype=np.float64)
qconv = xu.experiment.QConversion(['y-', 'z-'],
['z-', 'y-'],
[1, 0, 0])
# convention for coordinate system:
# x downstream
# z upwards
# y to the "outside"
# (righthanded)
hxrd = xu.HXRD([1, 0, 0],
[0, 0, 1],
en=beam_energy,
qconv=qconv)
hxrd.Ang2Q.init_area('z-',
'y+',
cch1=center_chan[0],
cch2=center_chan[1],
Nch1=img_size[0],
Nch2=img_size[1],
chpdeg1=chan_per_deg[0],
chpdeg2=chan_per_deg[1])
n_images = len(eta)
# shape of the array that will store the qx/qy/qz for all
# rocking angles
q_shape = n_images, img_size[0] * img_size[1], 3
# then the array
q_ar = np.zeros(q_shape, dtype=np.float64)
for index in range(n_images):
qx, qy, qz = hxrd.Ang2Q.area(
eta[index], phi[index], nu[index], delta[index])
q_ar[index, :, 0] = qx.reshape(-1)
q_ar[index, :, 1] = qy.reshape(-1)
q_ar[index, :, 2] = qz.reshape(-1)
return q_ar
class QSpaceConverter(object):
(READY, RUNNING, DONE,
ERROR, CANCELED, UNKNOWN) = __STATUSES = range(6)
......@@ -679,39 +746,13 @@ class QSpaceConverter(object):
center_chan = (center_chan[0] - image_roi_offset[0],
center_chan[1] - image_roi_offset[1])
# TODO : make this editable?
nx, ny, nz = qspace_dims
qconv = xu.experiment.QConversion(['y-', 'z-'],
['z-', 'y-'],
[1, 0, 0])
# convention for coordinate system:
# x downstream
# z upwards
# y to the "outside"
# (righthanded)
hxrd = xu.HXRD([1, 0, 0],
[0, 0, 1],
en=beam_energy,
qconv=qconv)
hxrd.Ang2Q.init_area('z-',
'y+',
cch1=center_chan[0],
cch2=center_chan[1],
Nch1=img_size[0],
Nch2=img_size[1],
chpdeg1=chan_per_deg[0],
chpdeg2=chan_per_deg[1])
# shape of the array that will store the qx/qy/qz for all
# rocking angles
q_shape = n_entries, img_size[0] * img_size[1], 3
# then the array
q_ar = np.zeros(q_shape, dtype=np.float64)
img_dtype = None
phi = []
eta = []
nu = []
delta = []
with XsocsH5.XsocsH5(xsocsH5_f, mode='r') as master_h5:
......@@ -735,15 +776,10 @@ class QSpaceConverter(object):
entry_file))
entry_files.append(entry_file)
phi = np.float64(master_h5.positioner(entry, 'phi'))
eta = np.float64(master_h5.positioner(entry, 'eta'))
nu = np.float64(master_h5.positioner(entry, 'nu'))
delta = np.float64(master_h5.positioner(entry, 'del'))
qx, qy, qz = hxrd.Ang2Q.area(eta, phi, nu, delta)
q_ar[entry_idx, :, 0] = qx.reshape(-1)
q_ar[entry_idx, :, 1] = qy.reshape(-1)
q_ar[entry_idx, :, 2] = qz.reshape(-1)
phi.append(master_h5.positioner(entry, 'phi'))
eta.append(master_h5.positioner(entry, 'eta'))
nu.append(master_h5.positioner(entry, 'nu'))
delta.append(master_h5.positioner(entry, 'del'))
entry_dtype = master_h5.image_dtype(entry=entry)
......@@ -755,6 +791,16 @@ class QSpaceConverter(object):
'be of the same type. Found {0} and {1}.'
''.format(img_dtype, entry_dtype))
# Convert image pixel positions to Q space
q_ar = qspace_conversion(img_size,
center_chan,
chan_per_deg,
beam_energy,
phi,
eta,
nu,
delta)
if mask is not None:
# Mark masked pixels (i.e., non zero in the mask) with NaN
q_ar[:, mask.reshape(-1) != 0, :] = np.nan
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment