- TODO: lire directement les fichiers NeXus pour la calibration.
- TODO: faire une calibration en utilisant toutes les images d’une série.
- Le bras détecteur bouge et il faut le prendre en compte.

voir ce qu’a fait marco voltolini

    [1]:

    # list of imports
    import os
    import sys
    import time
    import fabio
    import pyFAI
    import numpy
    import numpy.ma

    from typing import Iterator, List, NamedTuple, Optional, Text, Tuple, Union

    from collections import namedtuple
    from functools import partial
    from h5py import Dataset, File
    from numpy import ndarray
    from matplotlib import pyplot
    from scipy.ndimage import generic_filter

    from pyFAI.detectors import detector_factory

    # A local import with Soleil common methodes.
    from soleil import *

    [3]:

    # define a bunch of constants
    LAMBDA = 0.6887
    CALIBRANT = "CeO2"
    CEO2 = "XRD18keV_27.nxs"

    ROOT = os.path.join("/nfs", "ruche-diffabs", "diffabs-soleil", "com-diffabs", "2016", "Run2")
    PUBLISHED_DATA = os.path.join("/nfs", "ruche-diffabs", "diffabs-users", "99160066", "published-data")

    # temporary until the ruch is ON
    ROOT = os.path.join("/home", "experiences", "instrumentation", "picca", "data", "99160066", "2016", "Run2")
    PUBLISHED_DATA = os.path.join("/home", "experiences", "instrumentation", "picca", "data", "99160066", "published-data")

    [3]:

    # define a bunch of generic hdf5 structures.

    DatasetPathContains = NamedTuple("DatasetPathContains", [("path", Text)])

    DatasetPathWithAttribute = NamedTuple("DatasetPathWithAttribute", [('attribute', Text),
                                                                       ('value', bytes)])

    DatasetPath = Union[DatasetPathContains,
                        DatasetPathWithAttribute]

    def _v_attrs(attribute: Text, value: Text, _name: Text, obj) -> Dataset:
        """extract all the images and accumulate them in the acc variable"""
        if isinstance(obj, Dataset):
            if attribute in obj.attrs and obj.attrs[attribute] == value:
                return obj

    def _v_item(key: Text, name: Text, obj: Dataset) -> Dataset:
        if key in name:
            return obj

    def get_dataset(h5file: File, path: DatasetPath) -> Dataset:
        res = None
        if isinstance(path, DatasetPathContains):
            res = h5file.visititems(partial(_v_item, path.path))
        elif isinstance(path, DatasetPathWithAttribute):
            res = h5file.visititems(partial(_v_attrs,
                                            path.attribute, path.value))
        return res

    [4]:

    # Compute the mask and the dark
    # 3.5 could be simplify with 3.6

    DarkMaskSources = NamedTuple('DarkMaskSources', [('filenames', List[Text]),  # name of the files
                                                     ('threshold', int),  # mask all pixel above this treshold
                                                     ('detector', Optional[Text]),
                                                     ('images_path', DatasetPath)])

    def create_dark_and_mask(params: DarkMaskSources) -> Tuple[ndarray, ndarray]:
        """
        génère un dark moyen ainsi qu'un mask avec les dark mesurés. Sur
        un XPAD, le dark doit valoir zero (si le détecteur est bien
        calibré). Nous utilisons tous les pixels dont la valeur est
        supérieur à 2 pour définir un masque (mask ;). Pour la prochaine
        expérience penser à faire non pas un seul dark mais des séries de
        dark.
        """
        for i, filename in enumerate(params.filenames):
            with File(filename, mode='r') as f:
                img = get_dataset(f, params.images_path)[0]  # read the first image
                _mask = numpy.where(img > params.threshold, True, False)
                if i == 0:
                    dark = img.copy()
                    mask = _mask.copy()
                else:
                    dark += img
                    mask = numpy.logical_or(mask, _mask)

        if params.detector is not None:
            det = detector_factory(params.detector)
            mask = numpy.logical_or(mask, det.calc_mask())
            # on a repere des discidents
            mask[480:482, 480:482] = True

        dark = dark.astype(float)
        dark /= len(params.filenames)
        dark = dark.astype('uint32')

        return dark, mask

    dark, mask = create_dark_and_mask(DarkMaskSources([os.path.join(ROOT, "2016-03-26", "dark_%d.nxs" % n) for n in range(7, 12)],
                                      1, 'Xpad_flat',
                                      DatasetPathWithAttribute("interpretation", b"image")))

    [5]:

    # compute the flat

    FlatParams = NamedTuple('FlatParams', [('filename', Text),  # name of the file
                                           ('threshold', float),  # mask all pixel above this treshold
                                           ('dark', Optional[ndarray]),
                                           ('images_path', DatasetPath)])

    def get_flat(params: FlatParams) -> ndarray:
        """
        :param filename: name of the files
        :type filename: list(str)

        génère un flat corrigé du dark si dark is not None
        """
        with File(params.filename, mode='r') as f:
            images = get_dataset(f, params.images_path)[:]
            flat = images.mean(axis=0)
            if dark is not None:
                flat -= dark

            flat = numpy.where(numpy.abs(flat) <= params.threshold, params.threshold, flat)

            return flat

    flat = get_flat(FlatParams(os.path.join(ROOT, "2016-03-26", "flat_12.nxs"),
                               1, dark,
                               DatasetPathWithAttribute("interpretation", b"image")))

    [25]:

    # Unfold the images

    # defines the namedtuple used for the tomogaphie.
    # 3.5 could be simplify with 3.6

    TomoSources = NamedTuple('TomoSources', [('filename', Text),
                                             ('images_path', DatasetPath),
                                             ('rotations_path', DatasetPath),
                                             ('translations_path', DatasetPath)])

    TomoFrame = NamedTuple('TomoFrame', [("image", ndarray),
                                         ("shape", Tuple[int, int]),
                                         ("index", Tuple[int, int]),
                                         ("rotation", float),
                                         ("translation", float)]
                          )

    UnfoldFrame = NamedTuple('UnfoldFrame', [("tomoframe", TomoFrame),
                                             ("unfolded", ndarray)])

    Unfold = NamedTuple('Unfold', [('sources', TomoSources),  # sources of the data's
                                   ('poni', Text),  # the poni file use to do the unfold
                                   ('mask', Optional[ndarray]),  # name of the poni file used for the integration.
                                   ('dark', Optional[ndarray]),  # mask used for the computation.
                                   ('flat', Optional[ndarray]),  # number of bins used for the powder spectrum.
                                   ('npt_rad', int),  # he number of bins used for the powder spectrum.
                                   ('npt_azim', int),  # he number of bins used for the powder spectrum.
                                  ]
                       )

    TomoSave = NamedTuple('TomoSave', [('volume', ndarray),
                                       ('rotation', ndarray),
                                       ('translation', ndarray)])

    def read_multi(params: TomoSources) -> Iterator[TomoFrame]:
        """
        :param sources: list des fichiers à traiter.
        :type sources: list(str)

        :return: yield frames of data contain in all NxEntries located at the data_path location.
        :rtype: Frame
        """
        with File(params.filename, mode='r') as f:
            images = get_dataset(f, params.images_path)
            rotations = get_dataset(f, params.rotations_path)
            translations = get_dataset(f, params.translations_path)
            shape = images.shape[0], images.shape[1]
            for rotation_idx, rotation in enumerate(rotations):
                for translation_idx, translation in enumerate(translations):
                    yield TomoFrame(images[rotation_idx, translation_idx,:],
                                    shape,
                                    (rotation_idx, translation_idx),
                                    rotation, translation)

    def unfold(params: Unfold) -> Iterator[UnfoldFrame]:
        """
        :return: the tomography datas
        :rtype: numpy.ndarray

        the return data is a 5 dimensions numpy array (rotation,
        translation, 2, Nr, Na). the third dimension contain the x and y
        coordinates of the powder diffraction (0 for X and 1 for Y)
        """
        # load the Azimuthal Integration parameters
        ai = pyFAI.load(params.poni)

        for tomoframe in read_multi(params.sources):
            data = tomoframe.image

            """ mask all the non counting pixels """
            # TODO comment optimiser ce calcul du mask dynamic qui prend vraiment du temps
            # we will use dummy value (0) in order to define the dynamic mask
            data[data>500] = 0

            unfolded, r, a = ai.integrate2d(data, params.npt_rad, params.npt_azim,
                                            filename=None, dark=dark, flat=flat, dummy=0)

            yield UnfoldFrame(tomoframe, unfolded)

    def save_unfolded(filename: Text,
                      params: Unfold) -> None:
        with File(filename, mode="r+") as f:
            # now fill with the values.
            for i, frame in enumerate(unfold(params)):
                if i == 0:
                    # first create output volums at the right place
                    group = f.require_group("tomography_unfolded")
                    unfolded = group.require_dataset("unfolded", shape=frame.tomoframe.shape + (params.npt_azim, params.npt_rad),
                                                     dtype=frame.unfolded.dtype)
                    rotation = group.require_dataset("rotation", shape=(frame.tomoframe.shape[0],),
                                                     dtype='float32')
                    translation = group.require_dataset("translation", shape=(frame.tomoframe.shape[1],),
                                                        dtype='float32')
                    if dark is not None:
                        dark = group.require_dataset("dark", data=params.dark, dtype='uint32')
                        dark = params.dark
                    if flat is not None:
                        group.require_dataset("flat", data=params.flat, dtype='uint32')
                    if mask is not None:
                        group.require_dataset("mask", data=params.mask.astype('uint32'), dtype='uint32')

                unfolded[frame.tomoframe.index[0], frame.tomoframe.index[1], :] = frame.unfolded
                rotation[frame.tomoframe.index[0]] = frame.tomoframe.rotation
                translation[frame.tomoframe.index[1]] = frame.tomoframe.translation

    PONI = os.path.join(PUBLISHED_DATA, 'calibration', 'XRD18keV_27_0.poni')
    PONI = os.path.join(PUBLISHED_DATA, 'calibration', 'XRD18keV_26.nxs_03.poni')

    TOMOSOURCES = TomoSources(os.path.join(ROOT, "2016-03-27", 'P14_13_57.nxs'),
                              DatasetPathWithAttribute("interpretation", b"image"),
                              DatasetPathContains("scan_data/trajectory_2_1"),
                              DatasetPathContains("scan_data/trajectory_1_1"))

    t0 = time.time()
    params = Unfold(TOMOSOURCES, PONI, mask, dark, flat, 1000, 360)
    save_unfolded(os.path.join(PUBLISHED_DATA, "P14_13_57_unfold.h5"),
                  params)
    print("unfold time: ", time.time() - t0)

    ---------------------------------------------------------------------------
    TypeError                                 Traceback (most recent call last)
    <ipython-input-25-8186de126bd4> in <module>()
        118 params = Unfold(TOMOSOURCES, PONI, mask, dark, flat, 1000, 360)
        119 save_unfolded(os.path.join(PUBLISHED_DATA, "P14_13_57_unfold.h5"),
    --> 120               params)
        121 print("unfold time: ", time.time() - t0)

    <ipython-input-25-8186de126bd4> in save_unfolded(filename, params)
         96                                                   dtype='float32')
         97                 if dark is not None:
    ---> 98                     group.require_dataset("dark", data=params.dark)
         99                 if flat is not None:
        100                     group.require_dataset("flat", data=params.flat)

    TypeError: require_dataset() missing 2 required positional arguments: 'shape' and 'dtype'

    [ ]:

    # ici tout un tas de variables utilisées pour trouver les fichiers sur
    # le file system SOLEIL. Il y a de grand chance que cela ne
    # fonctionne pas pour vous :p donc veiller à adapter ces chemins.

    def convert(filename, numbers=None, data_path='scan_data/data_15'):
        """
        :param filename: name of the nexus file
        :type filename: str
        :param numbers: index of the image to extract if None extract all images
        :type numbers: list(int)

        conversion d'image nxs -> edf. L'image edf est ensuite utilisée
        avec pyFAI-calib pour générer le fichier .poni
        """
        imgs = read(filename, numbers, data_path)
        for index, img in enumerate(imgs):
            edf = fabio.edfimage.edfimage(img)
            name = os.path.splitext(os.path.basename(filename))[0]
            name += '_%d' % index + '.edf'
            saved_name = os.path.join(TREATMENT_PATH, name)
            edf.write(saved_name)

    def integrate2(sources, poni, mask=None, dark=None, flat=None, N=600,
                   data_path=None,
                   rotation_path=None,
                   translation_path=None,
                   monitor_path=None):
        """
        :param sources: list of nexus files to treat
        :type sources: list(str)
        :param poni: the name of the poni file used for the integration.
        :type poni: str
        :param mask: the mask used for the computation.
        :type mask: None or numpy.ndarray
        :param N: the number of bins used for the powder spectrum.
        :type N: int
        :param data_path: location of the images in the nexus file
        :type data_path: str
        :param rotation_path: location of the rotation coordinated in the nexus file
        :type rotation_path: str
        :param translation_path: location of the translation coordinates in the nexus file
        :type translation_path: str

        :return: the tomography datas
        :rtype: numpy.ndarray

        the return data is a 4 dimensions numpy array (rotation,
        translation, 2, N). the third dimension contain the x and y
        coordinates of the powder diffraction (0 for X and 1 for Y)
        """
        # load the Azimuthal Integration parameters
        ai = pyFAI.load(poni)

        volume = numpy.array([])
        rotation = numpy.array([])
        translation = numpy.array([])
        for i, frame in enumerate(read_multi(sources,
                                             data_path,
                                             rotation_path,
                                             translation_path,
                                             monitor_path)):
            if i == 0:
                volume = numpy.empty(frame.shape + (2, N))
                rotation = numpy.empty((frame.shape[0],))
                translation = numpy.empty((frame.shape[1],))

            data = frame.img

            """ mask all the non counting pixels """
            # TODO comment optimiser ce calcul du mask dynamic qui prend vraiment du temps

            # we will use dummy value (0) in order to define the dynamic mask
            data[data>500] = 0

            #mask_data = numpy.where(data == 0, True, False)
            #mask_data = numpy.logical_or(mask_data, numpy.where(data > 500, True, False))
            # mask_data = numpy.logical_or(mask_data,
            #                             numpy.where(data > 500, True, False))
            #if mask is not None:
            #    mask_data = numpy.logical_or(mask_data, mask)

            #filename = os.path.splitext(f)[0] + ".txt"
            t0 = time.time()
            spectrum_x_y = ai.integrate1d(data, N, filename=None, dark=dark, flat=flat, dummy=0, method="csr_ocl_0,0") #dark=dark, flat=flat) #, method="lut") #, method="csr_ocl_0,1")
            #spectrum_x_y = ai.integrate1d(data, N, mask=mask_data, method="csr_ocl_0,0") #, method="csr_ocl_0,1")
            print(i, frame.rotation, frame.translation, time.time() - t0)
            if frame.monitor:
                volume[frame.rotation, frame.translation] = spectrum_x_y[0], spectrum_x_y[1] / frame.monitor
            else:
                volume[frame.rotation, frame.translation] = spectrum_x_y

            rotation[frame.rotation] = frame.rotation_value
            translation[frame.translation] = frame.translation_value

        return volume, rotation, translation

    def from_hdf5(filename):
        """
        :param filename: name of the file where the tomographie was stored
        :type filename: str
        :return: volume, rotation, translation
        :rtype: truple(numpy.ndarray, numpy.ndarray, numpy.ndarray)

        read the datas from the filename
        """
        with tables.openFile(filename, mode="r") as h5file:
            volume = h5file.getNode("/tomography/diffraction")[:]
            rotation = h5file.getNode("/tomography/rotation")[:]
            translation = h5file.getNode("/tomography/translation")[:]

            return volume, rotation, translation

    def load_volume(filename, basedir):
        return from_hdf5(os.path.join(basedir, filename + ".h5"))

    # ugly but for now it is ok
    SIZE = 3

    def reject_outliers(data, m=5.):
        """
        :param data: the data to filter
        :type data: numpy.ndarray
        :param m: the filter rejection parameter
        :type m: float
        """
        center = SIZE / 2

        median = numpy.median(data)
        d = numpy.abs(data - median)
        mdev = numpy.median(d)
        return data[center] if (d[center] / mdev) < m else median

    def filter_sinogramme(data):
        """
        :param volume: the volume containing all the sinogram
        :param type: numpy.ndarray
        """
        return generic_filter(data, reject_outliers, size=SIZE)

    def all_conversions():
        convert(SI_0_100, data_path='scan_data/data_03')
        convert(SI_14, data_path='scan_data/data_15')

    def padd_sino(sino, offset=0):
        n = (256 - sino.shape[0]) / 2
        n1 = n + offset
        n2 = 256 - n1 - sino.shape[0]
        padd1 = numpy.tile(sino[0], (n1, 1))
        padd2 = numpy.tile(sino[-1], (n2, 1))
        tmp = numpy.vstack([padd1, sino, padd2])
        print(tmp.shape)
        return tmp

    ############
    # P7_19_12 #
    ############

    def p14_13_57():

        def create_unfold():
            t0 = time.time()
            params = Unfold(TOMOSOURCES, PONI, mask, dark, flat)
            save = unfold(params)
            print("unfold time: ", time.time() - t0)
            # TODO sauver la transposé de façon a ce que les données
            # soient organisée correctement pour les reconstructions et la
            # filtration look at
            to_hdf5(os.path.join(PUBLISHED_DATA, "P14_13_57_unfold.h5"),
                    save, params)

        def create_volume():
            t0 = time.time()
            dark, mask = create_mask_and_dark(MASKS_15S, detector='Xpad_flat', data_path=DATAPATH)
            #with h5py.File(P7_19_12, mode='r') as f:
            #    dark += f["scan_13/scan_data/data_58"][35, 65]
            flat = get_flat(FLAT, dark=dark, data_path=DATAPATH, threshold=1)  # TODO moyenner les flats. 12, 13
            volume, rotation, translation = \
                    integrate2(P14_13_57, PONI, mask=mask,
                               dark=dark, flat=flat,
                               data_path=DATAPATH,
                               rotation_path="scan_data/trajectory_2_1",
                               translation_path="scan_data/trajectory_1_1")
                               # monitor_path="scan_data/data_01")

            print("integration time: ", time.time() - t0)
            # TODO sauver la transposé de façon a ce que les données
            # soient organisée correctement pour les reconstructions et la
            # filtration look at
            save_volume("P14_13_57", volume[:,:, 1,:], rotation, translation, dark, flat, mask, basedir=PUBLISHED_DATA)

        def filter_volume():
            volume, rotation, translation = load_volume("P14_13_57", basedir=PUBLISHED_DATA)
            new_volume = volume.T.copy()
            for i, img in enumerate(new_volume):
                print("treat the %i sinogram" % i)
                new_volume[i] = filter_sinogramme(img)
                # new_volume[i] = medfilt2d(img)

            save_volume("P14_13_57_filtered", new_volume, rotation, translation, basedir=PUBLISHED_DATA)

        def _reconstruction(filename, transpose=False):
            CHANNEL = 366

            from skimage.transform import iradon
            volume, rotation, translation = load_volume(filename, basedir=PUBLISHED_DATA)
            if transpose:
                volume = volume.T.copy()
            print(rotation.shape)
            # sino = volume.mean(axis=0)
            # sino = volume[192:201].mean(axis=0)[:, :-1] # quartz
            sino = volume[CHANNEL][:, :-1]  # iron
            # sino = volume[52:63].mean(axis=0)[:,:-1]
            sino = padd_sino(sino, 3)
            pyplot.subplot(121)
            pyplot.title(filename + str(CHANNEL) + " channel")
            pyplot.imshow(sino, interpolation='nearest')
            pyplot.colorbar()
            rec = iradon(sino, theta=rotation[:-1], filter='shepp-logan')
            pyplot.subplot(122)
            pyplot.title(filename + str(CHANNEL) + " channel")
            pyplot.imshow(rec, interpolation='nearest')

        def reconstruction():
            _reconstruction("P14_13_57_filtered")
            pyplot.figure()
            _reconstruction("P14_13_57", transpose=True)
            pyplot.show()

        create_unfold()
        #create_volume()
        #filter_volume()
        #reconstruction()

    #p14_13_57()

    [2]:

    [ ]: