diff --git a/xsocs/process/qspace/QSpaceConverter.py b/xsocs/process/qspace/QSpaceConverter.py index d19a98e02a2182c9a5c85609b16abbac65dfd146..715208ada932fe2d06b4207d887215cb7f4009eb 100644 --- a/xsocs/process/qspace/QSpaceConverter.py +++ b/xsocs/process/qspace/QSpaceConverter.py @@ -30,6 +30,7 @@ __license__ = "MIT" __date__ = "01/03/2016" +import logging import os import time import ctypes @@ -47,6 +48,10 @@ from ...util.histogramnd_lut import histogramnd_get_lut, histogramnd_from_lut # from silx.math import histogramnd from ...io import XsocsH5, QSpaceH5, ShiftH5 + +logger = logging.getLogger(__name__) + + disp_times = False @@ -882,7 +887,8 @@ class QSpaceConverter(object): shiftH5=shiftH5, beam_energy=beam_energy, direct_beam=center_chan, - channels_per_degree=chan_per_deg) + channels_per_degree=chan_per_deg, + normalizer=normalizer) manager = mp.Manager() self.__term_evt = term_evt = manager.Event() @@ -967,6 +973,7 @@ class QSpaceConverter(object): entries, img_size, output_f, + normalizer, image_binning, medfilt_dims, img_dtype) @@ -1157,7 +1164,8 @@ def _create_result_file(h5_fn, shiftH5=None, beam_energy=None, direct_beam=None, - channels_per_degree=None): + channels_per_degree=None, + normalizer=''): """ Initializes the output file. :param h5_fn: name of the file to initialize @@ -1181,6 +1189,8 @@ def _create_result_file(h5_fn, :param beam_energy: Beam energy in eV :param direct_beam: Direct beam calibration position :param channels_per_degree: Channels per degree calibration + :param str normalizer: + Name of measurement group dataset used for normalization """ if not overwrite: @@ -1211,6 +1221,10 @@ def _create_result_file(h5_fn, qspace_h5.set_direct_beam(direct_beam) qspace_h5.set_channels_per_degree(channels_per_degree) + if normalizer is None: + normalizer = '' + qspace_h5.set_image_normalizer(normalizer) + if shiftH5: sample_shifts = [] grid_shifts = [] @@ -1240,6 +1254,7 @@ def _to_qspace(th_idx, entries, img_size, output_fn, + normalizer, image_binning, medfilt_dims, img_dtype): @@ -1250,6 +1265,8 @@ def _to_qspace(th_idx, :param entries: :param img_size: :param output_fn: + :param str normalizer: + Name of measurement group dataset to use for normalization :param image_binning: :param medfilt_dims: :param img_dtype: @@ -1292,6 +1309,11 @@ def _to_qspace(th_idx, # h_lut.shape = (n_xy, -1) h_lut = h_lut_shared + if normalizer and img_dtype.kind != 'f': + # Force the type to float64 + logger.info('Using float64 to perform normalization') + img_dtype = np.float64 + img = np.ascontiguousarray(np.zeros(img_size), dtype=img_dtype) # TODO : handle case when nav is not a multiple of img_size!! @@ -1333,6 +1355,7 @@ def _to_qspace(th_idx, image_indices = None for entry_idx, entry in enumerate(entries): + xsocsH5 = XsocsH5.XsocsH5(entry_files[entry_idx], mode='r') t0 = time.time() @@ -1345,9 +1368,7 @@ def _to_qspace(th_idx, # slows down things a big, not much tho) # TODO : add a lock on the files if there is no SWMR # test if it slows down things much - with XsocsH5.XsocsH5(entry_files[entry_idx], - mode='r').image_dset_ctx() \ - as img_data: # noqa + with xsocsH5.image_dset_ctx() as img_data: t1 = time.time() img_data.read_direct(img, source_sel=np.s_[image_idx], @@ -1363,11 +1384,29 @@ def _to_qspace(th_idx, t_read += time.time() - t0 t0 = time.time() + # Apply normalization + + if normalizer: + normalization = xsocsH5.measurement(entry, normalizer) + # Make sure to use float to perform division + assert img.dtype.kind == 'f' + img /= normalization[image_idx] + + # Perform binning + if image_binning[0] != 1 or image_binning[1] != 1: + # Increase type size for (u)int8|16 + if img.dtype.kind == 'u' and img.dtype.itemsize < 4: + binning_dtype = np.uint32 + elif img.dtype.kind == 'i' and img.dtype.itemsize < 4: + binning_dtype = np.int32 + else: # Keep same dtype for (u)int32|64 and float + binning_dtype = img.dtype + intensity = (img.reshape(img_shape_1). sum(axis=sum_axis_1, - dtype=np.uint32).reshape(img_shape_2). - sum(axis=sum_axis_2, dtype=np.uint32) * + dtype=binning_dtype).reshape(img_shape_2). + sum(axis=sum_axis_2, dtype=binning_dtype) * avg_weight) # intensity = xu.blockAverage2D(img, nav[0], # nav[1], roi=roi)