diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index dbfbf34df612d9923cfe668852796991adfac87d..68321144e096b13ec738d928d6a83a47d9e2383c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -67,7 +67,7 @@ doc: - install_silx - install_tomoscan - python -m pip install pytest-cov - - python -m pip install -e .[full] + - python -m pip install -e . script: - python -m pytest --cov-config=.coveragerc --cov=nxtomomill nxtomomill/ @@ -98,7 +98,7 @@ test:test-nxtomomill-tutorials: - python -m pip install tomoscan --upgrade --pre - python -m pip install .[doc] script: - - jupyter nbconvert --to notebook --execute doc/tutorials/create_nxtomo_from_scratch.ipynb + - jupyter nbconvert --to notebook --execute doc/tutorials/nx_tomo_tutorial.ipynb # deploy diff --git a/doc/design/hdf5_seq.txt b/doc/design/hdf5_converter_seq.txt similarity index 52% rename from doc/design/hdf5_seq.txt rename to doc/design/hdf5_converter_seq.txt index 1c8be2368cf1a6efcdf1c90e656c69840428412f..d01c830096c8994e4d951286ccba6c20c6bdeee2 100644 --- a/doc/design/hdf5_seq.txt +++ b/doc/design/hdf5_converter_seq.txt @@ -1,4 +1,5 @@ -title This is a title +# done from https://sequencediagram.org/ +title bliss hdf5 to nexus hdf5 converter user->ConfigHandler: provide configuration file\n and command options ConfigHandler->ConfigHandler: check options\nvalidity @@ -8,4 +9,9 @@ Configuration->user:return user->converter: request conversion from a configuration converter->converter: parse input files / entries converter->converter: create `BaseAcquisition` instances\n(from titles or from urls) -converter->converter: write NXTomo +converter->NXtomo: generate NXtomo instances from `BaseAcquisition` instances +NXtomo->converter:return +converter->NXtomo: Optional: apply modification to NXtomo instances from plugins +NXtomo->converter:return +converter->NXtomo: write NXtomo to disk +NXtomo->NXtomo: save \ No newline at end of file diff --git a/doc/design/img/hdf5_sequence_diagram.png b/doc/design/img/hdf5_sequence_diagram.png index da5f164eee73d7690fd153db7d42efa995e0d09a..bbc1257ee330fa154c806490224edc5be7721c58 100644 Binary files a/doc/design/img/hdf5_sequence_diagram.png and b/doc/design/img/hdf5_sequence_diagram.png differ diff --git a/doc/design/index.rst b/doc/design/index.rst index e6078a9af52514bf7e80fc95f4e0058d1e32e591..48df065011816d676933a3c6e590d24135a52199 100644 --- a/doc/design/index.rst +++ b/doc/design/index.rst @@ -7,8 +7,35 @@ HDF5Converter ''''''''''''' The more it goes the more use cases the hdf5 converter has to handle. -Today it handles "classical" tomography acquisition and zseries. -Tomorrow it will handle 3d-XRD and XRD-CT. +Today it handles "classical" tomography acquisition, zseries and pcotomo. +In the future it is expected to handle XRD-CT, XFD-CT and dual energy - CT. + +Behavior of the HDF5Converter +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The HDF5 is applying sequentially the following taks: + +1. preprocess raw data (During converter creation) + + browse raw data and collect information on Bliss scan sequence. + + During this step it will determine for each entry if the frame type are dark, flats, projections or alignment. + + Create an instance of Acquistion for each 'init' Bliss scan. + +2. HDF5Converter.convert function + + 2.1 create for Acquisition instance one or several instances of :class:`NXtomo`. + + 2.2 Optional split of NXtomo (use case of pcotomo) in order to obtain the final NXtomo. + + 2.3 apply plugin (allows user to modify NXtomo instances). + + 2.4 save data to disk. + + +Acquistion classes +^^^^^^^^^^^^^^^^^^ To manage all the use cases it has evolved. Each acquisition has now it's own acquisition class. All based on the :class:`nxtomomill.converter.hdf5.acquisition.baseacquisition.BaseAcquisition` diff --git a/doc/plugins.rst b/doc/plugins.rst index 8bffbd79f05cc6a8adba953282483da1a80c46b8..020c332d20c1b60a7e06814162ecea5c538e1369 100644 --- a/doc/plugins.rst +++ b/doc/plugins.rst @@ -3,23 +3,26 @@ Plugins ======= -In some cases all the information can be accessible but you might need to do some operations to get the valid NXTomo information. -This is why we have a simple plugin solution to define new 'arrays' and give you hand to create those from original bliss-hdf5 source file. +HDF5Plugin +"""""""""" -All the plugins will be discover during the conversion and will be apply as the file step of the output file writing. +In some cases you might want to modify information contained on an NXtomo. For example if you want to do some operations between raw dataset or append some information that does not appears on the raw file. -The process to create a plugin is: +This is the goal of the 'plugin mecanism'. + +All the plugins will be discover during the bliss-hdf5 to nexus_hdf5 conversion and will be apply just before saving NXtomo(s) to disk. + +In order to create a plugin you must: - create a directory that will contain the plugins. - define the environment variable 'NXTOMOMILL_PLUGINS_DIR' which register the plugin directory - inside the plugin directory create one or several python file which will implement the 'HDF5Plugin' class. -:warning: units are not managed by the plugin system. So you should know the units you are reading and writing. **Writing units is mandatory**. +You can retrieve data from bliss-hdf5 acquisition contained either in `instrument` or in `positioner` groups. -For now only the data contained in the 'positioners' groups are accessible. -For performance issues you also need to declare the positioners names you will access. +For performance issues, as we collect this information in the NXtomo creation preprocessing you must declare `dataset` name you want to retrieve from those groups. -The Plugin will look like: +The HDF5Plugin provide an API to modify the NXtomo before they will be saved on disk. This way you can modify them as you like using the NXtomo API by implementing the `update_nx_tomo` funciton. .. code-block:: python @@ -37,11 +40,9 @@ The Plugin will look like: positioners_names = 'diffty', 'd1ty' super().__init__(name='My plugin 1', positioner_keys=positioners_names) - def write(self, root_node, sample_node, detector_node, beam_node) -> list: + def update_nx_tomo(self, nx_tomo) -> nx_tomo: diffty = numpy.array(self.get_positioner_info('diffty')) d1ty = numpy.array(self.get_positioner_info('d1ty')) - if 'x_translation' in sample_node: - del sample_node['x_translation'] - sample_node['x_translation'] = diffty - d1ty - 14 - sample_node['x_translation'].attrs['units'] = 'mm' + nx_tomo.sample.x_translation = diffty - d1ty - 14 + nx_tomo.sample.x_translation.unit = "mm" diff --git a/doc/tutorials/create_nxtomo_from_scratch.ipynb b/doc/tutorials/create_nxtomo_from_scratch.ipynb deleted file mode 100644 index b08a7ce753873e30a9a3e6a64ed0553b2b465cff..0000000000000000000000000000000000000000 --- a/doc/tutorials/create_nxtomo_from_scratch.ipynb +++ /dev/null @@ -1,788 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# create an NXTomo from scratch\n", - "The goal of this tutorial is to create an [NXTomo](https://manual.nexusformat.org/classes/applications/NXtomo.html) from scratch. In the sense without converting **directly** from a bliss scan." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## description of the example" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let say we want to create an NXtomo matching the following sequence\n", - "\n", - "| frame index | Rotation Angle (in degree) | Frame Type | control Image Key | Image Key |\n", - "| ----------- | --------------------------- | ----------- | ------------------ | --------- |\n", - "| 0 | 0 | Dark | 2 | 2 |\n", - "| 1 | 0 | Flat | 1 | 1 | \n", - "| 2-201 | 0 - 89.9 | Projection | 0 | 0 | \n", - "| 202 | 90 | Flat | 1 | 1 | \n", - "| 203 - 402 | 90 - 180 | Projection | 0 | 0 | \n", - "| 403 | 180 | Flat | 1 | 1 | \n", - "| 404 | 90 | Alignment | -1 | 0 |" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## create dummy dataset" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Note: In order to siplify the setup we will take only create one dark and one flat frame that we will reuse. But keep in mind that those are raw data.\n", - "So in 'real' life this is expected that you will have several frame for dark and several frames for flats and there will differ of course depending on when you acquire them.*" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# %pylab inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from skimage.data import shepp_logan_phantom\n", - "from skimage.transform import radon\n", - "import numpy" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### create some projection" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "phantom = shepp_logan_phantom()\n", - "projections = {}\n", - "proj_rotation_angles = numpy.linspace(0., 180., max(phantom.shape), endpoint=False)\n", - "sinogram = radon(phantom, theta=proj_rotation_angles)\n", - "\n", - "sinograms = numpy.asarray([sinogram] * 20)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# imshow(phantom)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "radios = numpy.swapaxes(sinograms, 2, 0)\n", - "radios = numpy.swapaxes(radios, 1, 2)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# imshow(radios[0])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### create some dark" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "max_shape = max(phantom.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dark = numpy.zeros((20, max_shape))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# imshow(dark)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### create some flat" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "flat = numpy.ones((20, max_shape))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### add some noise to radios" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tmp_radios = []\n", - "for radio in radios:\n", - " tmp_radios.append(dark + radio * (flat - dark))\n", - "radios = numpy.asarray(tmp_radios)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# imshow(radios[0])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### create some alignment\n", - "In order to keep it simple we will pick one of the radio created" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "alignment = radios[200]\n", - "alignment_angle = proj_rotation_angles[200]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## create an NXtomo that fits the sequence we want" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from nxtomomill.nexus.nxtomo import NXtomo" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo = NXtomo()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### provide mandatory data for Contrast Tomography\n", - "Mandatory information for contrast tomography are:\n", - "* detector frames: raw data\n", - "* image-key (control): for each frame the \"type\" of frame (projections, flats, darks and alignment).\n", - "* rotation angles: for each frame the rotation angle in degree" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### detector frames\n", - "\n", - "to fit the sequence describe previously we need to create the following sequence: dark, flat, first half of the projections, flat, second half of the projections, flat and alignment frame.\n", - "\n", - "And we need to provide them as a numpy array (3d)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# reshape dark, flat and alignment that need to be 3d when numpy.concatenate is called\n", - "darks_stack = dark.reshape(1, dark.shape[0], dark.shape[1])\n", - "flats_stack = flat.reshape(1, flat.shape[0], flat.shape[1])\n", - "alignment_stack = alignment.reshape(1, alignment.shape[0], alignment.shape[1])\n", - "\n", - "assert darks_stack.ndim == 3\n", - "assert flats_stack.ndim == 3\n", - "assert alignment_stack.ndim == 3\n", - "assert radios.ndim == 3\n", - "print(\"radios shape is\", radios.shape)\n", - "# create the array\n", - "data = numpy.concatenate([\n", - " darks_stack,\n", - " flats_stack,\n", - " radios[:200],\n", - " flats_stack,\n", - " radios[200:],\n", - " flats_stack,\n", - " alignment_stack,\n", - "])\n", - "assert data.ndim == 3\n", - "print(data.shape)\n", - "# then register the data to the detector\n", - "my_nxtomo.detector.data = data" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### image key control" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from nxtomomill.utils import ImageKey\n", - "\n", - "image_key_control = numpy.concatenate([\n", - " [ImageKey.DARK_FIELD] * 1,\n", - " [ImageKey.FLAT_FIELD] * 1,\n", - " [ImageKey.PROJECTION] * 200,\n", - " [ImageKey.FLAT_FIELD] * 1,\n", - " [ImageKey.PROJECTION] * 200,\n", - " [ImageKey.FLAT_FIELD] * 1,\n", - " [ImageKey.ALIGNMENT] * 1,\n", - "])\n", - "\n", - "# insure with have the same number of frames and image key\n", - "assert len(image_key_control) == len(data)\n", - "# print position of flats in the sequence\n", - "print(\"flats are at indexes\", numpy.where(image_key_control == ImageKey.FLAT_FIELD))\n", - "# then register the image keys to the detector\n", - "my_nxtomo.detector.image_key_control = image_key_control" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### rotation angle" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "rotation_angle = numpy.concatenate([\n", - " [0.0, ], \n", - " [0.0, ], \n", - " proj_rotation_angles[:200],\n", - " [90.0, ],\n", - " proj_rotation_angles[200:],\n", - " [180.0, ],\n", - " [90.0, ],\n", - "])\n", - "assert len(rotation_angle) == len(data)\n", - "# register rotation angle to the sample\n", - "my_nxtomo.sample.rotation_angle = rotation_angle" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Field of view\n", - "field of view can either be `Half` or `Full`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo.detector.field_of_view = \"Full\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### pixel size" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo.detector.x_pixel_size = my_nxtomo.detector.y_pixel_size = 1e-7 # pixel size must be provided in SI: meter" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "but for attribute with a unit you can specify the unit the value should be \"converted to\" using the 'unit' attribute like:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo.detector.x_pixel_size = my_nxtomo.detector.y_pixel_size = 0.1\n", - "my_nxtomo.detector.x_pixel_size.unit = my_nxtomo.detector.y_pixel_size.unit = \"micrometer\"" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "When the unit is provided it will be stored as a property of the dataset. It must be interpreted by the software reading the NXtomo." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### save the nx to disk" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "nx_tomo_file_path = os.path.join(\"resources_cfg_create_nxtomo_from_scratch\", \"nxtomo.nx\")\n", - "my_nxtomo.save(file_path=nx_tomo_file_path, entry=\"entry\", overwrite=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### check data saved\n", - "We can use some validator to insure we have enought data to be treated by nabu" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from tomoscan.esrf import HDF5TomoScan\n", - "from tomoscan.validator import ReconstructionValidator\n", - "scan = HDF5TomoScan(nx_tomo_file_path, entry=\"entry\")\n", - "validator = ReconstructionValidator(scan, check_phase_retrieval=False, check_values=True)\n", - "assert validator.is_valid()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can check the layout of the file to insure it seems valid as well" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from h5glance import H5Glance\n", - "H5Glance(nx_tomo_file_path)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A good pratice is also to check frames, image_key and rotation angles to insure values seems valid." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ! silx view resources_cfg_create_nxtomo_from_scratch/nxtomo.nx" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### reconstruct using nabu\n", - "now that we have a valid nxtomo we are able to reconstruct it using [nabu](https://gitlab.esrf.fr/tomotools/nabu).\n", - "\n", - "We create a nabu configuration file for contrast tomography reconstruction name `nabu-ct.conf` in order to reconstruct one slice of the volume.\n", - "\n", - "*note: on the configuration you must disable take_logarithm due to the dataset.* \n", - "\n", - "if nabu is installed you can run it:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# ! nabu nabu-cf.conf" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### clean" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "if os.path.exists(nx_tomo_file_path):\n", - " os.remove(nx_tomo_file_path)\n", - "if os.path.exists(\"nxtomo_reconstruction.hdf5\"):\n", - " os.remove(\"nxtomo_reconstruction.hdf5\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### provide mandatory data for Phase Tomography\n", - "in order to compute Phase Tomography you must also register:\n", - "* incoming beam energy (in keV)\n", - "* sample / detector distance (in meter)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "we can take back the existing `my_nxtomo` and add it the missing elements" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo.energy = 12.5 # in keV by default\n", - "my_nxtomo.sample.distance = 0.2 # in meter" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And then you can reconstruct it with phase retrieval from modifing the nabu configuration file." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### provide more metadata" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "you can also provide x, y and z translation of the sample during the acquisition." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo.detector.x_translation = [0, 12]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "as a sample name, source information, start and end time" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_nxtomo.sample.name = \"my sample\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datetime import datetime\n", - "my_nxtomo.source.name = \"ESRF\" # default value\n", - "my_nxtomo.source.type = \"Synchrotron X-ray Source\" # default value\n", - "my_nxtomo.start_time = datetime.now()\n", - "my_nxtomo.stop_time = datetime(2022, 2, 27)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### tomoscan\n", - "once save on disk you can access it from a python API using the [tomoscan](https://tomotools.gitlab-pages.esrf.fr/tomoscan/) library like:\n", - "```\n", - "from tomoscan.esrf import HDF5TomoScan\n", - "scan = HDF5TomoScan(file_path, entry_name)\n", - "```" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Advance usage: using FrameAppender or how to handle a large amount of data in memory\n", - "The issue of using NXtomo as we presented above is that the memory to handle data can be used (if you have a large number of projections and / or large detector).\n", - "\n", - "The way to work around for now it to use the FrameAppender class.\n", - "\n", - "You will first save metadata to the hdf5 (and maybe some frame) then you can append to the dataset series of frames sequentially with there rotation angle, image key...\n", - "This way you have to load in memory only the frame to append.\n", - "\n", - "**FrameAppender only hande frames**\n", - "\n", - "Here is a quick example on how it will look like." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### create a simple nxtomo with all information at the exception of frames" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_large_nxtomo = NXtomo()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "provide all information at the exception of frames. Here lets say we will have a dataset with only 180 projections" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_large_nxtomo.detector.distance = 0.2\n", - "my_large_nxtomo.detector.x_pixel_size = my_large_nxtomo.detector.y_pixel_size = 1e-7\n", - "my_large_nxtomo.energy = 12.3\n", - "# ...\n", - "my_large_nxtomo.sample.rotation_angle = numpy.linspace(0, 180, 180, endpoint=False)\n", - "my_large_nxtomo.detector.image_key_control = [0] * 180 # 0 == Projection" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "now that all metadata are store we can save the file" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "my_large_nxtomo.save(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\", \"entry0000\", overwrite=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "once the data is on disk we can append missing frames (**and only frames**) to it.\n", - "\n", - "On this example we will append one frame at the time but of course this is up to you.\n", - "\n", - "In order to ease handling of virtual dataset you can provide data url to it or numpy array. For this example we will use numpy array directly." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from nxtomomill.utils.frameappender import FrameAppender\n", - "for i_frame in range(180):\n", - " frame_data = numpy.random.random(100 * 100).reshape([1, 100, 100]) # of course here this is most likely that you will load data from another file\n", - " frame_appender = FrameAppender(\n", - " data=frame_data,\n", - " file_path=\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\",\n", - " data_path=\"entry0000/instrument/detector/data\",\n", - " where=\"end\",\n", - " )\n", - " frame_appender.process()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Then you can see that the 'data' dataset now contains 180 frames (if you run several time the previous cell then it will continue appending data to it).\n", - "\n", - "If an url is provided instead of a numpy array then it will create be used to create a virtual dataset and avoid duplicating data. But be carreful in this case you must keep relative position of the two files.\n", - "\n", - "**append frames must have the same dimensions otherwise the operation will fail**" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "H5Glance(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "clean" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "if os.path.exists(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\"):\n", - " os.remove(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.3" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/doc/tutorials/nx_tomo_tutorial.ipynb b/doc/tutorials/nx_tomo_tutorial.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..fe6c3a28f1db8fa283f2fd8d19f8d4fb958a993a --- /dev/null +++ b/doc/tutorials/nx_tomo_tutorial.ipynb @@ -0,0 +1,1341 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# create an NXTomo (from scratch)\n", + "The goal of this tutorial is to create an [NXTomo](https://manual.nexusformat.org/classes/applications/NXtomo.html) from scratch. In the sense without converting **directly** from a bliss scan." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## description of the example" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let say we want to create an NXtomo matching the following sequence\n", + "\n", + "| frame index | Rotation Angle (in degree) | Frame Type | control Image Key | Image Key |\n", + "| ----------- | --------------------------- | ----------- | ------------------ | --------- |\n", + "| 0 | 0 | Dark | 2 | 2 |\n", + "| 1 | 0 | Flat | 1 | 1 | \n", + "| 2-201 | 0 - 89.9 | Projection | 0 | 0 | \n", + "| 202 | 90 | Flat | 1 | 1 | \n", + "| 203 - 402 | 90 - 180 | Projection | 0 | 0 | \n", + "| 403 | 180 | Flat | 1 | 1 | \n", + "| 404 | 90 | Alignment | -1 | 0 |" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## create dummy dataset" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Note: In order to siplify the setup we will take only create one dark and one flat frame that we will reuse. But keep in mind that those are raw data.\n", + "So in 'real' life this is expected that you will have several frame for dark and several frames for flats and there will differ of course depending on when you acquire them.*" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# %pylab inline" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from skimage.data import shepp_logan_phantom\n", + "from skimage.transform import radon\n", + "import numpy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### create some projection" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "phantom = shepp_logan_phantom()\n", + "projections = {}\n", + "proj_rotation_angles = numpy.linspace(0., 180., max(phantom.shape), endpoint=False)\n", + "sinogram = radon(phantom, theta=proj_rotation_angles)\n", + "\n", + "sinograms = numpy.asarray([sinogram] * 20)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# imshow(phantom)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "radios = numpy.swapaxes(sinograms, 2, 0)\n", + "radios = numpy.swapaxes(radios, 1, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# imshow(radios[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### create some dark" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "max_shape = max(phantom.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "dark = numpy.zeros((20, max_shape))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# imshow(dark)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### create some flat" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "flat = numpy.ones((20, max_shape))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### add some noise to radios" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "tmp_radios = []\n", + "for radio in radios:\n", + " tmp_radios.append(dark + radio * (flat - dark))\n", + "radios = numpy.asarray(tmp_radios)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "# imshow(radios[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### create some alignment\n", + "In order to keep it simple we will pick one of the radio created" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "alignment = radios[200]\n", + "alignment_angle = proj_rotation_angles[200]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "## create an NXtomo that fits the sequence we want" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "import tomoscan.nexus" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "from nxtomomill.nexus.nxtomo import NXtomo" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo = NXtomo(\"entry\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### provide mandatory data for Contrast Tomography\n", + "Mandatory information for contrast tomography are:\n", + "* detector frames: raw data\n", + "* image-key (control): for each frame the \"type\" of frame (projections, flats, darks and alignment).\n", + "* rotation angles: for each frame the rotation angle in degree" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### detector frames\n", + "\n", + "to fit the sequence describe previously we need to create the following sequence: dark, flat, first half of the projections, flat, second half of the projections, flat and alignment frame.\n", + "\n", + "And we need to provide them as a numpy array (3d)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "radios shape is (400, 20, 400)\n", + "(405, 20, 400)\n" + ] + } + ], + "source": [ + "# reshape dark, flat and alignment that need to be 3d when numpy.concatenate is called\n", + "darks_stack = dark.reshape(1, dark.shape[0], dark.shape[1])\n", + "flats_stack = flat.reshape(1, flat.shape[0], flat.shape[1])\n", + "alignment_stack = alignment.reshape(1, alignment.shape[0], alignment.shape[1])\n", + "\n", + "assert darks_stack.ndim == 3\n", + "assert flats_stack.ndim == 3\n", + "assert alignment_stack.ndim == 3\n", + "assert radios.ndim == 3\n", + "print(\"radios shape is\", radios.shape)\n", + "# create the array\n", + "data = numpy.concatenate([\n", + " darks_stack,\n", + " flats_stack,\n", + " radios[:200],\n", + " flats_stack,\n", + " radios[200:],\n", + " flats_stack,\n", + " alignment_stack,\n", + "])\n", + "assert data.ndim == 3\n", + "print(data.shape)\n", + "# then register the data to the detector\n", + "my_nxtomo.instrument.detector.data = data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### image key control" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "flats are at indexes (array([ 1, 202, 403]),)\n" + ] + } + ], + "source": [ + "from nxtomomill.utils import ImageKey\n", + "\n", + "image_key_control = numpy.concatenate([\n", + " [ImageKey.DARK_FIELD] * 1,\n", + " [ImageKey.FLAT_FIELD] * 1,\n", + " [ImageKey.PROJECTION] * 200,\n", + " [ImageKey.FLAT_FIELD] * 1,\n", + " [ImageKey.PROJECTION] * 200,\n", + " [ImageKey.FLAT_FIELD] * 1,\n", + " [ImageKey.ALIGNMENT] * 1,\n", + "])\n", + "\n", + "# insure with have the same number of frames and image key\n", + "assert len(image_key_control) == len(data)\n", + "# print position of flats in the sequence\n", + "print(\"flats are at indexes\", numpy.where(image_key_control == ImageKey.FLAT_FIELD))\n", + "# then register the image keys to the detector\n", + "my_nxtomo.instrument.detector.image_key_control = image_key_control" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### rotation angle" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "rotation_angle = numpy.concatenate([\n", + " [0.0, ], \n", + " [0.0, ], \n", + " proj_rotation_angles[:200],\n", + " [90.0, ],\n", + " proj_rotation_angles[200:],\n", + " [180.0, ],\n", + " [90.0, ],\n", + "])\n", + "assert len(rotation_angle) == len(data)\n", + "# register rotation angle to the sample\n", + "my_nxtomo.sample.rotation_angle = rotation_angle" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Field of view\n", + "field of view can either be `Half` or `Full`" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.instrument.detector.field_of_view = \"Full\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### pixel size" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.instrument.detector.x_pixel_size = my_nxtomo.instrument.detector.y_pixel_size = 1e-7 # pixel size must be provided in SI: meter" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "but for attribute with a unit you can specify the unit the value should be \"converted to\" using the 'unit' attribute like:" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.instrument.detector.x_pixel_size = my_nxtomo.instrument.detector.y_pixel_size = 0.1\n", + "my_nxtomo.instrument.detector.x_pixel_size.unit = my_nxtomo.instrument.detector.y_pixel_size.unit = \"micrometer\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When the unit is provided it will be stored as a property of the dataset. It must be interpreted by the software reading the NXtomo." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### save the nx to disk" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "nx_tomo_file_path = os.path.join(\"resources_cfg_create_nxtomo_from_scratch\", \"nxtomo.nx\")\n", + "my_nxtomo.save(file_path=nx_tomo_file_path, overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### check data saved\n", + "We can use some validator to insure we have enought data to be treated by nabu" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "from tomoscan.esrf import HDF5TomoScan\n", + "from tomoscan.validator import ReconstructionValidator\n", + "scan = HDF5TomoScan(nx_tomo_file_path, entry=\"entry\")\n", + "validator = ReconstructionValidator(scan, check_phase_retrieval=False, check_values=True)\n", + "assert validator.is_valid()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can check the layout of the file to insure it seems valid as well" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from h5glance import H5Glance\n", + "H5Glance(nx_tomo_file_path)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A good pratice is also to check frames, image_key and rotation angles to insure values seems valid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ! silx view resources_cfg_create_nxtomo_from_scratch/nxtomo.nx" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### reconstruct using nabu\n", + "now that we have a valid nxtomo we are able to reconstruct it using [nabu](https://gitlab.esrf.fr/tomotools/nabu).\n", + "\n", + "We create a nabu configuration file for contrast tomography reconstruction name `nabu-ct.conf` in order to reconstruct one slice of the volume.\n", + "\n", + "*note: on the configuration you must disable take_logarithm due to the dataset.* \n", + "\n", + "if nabu is installed you can run it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ! nabu nabu-cf.conf" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### provide mandatory data for Phase Tomography\n", + "in order to compute Phase Tomography you must also register:\n", + "* incoming beam energy (in keV)\n", + "* sample / detector distance (in meter)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "we can take back the existing `my_nxtomo` and add it the missing elements" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.energy = 12.5 # in keV by default\n", + "my_nxtomo.sample.distance = 0.2 # in meter" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And then you can reconstruct it with phase retrieval from modifing the nabu configuration file." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### provide more metadata" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "you can also provide x, y and z translation of the sample during the acquisition." + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.instrument.detector.x_translation = [0, 12]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "as a sample name, source information, start and end time" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.sample.name = \"my sample\"" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import datetime\n", + "my_nxtomo.instrument.source.name = \"ESRF\" # default value\n", + "my_nxtomo.instrument.source.type = \"Synchrotron X-ray Source\" # default value\n", + "my_nxtomo.start_time = datetime.now()\n", + "my_nxtomo.stop_time = datetime(2022, 2, 27)" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [], + "source": [ + "my_nxtomo.save(file_path=nx_tomo_file_path, overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# edit an NXTomo\n", + "\n", + "You can load an NXtomo from disk in order to edit it and save once done" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "nx_tomo type is \n", + "nx_tomo energy is 12.5 keV\n" + ] + } + ], + "source": [ + "nx_tomo_file_path = os.path.join(\"resources_cfg_create_nxtomo_from_scratch\", \"nxtomo.nx\")\n", + "nx_tomo = NXtomo(\"\").load(nx_tomo_file_path, \"entry\", detector_data_as=\"as_numpy_array\")\n", + "print(\"nx_tomo type is\", type(nx_tomo))\n", + "print(\"nx_tomo energy is\", nx_tomo.energy)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then you can modify your values as it was presented previously and overwrite the file." + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "new energy is 13.6 keV\n" + ] + } + ], + "source": [ + "my_nxtomo.energy = 13.6\n", + "my_nxtomo.save(file_path=nx_tomo_file_path, overwrite=True)\n", + "print(\"new energy is\", NXtomo(\"\").load(nx_tomo_file_path, \"entry\").energy)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Note:** the detector data is usually saved as a (h5py virtual dataset)[https://docs.h5py.org/en/stable/vds.html]. The amount of data assocaited can be heavy according to the acquisition and to the available memory. In order to allow a 'smooth' edition detector data can be load according to several strategies:\n", + "\n", + "* \"as_data_url\" (default): in this case each (Virtual Source)[https://docs.h5py.org/en/stable/vds.html#h5py.VirtualSource] will be saved as a DataUrl in order to ease it handling (see later on the tutorial)\n", + "* \"as_virtual_source\": retrieve original VirtualSource to allow edition of it\n", + "* \"as_numpy_array\": load all data in memory in order to modify it (and will dump the entire data). to avoid in case of \"real\" dataset. Can trigger huge IO." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### clean" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "if os.path.exists(nx_tomo_file_path):\n", + " os.remove(nx_tomo_file_path)\n", + "if os.path.exists(\"nxtomo_reconstruction.hdf5\"):\n", + " os.remove(\"nxtomo_reconstruction.hdf5\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Advance usage: Provide DataUrl to instrument.detector.data\n", + "The issue of using NXtomo we presented above is that the memory to handle data can be used (if you have a large number of projections and / or large detector).\n", + "\n", + "The way to work around for now it to use provide DataUrl (that can be pointing to an external file). Then this is the job to the FrameAppender to handle those.\n", + "\n", + "For this example we will first save metadata to the hdf5 (and maybe some frame) then you can append to the dataset series of frames sequentially with there rotation angle, image key..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### create some dataset to external files\n", + "Here we simply create some dataset on some external files and record the DataUrl\n", + "\n", + "**note:** **those datasets must be 3D otherwise virtual dataset creation will fail**" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "from nxtomomill.utils.frameappender import FrameAppender\n", + "from silx.io.url import DataUrl\n", + "import h5py\n", + "\n", + "detector_data_urls = []\n", + "for i_file in range(5):\n", + " os.makedirs(\"resources_cfg_create_nxtomo_from_scratch/external_files\", exist_ok=True)\n", + " external_file = os.path.join(f\"resources_cfg_create_nxtomo_from_scratch/external_files/file_{i_file}.nx\")\n", + " with h5py.File(external_file, mode=\"w\") as h5f:\n", + " h5f[\"data\"] = numpy.arange(\n", + " start=(5 * 100 * 100 * i_file),\n", + " stop=(5 * 100 * 100 * (i_file + 1))\n", + " ).reshape([5, 100, 100]) # of course here this is most likely that you will load data from another file\n", + "\n", + " detector_data_urls.append(\n", + " DataUrl(\n", + " file_path=external_file,\n", + " data_path=\"data\",\n", + " scheme=\"silx\",\n", + " )\n", + " )\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### create a simple nxtomo but this time provide the list of DataUrl to the instrument.detector.data attribute" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "my_large_nxtomo = NXtomo(\"entry0000\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "provide all information at the exception of frames. Here lets say we will have a dataset with only 180 projections" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "my_large_nxtomo.instrument.detector.distance = 0.2\n", + "my_large_nxtomo.instrument.detector.x_pixel_size = my_large_nxtomo.instrument.detector.y_pixel_size = 1e-7\n", + "my_large_nxtomo.energy = 12.3\n", + "# ...\n", + "my_large_nxtomo.sample.rotation_angle = numpy.linspace(0, 180, 180, endpoint=False)\n", + "my_large_nxtomo.instrument.detector.image_key_control = [0] * 180 # 0 == Projection" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "provide the list of DataUrl to `instrument.detector.data`" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "my_large_nxtomo.instrument.detector.data = detector_data_urls" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "my_large_nxtomo.save(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\", overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "note: this will create a virtual dataset under `instrument/detector/data` containing relative links from \"my_large_nxtomo.nx\" to other files." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then you can see that the 'data' dataset now contains 180 frames (if you run several time the previous cell then it will continue appending data to it).\n", + "\n", + "If an url is provided instead of a numpy array then it will create be used to create a virtual dataset and avoid duplicating data. But be carreful in this case you must keep relative position of the two files.\n", + "\n", + "**append frames must have the same dimensions otherwise the operation will fail**" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "H5Glance(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "check path of VirtualSources are relative (must start with './' string):" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dataset is virtual: True\n", + "file name is ./external_files/file_0.nx\n", + "file name is ./external_files/file_1.nx\n", + "file name is ./external_files/file_2.nx\n", + "file name is ./external_files/file_3.nx\n", + "file name is ./external_files/file_4.nx\n" + ] + } + ], + "source": [ + "with h5py.File(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\", mode=\"r\") as h5f:\n", + " dataset = h5f[\"entry0000/instrument/detector/data\"]\n", + " print(\"dataset is virtual:\", dataset.is_virtual)\n", + " for vs_info in dataset.virtual_sources():\n", + " print(\"file name is\", vs_info.file_name)\n", + " assert vs_info.file_name.startswith(\"./\")\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "clean" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [], + "source": [ + "import shutil\n", + "if os.path.exists(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\"):\n", + " os.remove(\"resources_cfg_create_nxtomo_from_scratch/my_large_nxtomo.nx\")\n", + "if os.path.exists(\"resources_cfg_create_nxtomo_from_scratch\"):\n", + " shutil.rmtree(\"resources_cfg_create_nxtomo_from_scratch\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**notes**:\n", + " * you can also provide a list of `h5py.VirtualSource` to the `detector.data` attribute.\n", + " * To append frame to an existing dataset you can also use the [`FrameAppender`](https://tomotools.gitlab-pages.esrf.fr/nxtomomill/_generated/nxtomomill.utils.frameappender.html) class directly*" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/nxtomomill/app/h52nx.py b/nxtomomill/app/h52nx.py index 7a21cc3f186b485fd6aa8ac9f60367030628a9cd..bac06ca6c83e1cc563ebd7bcf18d0efa4d792ec8 100644 --- a/nxtomomill/app/h52nx.py +++ b/nxtomomill/app/h52nx.py @@ -296,6 +296,12 @@ def main(argv): default=None, help="Titles corresponding to zserie init scans", ) + parser.add_argument( + "--init_pcotomo_titles", + "--init-pcotomo-titles", + default=None, + help="Titles corresponding to pcotomo init scans", + ) parser.add_argument( "--dark_titles", "--dark-titles", diff --git a/nxtomomill/converter/dxfile/dxfileconverter.py b/nxtomomill/converter/dxfile/dxfileconverter.py index ff743cb5d112427909d0f1d2e1fb300e2dce418f..bd960c2028776b3f8d4b4783b2762b654a32ef56 100644 --- a/nxtomomill/converter/dxfile/dxfileconverter.py +++ b/nxtomomill/converter/dxfile/dxfileconverter.py @@ -41,8 +41,8 @@ from silx.io.url import DataUrl from nxtomomill.utils import ImageKey from nxtomomill.converter.baseconverter import BaseConverter from silx.io.utils import h5py_read_dataset -from nxtomomill.converter.hdf5.acquisition.baseacquisition import EntryReader -from nxtomomill.converter.hdf5.acquisition.baseacquisition import DatasetReader +from nxtomomill.utils.hdf5 import EntryReader +from nxtomomill.utils.hdf5 import DatasetReader from nxtomomill.converter.version import version as converter_version from nxtomomill.nexus.nxdetector import FieldOfView from nxtomomill.utils import FrameAppender diff --git a/nxtomomill/converter/edf/edfconverter.py b/nxtomomill/converter/edf/edfconverter.py index 35af8c047d4e3136473a57c66e4149e667babe4a..58bfc8a29a499aeaf657e2fb01e6ea22d8f5f708 100644 --- a/nxtomomill/converter/edf/edfconverter.py +++ b/nxtomomill/converter/edf/edfconverter.py @@ -40,8 +40,8 @@ from nxtomomill.nexus.nxsource import SourceType from nxtomomill import utils from nxtomomill.utils import ImageKey from nxtomomill.converter.version import version as converter_version -from nxtomomill.converter.utils import create_nx_data_group -from nxtomomill.converter.utils import link_nxbeam_to_root +from nxtomomill.nexus.utils import create_nx_data_group +from nxtomomill.nexus.utils import link_nxbeam_to_root from nxtomomill.settings import Tomo try: diff --git a/nxtomomill/converter/hdf5/acquisition/baseacquisition.py b/nxtomomill/converter/hdf5/acquisition/baseacquisition.py index b46ad2cd4165610783cc714f206e21e7ac9bb9d0..10e1bd5d790caac20d44509979679f5586664109 100644 --- a/nxtomomill/converter/hdf5/acquisition/baseacquisition.py +++ b/nxtomomill/converter/hdf5/acquisition/baseacquisition.py @@ -33,7 +33,12 @@ __authors__ = [ __license__ = "MIT" __date__ = "27/11/2020" +import os import h5py +from nxtomomill.nexus.nxtomo import NXtomo +from nxtomomill.utils.hdf5 import EntryReader +from nxtomomill.utils.utils import embed_url, get_file_name # noqa F401 +from nxtomomill.utils.h5pyutils import from_data_url_to_virtual_source try: import hdf5plugin # noqa F401 @@ -47,62 +52,18 @@ from silx.io.url import DataUrl from nxtomomill.io.config import TomoHDF5Config from silx.io.utils import h5py_read_dataset from tomoscan.unitsystem import metricsystem -from tomoscan.io import HDF5File from collections import OrderedDict import numpy import typing import logging -import contextlib _logger = logging.getLogger(__name__) -class _BaseReader(contextlib.AbstractContextManager): - def __init__(self, url: DataUrl): - if not isinstance(url, DataUrl): - raise TypeError("url should be an instance of DataUrl") - if url.scheme() not in ("silx", "h5py"): - raise ValueError("Valid scheme are silx and h5py") - if url.data_slice() is not None: - raise ValueError( - "Data slices are not managed. Data path should " - "point to a bliss node (h5py.Group)" - ) - self._url = url - self._file_handler = None - - def __exit__(self, *exc): - return self._file_handler.close() - - -class EntryReader(_BaseReader): - """Context manager used to read a bliss node""" - - def __enter__(self): - self._file_handler = HDF5File(filename=self._url.file_path(), mode="r") - if self._url.data_path() == "": - entry = self._file_handler - else: - entry = self._file_handler[self._url.data_path()] - if not isinstance(entry, h5py.Group): - raise ValueError("Data path should point to a bliss node (h5py.Group)") - return entry - - -class DatasetReader(_BaseReader): - """Context manager used to read a bliss node""" - - def __enter__(self): - self._file_handler = HDF5File(filename=self._url.file_path(), mode="r") - entry = self._file_handler[self._url.data_path()] - if not isinstance(entry, h5py.Dataset): - raise ValueError( - "Data path ({}) should point to a dataset (h5py.Dataset)".format( - self._url.path() - ) - ) - return entry +def _ask_for_file_removal(file_path): + res = input("Overwrite %s ? (Y/n)" % file_path) + return res == "Y" class BaseAcquisition: @@ -127,6 +88,10 @@ class BaseAcquisition: _FOV_PATH = "technique/scan/field_of_view" + _NB_LOOP_PATH = "technique/scan/nb_loop" + + _NB_TOMO_PATH = "technique/scan/nb_tomo" + _TOMO_N_PATH = "technique/scan/tomo_n" _START_TIME_PATH = "start_time" @@ -166,14 +131,115 @@ class BaseAcquisition: def write_as_nxtomo( self, - output_file: str, - data_path: str, + shift_entry: int, input_file_path: str, request_input: bool, plugins, input_callback=None, - ) -> None: - raise NotImplementedError("Abstract class") + ) -> tuple: + self.set_plugins(plugins) + nx_tomos = self.to_NXtomos( + request_input=request_input, + input_callback=input_callback, + check_tomo_n=True, + ) + + # preprocessing to define output file name + possible_extensions = (".hdf5", ".h5", ".nx", ".nexus") + output_file_basename = os.path.basename(self.configuration.output_file) + file_extension_ = None + for possible_extension in possible_extensions: + if output_file_basename.endswith(possible_extension): + output_file_basename.rstrip(possible_extension) + file_extension_ = possible_extension + + def get_file_name_and_entry(index): + entry = "entry" + str(index).zfill(4) + + if self.configuration.single_file: + en_output_file = self.configuration.output_file + else: + ext = file_extension_ or self.configuration.file_extension + file_name = ( + os.path.splitext(output_file_basename)[0] + + "_" + + str(index).zfill(4) + + ext + ) + en_output_file = os.path.join( + os.path.dirname(self.configuration.output_file), file_name + ) + + if os.path.exists(en_output_file): + if self.configuration.overwrite is True: + _logger.warning(en_output_file + " will be removed") + _logger.info("remove " + en_output_file) + os.remove(en_output_file) + elif _ask_for_file_removal(en_output_file) is False: + _logger.info( + "unable to overwrite {}, exit".format(en_output_file) + ) + exit(0) + else: + os.remove(en_output_file) + return en_output_file, entry + + result = [] + for i_nx_tomo, nx_tomo in enumerate(nx_tomos): + output_file, data_path = get_file_name_and_entry((shift_entry + i_nx_tomo)) + output_file = os.path.abspath(os.path.relpath(output_file, os.getcwd())) + output_file = os.path.realpath(output_file) + + # For pcotomo for example the data path is modified in order to handle with the splitting + # that does not fit at 100 % with the current API. for now there is + # no convenient way to handle this in a better way + assert isinstance(nx_tomo, NXtomo) + self._update_nxtomo_from_plugins(nx_tomo=nx_tomo) + # embed data if requested + vs_0 = nx_tomo.instrument.detector.data[0] + if nx_tomo.detector_data_is_defined_by_url(): + new_urls = [] + for url in nx_tomo.instrument.detector.data: + if self._copy_frames[url.path()]: + created_url = embed_url( + url=url, + output_file=output_file, + ) + new_urls.append(created_url) + else: + new_urls.append(url) + nx_tomo.instrument.detector.data = new_urls + elif nx_tomo.detector_data_is_defined_by_virtual_source(): + new_vs = [] + for vs in nx_tomo.instrument.detector.data: + assert isinstance(vs, h5py.VirtualSource) + assert isinstance(vs_0, h5py.VirtualSource) + url = DataUrl(file_path=vs.path, data_path=vs.name, scheme="silx") + if self._copy_frames[url.path()]: + new_url = embed_url(url, output_file=output_file) + n_vs, _, _ = from_data_url_to_virtual_source(new_url) + new_vs.append(n_vs) + else: + new_vs.append(vs) + nx_tomo.instrument.detector.data = new_vs + # save data + nx_tomo.save(file_path=output_file, data_path=data_path) + result.append((output_file, data_path)) + return tuple(result) + + def _update_nxtomo_from_plugins(self, nx_tomo: NXtomo): + if len(self._plugins) > 0: + _logger.warning( + f"{len(self._plugins)} plugin found. Might modify NXtomo from them" + ) + for plugin in self._plugins: + _logger.info(f"Apply plugin {plugin.name}") + plugin.update_nx_tomo( + nx_tomo=nx_tomo, + ) + + def to_NXtomos(self, request_input, input_callback) -> tuple: + raise NotImplementedError("Base class") @property def raise_error_if_issue(self): @@ -332,17 +398,19 @@ class BaseAcquisition: with EntryReader(self._root_url) as entry: if self._TECHNIQUE_MOTOR_PATH in entry: try: - motor_and_aliases = h5py_read_dataset( - entry[self._TECHNIQUE_MOTOR_PATH] + rotation_motor = get_dataset_name_from_motor( + motors=h5py_read_dataset( + numpy.asarray(entry[self._TECHNIQUE_MOTOR_PATH]) + ), + motor_name="rotation", ) - rotation_motor = motor_and_aliases[0] except Exception as e: _logger.error(e) else: return rotation_motor else: _logger.warning( - "{} does not exists in {}".format( + "{} unable to find rotation motor from {}".format( self._TECHNIQUE_MOTOR_PATH, self._root_url ) ) @@ -536,30 +604,6 @@ class BaseAcquisition: nodes=(self._get_instrument_node(root_node), root_node), ) - def _ignore_sample_output(self, output_dataset_name): - """Should will ignore management of this dataset. Can be the case - if managed by a plugin""" - for plugin in self._plugins: - if output_dataset_name in plugin.sample_datasets_created: - return True - return False - - def _ignore_detector_output(self, output_dataset_name): - """Should will ignore management of this dataset. Can be the case - if managed by a plugin""" - for plugin in self._plugins: - if output_dataset_name in plugin.detector_datasets_created: - return True - return False - - def _ignore_beam_output(self, output_dataset_name): - """Should will ignore management of this dataset. Can be the case - if managed by a plugin""" - for plugin in self._plugins: - if output_dataset_name in plugin.beam_datasets_created: - return True - return False - def set_plugins(self, plugins): """ @@ -587,3 +631,21 @@ class BaseAcquisition: return "NXTomo" else: return self.root_url.path() + + +def get_dataset_name_from_motor(motors, motor_name): + motors = numpy.asarray(motors) + indexes = numpy.where(motors == motor_name)[0] + if len(indexes) == 0: + return None + elif len(indexes) == 1: + index = indexes[0] + index_dataset_id = index + 1 + if index_dataset_id < len(motors): + return motors[index_dataset_id] + else: + raise ValueError( + f"{motor_name} found but unable to find dataset name from {motors}" + ) + else: + raise ValueError(f"More than one instance of {motor_name} as been found.") diff --git a/nxtomomill/converter/hdf5/acquisition/pcotomoacquisition.py b/nxtomomill/converter/hdf5/acquisition/pcotomoacquisition.py new file mode 100644 index 0000000000000000000000000000000000000000..92ce033e45cc99021b5aa3bd03bd41eebb4ff90d --- /dev/null +++ b/nxtomomill/converter/hdf5/acquisition/pcotomoacquisition.py @@ -0,0 +1,182 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2015-2020 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +""" +module to define a pcotomo acquistion +""" + +__authors__ = [ + "H. Payno", +] +__license__ = "MIT" +__date__ = "08/02/2021" + + +from silx.utils.proxy import docstring +from nxtomomill.converter.hdf5.acquisition.baseacquisition import EntryReader +from nxtomomill.converter.hdf5.acquisition.standardacquisition import ( + StandardAcquisition, +) +from nxtomomill.io.acquisitionstep import AcquisitionStep +from nxtomomill.io.config import TomoHDF5Config +from silx.io.url import DataUrl +from typing import Optional, Union +import logging +from nxtomomill.nexus.nxtomo import NXtomo +from nxtomomill.utils import ImageKey +from nxtomomill.utils.nxsplitter import _NXtomoDetectorDataSplitter + +_logger = logging.getLogger(__name__) + + +class PCOTomoAcquisition(StandardAcquisition): + """ + A PCOTomo acquisition is an acquisition that can look like that at bliss side: + + * scan "descritpion", title can be tomo:pcotomo for example + * Optional scan dark + * Optional scan flat + * scan "projections" + * Optional scan flat + + The scan "projections" contains several "tomo" with a parameter evolving with time like heat or pressure. + The idea is that we want to split those into several NXtomo. + Those NXtomo must duplicate dark and flat scans + For example if scan "projections" contains nb_loop = 2 and nb_loop = 3 we must create nb_loop*nb_loop == 6 NXtomo as output. + + Split of those can are done in postprocessing on the to_NXtomos function + """ + + def __init__( + self, + root_url: Union[DataUrl, None], + configuration: TomoHDF5Config, + detector_sel_callback, + ): + super().__init__( + root_url=root_url, + configuration=configuration, + detector_sel_callback=detector_sel_callback, + ) + self._nb_loop = None + self._nb_tomo = None + self._entries_to_split = set() # set of URL as str + + def _preprocess_registered_entries(self): + # those must be defined before calling super because the super will call then `_preprocess_registered_entry` + self._nb_loop = None + self._nb_tomo = None + self._entries_to_split = set() + super()._preprocess_registered_entries() + + def _preprocess_registered_entry(self, entry_url, type_): + super()._preprocess_registered_entry(entry_url=entry_url, type_=type_) + if type_ is AcquisitionStep.PROJECTION: + # nb loop parameter must be present only on projection entries + nb_loop = self._get_nb_loop(entry_url) + if ( + nb_loop is not None + ): # at this moment 02/2022 nb_loop is only defined on projection type + if self._nb_loop is None or self._nb_loop == nb_loop: + self._nb_loop = nb_loop + self._entries_to_split.add(entry_url.path()) + else: + _logger.error( + f"Found entries with a different number of nb_loop: {entry_url.path()}" + ) + nb_tomo = self._get_nb_tomo(entry_url) + if ( + nb_tomo is not None + ): # at this moment 02/2022 nb_loop is only defined on projection type + if self._nb_tomo is None or self._nb_tomo == nb_tomo: + self._nb_tomo = nb_tomo + self._entries_to_split.add(entry_url.path()) + else: + _logger.error( + f"Found entries with a different number of _nb_tomo: {entry_url.path()}" + ) + + def _get_nb_loop(self, url) -> Optional[int]: + with EntryReader(url) as entry: + if self._NB_LOOP_PATH in entry: + return entry[self._NB_LOOP_PATH][()] + return None + + def _get_nb_tomo(self, url) -> Optional[int]: + with EntryReader(url) as entry: + if self._NB_LOOP_PATH in entry: + return entry[self._NB_TOMO_PATH][()] + return None + + @docstring(StandardAcquisition) + def to_NXtomos(self, request_input, input_callback, check_tomo_n: bool) -> tuple: + nx_tomos = super().to_NXtomos(request_input, input_callback, check_tomo_n=False) + results = [] + + for nx_tomo in nx_tomos: + splitter = _NXtomoDetectorDataSplitter(nx_tomo) + projections_slices = self._get_projections_slices(nx_tomo) + if len(projections_slices) > 1: + # insure projections are contiguous otherwise we don't know how to split it. + # not defined on the current design from bliss. should never happen + raise ValueError("Expect all projections to be contiguous") + elif len(projections_slices) == 0: + raise ValueError("No projection found") + else: + results.extend( + splitter.split( + projections_slices[0], + nb_part=(int(self._nb_loop * self._nb_tomo)), + ) + ) + + if check_tomo_n: + self.check_tomo_n() + + return tuple(results) + + @staticmethod + def _get_projections_slices(nx_tomo: NXtomo) -> tuple: + """Return a tuple of slices for each group of contiguous projections""" + if nx_tomo.instrument.detector.image_key_control is None: + return () + + res = [] + start_pos = -1 + browsing_projection = False + for i_frame, image_key in enumerate( + nx_tomo.instrument.detector.image_key_control + ): + image_key_value = ImageKey.from_value(image_key) + if image_key_value is ImageKey.PROJECTION and not browsing_projection: + browsing_projection = True + start_pos = i_frame + elif browsing_projection and image_key_value is not ImageKey.PROJECTION: + res.append(slice(start_pos, i_frame, 1)) + start_pos = -1 + browsing_projection = False + else: + if browsing_projection is True: + res.append(slice(start_pos, i_frame + 1, 1)) + return tuple(res) diff --git a/nxtomomill/converter/hdf5/acquisition/standardacquisition.py b/nxtomomill/converter/hdf5/acquisition/standardacquisition.py index 1262d584fafc11e281827e63b35d9f9558bd032a..3a63a50efd44daf80c7a5507cc74db28b3a7925d 100644 --- a/nxtomomill/converter/hdf5/acquisition/standardacquisition.py +++ b/nxtomomill/converter/hdf5/acquisition/standardacquisition.py @@ -31,25 +31,22 @@ __authors__ = [ "H. Payno", ] __license__ = "MIT" -__date__ = "27/11/2020" +__date__ = "14/02/2022" from .baseacquisition import BaseAcquisition from nxtomomill.nexus.nxsource import SourceType from nxtomomill.io.acquisitionstep import AcquisitionStep from .baseacquisition import EntryReader -from .baseacquisition import DatasetReader from .utils import get_entry_type from .utils import guess_nx_detector from .utils import get_nx_detectors from nxtomomill.utils import ImageKey -from nxtomomill.converter.version import version as converter_version from tomoscan.unitsystem import metricsystem -from tomoscan.io import HDF5File from silx.io.utils import h5py_read_dataset import h5py from silx.io.url import DataUrl -from typing import Union +from typing import Optional, Union try: import hdf5plugin # noqa F401 @@ -60,8 +57,7 @@ import fnmatch import numpy import os from nxtomomill.io.config import TomoHDF5Config -from nxtomomill.converter.utils import create_nx_data_group -from nxtomomill.converter.utils import link_nxbeam_to_root +from nxtomomill.nexus.nxtomo import NXtomo _logger = logging.getLogger(__name__) @@ -69,11 +65,15 @@ _logger = logging.getLogger(__name__) class StandardAcquisition(BaseAcquisition): """ - A standard acquisition where all registered scan are connected to - group an NXTomo entry + Class to collect information from a bliss - hdf scan (see https://bliss.gitlab-pages.esrf.fr/fscan). + Once all data is collected a set of NXtomo will be created. + Those NXtomo instances will then be pass to the plugins to allow users to modify. + Then NXtomo instances will be saved to disk. :param DataUrl root_url: url of the acquisition. Can be None if this is the initialization entry + :param TomoHDF5Config configuration: configuration to use to collect raw data and generate outputs + :param Optional[Function] detector_sel_callback: possible callback to retrieve missing information """ def __init__( @@ -87,13 +87,8 @@ class StandardAcquisition(BaseAcquisition): configuration=configuration, detector_sel_callback=detector_sel_callback, ) - # variables set by the `_preprocess_frames` function - self._data = None - """frames as a virtual dataset""" - self._image_key = None - """list of image keys""" + self._nx_tomos = [NXtomo("/")] self._image_key_control = None - """list of image keys""" self._rotation_angle = None """list of rotation angles""" self._x_translation = None @@ -101,20 +96,16 @@ class StandardAcquisition(BaseAcquisition): self._y_translation = None """y_translation""" self._z_translation = None - """z_translation""" - self._n_frames = None - self._dim_1 = None - self._dim_2 = None - self._data_type = None + + # self._n_frames = None + # self._dim_1 = None + # self._dim_2 = None + # self._data_type = None self._virtual_sources = None self._acq_expo_time = None self._copied_dataset = {} "register dataset copied. Key if the original location as" "DataUrl.path. Value is the DataUrl it has been moved to" - self._current_scan_n_frame = None - - @property - def image_key(self): - return self._image_key + # self._current_scan_n_frame = None @property def image_key_control(self): @@ -179,7 +170,12 @@ class StandardAcquisition(BaseAcquisition): def is_different_sequence(self, entry): return True - def register_step(self, url: DataUrl, entry_type=None, copy_frames=False) -> None: + def register_step( + self, + url: DataUrl, + entry_type: Optional[AcquisitionStep] = None, + copy_frames=False, + ) -> None: """ :param DataUrl url: entry to be registered and contained in the @@ -243,7 +239,6 @@ class StandardAcquisition(BaseAcquisition): data_dataset: h5py.Dataset, data_name, frame_type, - output_file, entry, entry_path, camera_dataset_url, @@ -281,19 +276,16 @@ class StandardAcquisition(BaseAcquisition): # Note: for now there is no image_key on the master file # should be added later. image_key_control = frame_type.to_image_key_control() - image_key = frame_type.to_image_key() self._image_key_control.extend([image_key_control.value] * n_frame) - self._image_key.extend([image_key.value] * n_frame) - # create virtual source (getting ready for writing) - rel_input = os.path.relpath( - camera_dataset_url.file_path(), os.path.dirname(output_file) - ) + data_dataset_path = data_dataset.name.replace(entry.name, entry_path, 1) # replace data_dataset name by the original entry_path. # this is a workaround to use the dataset path on the # "treated file". Because .name if the name on the 'target' # file of the virtual dataset - v_source = h5py.VirtualSource(rel_input, data_dataset_path, data_dataset.shape) + v_source = h5py.VirtualSource( + camera_dataset_url.file_path(), data_dataset_path, data_dataset.shape + ) self._virtual_sources.append(v_source) self._virtual_sources_len.append(n_frame) return n_frame @@ -304,7 +296,6 @@ class StandardAcquisition(BaseAcquisition): entry, frame_type, input_file_path, - output_file, entry_path, entry_url, ): @@ -323,63 +314,38 @@ class StandardAcquisition(BaseAcquisition): file_path=entry_url.file_path(), data_path=data_name, scheme="silx" ) - if self._copy_frames[entry_url.path()]: - from nxtomomill.utils import embed_url - - created_url = embed_url( - url=camera_dataset_url, - output_file=output_file, - ) - self._copied_dataset[camera_dataset_url.path()] = created_url - with DatasetReader(created_url) as copied_data_dataset: - n_frame = self.__get_data_from_camera( - copied_data_dataset, - data_name=data_name, - frame_type=frame_type, - output_file=output_file, - entry=entry, - entry_path=created_url.data_path(), - camera_dataset_url=created_url, - ) - else: - n_frame = self.__get_data_from_camera( - data_dataset, - data_name=data_name, - frame_type=frame_type, - output_file=output_file, - entry=entry, - entry_path=entry_path, - camera_dataset_url=camera_dataset_url, - ) + n_frame = self.__get_data_from_camera( + data_dataset, + data_name=data_name, + frame_type=frame_type, + entry=entry, + entry_path=entry_path, + camera_dataset_url=camera_dataset_url, + ) + # save information if this url must be embed / copy or not. Will be used later at nxtomo side + self._copy_frames[camera_dataset_url.path()] = self._copy_frames[ + entry_url.path() + ] # store rotation - if not self._ignore_sample_output("rotation_angle"): - rots = self._get_rotation_angle(root_node=entry, n_frame=n_frame)[0] - self._rotation_angle.extend(rots) - else: - self._rotation_angle = None + rots = self._get_rotation_angle(root_node=entry, n_frame=n_frame)[0] + self._rotation_angle.extend(rots) - if self.require_x_translation and not self._ignore_sample_output( - "x_translation" - ): + if self.require_x_translation: self._x_translation.extend( self._get_x_translation(root_node=entry, n_frame=n_frame)[0] ) else: self._x_translation = None - if self.require_y_translation and not self._ignore_sample_output( - "y_translation" - ): + if self.require_y_translation: self._y_translation.extend( self._get_y_translation(root_node=entry, n_frame=n_frame)[0] ) else: self._y_translation = None - if self.require_z_translation and not self._ignore_sample_output( - "z_translation" - ): + if self.require_z_translation: self._z_translation.extend( self._get_z_translation(root_node=entry, n_frame=n_frame)[0] ) @@ -424,7 +390,7 @@ class StandardAcquisition(BaseAcquisition): return True return False - def _preprocess_registered_entry(self, entry_url, type_, output_file): + def _preprocess_registered_entry(self, entry_url, type_): with EntryReader(entry_url) as entry: entry_path = self._entries_o_path[entry_url.path()] input_file_path = entry_url.file_path() @@ -471,6 +437,7 @@ class StandardAcquisition(BaseAcquisition): pass else: self._diode.extend(diode_vals) + self._diode_unit = diode_unit if not self.camera_is_valid(key): _logger.debug(f"ignore {key}, not a `valid` camera name") @@ -483,12 +450,11 @@ class StandardAcquisition(BaseAcquisition): entry=entry, frame_type=type_, input_file_path=input_file_path, - output_file=output_file, entry_path=entry_path, entry_url=entry_url, ) - def _preprocess_registered_entries(self, output_file): + def _preprocess_registered_entries(self): """parse all frames of the different steps and retrieve data, image_key...""" self._n_frames = 0 @@ -498,7 +464,6 @@ class StandardAcquisition(BaseAcquisition): self._x_translation = [] self._y_translation = [] self._z_translation = [] - self._image_key = [] self._image_key_control = [] self._rotation_angle = [] self._virtual_sources = [] @@ -519,7 +484,7 @@ class StandardAcquisition(BaseAcquisition): # list of data virtual source for the virtual dataset for entry_url, type_ in self._registered_entries.items(): url = DataUrl(path=entry_url) - self._preprocess_registered_entry(url, type_, output_file=output_file) + self._preprocess_registered_entry(url, type_) if len(self._diode) == 0: self._diode = None @@ -560,203 +525,6 @@ class StandardAcquisition(BaseAcquisition): res[alias] = self.configuration.param_already_defined[alias] return res - def _write_beam(self, root_node, request_input, input_callback): - instrument_node = root_node.require_group("instrument") - beam_node = instrument_node.require_group("beam") - - already_defined = self.get_already_defined_params( - key=TomoHDF5Config.EXTRA_PARAMS_ENERGY_DK - ) - if already_defined is not None: - energy = float(already_defined) - unit = "kev" - else: - energy, unit = self._get_energy( - ask_if_0=request_input, input_callback=input_callback - ) - if energy is not None: - beam_node["incident_energy"] = energy - beam_node["incident_energy"].attrs["unit"] = unit - - if "NX_class" not in beam_node.attrs: - beam_node.attrs["NX_class"] = "NXbeam" - - def _write_instrument(self, root_node): - # handle instrument - instrument_node = root_node.require_group("instrument") - instrument_node.attrs["NX_class"] = "NXinstrument" - instrument_title = self._get_instrument_name() - if instrument_title: - instrument_node["name"] = instrument_title - - detector_node = instrument_node.require_group("detector") - detector_node.attrs["NX_class"] = "NXdetector" - # write data - if self._virtual_sources is not None: - self._create_data_virtual_dataset(detector_node) - if self.image_key is not None: - detector_node["image_key"] = self.image_key - if self.image_key_control is not None: - detector_node["image_key_control"] = self.image_key_control - if self._acq_expo_time is not None: - detector_node["count_time"] = self._acq_expo_time - detector_node["count_time"].attrs["unit"] = "s" - # write distance - already_defined = self.get_already_defined_params( - key=TomoHDF5Config.EXTRA_PARAMS_DISTANCE - ) - if already_defined is not None: - distance = float(already_defined) - unit = "m" - else: - distance, unit = self._get_distance() - if distance is not None: - detector_node["distance"] = distance - detector_node["distance"].attrs["unit"] = unit - # write x and y pixel size - # if magnified pixel size is found then we right this value. - # otherwise will take pixel size (if found) - already_defined = self.get_already_defined_params( - key=TomoHDF5Config.EXTRA_PARAMS_X_PIXEL_SIZE_DK - ) - if already_defined is not None: - x_pixel_size = float(already_defined) - unit = "m" - else: - x_pixel_size, unit = self._get_pixel_size("x") - if x_pixel_size is not None: - detector_node["x_pixel_size"] = x_pixel_size - detector_node["x_pixel_size"].attrs["unit"] = unit - - already_defined = self.get_already_defined_params( - key=TomoHDF5Config.EXTRA_PARAMS_Y_PIXEL_SIZE_DK - ) - if already_defined is not None: - y_pixel_size = float(already_defined) - unit = "m" - else: - y_pixel_size, unit = self._get_pixel_size("y") - if y_pixel_size is not None: - detector_node["y_pixel_size"] = y_pixel_size - detector_node["y_pixel_size"].attrs["unit"] = unit - # write field of view - fov = self._get_field_of_fiew() - if fov is not None: - detector_node["field_of_view"] = fov - if fov.lower() == "half": - estimated_cor, unit = self._get_estimated_cor_from_motor( - pixel_size=y_pixel_size - ) - if estimated_cor is not None: - detector_node["estimated_cor_from_motor"] = estimated_cor - detector_node["estimated_cor_from_motor"].attrs["unit"] = unit - # write tomo_n - tomo_n = self._get_tomo_n() - if tomo_n is not None: - # save tomo n - detector_node["tomo_n"] = tomo_n - # write nx_source - source_grp = instrument_node.require_group("source") - source_name = self._get_source_name() - if source_name: - source_grp["name"] = source_name - source_type = self._get_source_type() - if source_type: - # as the value is not saved as a standard value at esrf adapt the value to a standard one - if "synchrotron" in source_type.lower(): - source_type = SourceType.SYNCHROTRON_X_RAY_SOURCE.value - # drop a warning if the source type is invalid - if source_type not in SourceType.values(): - _logger.warning( - "Source type ({}) is not a 'standard value'".format(source_type) - ) - source_grp["type"] = source_type - - if "NX_class" not in source_grp.attrs: - source_grp.attrs["NX_class"] = "NXsource" - - def _create_data_virtual_dataset(self, detector_node): - if ( - self.n_frames is None - or self.dim_1 is None - or self.dim_2 is None - or self.data_type is None - ): - if self.n_frames is None: - _logger.error("unable to get the number of frames") - if self.dim_1 is None: - _logger.error("unable to get frame dim_1") - if self.dim_2 is None: - _logger.error("unable to get frame dim_2") - if self.data_type is None: - _logger.error("unable to get data type") - raise ValueError( - "Preprocessing could not deduce all information " - "for creating the `data` virtual dataset" - ) - layout = h5py.VirtualLayout( - shape=(self.n_frames, self.dim_2, self.dim_1), dtype=self.data_type - ) - last = 0 - for v_source, vs_len in zip(self._virtual_sources, self._virtual_sources_len): - layout[last : vs_len + last] = v_source - last += vs_len - - detector_node.create_virtual_dataset("data", layout) - detector_node["data"].attrs["NX_class"] = "NXdata" - detector_node["data"].attrs[ - "SILX_style/axis_scale_types" - ] = self.get_axis_scale_types() - detector_node["data"].attrs["interpretation"] = "image" - - def _check_has_metadata(self): - if self.root_url is None: - raise ValueError( - "no initialization entry specify, unable to" "retrieve energy" - ) - - def _write_sample(self, root_node): - sample_node = root_node.create_group("sample") - sample_node.attrs["NX_class"] = "NXsample" - sample_name = self._get_sample_name() - if sample_name: - sample_node["name"] = sample_name - group_size = self._get_grp_size() - if group_size: - sample_node.attrs["group_size"] = group_size - if self.rotation_angle is not None: - sample_node["rotation_angle"] = self.rotation_angle - sample_node["rotation_angle"].attrs["unit"] = "degree" - if self.require_x_translation and self.x_translation is not None: - sample_node["x_translation"] = self.x_translation - sample_node["x_translation"].attrs["unit"] = "m" - if self.require_y_translation and self.y_translation is not None: - sample_node["y_translation"] = self.y_translation - sample_node["y_translation"].attrs["unit"] = "m" - if self.require_z_translation and self.z_translation is not None: - sample_node["z_translation"] = self.z_translation - sample_node["z_translation"].attrs["unit"] = "m" - - def _write_diode(self, root_node): - assert self.has_diode, "this acquisition does not expect diode" - diode_node = root_node.create_group("diode") - diode_node.attrs["NX_class"] = "NXdetector" - diode_node["data"] = self._diode - if self._diode_unit is not None: - diode_node["data"].attrs["unit"] = self._diode_unit - - def _write_plugins_output(self, root_node): - for plugin in self._plugins: - instrument_node = root_node["instrument"] - detector_node = instrument_node["detector"] - detector_node.attrs["NX_class"] = "NXdetector" - plugin.write( - root_node=root_node, - sample_node=root_node["sample"], - detector_node=detector_node, - beam_node=root_node["instrument/beam"], - ) - def _generic_path_getter(self, path, message, level="warning"): """ :param str level: level can be logging.level values : "warning", "error", "info" @@ -981,86 +749,117 @@ class StandardAcquisition(BaseAcquisition): else: return y_rot / pixel_size, "pixels" - def write_as_nxtomo( - self, - output_file: str, - data_path: str, - input_file_path: str, - request_input: bool, - plugins, - input_callback=None, - ) -> None: - """ - write the current sequence in an NXtomo like - - :param str output_file: destination file - :param str data_path: path to store the data in the destination file - :param str input_file_path: hdf5 source file - :param bool request_input: if some entries are missing should we ask - the user for input - :param list plugins: plugins to process - :param input_callback: if provided then will call this callback - function with (missing_entry, desc) instead of - input - """ - _logger.info( - "write data of {} to {}".format(self, output_file + "::/" + data_path) + def to_NXtomos(self, request_input, input_callback, check_tomo_n=True) -> tuple: + self._preprocess_registered_entries() + + nx_tomo = NXtomo("/") + + # 1. root level information + # start and end time + nx_tomo.start_time = self._get_start_time() + nx_tomo.end_time = self._get_end_time() + + # title + nx_tomo.title = self._get_dataset_name() + + # 2. define beam + energy, unit = self._get_user_settable_parameter( + param_key=TomoHDF5Config.EXTRA_PARAMS_ENERGY_DK, + fallback_fct=self._get_energy, + dtype=float, + input_callback=input_callback, + ask_if_0=request_input, ) - self.set_plugins(plugins) - - # work on absolute path. The conversion to relative path and - # then to absolute path is a trick in case there is some 'mounted' - # directory exposed differently. - # Like '/mnt/multipath-shares/tmp_14_days' - output_file = os.path.abspath(os.path.relpath(output_file, os.getcwd())) - output_file = os.path.realpath(output_file) - - # first retrieve the data and create some virtual dataset. - self._preprocess_registered_entries(output_file=output_file) - with HDF5File(output_file, "a") as h5_file: - entry = h5_file.require_group(data_path) - entry.attrs["NX_class"] = "NXentry" - entry.attrs["definition"] = "NXtomo" - entry.attrs["version"] = converter_version() - entry.attrs["default"] = "data" - start_time = self._get_start_time() - if start_time is not None: - entry["start_time"] = start_time - end_time = self._get_end_time() - if end_time is not None: - entry["end_time"] = end_time - dataset_name = self._get_dataset_name() - if dataset_name: - entry["title"] = dataset_name - entry["definition"] = "NXtomo" - - self._write_beam( - entry, request_input=request_input, input_callback=input_callback - ) - self._write_instrument(entry) - self._write_sample(entry) - self._write_plugins_output(entry) - if self.has_diode: - self._write_diode(entry) - - # create nxdata group - try: - create_nx_data_group( - file_path=output_file, - entry_path=data_path, - axis_scale=self.get_axis_scale_types(), - ) - except Exception as e: - if self.raise_error_if_issue: - raise e - else: - _logger.error( - "Fail to create NXdata group. Reason is {}".format(str(e)) + if energy is not None: + # TODO: better manamgent of energy ? might be energy.beam or energy.instrument.beam ? + nx_tomo.energy = energy + nx_tomo.energy.unit = unit + + # 3. define instrument + nx_tomo.instrument.detector.data = self._virtual_sources + nx_tomo.instrument.detector.image_key_control = self.image_key_control + nx_tomo.instrument.detector.count_time = self._acq_expo_time + nx_tomo.instrument.detector.count_time.unit = "s" + # distance + distance, unit = self._get_user_settable_parameter( + param_key=TomoHDF5Config.EXTRA_PARAMS_DISTANCE, + fallback_fct=self._get_distance, + dtype=float, + ) + if distance is not None: + nx_tomo.instrument.detector.distance = distance + if nx_tomo.instrument.detector.distance is not None: + nx_tomo.instrument.detector.distance.unit = unit + # x and y pixel size + x_pixel_size, unit = self._get_user_settable_parameter( + param_key=TomoHDF5Config.EXTRA_PARAMS_X_PIXEL_SIZE_DK, + fallback_fct=self._get_pixel_size, + dtype=float, + axis="x", + ) + nx_tomo.instrument.detector.x_pixel_size = x_pixel_size + if unit is not None: + nx_tomo.instrument.detector.x_pixel_size.unit = unit + + y_pixel_size, unit = self._get_user_settable_parameter( + param_key=TomoHDF5Config.EXTRA_PARAMS_Y_PIXEL_SIZE_DK, + fallback_fct=self._get_pixel_size, + dtype=float, + axis="y", + ) + nx_tomo.instrument.detector.y_pixel_size = y_pixel_size + if unit is not None: + nx_tomo.instrument.detector.y_pixel_size.unit = unit + + fov = self._get_field_of_fiew() + if fov is not None: + nx_tomo.instrument.field_of_view = fov + if fov.lower() == "half": + estimated_cor, unit = self._get_estimated_cor_from_motor( + pixel_size=y_pixel_size ) + if estimated_cor is not None: + nx_tomo.instrument.detector.estimated_cor_from_motor = estimated_cor + + # define tomo_n + nx_tomo.instrument.detector.tomo_n = self._get_tomo_n() - # create beam group at root for compatibility - link_nxbeam_to_root(file_path=output_file, entry_path=data_path) + # 3. define nx source + source_name = self._get_source_name() + nx_tomo.instrument.source.name = source_name + source_type = self._get_source_type() + if source_type is not None: + if "synchrotron" in source_type.lower(): + source_type = SourceType.SYNCHROTRON_X_RAY_SOURCE.value + # drop a warning if the source type is invalid + if source_type not in SourceType.values(): + _logger.warning( + "Source type ({}) is not a 'standard value'".format(source_type) + ) + nx_tomo.instrument.source.type = source_type + + # 4. define sample + nx_tomo.sample.name = self._get_sample_name() + nx_tomo.sample.group_size = self._get_grp_size() + nx_tomo.sample.rotation_angle = self.rotation_angle + nx_tomo.sample.x_translation.value = self.x_translation + nx_tomo.sample.x_translation.unit = "m" + nx_tomo.sample.y_translation.value = self.y_translation + nx_tomo.sample.y_translation.unit = "m" + nx_tomo.sample.z_translation.value = self.z_translation + nx_tomo.sample.z_translation.unit = "m" + + # 5. define diode + if self.has_diode: + nx_tomo.instrument.diode.data = self._diode + nx_tomo.instrument.diode.data.unit = self._diode_unit + + if check_tomo_n: + self.check_tomo_n() + return (nx_tomo,) + + def check_tomo_n(self): # check scan is complete tomo_n = self._get_tomo_n() if self.configuration.check_tomo_n and tomo_n is not None: @@ -1077,3 +876,31 @@ class StandardAcquisition(BaseAcquisition): raise ValueError(mess) else: _logger.error(mess) + + def _check_has_metadata(self): + if self.root_url is None: + raise ValueError( + "no initialization entry specify, unable to" "retrieve energy" + ) + + def _get_user_settable_parameter( + self, + param_key, + fallback_fct, + dtype: Optional[type] = None, + *fallback_args, + **fallback_kwargs, + ): + """ + return value, unit + """ + value = self.get_already_defined_params(param_key) + if value is not None: + unit = TomoHDF5Config.get_extra_params_default_unit(param_key) + else: + value, unit = fallback_fct(*fallback_args, **fallback_kwargs) + + if dtype is None or value is None: + return value, unit + else: + return dtype(value), unit diff --git a/nxtomomill/converter/hdf5/acquisition/test/test_acquisition.py b/nxtomomill/converter/hdf5/acquisition/test/test_acquisition.py index eb35e48cd705a8e8949473f51a9595e54720b70e..67f3964d6ca86488e6b55309ca973e3f933c87a0 100644 --- a/nxtomomill/converter/hdf5/acquisition/test/test_acquisition.py +++ b/nxtomomill/converter/hdf5/acquisition/test/test_acquisition.py @@ -27,7 +27,11 @@ __authors__ = ["H. Payno"] __license__ = "MIT" __date__ = "12/03/2021" -from nxtomomill.converter.hdf5.acquisition.baseacquisition import BaseAcquisition +import pytest +from nxtomomill.converter.hdf5.acquisition.baseacquisition import ( + BaseAcquisition, + get_dataset_name_from_motor, +) import os from silx.io.url import DataUrl from nxtomomill.io.config import TomoHDF5Config @@ -49,3 +53,18 @@ def test_BaseAquisition(): ) with std_acq.read_entry() as entry: assert "dataset" in entry + + +def test_get_dataset_name_from_motor(): + """test get_dataset_name_from_motor function""" + set_1 = ["rotation", "test1", "alias"] + assert get_dataset_name_from_motor(set_1, "rotation") == "test1" + assert get_dataset_name_from_motor(set_1, "my motor") is None + + set_2 = ["rotation", "test1", "alias", "x translation", "m2", "test1"] + assert get_dataset_name_from_motor(set_2, "rotation") == "test1" + assert get_dataset_name_from_motor(set_2, "x translation") == "m2" + + set_3 = ["rotation"] + with pytest.raises(ValueError): + get_dataset_name_from_motor(set_3, "rotation") diff --git a/nxtomomill/converter/hdf5/acquisition/test/test_pcotomoacquisition.py b/nxtomomill/converter/hdf5/acquisition/test/test_pcotomoacquisition.py new file mode 100644 index 0000000000000000000000000000000000000000..e51e8a98e14816a05d554131487a5fc3b957f64f --- /dev/null +++ b/nxtomomill/converter/hdf5/acquisition/test/test_pcotomoacquisition.py @@ -0,0 +1,126 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2017 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "08/02/2021" + + +import os +import pytest +from nxtomomill import converter +from nxtomomill.converter.hdf5.acquisition.pcotomoacquisition import PCOTomoAcquisition +from nxtomomill.io.config import TomoHDF5Config +from nxtomomill.nexus.nxtomo import NXtomo +from nxtomomill.test.utils.bliss import MockBlissAcquisition +from nxtomomill.utils import Format +import h5py +import numpy + +configs = ( + {"nb_tomo": 1, "nb_loop": 5}, + {"nb_tomo": 5, "nb_loop": 1}, + {"nb_tomo": 2, "nb_loop": 3}, +) + + +@pytest.mark.parametrize("config", configs) +def test_pcotomo_conversion(tmp_path, config): + n_darks = 10 + n_flats = 10 + bliss_scan_dir = str(tmp_path / "my_acquisition") + nb_tomo = config["nb_tomo"] + nb_loop = config["nb_loop"] + bliss_mock = MockBlissAcquisition( + n_sample=1, + n_sequence=1, + n_scan_per_sequence=1, + n_darks=n_darks, + n_flats=n_flats, + output_dir=bliss_scan_dir, + acqui_type="pcotomo", + with_rotation_motor_info=True, + nb_tomo=nb_tomo, + nb_loop=nb_loop, + ) + sample_file = bliss_mock.samples[0].sample_file + + configuration = TomoHDF5Config() + configuration.format = Format.STANDARD + + output_nx_file = os.path.join(str(tmp_path), "nexus_scans.nx") + assert not os.path.exists(output_nx_file) + configuration.output_file = output_nx_file + configuration.input_file = sample_file + configuration.request_input = False + configuration.raises_error = False + + result = converter.from_h5_to_nx(configuration=configuration) + assert os.path.exists( + output_nx_file + ) # TODO: must create 5 files with one master file + # check we have our 5 NXtomo converted + assert len(result) == nb_loop * nb_tomo + # check saved datasets + with h5py.File(sample_file, mode="r") as h5f: + darks = h5f["2.1/instrument/pcolinux/data"][()] + assert len(darks) == n_darks + flats_1 = h5f["3.1/instrument/pcolinux/data"][()] + assert len(flats_1) == n_flats + projections = h5f["4.1/instrument/pcolinux/data"][()] + assert len(projections) == 10 * nb_loop * nb_tomo + flats_2 = h5f["5.1/instrument/pcolinux/data"][()] + assert len(flats_2) == n_flats + + for i_nx_tomo, (file_path, data_path) in enumerate(result): + with h5py.File(file_path, mode="r") as h5f: + detector_path = "/".join([data_path, "instrument", "detector", "data"]) + detector_data = h5f[detector_path][()] + expected_data = numpy.concatenate( + [ + darks, + flats_1, + projections[(10 * i_nx_tomo) : (10 * (i_nx_tomo + 1))], + flats_2, + ] + ) + numpy.testing.assert_array_equal(detector_data, expected_data) + + +def test_get_projections_slices(): + """test the _get_projections_slices function""" + nx_tomo = NXtomo("test") + nx_tomo.instrument.detector.image_key_control = [2, 1, 0, 0, 0, 2] + assert PCOTomoAcquisition._get_projections_slices(nx_tomo) == (slice(2, 5, 1),) + nx_tomo.instrument.detector.image_key_control = [2, 1, 0, 0, 0, 2] + assert PCOTomoAcquisition._get_projections_slices(nx_tomo) == (slice(2, 5, 1),) + nx_tomo.instrument.detector.image_key_control = [2, 1] + assert PCOTomoAcquisition._get_projections_slices(nx_tomo) == () + nx_tomo.instrument.detector.image_key_control = [0, 0, 2, 1, 0, 0, 1, 0, 0] + assert PCOTomoAcquisition._get_projections_slices(nx_tomo) == ( + slice(0, 2, 1), + slice(4, 6, 1), + slice(7, 9, 1), + ) diff --git a/nxtomomill/converter/hdf5/acquisition/utils.py b/nxtomomill/converter/hdf5/acquisition/utils.py index 99e7417373a6a579f3245dc9831a10b2dbf43073..37e2b60dc305e73cfe2bdfd149accd2b00f2f94f 100644 --- a/nxtomomill/converter/hdf5/acquisition/utils.py +++ b/nxtomomill/converter/hdf5/acquisition/utils.py @@ -89,6 +89,7 @@ def get_entry_type( else: init_titles = list(configuration.init_titles) init_titles.extend(configuration.zserie_init_titles) + init_titles.extend(configuration.pcotomo_init_titles) step_titles = { AcquisitionStep.INITIALIZATION: init_titles, diff --git a/nxtomomill/converter/hdf5/acquisition/xrd3dacquisition.py b/nxtomomill/converter/hdf5/acquisition/xrd3dacquisition.py index 4a2d6eeb2441977409e58cc4eee5f2ace4a1bfe3..f63d959d6d31e1bb55be3bb84ec5c76772fd787d 100644 --- a/nxtomomill/converter/hdf5/acquisition/xrd3dacquisition.py +++ b/nxtomomill/converter/hdf5/acquisition/xrd3dacquisition.py @@ -131,26 +131,20 @@ class XRD3DAcquisition(StandardAcquisition): entry_url=entry_url, ) # store base tilt information - if not self._ignore_sample_output("base_tilt"): - base_tilt, _ = self._get_base_tilt_dataset( - entry=entry, n_frames=self._current_scan_n_frame - ) - self._base_tilt.extend(base_tilt) - else: - self._base_tilt = None + base_tilt, _ = self._get_base_tilt_dataset( + entry=entry, n_frames=self._current_scan_n_frame + ) + self._base_tilt.extend(base_tilt) # store rocking information - if not self._ignore_sample_output("rocking"): - rocking, _ = self._get_rocking_dataset( - entry=entry, n_frames=self._current_scan_n_frame - ) - self._rocking.extend(rocking) - else: - self._rocking = None + rocking, _ = self._get_rocking_dataset( + entry=entry, n_frames=self._current_scan_n_frame + ) + self._rocking.extend(rocking) - def _preprocess_registered_entries(self, output_file): + def _preprocess_registered_entries(self): self._base_tilt = [] self._rocking = [] - super()._preprocess_registered_entries(output_file=output_file) + super()._preprocess_registered_entries() def _write_sample(self, root_node): super()._write_sample(root_node) diff --git a/nxtomomill/converter/hdf5/acquisition/zseriesacquisition.py b/nxtomomill/converter/hdf5/acquisition/zseriesacquisition.py index fd8903755ff0b6a6fd3b531231a1938efd496ba5..2c5cb66e91e1ad0a11fca10a246e43511fabdc1a 100644 --- a/nxtomomill/converter/hdf5/acquisition/zseriesacquisition.py +++ b/nxtomomill/converter/hdf5/acquisition/zseriesacquisition.py @@ -65,6 +65,18 @@ def is_z_series_frm_titles(entry: h5py.Group, configuration: TomoHDF5Config) -> return title in configuration.zserie_init_titles +def is_pcotomo_frm_titles(entry: h5py.Group, configuration: TomoHDF5Config) -> bool: + """ + if the provided h5py.Group must be consider as an "initialization" entry/scan of a pcotomo acquistion + """ + try: + title = h5py_read_dataset(entry["title"]) + except Exception: + return False + else: + return title in configuration.pcotomo_init_titles + + def is_z_series_frm_z_translation( projection_urls: Iterable, configuration: TomoHDF5Config ): diff --git a/nxtomomill/converter/hdf5/hdf5converter.py b/nxtomomill/converter/hdf5/hdf5converter.py index 1b0b84f75f1ad4c67cc88eedf75c6052a8076dac..5fb75322799c0c2c2233f47deeb705e36bbae739 100644 --- a/nxtomomill/converter/hdf5/hdf5converter.py +++ b/nxtomomill/converter/hdf5/hdf5converter.py @@ -34,10 +34,14 @@ __license__ = "MIT" __date__ = "27/11/2020" +from nxtomomill.converter.hdf5.acquisition.pcotomoacquisition import PCOTomoAcquisition from .acquisition.utils import get_entry_type from nxtomomill.io.acquisitionstep import AcquisitionStep from .acquisition.standardacquisition import StandardAcquisition -from .acquisition.zseriesacquisition import ZSeriesBaseAcquisition +from .acquisition.zseriesacquisition import ( + ZSeriesBaseAcquisition, + is_pcotomo_frm_titles, +) from .acquisition.zseriesacquisition import is_z_series_frm_titles from .acquisition.zseriesacquisition import is_z_series_frm_z_translation from .acquisition.xrdctacquisition import XRDCTAcquisition @@ -47,6 +51,7 @@ from silx.utils.deprecation import deprecated from nxtomomill.io.config import TomoHDF5Config from nxtomomill.io.framegroup import FrameGroup from nxtomomill.converter.baseconverter import BaseConverter +from nxtomomill.converter.hdf5.acquisition.baseacquisition import _ask_for_file_removal from silx.io.url import DataUrl from tomoscan.io import HDF5File from typing import Union @@ -80,6 +85,7 @@ H5_X_PIXEL_SIZE = Tomo.H5.X_PIXEL_SIZE H5_Y_PIXEL_SIZE = Tomo.H5.Y_PIXEL_SIZE H5_DARK_TITLES = Tomo.H5.DARK_TITLES H5_INIT_TITLES = Tomo.H5.INIT_TITLES +H5_PCOTOMO_INIT_TITLES = Tomo.H5.PCOTOMO_INIT_TITLES H5_ZSERIE_INIT_TITLES = Tomo.H5.ZSERIE_INIT_TITLES H5_PROJ_TITLES = Tomo.H5.PROJ_TITLES H5_FLAT_TITLES = Tomo.H5.FLAT_TITLES @@ -90,6 +96,7 @@ H5_DIODE_KEYS = Tomo.H5.DIODE_KEYS DEFAULT_SCAN_TITLES = H5ScanTitles( H5_INIT_TITLES, H5_ZSERIE_INIT_TITLES, + H5_PCOTOMO_INIT_TITLES, H5_DARK_TITLES, H5_FLAT_TITLES, H5_PROJ_TITLES, @@ -114,11 +121,6 @@ DEFAULT_H5_KEYS = H5FileKeys( _logger = logging.getLogger(__name__) -def _ask_for_file_removal(file_path): - res = input("Overwrite %s ? (Y/n)" % file_path) - return res == "Y" - - class _H5ToNxConverter(BaseConverter): """ Class used to convert a HDF5Config to one or several NXTomoEntry. @@ -156,7 +158,6 @@ class _H5ToNxConverter(BaseConverter): "configuration should be an instance of HDFConfig" " not {}".format(type(configuration)) ) - self._configuration = configuration self._progress = progress self._input_callback = input_callback @@ -289,6 +290,12 @@ class _H5ToNxConverter(BaseConverter): configuration=self.configuration, detector_sel_callback=self.detector_sel_callback, ) + elif is_pcotomo_frm_titles(acqui_projs_urls, self.configuration): + root_acquisition = PCOTomoAcquisition( + root_url=frame_grp.url, + configuration=self.configuration, + detector_sel_callback=self.detector_sel_callback, + ) else: root_acquisition = StandardAcquisition( root_url=frame_grp.url, @@ -469,6 +476,14 @@ class _H5ToNxConverter(BaseConverter): configuration=self.configuration, detector_sel_callback=self.detector_sel_callback, ) + elif is_pcotomo_frm_titles( + entry=entry, configuration=self.configuration + ): + current_acquisition = PCOTomoAcquisition( + root_url=url, + configuration=self.configuration, + detector_sel_callback=self.detector_sel_callback, + ) else: current_acquisition = StandardAcquisition( root_url=url, @@ -528,13 +543,6 @@ class _H5ToNxConverter(BaseConverter): def write(self): res = [] - possible_extensions = (".hdf5", ".h5", ".nx", ".nexus") - output_file_basename = os.path.basename(self.configuration.output_file) - file_extension_ = None - for possible_extension in possible_extensions: - if output_file_basename.endswith(possible_extension): - output_file_basename.rstrip(possible_extension) - file_extension_ = possible_extension acq_str = [str(acq) for acq in self.acquisitions] acq_str.insert( @@ -552,62 +560,20 @@ class _H5ToNxConverter(BaseConverter): if self.progress is not None: self.progress.set_name("write sequences") self.progress.set_max_advancement(len(self.acquisitions)) - for i_acquisition, acquisition in enumerate(self.acquisitions): + shift_entry = 0 + for acquisition in self.acquisitions: if self._ignore_sub_entry(acquisition.root_url): continue - entry = "entry" + str(i_acquisition).zfill(4) - if self.configuration.single_file or len(self.acquisitions) == 1: - en_output_file = self.configuration.output_file - else: - ext = file_extension_ or self.configuration.file_extension - file_name = ( - os.path.splitext(output_file_basename)[0] - + "_" - + str(i_acquisition).zfill(4) - + ext - ) - en_output_file = os.path.join( - os.path.dirname(self.configuration.output_file), file_name - ) - if os.path.exists(en_output_file): - if self.configuration.overwrite is False: - _logger.warning(en_output_file + " will be removed") - _logger.info("remove " + en_output_file) - os.remove(en_output_file) - elif _ask_for_file_removal(en_output_file) is False: - _logger.info( - "unable to overwrite {}, exit".format(en_output_file) - ) - exit(0) - else: - os.remove(en_output_file) - # if acquisition.root_url is None: - # continue try: - acquisition.write_as_nxtomo( - output_file=en_output_file, - data_path=entry, + new_entries = acquisition.write_as_nxtomo( + shift_entry=shift_entry, input_file_path=self.configuration.input_file, request_input=self.configuration.request_input, input_callback=self.input_callback, plugins=self.plugins, ) - # if split files create a master file with link to those entries - if ( - self.configuration.single_file is False - and len(self.acquisitions) > 1 - ): - _logger.info("create link in %s" % self.configuration.output_file) - with HDF5File(self.configuration.output_file, "a") as master_file: - link_file = os.path.relpath( - en_output_file, - os.path.dirname(self.configuration.output_file), - ) - master_file[entry] = h5py.ExternalLink(link_file, entry) - res.append((en_output_file, entry)) - else: - res.append((en_output_file, entry)) + shift_entry += len(new_entries) except Exception as e: if self.configuration.raises_error: raise e @@ -616,9 +582,22 @@ class _H5ToNxConverter(BaseConverter): "Fails to write %s. Error is %s" % (str(acquisition.root_url), str(e)) ) + else: + res.extend(new_entries) if self.progress is not None: self.progress.increase_advancement() + # if we created one file per entry then create a master file with link to those entries + if self.configuration.single_file is False: + _logger.info("create link in %s" % self.configuration.output_file) + for (en_output_file, entry) in res: + with HDF5File(self.configuration.output_file, "a") as master_file: + link_file = os.path.relpath( + en_output_file, + os.path.dirname(self.configuration.output_file), + ) + master_file[entry] = h5py.ExternalLink(link_file, entry) + return tuple(res) def _check_conversion_is_possible(self): @@ -770,6 +749,7 @@ def h5_to_nx( # handle titles keys configuration.init_titles = scan_titles.init_titles configuration.zserie_init_titles = scan_titles.init_zserie_titles + configuration.pcotomo_init_titles = scan_titles.init_pcotomo_titles configuration.dark_titles = scan_titles.dark_titles configuration.flat_titles = scan_titles.flat_titles configuration.projections_titles = scan_titles.proj_titles diff --git a/nxtomomill/converter/hdf5/test/test_hdf5converter.py b/nxtomomill/converter/hdf5/test/test_hdf5converter.py index 17398fe470a67b5af608270f957fce1f3568266d..d722e026e33298cc38be4346b72c451fc1ed42de 100644 --- a/nxtomomill/converter/hdf5/test/test_hdf5converter.py +++ b/nxtomomill/converter/hdf5/test/test_hdf5converter.py @@ -37,8 +37,8 @@ import os from nxtomomill import converter from nxtomomill.converter.hdf5.acquisition.utils import get_nx_detectors from nxtomomill.converter.hdf5.acquisition.utils import guess_nx_detector -from nxtomomill.converter.hdf5.acquisition.baseacquisition import EntryReader -from nxtomomill.converter.hdf5.acquisition.baseacquisition import DatasetReader +from nxtomomill.utils.hdf5 import EntryReader +from nxtomomill.utils.hdf5 import DatasetReader from nxtomomill.io.config import TomoHDF5Config try: @@ -53,6 +53,7 @@ from glob import glob from nxtomomill.io.framegroup import FrameGroup from nxtomomill.utils import Format import subprocess +from nxtomomill.nexus.nxtomo import NXtomo def url_has_been_copied(file_path: str, url: DataUrl): @@ -101,6 +102,7 @@ class TestH5ToNxConverter(unittest.TestCase): self.assertTrue(os.path.exists(sample.sample_file)) self.config.output_file = sample.sample_file.replace(".h5", ".nx") self.config.input_file = sample.sample_file + self.config.raises_error = True assert ( len(converter.get_bliss_tomo_entries(sample.sample_file, self.config)) == 1 @@ -598,9 +600,15 @@ class TestStandardAcqConversionWithExternalUrls(unittest.TestCase): self.assertTrue(os.path.exists(self.config.output_file)) with h5py.File(self.config.output_file, mode="r") as h5s: - self.assertEqual(len(h5s.items()), 2) + self.assertTrue("entry0000" in h5s) + self.assertEqual(len(h5s.items()), 1) + + with h5py.File( + self.config.output_file.replace(".nx", "_0000.nx"), mode="r" + ) as h5s: self.assertTrue("entry0000" in h5s) self.assertTrue("duplicate_data" in h5s) + self.assertEqual(len(h5s.items()), 2) detector_url = DataUrl( file_path=self.config.output_file, @@ -612,9 +620,9 @@ class TestStandardAcqConversionWithExternalUrls(unittest.TestCase): for i_vs, vs in enumerate(detector_dataset.virtual_sources()): self.assertFalse(os.path.isabs(vs.file_name)) if i_vs in (0, 1, 4, 5): - self.assertEqual(vs.file_name, "output.nx") + self.assertEqual(vs.file_name, "./output_0000.nx") else: - self.assertEqual(vs.file_name, "acqui_1/sample_0/sample_0.h5") + self.assertEqual(vs.file_name, "./acqui_1/sample_0/sample_0.h5") scan = HDF5TomoScan(scan=self.config.output_file, entry="entry0000") self.assertTrue(is_valid_for_reconstruction(scan)) @@ -668,12 +676,18 @@ class TestStandardAcqConversionWithExternalUrls(unittest.TestCase): for url in urls_copied: self.assertTrue( - url_has_been_copied(file_path=self.config.output_file, url=url) + url_has_been_copied( + file_path=self.config.output_file.replace(".nx", "_0000.nx"), + url=url, + ) ) for url in urls_not_copied: self.assertFalse( - url_has_been_copied(file_path=self.config.output_file, url=url) + url_has_been_copied( + file_path=self.config.output_file.replace(".nx", "_0000.nx"), + url=url, + ) ) # test with some extra parameters @@ -702,14 +716,18 @@ class TestStandardAcqConversionWithExternalUrls(unittest.TestCase): configuration=self.config, ) scan.clear_caches() - assert numpy.isclose(scan.energy, 12.2) + energy = scan.energy + assert numpy.isclose(energy, 12.2) + assert scan.x_pixel_size is not None assert numpy.isclose(scan.x_pixel_size, 2.6 * 10e-6) + assert scan.y_pixel_size is not None assert numpy.isclose(scan.y_pixel_size, 2.7 * 10e-6) with EntryReader( DataUrl(file_path=scan.master_file, data_path=scan.entry, scheme="h5py") ) as entry: assert "instrument/detector" in entry assert "instrument/detector/estimated_cor_from_motor" in entry + assert "instrument/diode" not in entry class TestZSeriesConversionWithExternalUrls(unittest.TestCase): @@ -895,6 +913,8 @@ class TestH5ToNxFrmCommandLine(unittest.TestCase): input_file, output_file, "--copy-data", + "--raises-error", + "--single-file", ) ) subprocess.call(cmd, shell=True, cwd=os.path.dirname(self.input_file)) @@ -913,3 +933,10 @@ class TestH5ToNxFrmCommandLine(unittest.TestCase): for vs_info in dataset.virtual_sources(): self.assertTrue(dataset.is_virtual) self.assertEqual(os.path.realpath(vs_info.file_name), self.output_file) + + with h5py.File(output_file, "r") as h5f: + assert "/entry0000/instrument/diode" not in h5f + + # insure an nxtomo can be created from it + nx_tomo = NXtomo("").load(output_file, "entry0000") + assert nx_tomo.energy is not None diff --git a/nxtomomill/converter/hdf5/utils.py b/nxtomomill/converter/hdf5/utils.py index fa890d35858b8497de17acf29b5ef53ffa59ffcf..bb33669e5a65518333daafcb774710e44e2c8e17 100644 --- a/nxtomomill/converter/hdf5/utils.py +++ b/nxtomomill/converter/hdf5/utils.py @@ -52,6 +52,7 @@ H5ScanTitles = namedtuple( [ "init_titles", "init_zserie_titles", + "init_pcotomo_titles", "dark_titles", "flat_titles", "proj_titles", diff --git a/nxtomomill/converter/utils.py b/nxtomomill/converter/utils.py deleted file mode 100644 index 7ec3bfc4a265b601ea6600a1526c3632ca7637bb..0000000000000000000000000000000000000000 --- a/nxtomomill/converter/utils.py +++ /dev/null @@ -1,101 +0,0 @@ -# coding: utf-8 -# /*########################################################################## -# -# Copyright (c) 2015-2020 European Synchrotron Radiation Facility -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in -# all copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -# THE SOFTWARE. -# -# ###########################################################################*/ - -""" -module to define some converter utils function -""" - -__authors__ = [ - "H. Payno", -] -__license__ = "MIT" -__date__ = "02/08/2021" - - -from tomoscan.io import HDF5File -import h5py -from typing import Iterable - - -def create_nx_data_group(file_path: str, entry_path: str, axis_scale: list): - """ - Create the 'Nxdata' group at entry level with soft links on the NXDetector - and NXsample. - - :param file_path: - :param entry_path: - :param axis_scale: - :return: - """ - if not isinstance(file_path, str): - raise TypeError("file_path is expected to be a file") - if not isinstance(entry_path, str): - raise TypeError("entry_path is expected to be a file") - if not isinstance(axis_scale, Iterable): - raise TypeError("axis_scale is expected to be an Iterable") - - with HDF5File(filename=file_path, mode="a") as h5f: - entry_group = h5f[entry_path] - - nx_data_grp = entry_group.require_group("data") - # link detector datasets: - if not entry_path.startswith("/"): - entry_path = "/" + entry_path - for dataset in ("data", "image_key", "image_key_control"): - dataset_path = "/".join((entry_path, "instrument", "detector", dataset)) - nx_data_grp[dataset] = h5py.SoftLink(dataset_path) - # link rotation angle - nx_data_grp["rotation_angle"] = h5py.SoftLink( - "/".join((entry_path, "sample", "rotation_angle")) - ) - - # write NX attributes - nx_data_grp.attrs["NX_class"] = "NXdata" - nx_data_grp.attrs["signal"] = "data" - nx_data_grp.attrs["SILX_style/axis_scale_types"] = axis_scale - nx_data_grp["data"].attrs["interpretation"] = "image" - - -def link_nxbeam_to_root(file_path, entry_path): - """ - Create the 'Nxdata' group at entry level with soft links on the NXDetector - and NXsample. - - :param file_path: - :param entry_path: - :return: - """ - if not isinstance(file_path, str): - raise TypeError("file_path is expected to be a file") - if not isinstance(entry_path, str): - raise TypeError("entry_path is expected to be a file") - - if not entry_path.startswith("/"): - entry_path = "/" + entry_path - with HDF5File(filename=file_path, mode="a") as h5f: - entry_group = h5f[entry_path] - entry_group["beam"] = h5py.SoftLink( - "/".join((entry_path, "instrument", "beam")) - ) diff --git a/nxtomomill/io/config.py b/nxtomomill/io/config.py index 08f0b276014c10b5c85f63d72451c87b1358879a..fb811b4f2dc385aba1165897da58324da04aba08 100644 --- a/nxtomomill/io/config.py +++ b/nxtomomill/io/config.py @@ -299,6 +299,20 @@ class TomoHDF5Config: COMMENTS.update(COMMENTS_FRAME_TYPE_SECTION) COMMENTS.update(COMMENTS_EXTRA_PARAMS_SECTION) + @staticmethod + def get_extra_params_default_unit(key) -> str: + """return the default unit for the extra parameters that can be defined by the user""" + if key in ( + TomoHDF5Config.EXTRA_PARAMS_DISTANCE, + TomoHDF5Config.EXTRA_PARAMS_X_PIXEL_SIZE_DK, + TomoHDF5Config.EXTRA_PARAMS_Y_PIXEL_SIZE_DK, + ): + return "m" + elif key in (TomoHDF5Config.EXTRA_PARAMS_ENERGY_DK): + return "keV" + else: + raise ValueError(f"No default unit for {key}") + def __init__(self): # general information self._output_file = None @@ -328,6 +342,7 @@ class TomoHDF5Config: self._sub_entries_to_ignore = None self._init_titles = settings.Tomo.H5.INIT_TITLES self._zserie_init_titles = settings.Tomo.H5.ZSERIE_INIT_TITLES + self._pcotomo_init_titles = settings.Tomo.H5.PCOTOMO_INIT_TITLES self._dark_titles = settings.Tomo.H5.DARK_TITLES self._flat_titles = settings.Tomo.H5.FLAT_TITLES self._projection_titles = settings.Tomo.H5.PROJ_TITLES @@ -651,6 +666,19 @@ class TomoHDF5Config: else: self._zserie_init_titles = titles + @property + def pcotomo_init_titles(self) -> Union[None, Iterable]: + return self._pcotomo_init_titles + + @pcotomo_init_titles.setter + def pcotomo_init_titles(self, titles: Union[None, Iterable]) -> None: + if titles is None: + self._pcotomo_init_titles = None + elif not isinstance(titles, Iterable): + raise TypeError("'titles' should be None or an Iterable") + else: + self._pcotomo_init_titles = titles + @property def dark_titles(self) -> Union[None, Iterable]: return self._dark_titles diff --git a/nxtomomill/io/confighandler.py b/nxtomomill/io/confighandler.py index 01ff86beccec4ffc6a1543b379ddff9723e6c50c..13ee60e9cb5716ea55da9808002db76bbb1fd1dc 100644 --- a/nxtomomill/io/confighandler.py +++ b/nxtomomill/io/confighandler.py @@ -249,6 +249,10 @@ class TomoHDF5ConfigHandler(BaseHDF5ConfigHandler): "zserie_init_titles", self.get_tuple_of_keys_from_cmd, ), + "init_pcotomo_titles": ( + "pcotomo_init_titles", + self.get_tuple_of_keys_from_cmd, + ), "dark_titles": ("dark_titles", self.get_tuple_of_keys_from_cmd), "flat_titles": ("flat_titles", self.get_tuple_of_keys_from_cmd), "proj_titles": ("projections_titles", self.get_tuple_of_keys_from_cmd), @@ -283,7 +287,7 @@ class TomoHDF5ConfigHandler(BaseHDF5ConfigHandler): "dark_titles", "init_zserie_titles", "init_titles", - "init_zserie_titles", + "init_pcotomo_titles", "x_pixel_size_key", "y_pixel_size_key", "acq_expo_time_keys", @@ -360,6 +364,10 @@ class XRD3DHDF5ConfigHandler(BaseHDF5ConfigHandler): "zserie_init_titles", self.get_tuple_of_keys_from_cmd, ), + "init_pcotomo_titles": ( + "pcotomo_init_titles", + self.get_tuple_of_keys_from_cmd, + ), "dark_titles": ("dark_titles", self.get_tuple_of_keys_from_cmd), "flat_titles": ("flat_titles", self.get_tuple_of_keys_from_cmd), "proj_titles": ("projections_titles", self.get_tuple_of_keys_from_cmd), @@ -394,7 +402,7 @@ class XRD3DHDF5ConfigHandler(BaseHDF5ConfigHandler): "dark_titles", "init_zserie_titles", "init_titles", - "init_zserie_titles", + "init_pcotomo_titles", "x_pixel_size_key", "y_pixel_size_key", "acq_expo_time_keys", diff --git a/nxtomomill/io/test/test_confighandler.py b/nxtomomill/io/test/test_confighandler.py index 7f4389abb8347e2ed57787f262e807d4063e4b23..e8ad444a0fb1e8afe9b16c221c7cea3fde00423c 100644 --- a/nxtomomill/io/test/test_confighandler.py +++ b/nxtomomill/io/test/test_confighandler.py @@ -48,6 +48,7 @@ class _ARParseMock(object): self.flat_titles = None self.dark_titles = None self.init_zserie_titles = None + self.init_pcotomo_titles = None self.init_titles = None self.x_pixel_size_key = None self.y_pixel_size_key = None diff --git a/nxtomomill/nexus/nxdetector.py b/nxtomomill/nexus/nxdetector.py index 347ab0526dce17ed4842ea010831043373a107f6..7a1c5a40ef9f86b7c9f7d32d9e4fcf5ae8ccd247 100644 --- a/nxtomomill/nexus/nxdetector.py +++ b/nxtomomill/nexus/nxdetector.py @@ -30,25 +30,42 @@ __date__ = "03/02/2022" from silx.utils.proxy import docstring +from silx.io.url import DataUrl +from h5py import VirtualSource from typing import Iterable, Optional, Union import numpy +import h5py +from nxtomomill.utils.h5pyutils import from_virtual_source_to_data_url try: from tomoscan.esrf.scan.hdf5scan import ImageKey - from tomoscan.esrf.scan.hdf5scan import get_nexus_paths except ImportError: from tomoscan.esrf.hdf5scan import ImageKey - from tomoscan.esrf.hdf5scan import get_nexus_paths +from tomoscan.nexus.paths.nxtomo import get_paths as get_nexus_path from tomoscan.unitsystem.metricsystem import MetricSystem -from nxtomomill.nexus.utils import cast_and_check_array_1D +from nxtomomill.nexus.utils import cast_and_check_array_1D, get_data, get_data_and_unit from .nxobject import ElementWithUnit, NXobject from tomoscan.scanbase import FOV as FieldOfView from tomoscan.unitsystem import TimeSystem +from tomoscan.unitsystem import Unit +from tomoscan.io import HDF5File +from h5py import h5s as h5py_h5s +import h5py._hl.selections as selection class NXdetector(NXobject): - def __init__(self) -> None: - super().__init__() + def __init__( + self, + node_name="detector", + parent=None, + field_of_view=None, + expected_dim: Optional[tuple] = None, + ) -> None: + """ + :param Optional[tuple] expected_dim: user can provide expected dimesions as a tuple of int to be checked when data is set + """ + super().__init__(node_name=node_name, parent=parent) + self._expected_dim = expected_dim self._data = None self._image_key_control = None @@ -57,22 +74,39 @@ class NXdetector(NXobject): self._distance = ElementWithUnit( default_unit=MetricSystem.METER ) # detector / sample distance - self._field_of_view = FieldOfView.FULL + self.field_of_view = field_of_view self._count_time = ElementWithUnit(default_unit=TimeSystem.SECOND) self._estimated_cor_from_motor = None @property - def data(self) -> numpy.ndarray: + def data(self) -> Union[numpy.ndarray, tuple]: + """data can be None, a numpy array or a list of DataUrl xor h5py Virtual Source""" return self._data @data.setter - def data(self, data: Optional[numpy.ndarray]): - if not isinstance(data, (type(None), numpy.ndarray)): + def data(self, data: Optional[Union[numpy.ndarray, tuple]]): + if isinstance(data, numpy.ndarray): + if ( + self._expected_dim is not None + and data is not None + and data.ndim not in self._expected_dim + ): + raise ValueError( + f"data is expected to be {self._expected_dim}d not {data.ndim}d" + ) + elif isinstance(data, (tuple, list)): + for elmt in data: + if not isinstance(elmt, (DataUrl, VirtualSource)): + raise TypeError( + f"element of 'data' are expected to be silx DataUrl or h5py virtualSource. Not {type(elmt)}" + ) + data = tuple(data) + elif data is None: + pass + else: raise TypeError( - f"data is expected ot be an instance of {numpy.ndarray} or None. Not {type(data)}" + f"data is expected to be an instance of {numpy.ndarray}, None or a list of silx DataUrl or h5py Virtual Source. Not {type(data)}" ) - if data.ndim not in (2, 3): - raise ValueError(f"data is expected to be 2 or 3d not {data.ndim}d") self._data = data @property @@ -121,7 +155,9 @@ class NXdetector(NXobject): return self._field_of_view @field_of_view.setter - def field_of_view(self, field_of_view: Optional[Union[FieldOfView, str]]) -> None: + def field_of_view( + self, field_of_view: Optional[Union[FieldOfView, str, None]] + ) -> None: if field_of_view is not None: field_of_view = FieldOfView.from_value(field_of_view) self._field_of_view = field_of_view @@ -181,54 +217,261 @@ class NXdetector(NXobject): return control_image_key @docstring(NXobject) - def to_nx_dict(self, nexus_path_version: Optional[float] = None) -> dict: - nexus_paths = get_nexus_paths(nexus_path_version) + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: + nexus_paths = get_nexus_path(nexus_path_version) + nexus_detector_paths = nexus_paths.nx_detector_paths + nx_dict = {} - if self.data is not None: - nx_dict[nexus_paths.PROJ_PATH] = self.data - nx_dict["@".join([nexus_paths.PROJ_PATH, "interpretation"])] = "image" if self.image_key_control is not None: - nx_dict[nexus_paths.IMG_KEY_PATH] = [ - img_key.value for img_key in self.image_key - ] - nx_dict[nexus_paths.IMG_KEY_CONTROL_PATH] = [ + path_img_key = f"{self.path}/{nexus_detector_paths.IMAGE_KEY}" + nx_dict[path_img_key] = [img_key.value for img_key in self.image_key] + path_img_key_ctrl = f"{self.path}/{nexus_detector_paths.IMAGE_KEY_CONTROL}" + nx_dict[path_img_key_ctrl] = [ img_key.value for img_key in self.image_key_control ] - if self.x_pixel_size is not None: - nx_dict[nexus_paths.X_PIXEL_SIZE_PATH] = self.x_pixel_size.value - nx_dict["@".join([nexus_paths.X_PIXEL_SIZE_PATH, "unit"])] = str( - self.distance.unit + if self.x_pixel_size.value is not None: + path_x_pixel_size = f"{self.path}/{nexus_detector_paths.X_PIXEL_SIZE}" + nx_dict[path_x_pixel_size] = self.x_pixel_size.value + nx_dict["@".join([path_x_pixel_size, "unit"])] = str(self.distance.unit) + if self.y_pixel_size.value is not None: + path_y_pixel_size = f"{self.path}/{nexus_detector_paths.Y_PIXEL_SIZE}" + nx_dict[path_y_pixel_size] = self.y_pixel_size.value + nx_dict["@".join([path_y_pixel_size, "unit"])] = str(self.distance.unit) + if self.distance.value is not None: + path_distance = f"{self.path}/{nexus_detector_paths.DISTANCE}" + nx_dict[path_distance] = self.distance.value + nx_dict["@".join([path_distance, "unit"])] = str(self.distance.unit) + if self.field_of_view is not None: + path_fov = f"{self.path}/{nexus_detector_paths.FOV}" + nx_dict[path_fov] = self.field_of_view.value + if self.count_time.value is not None: + path_count_time = f"{self.path}/{nexus_detector_paths.EXPOSURE_TIME}" + nx_dict[path_count_time] = self.count_time.value + nx_dict["@".join([path_count_time, "unit"])] = str(self.count_time.unit) + if self.estimated_cor_from_motor is not None: + path_estimated_cor = ( + f"{self.path}/{nexus_detector_paths.ESTIMATED_COR_FRM_MOTOR}" ) - if self.y_pixel_size is not None: - nx_dict[nexus_paths.Y_PIXEL_SIZE_PATH] = self.y_pixel_size.value - nx_dict["@".join([nexus_paths.Y_PIXEL_SIZE_PATH, "unit"])] = str( - self.distance.unit + nx_dict[path_estimated_cor] = self.estimated_cor_from_motor + nx_dict["@".join([path_estimated_cor, "unit"])] = "pixel" + + nx_dict.update( + self._data_to_nx_dict( + nexus_path_version=nexus_path_version, + data_path=data_path, ) - if self.distance is not None: - nx_dict[nexus_paths.DISTANCE_PATH] = self.distance.value - nx_dict["@".join([nexus_paths.DISTANCE_PATH, "unit"])] = str( - self.distance.unit + ) + return nx_dict + + def _data_to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: + nexus_paths = get_nexus_path(nexus_path_version) + nexus_detector_paths = nexus_paths.nx_detector_paths + + nx_dict = {} + if self.data is not None: + # add data + path_data = f"{self.path}/{nexus_detector_paths.DATA}" + nx_dict[path_data] = self.data + nx_dict["@".join([path_data, "interpretation"])] = "image" + # add attributes to data + nx_dict[f"{self.path}@NX_class"] = "NXdetector" + nx_dict[f"{self.path}@signal"] = nexus_detector_paths.DATA + nx_dict[f"{self.path}@SILX_style/axis_scale_types"] = [ + "linear", + "linear", + ] + return nx_dict + + def _load( + self, file_path: str, data_path: str, nexus_version: float, load_data_as: str + ) -> None: + possible_as_values = ("as_virtual_source", "as_data_url", "as_numpy_array") + if load_data_as not in possible_as_values: + raise ValueError( + f"load_data_as is expected to be in {possible_as_values} and not {load_data_as}" ) - if self.field_of_view is not None: - nx_dict[nexus_paths.FOV_PATH] = self.field_of_view.value - if self.count_time is not None: - nx_dict[nexus_paths.EXPOSURE_TIME_PATH] = self.count_time.value - nx_dict["@".join([nexus_paths.EXPOSURE_TIME_PATH, "unit"])] = str( - self.count_time.unit + + nexus_paths = get_nexus_path(nexus_version) + nexus_detector_paths = nexus_paths.nx_detector_paths + + data_dataset_path = f"{data_path}/{nexus_detector_paths.DATA}" + with HDF5File(file_path, mode="r") as h5f: + if data_dataset_path in h5f: + dataset = h5f[data_dataset_path] + if load_data_as == "as_numpy_array": + self.data = dataset[()] + elif load_data_as == "as_data_url": + if dataset.is_virtual: + urls = [] + for vs_info in dataset.virtual_sources(): + dataset = h5f[vs_info.dset_name] + select_bounds = vs_info.vspace.get_select_bounds() + left_bound = select_bounds[0] + right_bound = select_bounds[1] + # warning: for now step is not managed with virtual + # dataset + + length = right_bound[0] - left_bound[0] + 1 + # warning: for now step is not managed with virtual + # dataset + + virtual_source = h5py.VirtualSource( + vs_info.file_name, + vs_info.dset_name, + shape=dataset.shape, + ) + # here we could provide dataset but we won't to + # insure file path will be relative. + type_code = vs_info.src_space.get_select_type() + # check for unlimited selections in case where selection is regular + # hyperslab, which is the only allowed case for h5s.UNLIMITED to be + # in the selection + if ( + type_code == h5py_h5s.SEL_HYPERSLABS + and vs_info.src_space.is_regular_hyperslab() + ): + ( + source_start, + stride, + count, + block, + ) = vs_info.src_space.get_regular_hyperslab() + source_end = source_start[0] + length + + sel = selection.select( + dataset.shape, + slice(source_start[0], source_end), + dataset=dataset, + ) + virtual_source.sel = sel + + urls.append(from_virtual_source_to_data_url(virtual_source)) + else: + urls = [ + DataUrl( + file_path=file_path, + data_path=data_dataset_path, + scheme="silx", + ) + ] + self.data = urls + elif load_data_as == "as_virtual_source": + if dataset.is_virtual: + self.data = dataset.virtual_sources() + else: + raise ValueError(f"{data_dataset_path} is not virtual") + + self.x_pixel_size, self.x_pixel_size.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.X_PIXEL_SIZE]), + default_unit=MetricSystem.METER, + ) + self.y_pixel_size, self.y_pixel_size.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.Y_PIXEL_SIZE]), + default_unit=MetricSystem.METER, + ) + self.distance, self.distance.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.DISTANCE]), + default_unit=MetricSystem.METER, + ) + self.distance, self.distance.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.DISTANCE]), + default_unit=MetricSystem.METER, + ) + self.field_of_view = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.FOV]), + ) + self.count_time, self.count_time.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.EXPOSURE_TIME]), + default_unit=TimeSystem.SECOND, + ) + self.estimated_cor_from_motor = get_data( + file_path=file_path, + data_path="/".join( + [data_path, nexus_detector_paths.ESTIMATED_COR_FRM_MOTOR] + ), + ) + self.image_key_control = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_detector_paths.IMAGE_KEY_CONTROL]), + ) + + +class NXdetectorWithUnit(NXdetector): + def __init__( + self, + default_unit: Unit, + node_name="detector", + parent=None, + field_of_view=None, + expected_dim: Optional[tuple] = None, + ) -> None: + super().__init__(node_name, parent, field_of_view, expected_dim) + self._data = ElementWithUnit(default_unit=default_unit) + + @property + def data(self) -> Union[numpy.ndarray, tuple]: + """data can be None, a numpy array or a list of DataUrl xor h5py Virtual Source""" + return self._data + + @data.setter + def data(self, data: Optional[Union[numpy.ndarray, tuple]]): + if isinstance(data, numpy.ndarray): + if ( + self._expected_dim is not None + and data is not None + and data.ndim not in self._expected_dim + ): + raise ValueError( + f"data is expected to be {self._expected_dim}d not {data.ndim}d" + ) + elif isinstance(data, (tuple, list)): + for elmt in data: + if not isinstance(elmt, (DataUrl, VirtualSource)): + raise TypeError( + f"'data' is expected to be a numpy array or a list/tuple composed of DataUrl or h5py virtualSource. Not {type(elmt)}" + ) + data = tuple(data) + elif data is None: + pass + else: + raise TypeError( + f"data is expected to be an instance of {numpy.ndarray}, None or a list of silx DataUrl or h5py Virtual Source. Not {type(data)}" ) - if self.estimated_cor_from_motor is not None: - nx_dict[ - nexus_paths.ESTIMATED_COR_FRM_MOTOR_PATH - ] = self.estimated_cor_from_motor + self._data.value = data - if self.data is not None: - nx_dict[f"{self._entry}@default"] = "instrument/detector" - nx_dict["instrument@NX_class"] = "NXdata" - nx_dict["instrument@default"] = "detector" - nx_dict["instrument/detector@NX_class"] = "NXdata" - nx_dict["instrument/detector@signal"] = "data" - nx_dict["instrument/detector@SILX_style/axis_scale_types"] = [ + def _data_to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: + nexus_paths = get_nexus_path(nexus_path_version) + nexus_detector_paths = nexus_paths.nx_detector_paths + + nx_dict = {} + if self.data.value is not None: + # add data + path_data = f"{self.path}/{nexus_detector_paths.DATA}" + nx_dict[path_data] = self.data.value + nx_dict["@".join([path_data, "interpretation"])] = "image" + # add attributes to data + nx_dict[f"{self.path}@NX_class"] = "NXdetector" + nx_dict[f"{self.path}@signal"] = nexus_detector_paths.DATA + nx_dict[f"{self.path}@SILX_style/axis_scale_types"] = [ "linear", "linear", ] diff --git a/nxtomomill/nexus/nxinstrument.py b/nxtomomill/nexus/nxinstrument.py new file mode 100644 index 0000000000000000000000000000000000000000..ef003963a279787409deff645fcd1ed7b983ed69 --- /dev/null +++ b/nxtomomill/nexus/nxinstrument.py @@ -0,0 +1,205 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2022 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "10/02/2022" + + +from silx.utils.proxy import docstring +from typing import Optional +from nxtomomill.nexus.nxdetector import NXdetector, NXdetectorWithUnit +from nxtomomill.nexus.nxsource import NXsource +from nxtomomill.nexus.nxsource import DefaultESRFSource +from tomoscan.nexus.paths.nxtomo import get_paths as get_nexus_paths +from tomoscan.unitsystem.voltagesystem import VoltageSystem +from .utils import get_data +from tomoscan.io import HDF5File +from .nxobject import NXobject + + +class NXinstrument(NXobject): + def __init__( + self, node_name: str = "instrument", parent: Optional[NXobject] = None + ) -> None: + super().__init__(node_name=node_name, parent=parent) + self._detector = NXdetector( + node_name="detector", + parent=self, + field_of_view="Full", + expected_dim=(2, 3), + ) + self._diode = NXdetectorWithUnit( + node_name="diode", + parent=self, + expected_dim=(1,), + default_unit=VoltageSystem.VOLT, + ) + self._source = DefaultESRFSource(node_name="source", parent=self) + self._name = None + + @property + def name(self) -> Optional[str]: + return self._name + + @name.setter + def name(self, name: Optional[str]): + if not isinstance(name, (str, type(None))): + raise TypeError(f"name is expected to be None or a str. Not {type(name)}") + self._name = name + + @property + def detector(self) -> Optional[NXdetector]: + return self._detector + + @detector.setter + def detector(self, detector: Optional[NXdetector]): + if not isinstance(detector, (NXdetector, type(None))): + raise TypeError( + f"detector is expected to be None or an instance of NXdetecetor. Not {type(detector)}" + ) + self._detector = detector + + @property + def diode(self) -> Optional[NXdetector]: + return self._diode + + @diode.setter + def diode(self, diode: Optional[NXdetector]): + if not isinstance(diode, (NXdetector, type(None))): + raise TypeError( + f"diode is expected to be None or an instance of NXdetecetor. Not {type(diode)}" + ) + self._diode = diode + + @property + def source(self) -> Optional[NXsource]: + return self._source + + @source.setter + def source(self, source: Optional[NXsource]) -> None: + if not isinstance(source, (NXsource, type(None))): + raise TypeError( + f"source is expected to be None or an instance of NXsource. Not {type(source)}" + ) + self._source = source + + @property + def name(self) -> Optional[str]: + return self._name + + @name.setter + def name(self, name: Optional[str]) -> None: + if not isinstance(name, (str, type(None))): + raise TypeError( + f"name is expected to be None or an instance of str. Not {type(name)}" + ) + self._name = name + + @docstring(NXobject) + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: + nexus_paths = get_nexus_paths(nexus_path_version) + nexus_instrument_paths = nexus_paths.nx_instrument_paths + nx_dict = {} + + if self._detector is not None: + nx_dict.update( + self._detector.to_nx_dict(nexus_path_version=nexus_path_version) + ) + + if self._diode is not None: + nx_dict.update( + self._diode.to_nx_dict(nexus_path_version=nexus_path_version) + ) + + if self._source is not None: + nx_dict.update( + self.source.to_nx_dict(nexus_path_version=nexus_path_version) + ) + + if self.name is not None: + nx_dict[f"{self.name}/{nexus_instrument_paths.NAME}"] = self.name + if nx_dict != {}: + nx_dict[f"{self.path}@NX_class"] = "NXinstrument" + + return nx_dict + + def _load( + self, + file_path: str, + data_path: str, + nexus_version: float, + detector_data_as: str, + ) -> NXobject: + """ + Create and load an NXsample from data on disk + """ + nexus_paths = get_nexus_paths(nexus_version) + nexus_instrument_paths = nexus_paths.nx_instrument_paths + + with HDF5File(file_path, mode="r") as h5f: + if data_path in h5f: + has_detector = "detector" in h5f[data_path] + has_diode = "diode" in h5f[data_path] + has_source = "source" in h5f[data_path] + else: + has_detector = False + has_diode = False + has_source = False + # TODO: loading detector might be done using the NXclass instead of some hard coded names + if has_detector: + self.detector._load( + file_path=file_path, + data_path="/".join( + [data_path, "detector"], + ), + nexus_version=nexus_version, + load_data_as=detector_data_as, + ) + if has_diode: + self.diode._load( + file_path=file_path, + data_path="/".join( + [data_path, "diode"], + ), + nexus_version=nexus_version, + load_data_as="as_numpy_array", + ) + if has_source: + self.source._load( + file_path=file_path, + data_path="/".join([data_path, "source"]), + nexus_version=nexus_version, + ) + if nexus_instrument_paths.NAME is not None: + self.name = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_instrument_paths.NAME]), + ) diff --git a/nxtomomill/nexus/nxobject.py b/nxtomomill/nexus/nxobject.py index 9395c9a64e089bd352e94cd1f9e2beb65ec49780..b59bbfeac64a7d5b45e73da35b332f87f1c21d9b 100644 --- a/nxtomomill/nexus/nxobject.py +++ b/nxtomomill/nexus/nxobject.py @@ -30,7 +30,9 @@ __date__ = "03/02/2022" from typing import Optional +from silx.io.url import DataUrl from tomoscan.unitsystem import Unit +from tomoscan.nexus.paths.nxtomo import LATEST_VERSION as LATEST_NXTOMO_VERSION from silx.io.dictdump import dicttonx import h5py import os @@ -57,9 +59,12 @@ class ElementWithUnit: except Exception: pass if not isinstance(unit, self._unit_type): - raise TypeError( - f"invalid unit type. {type(unit)} provided when {type(self._unit_type)} expected" - ) + if isinstance(unit, str): + raise ValueError(f"Unable to cast {unit} to a {type(self._unit_type)}") + else: + raise TypeError( + f"invalid unit type. {type(unit)} provided when {type(self._unit_type)} expected" + ) self._unit = unit @property @@ -79,13 +84,64 @@ class ElementWithUnit: class NXobject: - def __init__(self) -> None: - self._entry = "./" # set when save it + def __init__(self, node_name: str, parent=None) -> None: + if not isinstance(node_name, str): + raise TypeError( + f"name is expected to be an instance of str. Not {type(node_name)}" + ) + self.node_name = node_name + self.parent = parent + + @property + def parent(self): # -> Optional[NXobject]: + return self._parent + + @parent.setter + def parent(self, parent) -> None: + if not isinstance(parent, (type(None), NXobject)): + raise TypeError( + f"parent is expected to be None or an instance of {NXobject}" + ) + self._parent = parent + + @property + def is_root(self) -> bool: + return self.parent is None + + @property + def root_path(self) -> str: + """return path of the root NXobject""" + if self.is_root: + return self.path + else: + return self.parent.root_path + + @property + def path(self): + if self.parent is not None: + path = "/".join([self.parent.path, self.node_name]) + else: + path = self.node_name + # clean some possible issues with "//" + path = path.replace("//", "/") + return path + + @property + def node_name(self) -> str: + return self._node_name + + @node_name.setter + def node_name(self, node_name: str): + if not isinstance(node_name, str): + raise TypeError( + f"nexus_name should be an instance of str and not {type(node_name)}" + ) + self._node_name = node_name def save( self, file_path: str, - entry: str, + data_path: str = "/", nexus_path_version: Optional[float] = None, overwrite: bool = False, ) -> None: @@ -95,7 +151,9 @@ class NXobject: :param str file_path: hdf5 file :param Optional[str] entry: HDF5 group to save the data. If entry is already provided in """ - for key, value in dict([("file_path", file_path), ("entry", entry)]).items(): + for key, value in dict( + [("file_path", file_path), ("entry", data_path)] + ).items(): if not isinstance(value, (type(None), str)): raise TypeError( f"{key} is expected to be None or an instance of str not {type(value)}" @@ -103,27 +161,136 @@ class NXobject: if not isinstance(overwrite, bool): raise TypeError + entry_path = "/".join([data_path, self.path]) # not fully sure about the dicttoh5 "add" behavior if os.path.exists(file_path): with h5py.File(file_path, mode="a") as h5f: - if entry in h5f: + if entry_path in h5f: if overwrite: - del h5f[entry] + del h5f[entry_path] else: - raise KeyError(f"{entry} already exists") + raise KeyError(f"{entry_path} already exists") + if nexus_path_version is None: + nexus_path_version = LATEST_NXTOMO_VERSION + + nx_dict = self.to_nx_dict( + nexus_path_version=nexus_path_version, data_path=data_path + ) + # retrieve virtual sources and DataUrl + datasets_to_handle_in_postprocessing = {} + for key in self._get_virtual_sources(nx_dict): + datasets_to_handle_in_postprocessing[key] = nx_dict.pop(key) + for key in self._get_data_urls(nx_dict): + datasets_to_handle_in_postprocessing[key] = nx_dict.pop(key) + + # retrieve attributes + attributes = {} + + dataset_to_postpone = tuple(datasets_to_handle_in_postprocessing.keys()) + for key, value in nx_dict.items(): + if key.startswith(dataset_to_postpone): + attributes[key] = value + # clean attributes + for key in attributes: + del nx_dict[key] - self._entry = entry dicttonx( - self.to_nx_dict(nexus_path_version=nexus_path_version), + nx_dict, h5file=file_path, - h5path=entry, + h5path=data_path, update_mode="replace", + mode="a", + ) + assert os.path.exists(file_path) + # now handle nx_dict containing h5py.virtualSource or DataUrl + # this cannot be handled from the nxdetector class because not aware about + # the output file. + for ( + dataset_path, + v_sources_or_data_urls, + ) in datasets_to_handle_in_postprocessing.items(): + for v_source_or_data_url in v_sources_or_data_urls: + if isinstance(v_source_or_data_url, DataUrl): + data = DataUrl( + file_path=v_source_or_data_url.file_path(), + data_path=v_source_or_data_url.data_path(), + scheme=v_source_or_data_url.scheme(), + data_slice=v_source_or_data_url.data_slice(), + ) + else: + data = v_source_or_data_url + from nxtomomill.utils.frameappender import ( + FrameAppender, + ) # avoid cyclic import + + frame_appender = FrameAppender( + data=data, + file_path=file_path, + data_path="/".join([data_path, dataset_path]), + where="end", + ) + frame_appender.process() + + # write attributes of dataset defined from a list of DataUrl or VirtualSource + assert os.path.exists(file_path) + dicttonx( + attributes, + h5file=file_path, + h5path=data_path, + update_mode="add", + mode="a", ) - def to_nx_dict(self, nexus_path_version: Optional[float] = None) -> dict: + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: """ convert the NXobject to an nx dict. Dictionnary that we can dump to hdf5 file :param Optional[float] nexus_path_version: version of the nexus path version to use + :param Optional[str] data_path: can be provided to create some link in the file """ raise NotImplementedError("Base class") + + def __str__(self): + f"{type(self)}: {self.path}" + + @staticmethod + def _get_virtual_sources(ddict) -> tuple: + """Return key / path containing a list or a tuple of h5py.VirtualSource""" + + def has_virtual_sources(value): + if isinstance(value, h5py.VirtualSource): + return True + elif isinstance(value, (list, tuple)): + for v in value: + if has_virtual_sources(v): + return True + return False + + keys = [] + for key, value in ddict.items(): + if has_virtual_sources(value): + keys.append(key) + return tuple(keys) + + @staticmethod + def _get_data_urls(ddict) -> tuple: + """Return key / path containing a list or a tuple of silx.io.url.DataUrl""" + + def has_data_url(value): + if isinstance(value, DataUrl): + return True + elif isinstance(value, (list, tuple)): + for v in value: + if has_data_url(v): + return True + return False + + keys = [] + for key, value in ddict.items(): + if has_data_url(value): + keys.append(key) + return tuple(keys) diff --git a/nxtomomill/nexus/nxsample.py b/nxtomomill/nexus/nxsample.py index 74ecf9702a4eea94efbbf5062569d247fe5b3faa..a8fb86d1cb967d96af86f8c000bc9d47ff3bbf14 100644 --- a/nxtomomill/nexus/nxsample.py +++ b/nxtomomill/nexus/nxsample.py @@ -33,19 +33,15 @@ from typing import Iterable, Optional import numpy from .nxobject import NXobject, ElementWithUnit -from .utils import cast_and_check_array_1D +from .utils import cast_and_check_array_1D, get_data_and_unit, get_data from silx.utils.proxy import docstring - -try: - from tomoscan.esrf.scan.hdf5scan import get_nexus_paths -except ImportError: - from tomoscan.esrf.hdf5scan import get_nexus_paths +from tomoscan.nexus.paths.nxtomo import get_paths as get_nexus_paths from tomoscan.unitsystem.metricsystem import MetricSystem class NXsample(NXobject): - def __init__(self) -> None: - super().__init__() + def __init__(self, node_name="sample", parent: Optional[NXobject] = None) -> None: + super().__init__(node_name=node_name, parent=parent) self._name = None self._rotation_angle = None self._x_translation = ElementWithUnit(default_unit=MetricSystem.METER) @@ -101,31 +97,82 @@ class NXsample(NXobject): ) @docstring(NXobject) - def to_nx_dict(self, nexus_path_version: Optional[float] = None) -> dict: + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: nexus_paths = get_nexus_paths(nexus_path_version) + nexus_sample_paths = nexus_paths.nx_sample_paths nx_dict = {} if self.name is not None: - nx_dict[nexus_paths.SAMPLE_NAME_PATH] = self.name + path_name = f"{self.path}/{nexus_sample_paths.NAME}" + nx_dict[path_name] = self.name if self.rotation_angle is not None: - nx_dict[nexus_paths.ROTATION_ANGLE_PATH] = self.rotation_angle - if self.x_translation is not None: - nx_dict[nexus_paths.X_TRANS_PATH] = self.x_translation.value - nx_dict["@".join([nexus_paths.X_TRANS_PATH, "unit"])] = str( + path_rotation_angle = f"{self.path}/{nexus_sample_paths.ROTATION_ANGLE}" + nx_dict[path_rotation_angle] = self.rotation_angle + nx_dict["@".join([path_rotation_angle, "unit"])] = "degree" + if self.x_translation.value is not None: + path_x_translation = f"{self.path}/{nexus_sample_paths.X_TRANSLATION}" + nx_dict[path_x_translation] = self.x_translation.value + nx_dict["@".join([path_x_translation, "unit"])] = str( self.x_translation.unit ) - if self.y_translation is not None: - nx_dict[nexus_paths.Y_TRANS_PATH] = self.y_translation.value - nx_dict["@".join([nexus_paths.Y_TRANS_PATH, "unit"])] = str( + if self.y_translation.value is not None: + path_y_translation = f"{self.path}/{nexus_sample_paths.Y_TRANSLATION}" + nx_dict[path_y_translation] = self.y_translation.value + nx_dict["@".join([path_y_translation, "unit"])] = str( self.y_translation.unit ) - if self.z_translation is not None: - nx_dict[nexus_paths.Z_TRANS_PATH] = self.z_translation.value - nx_dict["@".join([nexus_paths.Z_TRANS_PATH, "unit"])] = str( + if self.z_translation.value is not None: + path_z_translation = f"{self.path}/{nexus_sample_paths.Z_TRANSLATION}" + nx_dict[path_z_translation] = self.z_translation.value + nx_dict["@".join([path_z_translation, "unit"])] = str( self.z_translation.unit ) - nx_dict["sample@NX_class"] = "NXsample" + if nx_dict != {}: + nx_dict[f"{self.path}@NX_class"] = "NXsample" return nx_dict + + def _load(self, file_path: str, data_path: str, nexus_version: float) -> NXobject: + """ + Create and load an NXsample from data on disk + """ + nexus_paths = get_nexus_paths(nexus_version) + nexus_sample_paths = nexus_paths.nx_sample_paths + + self.name = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_sample_paths.NAME]), + ) + self.rotation_angle, angle_unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_sample_paths.ROTATION_ANGLE]), + default_unit="degree", + ) + if angle_unit == "degree": + pass + elif isinstance(angle_unit, str) and angle_unit.lower() in ("rad", "radian"): + self.rotation_angle = numpy.rad2deg(self.rotation_angle) + elif angle_unit is not None: + raise ValueError(f"rotation angle unit not recognized: {angle_unit}") + + self.x_translation, self.x_translation.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_sample_paths.X_TRANSLATION]), + default_unit=MetricSystem.METER, + ) + self.y_translation, self.y_translation.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_sample_paths.Y_TRANSLATION]), + default_unit=MetricSystem.METER, + ) + self.z_translation, self.z_translation.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_sample_paths.Z_TRANSLATION]), + default_unit=MetricSystem.METER, + ) diff --git a/nxtomomill/nexus/nxsource.py b/nxtomomill/nexus/nxsource.py index b932cedbac3fd6d9be243e86c5172a635816502a..bd914d6f646d05483018ae8d2f9f1ba37047b8c9 100644 --- a/nxtomomill/nexus/nxsource.py +++ b/nxtomomill/nexus/nxsource.py @@ -32,11 +32,8 @@ from typing import Optional, Union from silx.utils.proxy import docstring from .nxobject import NXobject from silx.utils.enum import Enum as _Enum - -try: - from tomoscan.esrf.scan.hdf5scan import get_nexus_paths, _NEXUS_PATHS_V_1_1 -except ImportError: - from tomoscan.esrf.hdf5scan import get_nexus_paths, _NEXUS_PATHS_V_1_1 +from tomoscan.nexus.paths.nxtomo import get_paths as get_nexus_paths +from .utils import get_data class SourceType(_Enum): @@ -58,46 +55,80 @@ class SourceType(_Enum): class NXsource(NXobject): """Information regarding the x-ray storage ring/facility""" - def __init__(self, name=None, type=None): - self._name = name - self._type = type + def __init__( + self, node_name="source", parent=None, source_name=None, source_type=None + ): + super().__init__(node_name=node_name, parent=parent) + self._source_name = source_name + self._source_type = source_type @property - def name(self) -> Union[None, str]: - return self._name + def source_name(self) -> Union[None, str]: + return self._source_name - @name.setter - def name(self, name: Union[str, None]): - if not isinstance(name, (str, type(None))): - raise TypeError("name is expected to be None or a str") - self._name = name + @source_name.setter + def source_name(self, source_name: Union[str, None]): + if not isinstance(source_name, (str, type(None))): + raise TypeError( + f"source_name is expected to be None or a str not {type(source_name)}" + ) + self._source_name = source_name @property - def type(self) -> Union[None, SourceType]: - return self._type + def source_type(self) -> Union[None, SourceType]: + return self._source_type - @type.setter - def type(self, type_: Union[None, str, SourceType]): + @source_type.setter + def source_type(self, type_: Union[None, str, SourceType]): if type_ is None: self._type = None else: type_ = SourceType.from_value(type_) - self._type = type_ + self._source_type = type_ def __str__(self): - return f"source (name: {self.name}, type: {self.type})" + return f"{super().__str__}, (source name: {self.source_name}, source type: {self.source_type})" @docstring(NXobject) - def to_nx_dict(self, nexus_path_version: Optional[float] = None) -> dict: - nx_dict = {} + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: nexus_paths = get_nexus_paths(nexus_path_version) + nexus_source_paths = nexus_paths.nx_source_paths + nx_dict = {} # warning: source is integrated only since 1.1 version of the nexus path - if self.name is not None and nexus_paths.SOURCE_NAME is not None: - nx_dict[nexus_paths.SOURCE_NAME] = self.name - if self.type is not None and nexus_paths.SOURCE_TYPE is not None: - nx_dict[nexus_paths.SOURCE_TYPE] = self.type.value - if isinstance(nexus_paths, _NEXUS_PATHS_V_1_1): - nx_dict["source@NX_class"] = "NXsource" + if self.source_name is not None and nexus_paths.SOURCE_NAME is not None: + path_name = f"{self.path}/{nexus_source_paths.NAME}" + nx_dict[path_name] = self.source_name + if self.source_type is not None and nexus_paths.SOURCE_TYPE is not None: + path_type = f"{self.path}/{nexus_source_paths.TYPE}" + nx_dict[path_type] = self.source_type.value + if nx_dict != {}: + nx_dict[f"{self.path}@NX_class"] = "NXsource" return nx_dict + + def _load(self, file_path: str, data_path: str, nexus_version: float) -> None: + nexus_paths = get_nexus_paths(nexus_version) + nexus_source_paths = nexus_paths.nx_source_paths + self.source_name = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_source_paths.NAME]), + ) + self.source_type = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_source_paths.TYPE]), + ) + + +class DefaultESRFSource(NXsource): + def __init__(self, node_name="source", parent=None) -> None: + super().__init__( + node_name=node_name, + parent=parent, + source_name="ESRF", + source_type=SourceType.SYNCHROTRON_X_RAY_SOURCE, + ) diff --git a/nxtomomill/nexus/nxtomo.py b/nxtomomill/nexus/nxtomo.py index c3b02f6be0d3b572332372bebabe98cfeff5d74b..5724891c1edd53ffd7cf1e3fbc75065a33d345be 100644 --- a/nxtomomill/nexus/nxtomo.py +++ b/nxtomomill/nexus/nxtomo.py @@ -28,95 +28,92 @@ __license__ = "MIT" __date__ = "03/02/2022" -from typing import Optional +import os +from typing import Optional, Union +from tomoscan.io import HDF5File +from tomoscan.nexus.paths.nxtomo import LATEST_VERSION as LATEST_NXTOMO_VERSION from .nxobject import NXobject, ElementWithUnit -from .nxdetector import NXdetector +from .nxinstrument import NXinstrument from .nxsample import NXsample +from .utils import get_data_and_unit, get_data from datetime import datetime -from nxtomomill.nexus.nxsource import NXsource, SourceType from silx.utils.proxy import docstring from tomoscan.unitsystem.metricsystem import EnergySI - -try: - from tomoscan.esrf.scan.hdf5scan import get_nexus_paths -except ImportError: - from tomoscan.esrf.hdf5scan import get_nexus_paths -from nxtomomill.converter.version import LATEST_VERSION as LATEST_NX_PATH_VERSION +from tomoscan.nexus.paths.nxtomo import get_paths as get_nexus_paths +from silx.io.url import DataUrl import logging +import h5py _logger = logging.getLogger(__name__) -class DefaultESRFSource(NXsource): - def __init__(self) -> None: - super().__init__(name="ESRF", type=SourceType.SYNCHROTRON_X_RAY_SOURCE) - - class NXtomo(NXobject): """ Class defining an NXTomo. His final goal is to save data to disk. """ - def __init__(self) -> None: - super().__init__() + def __init__( + self, node_name: Optional[str] = None, parent: Optional[NXobject] = None + ) -> None: + super().__init__(node_name=node_name, parent=parent) self._start_time = None self._end_time = None - self._source = DefaultESRFSource() - self._sample = NXsample() - self._detector = NXdetector() + self._instrument = NXinstrument(node_name="instrument", parent=self) + self._sample = NXsample(node_name="sample", parent=self) self._group_size = None self._energy = ElementWithUnit( default_unit=EnergySI.KILOELECTRONVOLT ) # energy in kev + self._title = None @property - def start_time(self) -> Optional[datetime]: + def start_time(self) -> Optional[Union[datetime, str]]: return self._start_time @start_time.setter - def start_time(self, start_time: Optional[datetime]): - if not isinstance(start_time, (type(None), datetime)): + def start_time(self, start_time: Optional[Union[datetime, str]]): + if not isinstance(start_time, (type(None), datetime, str)): raise TypeError( f"start_time is expected ot be an instance of datetime or None. Not {type(start_time)}" ) self._start_time = start_time @property - def end_time(self) -> Optional[datetime]: + def end_time(self) -> Optional[Union[datetime, str]]: return self._end_time @end_time.setter - def end_time(self, end_time: Optional[datetime]): - if not isinstance(end_time, (type(None), datetime)): + def end_time(self, end_time: Optional[Union[datetime, str]]): + if not isinstance(end_time, (type(None), datetime, str)): raise TypeError( f"end_time is expected ot be an instance of datetime or None. Not {type(end_time)}" ) self._end_time = end_time @property - def source(self) -> Optional[NXsource]: - return self._source + def title(self) -> Optional[str]: + return self._title - @source.setter - def source(self, source: Optional[NXsource]): - if not isinstance(source, (type(None), NXsource)): + @title.setter + def title(self, title: Optional[str]): + if not isinstance(title, (type(None), str)): raise TypeError( - f"source is expected ot be an instance of {NXsource} or None. Not {type(source)}" + f"title is expected ot be an instance of str or None. Not {type(title)}" ) - self._source = source + self._title = title @property - def detector(self) -> Optional[NXdetector]: - return self._detector + def instrument(self) -> Optional[NXinstrument]: + return self._instrument - @detector.setter - def detector(self, detector: Optional[NXdetector]) -> None: - if not isinstance(detector, (type(None), NXdetector)): + @instrument.setter + def instrument(self, instrument: Optional[NXinstrument]) -> None: + if not isinstance(instrument, (type(None), NXinstrument)): raise TypeError( - f"detector is expected ot be an instance of {NXdetector} or None. Not {type(detector)}" + f"instrument is expected ot be an instance of {NXinstrument} or None. Not {type(instrument)}" ) - self._detector = detector + self._instrument = instrument @property def sample(self) -> Optional[NXsample]: @@ -157,9 +154,14 @@ class NXtomo(NXobject): self._group_size = group_size @docstring(NXobject) - def to_nx_dict(self, nexus_path_version: Optional[float] = None) -> dict: - if nexus_path_version is None: - nexus_path_version = LATEST_NX_PATH_VERSION + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: + if data_path is None: + data_path = "" + nexus_paths = get_nexus_paths(nexus_path_version) nx_dict = {} @@ -170,47 +172,159 @@ class NXtomo(NXobject): else: _logger.info("no sample found. Won't be saved") - if self.detector is not None: - nx_dict.update( - self.detector.to_nx_dict(nexus_path_version=nexus_path_version) - ) - if self.detector.image_key is not None: - nx_dict[">data/data"] = f"/{self._entry}/instrument/detector/data" - nx_dict[ - ">data/image_key" - ] = f"/{self._entry}/instrument/detector/image_key" - nx_dict[ - ">data/image_key_control" - ] = f"/{self._entry}/instrument/detector/image_key_control" - nx_dict["instrument@NX_class"] = "NXinstrument" - - else: - _logger.info("no detector found. Won't be saved") - - if self.source is not None: + if self.instrument is not None: nx_dict.update( - self.source.to_nx_dict(nexus_path_version=nexus_path_version) + self.instrument.to_nx_dict(nexus_path_version=nexus_path_version) ) else: - _logger.info("no source found. Won't be saved") + _logger.info("no instrument found. Won't be saved") if self.start_time is not None: - nx_dict[nexus_paths.START_TIME_PATH] = self.start_time.isoformat() + path_start_time = f"{self.path}/{nexus_paths.START_TIME_PATH}" + if isinstance(self.start_time, datetime): + start_time = self.start_time.isoformat() + else: + start_time = self.start_time + nx_dict[path_start_time] = start_time if self.end_time is not None: - nx_dict[nexus_paths.END_TIME_PATH] = self.end_time.isoformat() + path_end_time = f"{self.path}/{nexus_paths.END_TIME_PATH}" + if isinstance(self.end_time, datetime): + end_time = self.end_time.isoformat() + else: + end_time = self.end_time + nx_dict[path_end_time] = end_time if self.group_size is not None: - nx_dict[nexus_paths.GRP_SIZE_ATTR] = self.group_size - if self.energy is not None: - nx_dict[nexus_paths.ENERGY_PATH] = self.energy.value - nx_dict["@".join([nexus_paths.ENERGY_PATH, "unit"])] = str(self.energy.unit) - if nexus_path_version > 1.0: - nx_dict[">beam"] = f"{self._entry}/{nexus_paths.ENERGY_PATH}" + path_grp_size = f"{self.path}/{nexus_paths.GRP_SIZE_ATTR}" + nx_dict[path_grp_size] = self.group_size + if self.energy.value is not None: + path_energy = f"{self.path}/{nexus_paths.ENERGY_PATH}" - if self._entry is not None: - nx_dict["@NX_class"] = "NXentry" + nx_dict[path_energy] = self.energy.value + nx_dict["@".join([path_energy, "unit"])] = str(self.energy.unit) + path_beam = f"{self.path}/{nexus_paths.BEAM_PATH}" + nx_dict["@".join([path_beam, "NX_class"])] = "NXbeam" - # save version used - if nexus_path_version is not None: - nx_dict["@version"] = nexus_path_version + if nexus_paths.VERSION > 1.0: + nx_dict[ + f">/{self.path}/beam/incident_energy" + ] = f"/{data_path}/{self.path}/{nexus_paths.ENERGY_PATH}" + if self.title is not None: + path_title = f"{self.path}/{nexus_paths.NAME_PATH}" + nx_dict[path_title] = self.title + + # create data group from symbolic links + + if self.instrument.detector.image_key is not None: + nx_dict[ + f">/{self.path}/data/image_key" + ] = f"/{data_path}/{self.instrument.detector.path}/{nexus_paths.nx_detector_paths.IMAGE_KEY}" + nx_dict[ + f">/{self.path}/data/image_key_control" + ] = f"/{data_path}/{self.instrument.detector.path}/{nexus_paths.nx_detector_paths.IMAGE_KEY_CONTROL}" + if self.instrument.detector.data is not None: + nx_dict[ + f">/{self.path}/data/data" + ] = f"/{data_path}/{self.instrument.detector.path}/{nexus_paths.nx_detector_paths.DATA}" + nx_dict[f"/{self.path}/data@NX_class"] = "NXdata" + nx_dict[f"/{self.path}/data@signal"] = "data" + nx_dict[f"/{self.path}@default"] = "data" + nx_dict[f"{self.path}/data@SILX_style/axis_scale_types"] = [ + "linear", + "linear", + ] + if self.sample.rotation_angle is not None: + nx_dict[ + f">/{self.path}/data/rotation_angle" + ] = f"/{data_path}/{self.sample.path}/{nexus_paths.nx_sample_paths.ROTATION_ANGLE}" + + if nx_dict != {}: + nx_dict[f"{self.path}@NX_class"] = "NXentry" + + nx_dict[f"{self.path}@version"] = nexus_paths.VERSION return nx_dict + + def detector_data_is_defined_by_url(self) -> bool: + return self._detector_data_is_defined_by_type(DataUrl) + + def detector_data_is_defined_by_virtual_source(self) -> bool: + return self._detector_data_is_defined_by_type(h5py.VirtualSource) + + def _detector_data_is_defined_by_type(self, type_): + return ( + self.instrument is not None + and self.instrument.detector is not None + and self.instrument.detector.data is not None + and isinstance(self.instrument.detector.data, (str, tuple)) + and isinstance(self.instrument.detector.data[0], type_) + ) + + def load( + self, file_path: str, data_path: str, detector_data_as="as_data_url" + ) -> NXobject: + """ + Load NXtomo instance from file_path and data_path + + :param str file_path: hdf5 file path containing the NXtomo + :param str data_path: location of the NXtomo + """ + possible_as_values = ("as_virtual_source", "as_data_url", "as_numpy_array") + if detector_data_as not in possible_as_values: + raise ValueError( + f"detector_data_as is expected to be in {possible_as_values} and not {detector_data_as}" + ) + + if not os.path.exists(file_path): + raise IOError(f"{file_path} does not exists") + with HDF5File(file_path, mode="r") as h5f: + if data_path not in h5f: + raise ValueError(f"{data_path} cannot be find in {file_path}") + root_node = h5f[data_path] + + if "version" in root_node.attrs: + nexus_version = root_node.attrs["version"] + else: + _logger.warning( + f"Unable to find nexus version associated with {data_path}@{file_path}" + ) + nexus_version = LATEST_NXTOMO_VERSION + nexus_paths = get_nexus_paths(nexus_version) + self.energy, self.energy.unit = get_data_and_unit( + file_path=file_path, + data_path="/".join([data_path, nexus_paths.ENERGY_PATH]), + default_unit="kev", + ) + start_time = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_paths.START_TIME_PATH]), + ) + try: + start_time = datetime.fromisoformat(start_time) + except Exception: + start_time = str(start_time) if start_time is not None else None + self.start_time = start_time + + end_time = get_data( + file_path=file_path, + data_path="/".join([data_path, nexus_paths.END_TIME_PATH]), + ) + try: + end_time = datetime.fromisoformat(end_time) + except Exception: + end_time = str(end_time) if end_time is not None else None + self.end_time = end_time + + self.title = get_data( + file_path=file_path, data_path="/".join([data_path, nexus_paths.NAME_PATH]) + ) + + self.sample._load( + file_path, "/".join([data_path, "sample"]), nexus_version=nexus_version + ) + self.instrument._load( + file_path, + "/".join([data_path, "instrument"]), + nexus_version=nexus_version, + detector_data_as=detector_data_as, + ) + return self diff --git a/nxtomomill/nexus/test/test_nxdetector.py b/nxtomomill/nexus/test/test_nxdetector.py index 0e39188c468f293f45720ad43a1496765fcb33a4..0d670d96f105e52e55443039448d065826e44ab9 100644 --- a/nxtomomill/nexus/test/test_nxdetector.py +++ b/nxtomomill/nexus/test/test_nxdetector.py @@ -29,19 +29,33 @@ __license__ = "MIT" __date__ = "04/02/2022" -from nxtomomill.nexus.nxdetector import NXdetector, FieldOfView +import tempfile +from nxtomomill.nexus.nxdetector import NXdetector, NXdetectorWithUnit, FieldOfView +from tomoscan.unitsystem.voltagesystem import VoltageSystem from nxtomomill.utils import ImageKey +from silx.io.url import DataUrl import pytest import numpy.random +import os +import h5py +from nxtomomill.utils.context import cwd_context def test_nx_detector(): """test creation and saving of an nxdetector""" - nx_detector = NXdetector() + nx_detector = NXdetector(expected_dim=(2, 3)) # check data with pytest.raises(TypeError): nx_detector.data = 12 + # if expected dims is not fulfill + with pytest.raises(ValueError): + nx_detector.data = numpy.random.random(100 * 100 * 5) + with pytest.raises(TypeError): + nx_detector.data = ( + 12, + 13, + ) nx_detector.data = numpy.random.random(100 * 100 * 5).reshape(5, 100, 100) # check image key control @@ -80,3 +94,177 @@ def test_nx_detector(): nx_detector.estimated_cor_from_motor = 0.5 assert isinstance(nx_detector.to_nx_dict(), dict) + + +def test_nx_detector_with_unit(): + diode = NXdetectorWithUnit( + node_name="diode", + expected_dim=(1,), + default_unit=VoltageSystem.VOLT, + ) + with pytest.raises(ValueError): + diode.data = numpy.arange(10 * 10).reshape([10, 10]) + with pytest.raises(TypeError): + diode.data = [10, 12] + with pytest.raises(TypeError): + diode.data = "test" + diode.data = None + diode.data = numpy.random.random(12) + diode.data = (DataUrl(),) + + +def test_nx_detector_with_virtual_source(): + """Insure detector data can be write from Virtual sources""" + cwd = os.getcwd() + with tempfile.TemporaryDirectory() as tmp_folder: + # create virtual dataset + n_base_raw_dataset = 5 + n_z, n_y, n_x = 4, 100, 100 + base_raw_dataset_shape = (n_z, n_y, n_x) + n_base_raw_dataset_elmts = n_z * n_y * n_x + + v_sources = [] + + raw_files = [ + os.path.join(tmp_folder, f"raw_file_{i_file}.hdf5") + for i_file in range(n_base_raw_dataset) + ] + for i_raw_file, raw_file in enumerate(raw_files): + with h5py.File(raw_file, mode="w") as h5f: + h5f["data"] = numpy.arange( + start=n_base_raw_dataset_elmts * i_raw_file, + stop=n_base_raw_dataset_elmts * (i_raw_file + 1), + ).reshape(base_raw_dataset_shape) + v_sources.append(h5py.VirtualSource(h5f["data"])) + + nx_detecteur = NXdetector() + nx_detecteur.data = v_sources + + detector_file = os.path.join(tmp_folder, "detector_file.hdf5") + nx_detecteur.save(file_path=detector_file, data_path="/") + + # check the virtual dataset has been properly createde and linked + with h5py.File(detector_file, mode="r") as h5f_master: + dataset = h5f_master["/detector/data"] + assert dataset.is_virtual + for i_raw_file, raw_file in enumerate(raw_files): + with h5py.File(raw_file, mode="r") as h5f_raw: + numpy.testing.assert_array_equal( + dataset[i_raw_file * n_z : (i_raw_file + 1) * n_z], + h5f_raw["data"], + ) + # check attributes have beem rewrite as expected + assert "interpretation" in dataset.attrs + + # check virtual dataset is composed of relative links + for vs_info in dataset.virtual_sources(): + assert vs_info.file_name.startswith("./") + assert cwd == os.getcwd() + + +def test_nx_detector_with_local_urls(): + """Insure detector data can be write from DataUrl linking to local dataset (in the same file)""" + + cwd = os.getcwd() + n_base_dataset = 3 + n_z, n_y, n_x = 2, 10, 20 + base_dataset_shape = (n_z, n_y, n_x) + n_base_dataset_elmts = n_z * n_y * n_x + urls = [] + + with tempfile.TemporaryDirectory() as tmp_folder: + master_file = os.path.join(tmp_folder, "master_file.hdf5") + with h5py.File(master_file, mode="a") as h5f: + for i in range(n_base_dataset): + data_path = f"/data_{i}" + h5f[data_path] = numpy.arange( + start=n_base_dataset_elmts * i, + stop=n_base_dataset_elmts * (i + 1), + ).reshape(base_dataset_shape) + urls.append( + DataUrl( + file_path=master_file, + data_path=data_path, + scheme="silx", + ) + ) + nx_detecteur = NXdetector() + nx_detecteur.data = urls + nx_detecteur.save(file_path=master_file, data_path="/") + + # check the virtual dataset has been properly createde and linked + with h5py.File(master_file, mode="r") as h5f_master: + dataset = h5f_master["/detector/data"] + assert dataset.is_virtual + for i in range(n_base_dataset): + numpy.testing.assert_array_equal( + dataset[i * n_z : (i + 1) * n_z], + numpy.arange( + start=n_base_dataset_elmts * i, + stop=n_base_dataset_elmts * (i + 1), + ).reshape(base_dataset_shape), + ) + # check virtual dataset is composed of relative links + for vs_info in dataset.virtual_sources(): + assert vs_info.file_name.startswith("./") + assert cwd == os.getcwd() + + +def test_nx_detector_with_external_urls(): + """Insure detector data can be write from DataUrl linking to external dataset""" + cwd = os.getcwd() + with tempfile.TemporaryDirectory() as tmp_folder: + # create virtual dataset + n_base_raw_dataset = 5 + n_z, n_y, n_x = 4, 100, 100 + base_raw_dataset_shape = (n_z, n_y, n_x) + n_base_raw_dataset_elmts = n_z * n_y * n_x + + urls = [] + + raw_files = [ + os.path.join(tmp_folder, f"raw_file_{i_file}.hdf5") + for i_file in range(n_base_raw_dataset) + ] + for i_raw_file, raw_file in enumerate(raw_files): + with h5py.File(raw_file, mode="w") as h5f: + h5f["data"] = numpy.arange( + start=n_base_raw_dataset_elmts * i_raw_file, + stop=n_base_raw_dataset_elmts * (i_raw_file + 1), + ).reshape(base_raw_dataset_shape) + # provide one file path each two as an absolue path + if i_raw_file % 2 == 0: + file_path = os.path.abspath(raw_file) + else: + file_path = os.path.relpath(raw_file, tmp_folder) + urls.append( + DataUrl( + file_path=file_path, + data_path="data", + scheme="silx", + ) + ) + + nx_detecteur = NXdetector() + nx_detecteur.data = urls + + detector_file = os.path.join(tmp_folder, "detector_file.hdf5") + # needed as we provide some link with relative path + with cwd_context(tmp_folder): + nx_detecteur.save(file_path=detector_file) + + # check the virtual dataset has been properly createde and linked + with h5py.File(detector_file, mode="r") as h5f_master: + dataset = h5f_master["/detector/data"] + assert dataset.is_virtual + for i_raw_file, raw_file in enumerate(raw_files): + with h5py.File(raw_file, mode="r") as h5f_raw: + numpy.testing.assert_array_equal( + dataset[i_raw_file * n_z : (i_raw_file + 1) * n_z], + h5f_raw["data"], + ) + # check virtual dataset is composed of relative links + for vs_info in dataset.virtual_sources(): + assert vs_info.file_name.startswith("./") + + assert cwd == os.getcwd() diff --git a/nxtomomill/nexus/test/test_nxinstrument.py b/nxtomomill/nexus/test/test_nxinstrument.py new file mode 100644 index 0000000000000000000000000000000000000000..9416d3d176d3aae41d55a8cbf5e6a05faadc7971 --- /dev/null +++ b/nxtomomill/nexus/test/test_nxinstrument.py @@ -0,0 +1,64 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2022 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "10/02/2022" + + +from nxtomomill.nexus.nxdetector import NXdetector +from nxtomomill.nexus.nxinstrument import NXinstrument +from nxtomomill.nexus.nxsource import NXsource, DefaultESRFSource +import pytest + + +def test_nx_instrument(): + """test creation and saving of an nxinstrument""" + nx_instrument = NXinstrument() + + # check data + with pytest.raises(TypeError): + nx_instrument.detector = 12 + nx_instrument.detector = NXdetector(node_name="test") + + with pytest.raises(TypeError): + nx_instrument.diode = 12 + nx_instrument.diode = NXdetector(node_name="test 2") + + with pytest.raises(TypeError): + nx_instrument.source = 12 + nx_instrument.source = DefaultESRFSource() + + with pytest.raises(TypeError): + nx_instrument.diode = NXsource(node_name="my source") + nx_instrument.diode = NXdetector(node_name="det34") + + assert isinstance(nx_instrument.to_nx_dict(), dict) + + with pytest.raises(TypeError): + nx_instrument.name = 12 + nx_instrument.name = "test name" + assert nx_instrument.name == "test name" diff --git a/nxtomomill/nexus/test/test_nxobject.py b/nxtomomill/nexus/test/test_nxobject.py index 6ce10a25cd86fc8d90f6025527462c8c24863e3d..eacccdd2af950eaa39951df14fa790782861991f 100644 --- a/nxtomomill/nexus/test/test_nxobject.py +++ b/nxtomomill/nexus/test/test_nxobject.py @@ -42,34 +42,47 @@ import os class test_nx_object: """Tets API of the nx object""" - nx_object = NXobject() + with pytest.raises(TypeError): + NXobject(node_name=12) + with pytest.raises(TypeError): + NXobject(node_name="test", parent=12) + + nx_object = NXobject(node_name="NXobject") with pytest.raises(NotImplementedError): nx_object.to_nx_dict(nexus_path_version=1.0) + assert nx_object.is_root is True + + with pytest.raises(TypeError): + nx_object.node_name = 12 class MyNXObject(NXobject): - def to_nx_dict(self, nexus_path_version: Optional[float] = None) -> dict: + def to_nx_dict( + self, + nexus_path_version: Optional[float] = None, + data_path: Optional[str] = None, + ) -> dict: return { - "test": "toto", + f"{self.path}/test": "toto", } - my_nx_object = MyNXObject() + my_nx_object = MyNXObject(node_name="NxObject2") with TemporaryDirectory() as folder: file_path = os.path.join(folder, "my_nexus.nx") assert not os.path.exists(file_path) - my_nx_object.save(file_path=file_path, entry="entry", nexus_path_version=1.0) + my_nx_object.save(file_path=file_path, data_path="/", nexus_path_version=1.0) assert os.path.exists(file_path) with pytest.raises(KeyError): my_nx_object.save( file_path=file_path, - entry="entry", + data_path="/", nexus_path_version=1.0, overwrite=False, ) my_nx_object.save( - file_path=file_path, entry="entry", nexus_path_version=1.0, overwrite=True + file_path=file_path, data_path="/", nexus_path_version=1.0, overwrite=True ) @@ -98,3 +111,8 @@ def test_ElementWithUnit(): elmt.unit = "minute" elmt.si_value == 8.0 / 60.0 str(elmt) + + with pytest.raises(ValueError): + elmt.unit = "not minute" + with pytest.raises(TypeError): + elmt.unit = 123 diff --git a/nxtomomill/nexus/test/test_nxsource.py b/nxtomomill/nexus/test/test_nxsource.py index a5b35b0811c8ffc44ee1a15defb41628e2490683..344eb03ef3be38accc56f4deb4ceea628856916b 100644 --- a/nxtomomill/nexus/test/test_nxsource.py +++ b/nxtomomill/nexus/test/test_nxsource.py @@ -37,11 +37,14 @@ def test_nx_source(): """test creation and saving of an nxsource""" nx_source = NXsource() with pytest.raises(TypeError): - nx_source.name = 12 - nx_source.name = "my source" + nx_source.source_name = 12 + nx_source.source_name = "my source" with pytest.raises(ValueError): - nx_source.type = "toto" - nx_source.type = "Synchrotron X-ray Source" + nx_source.source_type = "toto" + nx_source.source_type = "Synchrotron X-ray Source" + str(nx_source) + nx_source.source_type = None + str(nx_source) assert isinstance(nx_source.to_nx_dict(), dict) diff --git a/nxtomomill/nexus/test/test_nxtomo.py b/nxtomomill/nexus/test/test_nxtomo.py index f4ec2ea47fad96d1d770adaf4e669a2f644d8a30..8a35b42f9a3b6fe4cf9531a72ab06129e2156fcc 100644 --- a/nxtomomill/nexus/test/test_nxtomo.py +++ b/nxtomomill/nexus/test/test_nxtomo.py @@ -36,6 +36,7 @@ from nxtomomill.nexus.nxtomo import NXtomo from datetime import datetime from tomoscan.esrf import HDF5TomoScan from tomoscan.validator import ReconstructionValidator +from silx.io.url import DataUrl try: from tomoscan.esrf.scan.hdf5scan import ImageKey @@ -49,7 +50,7 @@ nexus_path_versions = (1.1, 1.0, None) @pytest.mark.parametrize("nexus_path_version", nexus_path_versions) def test_nx_tomo(nexus_path_version): - nx_tomo = NXtomo() + nx_tomo = NXtomo(node_name="") # check start time with pytest.raises(TypeError): @@ -61,17 +62,13 @@ def test_nx_tomo(nexus_path_version): nx_tomo.end_time = 12 nx_tomo.end_time = datetime.now() - # check source - with pytest.raises(TypeError): - nx_tomo.source = "tata" - # check sample with pytest.raises(TypeError): nx_tomo.sample = "tata" # check detector with pytest.raises(TypeError): - nx_tomo.detector = "tata" + nx_tomo.instrument.detector = "tata" # check energy with pytest.raises(TypeError): @@ -83,13 +80,22 @@ def test_nx_tomo(nexus_path_version): nx_tomo.group_size = "tata" nx_tomo.group_size = 3 + # check title + with pytest.raises(TypeError): + nx_tomo.title = 12 + nx_tomo.title = "title" + + # check instrument + with pytest.raises(TypeError): + nx_tomo.instrument = "test" + # create detector for test projections = numpy.random.random(100 * 100 * 8).reshape([8, 100, 100]) flats_1 = numpy.random.random(100 * 100 * 2).reshape([2, 100, 100]) darks = numpy.random.random(100 * 100 * 3).reshape([3, 100, 100]) flats_2 = numpy.random.random(100 * 100 * 2).reshape([2, 100, 100]) alignment = numpy.random.random(100 * 100 * 1).reshape([1, 100, 100]) - nx_tomo.detector.data = numpy.concatenate( + nx_tomo.instrument.detector.data = numpy.concatenate( [ darks, flats_1, @@ -99,7 +105,7 @@ def test_nx_tomo(nexus_path_version): ] ) - nx_tomo.detector.image_key_control = numpy.concatenate( + nx_tomo.instrument.detector.image_key_control = numpy.concatenate( [ [ImageKey.DARK_FIELD] * 3, [ImageKey.FLAT_FIELD] * 2, @@ -109,10 +115,12 @@ def test_nx_tomo(nexus_path_version): ] ) - nx_tomo.detector.x_pixel_size = nx_tomo.detector.y_pixel_size = 1e-7 - nx_tomo.detector.distance = 0.2 - nx_tomo.detector.field_of_view = FieldOfView.HALF - nx_tomo.detector.count_time = numpy.concatenate( + nx_tomo.instrument.detector.x_pixel_size = ( + nx_tomo.instrument.detector.y_pixel_size + ) = 1e-7 + nx_tomo.instrument.detector.distance = 0.2 + nx_tomo.instrument.detector.field_of_view = FieldOfView.HALF + nx_tomo.instrument.detector.count_time = numpy.concatenate( [ [0.2] * 3, # darks [0.1] * 2, # flats 1 @@ -138,12 +146,21 @@ def test_nx_tomo(nexus_path_version): nx_tomo.sample.x_translation = [0.6] * n_frames nx_tomo.sample.y_translation = [0.2] * n_frames nx_tomo.sample.z_translation = [0.1] * n_frames + assert nx_tomo.is_root is True + assert nx_tomo.instrument.is_root is False + assert ( + nx_tomo.root_path + == nx_tomo.instrument.root_path + == nx_tomo.instrument.detector.root_path + ) with TemporaryDirectory() as folder: file_path = os.path.join(folder, "nexus_file.hdf5") nx_tomo.save( - file_path=file_path, entry="entry", nexus_path_version=nexus_path_version + file_path=file_path, + data_path="entry", + nexus_path_version=nexus_path_version, ) assert os.path.exists(file_path) @@ -174,3 +191,93 @@ def test_nx_tomo(nexus_path_version): # check values validity using the validator validator = ReconstructionValidator(scan=scan) assert validator.is_valid() + + # try to load it from the disk + loaded_nx_tomo = NXtomo("test").load(file_path=file_path, data_path="entry") + assert isinstance(loaded_nx_tomo, NXtomo) + assert loaded_nx_tomo.energy.value == nx_tomo.energy.value + assert loaded_nx_tomo.energy.unit == nx_tomo.energy.unit + assert loaded_nx_tomo.start_time == nx_tomo.start_time + assert loaded_nx_tomo.end_time == nx_tomo.end_time + assert ( + loaded_nx_tomo.instrument.detector.x_pixel_size.value + == nx_tomo.instrument.detector.x_pixel_size.value + ) + assert ( + loaded_nx_tomo.instrument.detector.x_pixel_size.unit + == nx_tomo.instrument.detector.x_pixel_size.unit + ) + assert ( + loaded_nx_tomo.instrument.detector.y_pixel_size.value + == nx_tomo.instrument.detector.y_pixel_size.value + ) + assert ( + loaded_nx_tomo.instrument.detector.field_of_view + == nx_tomo.instrument.detector.field_of_view + ) + numpy.testing.assert_array_equal( + loaded_nx_tomo.instrument.detector.count_time.value, + nx_tomo.instrument.detector.count_time.value, + ) + assert ( + loaded_nx_tomo.instrument.detector.count_time.unit + == nx_tomo.instrument.detector.count_time.unit + ) + assert ( + loaded_nx_tomo.instrument.detector.distance.value + == nx_tomo.instrument.detector.distance.value + ) + assert ( + loaded_nx_tomo.instrument.detector.distance.unit + == nx_tomo.instrument.detector.distance.unit + ) + assert ( + loaded_nx_tomo.instrument.detector.estimated_cor_from_motor + == nx_tomo.instrument.detector.estimated_cor_from_motor + ) + + numpy.testing.assert_array_equal( + loaded_nx_tomo.instrument.detector.image_key_control, + nx_tomo.instrument.detector.image_key_control, + ) + numpy.testing.assert_array_equal( + loaded_nx_tomo.instrument.detector.image_key, + nx_tomo.instrument.detector.image_key, + ) + + assert loaded_nx_tomo.sample.name == nx_tomo.sample.name + assert loaded_nx_tomo.sample.rotation_angle is not None + numpy.testing.assert_array_almost_equal( + loaded_nx_tomo.sample.rotation_angle, nx_tomo.sample.rotation_angle + ) + numpy.testing.assert_array_almost_equal( + loaded_nx_tomo.sample.x_translation.value, + nx_tomo.sample.x_translation.value, + ) + numpy.testing.assert_array_almost_equal( + loaded_nx_tomo.sample.y_translation.value, + nx_tomo.sample.y_translation.value, + ) + numpy.testing.assert_array_almost_equal( + loaded_nx_tomo.sample.z_translation.value, + nx_tomo.sample.z_translation.value, + ) + + loaded_nx_tomo = NXtomo("test").load( + file_path=file_path, data_path="entry", detector_data_as="as_numpy_array" + ) + numpy.testing.assert_array_almost_equal( + loaded_nx_tomo.instrument.detector.data, + nx_tomo.instrument.detector.data, + ) + loaded_nx_tomo = NXtomo("test").load( + file_path=file_path, data_path="entry", detector_data_as="as_data_url" + ) + assert isinstance(loaded_nx_tomo.instrument.detector.data[0], DataUrl) + with pytest.raises(ValueError): + # check an error is raise because the dataset is not virtual + loaded_nx_tomo = NXtomo("test").load( + file_path=file_path, + data_path="entry", + detector_data_as="as_virtual_source", + ) diff --git a/nxtomomill/nexus/utils.py b/nxtomomill/nexus/utils.py index 65877f92ea49c11d5372ee036e753752a5c2ec96..03a069118869e4eaf8849631f3ab0fd87ea5221b 100644 --- a/nxtomomill/nexus/utils.py +++ b/nxtomomill/nexus/utils.py @@ -28,8 +28,12 @@ __authors__ = ["H. Payno"] __license__ = "MIT" __date__ = "03/02/2022" + from typing import Iterable +from tomoscan.io import HDF5File import numpy +import h5py +from silx.io.utils import h5py_read_dataset def cast_and_check_array_1D(array, array_name): @@ -42,3 +46,89 @@ def cast_and_check_array_1D(array, array_name): if array is not None and array.ndim > 1: raise ValueError(f"{array_name} is expected to be 0 or 1d not {array.ndim}") return array + + +def create_nx_data_group(file_path: str, entry_path: str, axis_scale: list): + """ + Create the 'Nxdata' group at entry level with soft links on the NXDetector + and NXsample. + + :param file_path: + :param entry_path: + :param axis_scale: + :return: + """ + if not isinstance(file_path, str): + raise TypeError("file_path is expected to be a file") + if not isinstance(entry_path, str): + raise TypeError("entry_path is expected to be a file") + if not isinstance(axis_scale, Iterable): + raise TypeError("axis_scale is expected to be an Iterable") + + with HDF5File(filename=file_path, mode="a") as h5f: + entry_group = h5f[entry_path] + + nx_data_grp = entry_group.require_group("data") + # link detector datasets: + if not entry_path.startswith("/"): + entry_path = "/" + entry_path + for dataset in ("data", "image_key", "image_key_control"): + dataset_path = "/".join((entry_path, "instrument", "detector", dataset)) + nx_data_grp[dataset] = h5py.SoftLink(dataset_path) + # link rotation angle + nx_data_grp["rotation_angle"] = h5py.SoftLink( + "/".join((entry_path, "sample", "rotation_angle")) + ) + + # write NX attributes + nx_data_grp.attrs["NX_class"] = "NXdata" + nx_data_grp.attrs["signal"] = "data" + nx_data_grp.attrs["SILX_style/axis_scale_types"] = axis_scale + nx_data_grp["data"].attrs["interpretation"] = "image" + + +def link_nxbeam_to_root(file_path, entry_path): + """ + Create the 'Nxdata' group at entry level with soft links on the NXDetector + and NXsample. + + :param file_path: + :param entry_path: + :return: + """ + if not isinstance(file_path, str): + raise TypeError("file_path is expected to be a file") + if not isinstance(entry_path, str): + raise TypeError("entry_path is expected to be a file") + + if not entry_path.startswith("/"): + entry_path = "/" + entry_path + with HDF5File(filename=file_path, mode="a") as h5f: + entry_group = h5f[entry_path] + entry_group["beam"] = h5py.SoftLink( + "/".join((entry_path, "instrument", "beam")) + ) + + +def get_data_and_unit(file_path, data_path, default_unit): + with HDF5File(file_path, mode="r") as h5f: + if data_path in h5f and isinstance(h5f[data_path], h5py.Dataset): + dataset = h5f[data_path] + unit = None + if "unit" in dataset.attrs: + unit = dataset.attrs["unit"] + elif "units" in dataset.attrs: + unit = dataset.attrs["units"] + else: + unit = default_unit + return h5py_read_dataset(dataset), unit + else: + return None, default_unit + + +def get_data(file_path, data_path): + with HDF5File(file_path, mode="r") as h5f: + if data_path in h5f: + return h5py_read_dataset(h5f[data_path]) + else: + return None diff --git a/nxtomomill/plugins.py b/nxtomomill/plugins.py index c495b321dfefb6c03c7441fd2fa956260db677fc..54372c1eecb79743a4c96e2ce2487f1b46378fe4 100644 --- a/nxtomomill/plugins.py +++ b/nxtomomill/plugins.py @@ -35,6 +35,8 @@ from importlib.machinery import SourceFileLoader import logging import inspect +from nxtomomill.nexus.nxtomo import NXtomo + _logger = logging.getLogger(__name__) @@ -150,22 +152,6 @@ class HDF5Plugin(_PluginBase): :param positioner_keys: which information you will need for the plugin. :param instrument_keys: which information to be read from the "instrument" dataset of the input bliss scan. - :param sample_datasets_created: name of the field created by the plugin in - the sample group. Allows those to be skip - by nxtomomill and ignore some warnings. - :type sample_datasets_created: Union[None, Iterable] - :param detector_datasets_created: name of the field created by the plugin - in the detector group. Allows those to be - skip by nxtomomill and ignore some - warnings. - :type detector_datasets_created: Union[None, Iterable] - :param beam_datasets_created: name of the field created by the plugin in - the beam group. Allows those to be skip by - nxtomomill and ignore some warnings. - :type beam_datasets_created: Union[None, Iterable] - :warning: unit are not managed by the plugin. So you will retrieve - raw data and should know the unit of the data you are reading - """ def __init__( @@ -173,9 +159,6 @@ class HDF5Plugin(_PluginBase): name, positioner_keys=tuple(), instrument_keys=tuple(), - sample_datasets_created=None, - detector_datasets_created=None, - beam_datasets_created=None, ): super().__init__(name) """As the sequence is splitted into several scan we need to parse @@ -194,9 +177,6 @@ class HDF5Plugin(_PluginBase): if not isinstance(key, str): raise TypeError("resource names should be strings") self.__instrument_info[key] = None - self._sample_datasets_created = sample_datasets_created - self._detector_datasets_created = detector_datasets_created - self._beam_datasets_created = beam_datasets_created def process(self) -> list: """ @@ -243,14 +223,5 @@ class HDF5Plugin(_PluginBase): """ return self.__instrument_info[name] - @property - def sample_datasets_created(self): - return self._sample_datasets_created - - @property - def detector_datasets_created(self): - return self._detector_datasets_created - - @property - def beam_datasets_created(self): - return self._beam_datasets_created + def update_nx_tomo(self, nx_tomo: NXtomo): + raise NotImplementedError("Base class") diff --git a/nxtomomill/settings.py b/nxtomomill/settings.py index 8928d0424736a1578edfacfb5f54ad68037a132a..efb097cb780346d2a4ebe8326e3b051db389ae94 100644 --- a/nxtomomill/settings.py +++ b/nxtomomill/settings.py @@ -74,7 +74,7 @@ class Tomo: ACQ_EXPO_TIME_KEYS = ("acq_expo_time",) INIT_TITLES = ( - "pcotomo:basic", + # "pcotomo:basic", seems to not be used anymore "tomo:basic", "tomo:fullturn", "sequence_of_scans", @@ -85,7 +85,10 @@ class Tomo: initialization scan""" ZSERIE_INIT_TITLES = ("tomo:zseries",) - """specific titles for zserie. Will extend H5_INIT_TITLES""" + """specific titles for zserie. Will extend DEFAULT_SCAN_TITLES""" + + PCOTOMO_INIT_TITLES = ("tomo:pcotomo",) + """specific titles for pcotomo. Will extend DEFAULT_SCAN_TITLES""" DARK_TITLES = ("dark images", "dark") """if a scan starts by one of those string then is considered as diff --git a/nxtomomill/test/test_plugins.py b/nxtomomill/test/test_plugins.py index e62e62ed9e316374e7d3437060969bb840deccb9..6b67b487656b54ba7790f0f03ddf28d707192bc3 100644 --- a/nxtomomill/test/test_plugins.py +++ b/nxtomomill/test/test_plugins.py @@ -88,25 +88,21 @@ class PluginTest(HDF5Plugin): super().__init__( name="plugin test", positioner_keys=("sy", ), - sample_datasets_created=("y_translation",), ) - def write(self, root_node, sample_node, detector_node, beam_node): + def update_nx_tomo(self, nx_tomo): try: sy = self.get_positioner_info("sy") except Exception as e: print(e) else: try: - sample_node["y_translation"] = numpy.ones(numpy.array(sy).shape) + nx_tomo.sample.y_translation = numpy.ones(numpy.array(sy).shape) except Exception as e: print(e) """ ) with EnvContext(NXTOMOMILL_PLUGINS_DIR=plugin_folder): - import shutil - - shutil.copyfile(my_plug_file, "plug_test.py") assert len(plugins.get_plugins_instances_frm_env_var()) == 1 with tempfile.TemporaryDirectory() as folder: nx_file_path = os.path.join(folder, "acquisition.nx") @@ -123,11 +119,9 @@ class PluginTest(HDF5Plugin): ) h5_file_path = bliss_mock.samples[0].sample_file - shutil.copyfile(h5_file_path, "raw.h5") assert not os.path.exists(nx_file_path), "outputfile exists already" main(["h52nx", h5_file_path, nx_file_path, "--single-file"]) assert os.path.exists(nx_file_path), "outputfile doesn't exists" - # TODO: insure y_translation has been correctly saved with h5py.File(nx_file_path, mode="r") as h5f: assert h5f["entry0000/sample/y_translation"][()].min() == 1 assert h5f["entry0000/sample/y_translation"][()].max() == 1 diff --git a/nxtomomill/test/utils/bliss.py b/nxtomomill/test/utils/bliss.py index 3729626e55b4b25c84a9e2203b9b5b138a4646e2..2f1800c7fa4198ea5b513d4773718ff87c283ecf 100644 --- a/nxtomomill/test/utils/bliss.py +++ b/nxtomomill/test/utils/bliss.py @@ -65,11 +65,22 @@ class MockBlissAcquisition: acqui_type="basic", z_values=None, y_rot=False, + nb_loop=None, + nb_tomo=None, + with_rotation_motor_info=True, ): self.__folder = output_dir if not os.path.exists(output_dir): os.makedirs(output_dir) self.__proposal_file = os.path.join(self.__folder, "ihproposal_file.h5") + if acqui_type not in ("pcotomo"): + if nb_loop is not None or nb_tomo is not None: + raise ValueError( + "nb_loop and nb_tomo are only handled by acqui_type: `pcotomo`" + ) + else: + if nb_loop is None or nb_tomo is None: + raise ValueError("nb_loop and nb_tomo must be provided") # create sample self.__samples = [] @@ -89,6 +100,22 @@ class MockBlissAcquisition: with_nx_detector_attr=with_nx_detector_attr, detector_name=detector_name, y_rot=y_rot, + with_rotation_motor_info=with_rotation_motor_info, + ) + elif acqui_type == "pcotomo": + acqui_tomo = _BlissPCOTomo( + sample_dir=sample_dir, + sample_file=sample_file, + n_sequence=n_sequence, + n_scan_per_sequence=n_scan_per_sequence, + n_darks=n_darks, + n_flats=n_flats, + with_nx_detector_attr=with_nx_detector_attr, + detector_name=detector_name, + y_rot=y_rot, + with_rotation_motor_info=with_rotation_motor_info, + nb_loop=nb_loop, + nb_tomo=nb_tomo, ) elif acqui_type == "zseries": if z_values is None: @@ -104,6 +131,7 @@ class MockBlissAcquisition: detector_name=detector_name, z_values=z_values, y_rot=y_rot, + with_rotation_motor_info=with_rotation_motor_info, ) elif acqui_type == "xrd-ct": acqui_tomo = _BlissXRD_CT( @@ -116,6 +144,7 @@ class MockBlissAcquisition: with_nx_detector_attr=with_nx_detector_attr, detector_name=detector_name, y_rot=y_rot, + with_rotation_motor_info=with_rotation_motor_info, ) else: raise NotImplementedError("") @@ -152,6 +181,7 @@ class _BlissSample: detector_name, with_nx_detector_attr=True, y_rot=False, + with_rotation_motor_info=True, ): self._with_nx_detector_attr = with_nx_detector_attr self._sample_dir = sample_dir @@ -165,11 +195,12 @@ class _BlissSample: self._detector_name = detector_name self._det_width = 64 self._det_height = 64 - self._n_frame_per_scan = 10 + self._tomo_n = 10 self._energy = 19.0 self._distance = 1.0 self._pixel_size = (0.0065, 0.0066) self._y_rot = y_rot + self._with_rotation_motor_info = with_rotation_motor_info for i_sequence in range(n_sequence): self.add_sequence() @@ -201,9 +232,7 @@ class _BlissSample: seq_node.require_group("instrument/positioners") # write energy seq_node["technique/scan/energy"] = self._energy - seq_node["technique/scan/tomo_n"] = ( - self._n_frame_per_scan * self._n_scan_per_seq - ) + seq_node["technique/scan/tomo_n"] = self._tomo_n * self._n_scan_per_seq seq_node["technique/scan/sample_detector_distance"] = self._distance seq_node["technique/optic/sample_pixel_size"] = numpy.asarray( self._pixel_size @@ -232,7 +261,13 @@ class _BlissSample: def get_next_group_name(seq_ini_index, scan_idx): return str(scan_idx) + ".1" - def add_scan(self, scan_type, seq_ini_index, z_value, skip_title=False): + def add_scan( + self, scan_type, seq_ini_index, z_value, skip_title=False, nb_loop=1, nb_tomo=1 + ): + """ + :param int nb_loop: number of loop in pcotomo use case. Else must be 1 + :param int nb_tomo: number of tomography done in pcotomo 'per iteration' use case. Else must be 1 + """ scan_idx = self.get_next_free_index() scan_name = str(scan_idx).zfill(4) scan_path = os.path.join(self.path, scan_name) @@ -250,12 +285,16 @@ class _BlissSample: # write data data = ( numpy.random.random( - self._det_height * self._det_width * self._n_frame_per_scan + self._det_height + * self._det_width + * self._tomo_n + * nb_loop + * nb_tomo ) * 256 ) data = data.reshape( - self._n_frame_per_scan, self._det_height, self._det_width + (self._tomo_n * nb_tomo * nb_loop), self._det_height, self._det_width ) data = data.astype(numpy.uint16) det_path_1 = "/".join(("instrument", self._detector_name)) @@ -271,13 +310,15 @@ class _BlissSample: # write rotation angle value and translations hrsrot_pos = seq_node.require_group("/".join(("instrument", "positioners"))) hrsrot_pos["hrsrot"] = numpy.random.randint( - low=0.0, high=360, size=self._n_frame_per_scan - ) - hrsrot_pos["sx"] = numpy.array( - numpy.random.random(size=self._n_frame_per_scan) + low=0.0, high=360, size=self._tomo_n ) - hrsrot_pos["sy"] = numpy.random.random(size=self._n_frame_per_scan) - hrsrot_pos["sz"] = numpy.asarray([z_value] * self._n_frame_per_scan) + hrsrot_pos["sx"] = numpy.array(numpy.random.random(size=self._tomo_n)) + hrsrot_pos["sy"] = numpy.random.random(size=self._tomo_n) + hrsrot_pos["sz"] = numpy.asarray([z_value] * self._tomo_n) + if self._with_rotation_motor_info: + scan_node = seq_node.require_group("technique/scan") + scan_node["motor"] = ("rotation", "hrsrot", "srot") + # motors are saved as (["motor_name_1", "dataset_name", "alias_name", "motor_name_2", ...]) def add_sequence(self): """Add a sequence to the bliss file""" @@ -302,6 +343,10 @@ class _BlissSample: def n_darks(self): return self._n_darks + @property + def with_rotation_motor_info(self): + return self._with_rotation_motor_info + class _BlissScan: """ @@ -339,6 +384,87 @@ class _BlissBasicTomo(_BlissSample): ) +class _BlissPCOTomo(_BlissSample): + def __init__( + self, + sample_dir, + sample_file, + n_sequence, + n_scan_per_sequence, + n_darks, + n_flats, + detector_name, + with_nx_detector_attr=True, + y_rot=False, + with_rotation_motor_info=True, + nb_loop=1, + nb_tomo=1, + ): + self.nb_loop = nb_loop + self.nb_tomo = nb_tomo + super().__init__( + sample_dir, + sample_file, + n_sequence, + n_scan_per_sequence, + n_darks, + n_flats, + detector_name, + with_nx_detector_attr, + y_rot, + with_rotation_motor_info=with_rotation_motor_info, + ) + + def get_main_entry_title(self): + return "tomo:pcotomo" + + def add_sequence(self): + # reserve the index for the 'initialization' sequence. No scan folder + # will be created for this one. + seq_ini_index = self.get_next_free_index() + + self.create_entry_and_technique(seq_ini_index=seq_ini_index) + + # start dark + if self.n_darks > 0: + self.add_scan(scan_type="dark", seq_ini_index=seq_ini_index, z_value=1) + + # start flat + if self._n_flats > 0: + self.add_scan(scan_type="flat", seq_ini_index=seq_ini_index, z_value=1) + + for i_proj_seq in range(self._n_scan_per_seq): + self.add_scan( + scan_type="projection", + seq_ini_index=seq_ini_index, + z_value=1, + nb_loop=self.nb_loop, + nb_tomo=self.nb_tomo, + ) + + # end flat + if self._n_flats > 0: + self.add_scan(scan_type="flat", seq_ini_index=seq_ini_index, z_value=1) + + def add_scan( + self, scan_type, seq_ini_index, z_value, skip_title=False, nb_loop=1, nb_tomo=1 + ): + super().add_scan( + scan_type, seq_ini_index, z_value, skip_title, nb_loop, nb_tomo + ) + if scan_type == "projection": + # register pcotomo specific informations (only for projections) + with h5py.File(self.sample_file, mode="a") as h5f: + seq_node = h5f.require_group(str(self._index - 1) + ".1") + scan_grp = seq_node.require_group("technique/scan") + if "nb_loop" not in scan_grp: + scan_grp["nb_loop"] = nb_loop + if "nb_tomo" not in scan_grp: + scan_grp["nb_tomo"] = nb_tomo + if "tomo_n" not in scan_grp: + scan_grp["tomo_n"] = self._tomo_n + + class _BlissZseriesTomo(_BlissSample): def __init__( self, @@ -352,6 +478,7 @@ class _BlissZseriesTomo(_BlissSample): z_values, with_nx_detector_attr=True, y_rot=False, + with_rotation_motor_info=True, ): self._z_values = z_values super().__init__( @@ -364,6 +491,7 @@ class _BlissZseriesTomo(_BlissSample): detector_name=detector_name, with_nx_detector_attr=with_nx_detector_attr, y_rot=y_rot, + with_rotation_motor_info=with_rotation_motor_info, ) def get_main_entry_title(self): @@ -405,6 +533,7 @@ class _BlissXRD_CT(_BlissSample): detector_name, with_nx_detector_attr=True, y_rot=False, + with_rotation_motor_info=True, ): """XRD-CT are scan with only projection. Projections are store in the init group. Data can be store under 1.1 but also 1.2...""" @@ -418,6 +547,7 @@ class _BlissXRD_CT(_BlissSample): detector_name=detector_name, with_nx_detector_attr=with_nx_detector_attr, y_rot=y_rot, + with_rotation_motor_info=with_rotation_motor_info, ) def get_main_entry_title(self): @@ -454,7 +584,7 @@ class _BlissXRD_CT(_BlissSample): scan_node.attrs["NX_class"] = "NXentry" measrurement_node = h5f.require_group("/".join((scan_name, "measurement"))) - data = numpy.random.random(self._n_frame_per_scan) * 256.2 + data = numpy.random.random(self._tomo_n) * 256.2 measrurement_node["fpico3"] = data h5f.require_group(scan_name)["title"] = self.get_main_entry_title() instrument_node = h5f.require_group("/".join((scan_name, "instrument"))) diff --git a/nxtomomill/utils/__init__.py b/nxtomomill/utils/__init__.py index c742905bf8a201b5dab8a171e85c10b719c88855..856b921421361e309d43cc669c6c71dbb9056b8d 100644 --- a/nxtomomill/utils/__init__.py +++ b/nxtomomill/utils/__init__.py @@ -1,4 +1,17 @@ from .frameappender import FrameAppender # noqa F401 -from .utils import * # noqa F401 +from .utils import embed_url # noqa F401 +from .utils import Format # noqa F401 +from .utils import FileExtension # noqa F401 +from .utils import get_file_name # noqa F401 +from .utils import get_tuple_of_keys_from_cmd # noqa F401 +from .utils import is_nx_tomo_entry # noqa F401 +from .utils import add_dark_flat_nx_file # noqa F401 +from .utils import change_image_key_control # noqa F401 from .progress import Progress # noqa F401 from .context import cwd_context # noqa F401 + +# expose ImageKey from tomoscan +try: + from tomoscan.esrf.scan.hdf5scan import ImageKey # noqa F401 +except ImportError: + from tomoscan.esrf.hdf5scan import ImageKey # noqa F401 diff --git a/nxtomomill/utils/context.py b/nxtomomill/utils/context.py index 4e3f06559f101b04cfb0075318cc8138a5d2fc8c..3684e728ef9d86e1862f69efbd7f1e5e0ac01625 100644 --- a/nxtomomill/utils/context.py +++ b/nxtomomill/utils/context.py @@ -33,9 +33,12 @@ import contextlib @contextlib.contextmanager -def cwd_context(): +def cwd_context(new_cwd=None): + curdir = os.getcwd() try: + if new_cwd not in (None, ""): + os.chdir(new_cwd) yield finally: os.chdir(curdir) diff --git a/nxtomomill/utils/file_path.py b/nxtomomill/utils/file_path.py new file mode 100644 index 0000000000000000000000000000000000000000..b0d084f30ac0a56bb255c93a710885a07296e253 --- /dev/null +++ b/nxtomomill/utils/file_path.py @@ -0,0 +1,43 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2022 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "15/02/2022" + + +import os + + +def to_target_rel_path(file_path, target_path): + """cast file_path to a relative path according to target_path. + This is used to deduce h5py.VirtualSource path + """ + file_path = os.path.realpath(file_path) + target_path = os.path.realpath(target_path) + path = os.path.relpath(file_path, os.path.dirname(target_path)) + if not path.startswith("./"): + path = "./" + path + return path diff --git a/nxtomomill/utils/frameappender.py b/nxtomomill/utils/frameappender.py index 47223592aa145e62e79b051f97150d256355e022..eedaef1662fab2e1d107a312e0b7051788dca435 100644 --- a/nxtomomill/utils/frameappender.py +++ b/nxtomomill/utils/frameappender.py @@ -36,6 +36,9 @@ import os import h5py from silx.io.utils import get_data import h5py._hl.selections as selection +from nxtomomill.utils.h5pyutils import from_data_url_to_virtual_source +from nxtomomill.utils.hdf5 import DatasetReader +from nxtomomill.utils.file_path import to_target_rel_path from .context import cwd_context from silx.io.utils import h5py_read_dataset from h5py import h5s as h5py_h5s @@ -57,7 +60,9 @@ class FrameAppender: if where not in ("start", "end"): raise ValueError("`where` should be `start` or `end`") - if not isinstance(data, (DataUrl, numpy.ndarray, list, tuple)): + if not isinstance( + data, (DataUrl, numpy.ndarray, list, tuple, h5py.VirtualSource) + ): raise TypeError( "data should be an instance of DataUrl or a numpy " "array not {}".format(type(data)) @@ -98,14 +103,10 @@ class FrameAppender: "directory without notifying." "See https://github.com/silx-kit/silx/issues/3277" ) - if isinstance(self.data, (numpy.ndarray, list, tuple)): - raise TypeError( - "Provided data is a numpy array when given" - "dataset path is a virtual dataset. " - "You must store the data somewhere else " - "and provide a DataUrl" - ) - else: + if isinstance(self.data, h5py.VirtualSource): + self.__insert_virtual_source_in_vds(h5s=h5s, new_virtual_source=self.data) + + elif isinstance(self.data, DataUrl): if self.logger is not None: self.logger.debug( "Update virtual dataset: " @@ -115,10 +116,6 @@ class FrameAppender: # store DataUrl in the current virtual dataset url = self.data - from nxtomomill.converter.hdf5.acquisition.baseacquisition import ( - DatasetReader, - ) - def check_dataset(dataset_frm_url): data_need_reshape = False """check if the dataset is valid or might need a reshape""" @@ -149,41 +146,32 @@ class FrameAppender: if url.data_slice() is None and not data_need_reshape: # case we can avoid to load the data in memory with DatasetReader(self.data) as data_frm_url: - self.__insert_in_vds(h5s, url, data_frm_url) + self.__insert_url_in_vds(h5s, url, data_frm_url) else: if loaded_dataset is None: data_frm_url = get_data(url) else: data_frm_url = loaded_dataset - self.__insert_in_vds(h5s, url, data_frm_url) + self.__insert_url_in_vds(h5s, url, data_frm_url) + else: + raise TypeError( + "Provided data is a numpy array when given" + "dataset path is a virtual dataset. " + "You must store the data somewhere else " + "and provide a DataUrl" + ) - def __insert_in_vds(self, h5s, url, data_frm_url): + def __insert_url_in_vds(self, h5s, url, data_frm_url): if data_frm_url.ndim == 2: - n_frames = 1 dim_2, dim_1 = data_frm_url.shape data_frm_url = data_frm_url.reshape(1, dim_2, dim_1) elif data_frm_url.ndim == 3: - n_frames, dim_2, dim_1 = data_frm_url.shape + _, dim_2, dim_1 = data_frm_url.shape else: raise ValueError("data to had is expected to be 2 or 3 d") - virtual_sources_len = [] - virtual_sources = [] - # we need to recreate the VirtualSource they are not - # store or available from the API - for vs_info in h5s[self.data_path].virtual_sources(): - length, vs = self._recreate_vs(vs_info=vs_info, vds_file=self.file_path) - virtual_sources.append(vs) - virtual_sources_len.append(length) - - vds_file_path = os.path.abspath(os.path.relpath(url.file_path(), os.getcwd())) - vds_file_path = os.path.realpath(vds_file_path) - vds_file_path = os.path.relpath(vds_file_path, os.path.dirname(self.file_path)) - if not vds_file_path.startswith("./"): - vds_file_path = "./" + vds_file_path - new_virtual_source = h5py.VirtualSource( - path_or_dataset=vds_file_path, + path_or_dataset=url.file_path(), name=url.data_path(), shape=data_frm_url.shape, ) @@ -196,19 +184,60 @@ class FrameAppender: h5sd[url.data_path()].shape, url.data_slice(), dst ) new_virtual_source.sel = sel + self.__insert_virtual_source_in_vds( + h5s=h5s, new_virtual_source=new_virtual_source, relative_path=True + ) - n_frames += h5s[self.data_path].shape[0] + def __insert_virtual_source_in_vds( + self, h5s, new_virtual_source: h5py.VirtualSource, relative_path=True + ): + if not isinstance(new_virtual_source, h5py.VirtualSource): + raise TypeError( + f"{new_virtual_source} is expected to be an instance of h5py.VirtualSource and not {type(new_virtual_source)}" + ) + if not len(new_virtual_source.shape) == 3: + raise ValueError( + f"virtual source shape is expected to be 3D and not {len(new_virtual_source.shape)}D." + ) + # preprocess virtualSource to insure having a relative path + if relative_path: + vds_file_path = to_target_rel_path(new_virtual_source.path, self.file_path) + new_virtual_source_sel = new_virtual_source.sel + new_virtual_source = h5py.VirtualSource( + path_or_dataset=vds_file_path, + name=new_virtual_source.name, + shape=new_virtual_source.shape, + ) + new_virtual_source.sel = new_virtual_source_sel + + virtual_sources_len = [] + virtual_sources = [] + # we need to recreate the VirtualSource they are not + # store or available from the API + for vs_info in h5s[self.data_path].virtual_sources(): + length, vs = self._recreate_vs(vs_info=vs_info, vds_file=self.file_path) + virtual_sources.append(vs) + virtual_sources_len.append(length) + + n_frames = h5s[self.data_path].shape[0] + new_virtual_source.shape[0] data_type = h5s[self.data_path].dtype if self.where == "start": virtual_sources.insert(0, new_virtual_source) - virtual_sources_len.insert(0, data_frm_url.shape[0]) + virtual_sources_len.insert(0, new_virtual_source.shape[0]) else: virtual_sources.append(new_virtual_source) - virtual_sources_len.append(data_frm_url.shape[0]) + virtual_sources_len.append(new_virtual_source.shape[0]) # create the new virtual dataset - layout = h5py.VirtualLayout(shape=(n_frames, dim_2, dim_1), dtype=data_type) + layout = h5py.VirtualLayout( + shape=( + n_frames, + new_virtual_source.shape[-2], + new_virtual_source.shape[-1], + ), + dtype=data_type, + ) last = 0 for v_source, vs_len in zip(virtual_sources, virtual_sources_len): layout[last : vs_len + last] = v_source @@ -291,20 +320,41 @@ class FrameAppender: if isinstance(self.data, DataUrl): url = self.data - with HDF5File(url.file_path(), mode="r") as o_h5s: - if not o_h5s[url.data_path()].ndim in (2, 3): - raise ValueError( - "{} should point to 2D or 3D dataset ".format(url.path()) - ) - data_shape = o_h5s[url.data_path()].shape - data_type = o_h5s[url.data_path()].dtype - if len(data_shape) == 2: - data_shape = (1, data_shape[0], data_shape[1]) - layout = h5py.VirtualLayout(shape=data_shape, dtype=data_type) - layout[:] = h5py.VirtualSource( - url.file_path(), url.data_path(), shape=data_shape - ) + url_file_path = to_target_rel_path(url.file_path(), self.file_path) + url = DataUrl( + file_path=url_file_path, + data_path=url.data_path(), + scheme=url.scheme(), + data_slice=url.data_slice(), + ) + + with cwd_context(os.path.dirname(self.file_path)): + vs, vs_shape, data_type = from_data_url_to_virtual_source(url) + layout = h5py.VirtualLayout(shape=vs_shape, dtype=data_type) + layout[:] = vs h5s.create_virtual_dataset(self.data_path, layout) + + elif isinstance(self.data, h5py.VirtualSource): + virtual_source = self.data + layout = h5py.VirtualLayout( + shape=virtual_source.shape, + dtype=virtual_source.dtype, + ) + + vds_file_path = to_target_rel_path(virtual_source.path, self.file_path) + virtual_source_rel_path = h5py.VirtualSource( + path_or_dataset=vds_file_path, + name=virtual_source.name, + shape=virtual_source.shape, + ) + virtual_source_rel_path.sel = virtual_source.sel + layout[:] = virtual_source_rel_path + # convert path to relative + h5s.create_virtual_dataset(self.data_path, layout) + elif not isinstance(self.data, numpy.ndarray): + raise TypeError( + f"self.data should be an instance of DataUrl, a numpy array or a VirtualSource. Not {type(self.data)}" + ) else: h5s[self.data_path] = self.data @@ -317,10 +367,13 @@ class FrameAppender: the use case exposed in issue: https://gitlab.esrf.fr/tomotools/nxtomomill/-/issues/40 """ - with cwd_context(): - if os.path.dirname(vds_file) not in ("", None): - os.chdir(os.path.dirname(vds_file)) - with HDF5File(vs_info.file_name, mode="r") as vs_node: + with cwd_context(os.path.dirname(vds_file)): + dataset_file_path = vs_info.file_name + # in case the virtual source is in the same file + if dataset_file_path == ".": + dataset_file_path = vds_file + + with HDF5File(dataset_file_path, mode="r") as vs_node: dataset = vs_node[vs_info.dset_name] select_bounds = vs_info.vspace.get_select_bounds() left_bound = select_bounds[0] diff --git a/nxtomomill/utils/h5pyutils.py b/nxtomomill/utils/h5pyutils.py new file mode 100644 index 0000000000000000000000000000000000000000..6914c83d574aeaa4b6f58f77fa18454907285980 --- /dev/null +++ b/nxtomomill/utils/h5pyutils.py @@ -0,0 +1,86 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2015-2020 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + +""" +module to define some converter utils function +""" + +__authors__ = [ + "H. Payno", +] +__license__ = "MIT" +__date__ = "08/03/2022" + + +from tomoscan.io import HDF5File +import h5py +from silx.io.url import DataUrl +import h5py._hl.selections as selection + + +def from_data_url_to_virtual_source(url: DataUrl) -> tuple: + """ + :param DataUrl url: url to be converted to a virtual source. It must target a 2D detector + :return: (h5py.VirtualSource, tuple(shape of the virtual source), numpy.drype: type of the dataset associated with the virtual source) + :rtype: tuple + """ + if not isinstance(url, DataUrl): + raise TypeError( + f"url is expected to be an instance of DataUrl and not {type(url)}" + ) + + with HDF5File(url.file_path(), mode="r") as o_h5s: + + original_data_shape = o_h5s[url.data_path()].shape + data_type = o_h5s[url.data_path()].dtype + if len(original_data_shape) == 2: + original_data_shape = ( + 1, + original_data_shape[0], + original_data_shape[1], + ) + + vs_shape = original_data_shape + if url.data_slice() is not None: + vs_shape = ( + url.data_slice().stop - url.data_slice().start, + original_data_shape[-2], + original_data_shape[-1], + ) + + vs = h5py.VirtualSource(url.file_path(), url.data_path(), shape=vs_shape) + + if url.data_slice() is not None: + vs.sel = selection.select(original_data_shape, url.data_slice()) + return vs, vs_shape, data_type + + +def from_virtual_source_to_data_url(vs: h5py.VirtualSource) -> DataUrl: + if not isinstance(vs, h5py.VirtualSource): + raise TypeError( + f"vs is expected to be an instance of h5py.VirtualSorce and not {type(vs)}" + ) + url = DataUrl(file_path=vs.path, data_path=vs.name, scheme="silx") + return url diff --git a/nxtomomill/utils/hdf5.py b/nxtomomill/utils/hdf5.py new file mode 100644 index 0000000000000000000000000000000000000000..d482bb649b5a6f4e12ffc0c699ffd1576fe957ec --- /dev/null +++ b/nxtomomill/utils/hdf5.py @@ -0,0 +1,86 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2022 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "23/02/2022" + + +import contextlib +import h5py +from tomoscan.io import HDF5File + +try: + import hdf5plugin # noqa F401 +except ImportError: + pass +from silx.io.url import DataUrl + + +class _BaseReader(contextlib.AbstractContextManager): + def __init__(self, url: DataUrl): + if not isinstance(url, DataUrl): + raise TypeError("url should be an instance of DataUrl") + if url.scheme() not in ("silx", "h5py"): + raise ValueError("Valid scheme are silx and h5py") + if url.data_slice() is not None: + raise ValueError( + "Data slices are not managed. Data path should " + "point to a bliss node (h5py.Group)" + ) + self._url = url + self._file_handler = None + + def __exit__(self, *exc): + return self._file_handler.close() + + +class EntryReader(_BaseReader): + """Context manager used to read a bliss node""" + + def __enter__(self): + self._file_handler = HDF5File(filename=self._url.file_path(), mode="r") + if self._url.data_path() == "": + entry = self._file_handler + else: + entry = self._file_handler[self._url.data_path()] + if not isinstance(entry, h5py.Group): + raise ValueError("Data path should point to a bliss node (h5py.Group)") + return entry + + +class DatasetReader(_BaseReader): + """Context manager used to read a bliss node""" + + def __enter__(self): + self._file_handler = HDF5File(filename=self._url.file_path(), mode="r") + entry = self._file_handler[self._url.data_path()] + if not isinstance(entry, h5py.Dataset): + raise ValueError( + "Data path ({}) should point to a dataset (h5py.Dataset)".format( + self._url.path() + ) + ) + return entry diff --git a/nxtomomill/utils/nxsplitter.py b/nxtomomill/utils/nxsplitter.py new file mode 100644 index 0000000000000000000000000000000000000000..f3e63bb0e31794b93d88401164c1c71bc21dcc8b --- /dev/null +++ b/nxtomomill/utils/nxsplitter.py @@ -0,0 +1,417 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2022 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "17/02/2022" + + +from typing import Optional, Union +from nxtomomill.utils.hdf5 import DatasetReader +from silx.io.url import DataUrl +from nxtomomill.nexus.nxtomo import NXtomo +import numpy +import h5py +from silx.io.utils import get_data +import copy +import h5py._hl.selections as selection + + +class _NXtomoDetectorDataSplitter: + """Splitter to split the dataset nxtomo.instrument.detector.data + into several NXtomo. + This will also keep up to date rotation_angle, image_key, x_translation... datasets. + In order to start the processing it requires a correctly formed NXtomo (same number of image_key, rotation_angle...) + This is required for the pcotomo acquisition. + """ + + def __init__(self, nx_tomo: NXtomo) -> None: + if not isinstance(nx_tomo, NXtomo): + raise TypeError( + f"nxtomo is expected to be an instance of {NXtomo} and not {type(nx_tomo)}" + ) + self._nx_tomo = nx_tomo + + @property + def nx_tomo(self) -> NXtomo: + return self._nx_tomo + + def split(self, data_slice: slice, nb_part: int) -> tuple: + """ + split the dataset targetted to have a set of h5py.VirtualSource. + + :param int nb_part: in how many contiguous dataset the instruement.detector.data must be splitted. + :raises: ValueError if the number of frame, image_key, x_translation... is incoherent. + """ + if not isinstance(nb_part, int): + raise TypeError(f"nb_part is expected to be an int not {type(nb_part)}") + invalid_datasets = self.get_invalid_datasets() + if len(invalid_datasets) > 0: + raise ValueError( + f"Some datasets have incoherent length compared to nx_tomo.instrument.detector.data length: {invalid_datasets}" + ) + elif nb_part <= 0: + raise ValueError(f"nb_part is expected to be >=1 not {nb_part}") + elif nb_part == 1: + return [ + self.nx_tomo, + ] + elif data_slice.step not in (1, None): + raise ValueError("slice step must be one.") + elif (data_slice.stop - data_slice.start) % nb_part != 0: + raise ValueError( + f"incoherent split requested. Request to spluit {(data_slice.stop - data_slice.start)} slices into {nb_part} parts" + ) + else: + parts = [] + for i_part in range(nb_part): + current_slice = data_slice + new_slice_size = (current_slice.stop - current_slice.start) // nb_part + new_slice = slice( + current_slice.start + new_slice_size * i_part, + current_slice.start + new_slice_size * (i_part + 1), + 1, + ) + nx_tomo_part = self.replace(old_slice=data_slice, new_slice=new_slice) + parts.append(nx_tomo_part) + return parts + + def replace(self, old_slice: slice, new_slice: slice) -> NXtomo: + """ + replace a section of the NXtomo instrument.detector.data by a subsection of it + """ + if not isinstance(old_slice, slice): + raise TypeError("old_slice is expected to be a slice") + if not isinstance(new_slice, slice): + raise TypeError("new_slice is expected to be a slice") + if old_slice.step not in (None, 1): + raise ValueError("old_slice step is expected to be one") + if new_slice.step not in (None, 1): + raise ValueError("new_slice step is expected to be one") + + if new_slice.start < old_slice.start or new_slice.stop > old_slice.stop: + raise ValueError( + f"new_slice ({new_slice}) must be contained in old_slice ({old_slice})" + ) + + if old_slice.start < 0: + raise ValueError( + f"old_slice.start must be at least 0 not {old_slice.start}" + ) + n_frames = self._get_n_frames() + if n_frames is not None and old_slice.stop > n_frames: + raise ValueError( + f"old_slice.start must be at most {n_frames} not {old_slice.stop}" + ) + + # handles datasets other than instrument.detector.data + result_nx_tomo = copy.deepcopy(self.nx_tomo) + if result_nx_tomo.sample.rotation_angle is not None: + result_nx_tomo.sample.rotation_angle = numpy.concatenate( + [ + self.nx_tomo.sample.rotation_angle[: old_slice.start], + self.nx_tomo.sample.rotation_angle[new_slice], + self.nx_tomo.sample.rotation_angle[old_slice.stop :], + ] + ) + + if result_nx_tomo.sample.x_translation.value is not None: + result_nx_tomo.sample.x_translation.value = numpy.concatenate( + [ + self.nx_tomo.sample.x_translation.value[: old_slice.start], + self.nx_tomo.sample.x_translation.value[new_slice], + self.nx_tomo.sample.x_translation.value[old_slice.stop :], + ] + ) + + if result_nx_tomo.sample.y_translation.value is not None: + result_nx_tomo.sample.y_translation.value = numpy.concatenate( + [ + self.nx_tomo.sample.y_translation.value[: old_slice.start], + self.nx_tomo.sample.y_translation.value[new_slice], + self.nx_tomo.sample.y_translation.value[old_slice.stop :], + ] + ) + + if result_nx_tomo.sample.z_translation.value is not None: + result_nx_tomo.sample.z_translation.value = numpy.concatenate( + [ + self.nx_tomo.sample.z_translation.value[: old_slice.start], + self.nx_tomo.sample.z_translation.value[new_slice], + self.nx_tomo.sample.z_translation.value[old_slice.stop :], + ] + ) + + if result_nx_tomo.instrument.detector.image_key_control is not None: + result_nx_tomo.instrument.detector.image_key_control = numpy.concatenate( + [ + self.nx_tomo.instrument.detector.image_key_control[ + : old_slice.start + ], + self.nx_tomo.instrument.detector.image_key_control[new_slice], + self.nx_tomo.instrument.detector.image_key_control[ + old_slice.stop : + ], + ] + ) + + # handles detector.data dataset. This one is special because it can contains + # numpy arrays (raw data), h5py.VirtualSource or DataUrl (or be None) + det_data = self.nx_tomo.instrument.detector.data + if det_data is None: + pass + elif isinstance(det_data, numpy.ndarray): + result_nx_tomo.instrument.detector.data = numpy.concatenate( + [ + det_data[: old_slice.start], + det_data[new_slice], + det_data[old_slice.stop :], + ] + ) + elif isinstance(det_data, (tuple, list)): + result_nx_tomo.instrument.detector.data = numpy.concatenate( + [ + self._get_detector_data_sub_section(slice(0, old_slice.start, 1)), + self._get_detector_data_sub_section(new_slice), + self._get_detector_data_sub_section( + slice(old_slice.stop, n_frames + 1, 1) + ), + ] + ).tolist() + else: + raise TypeError( + f"instrument.detector.data must be a numpy array or a VirtualSource or a DataUrl. Not {type(det_data)}" + ) + return result_nx_tomo + + def _get_detector_data_sub_section(self, section: slice) -> tuple: + """ + return a tuple of DataUrl or h5py.VirtualSource fitting the slice requested + """ + det_data = self.nx_tomo.instrument.detector.data + res = [] + if section.start == section.stop: + return () + + def get_elmt_shape(elmt: Union[h5py.VirtualSource, DataUrl]) -> tuple: + if isinstance(elmt, h5py.VirtualSource): + return elmt.shape + elif isinstance(elmt, DataUrl): + with DatasetReader(elmt) as dataset: + return dataset.shape + else: + raise TypeError( + f"elmt must be a DataUrl or h5py.VirtualSource. Not {type(elmt)}" + ) + + def get_elmt_nb_frame(elmt: Union[h5py.VirtualSource, DataUrl]) -> int: + shape = get_elmt_shape(elmt) + if len(shape) == 3: + return shape[0] + elif len(shape) == 2: + return 1 + else: + raise ValueError(f"virtualSource: {elmt} is not 2D or 3D") + + def construct_slices_elmt_list() -> dict: + "create a dictionary with slice as key and DataUrl or h5py.VirtualSource as value" + slices_elmts = [] + + current_index = 0 + for elmt in det_data: + n_frame = get_elmt_nb_frame(elmt) + slice_ = slice(current_index, current_index + n_frame, 1) + slices_elmts.append([slice_, elmt]) + current_index += n_frame + return slices_elmts + + def intersect(slice_1, slice_2): + """check if the two slices intersect""" + assert isinstance(slice_1, slice) and slice_1.step == 1 + assert isinstance(slice_2, slice) and slice_2.step == 1 + return slice_1.start < slice_2.stop and slice_1.stop > slice_2.start + + def select( + elmt: Union[h5py.VirtualSource, DataUrl], region: slice + ) -> Union[h5py.VirtualSource, DataUrl]: + """select a region on the elmt. + Can return at most the elmt itself or a region of it""" + elmt_n_frame = get_elmt_nb_frame(elmt) + assert elmt_n_frame != 0 + clamp_region = slice( + max(0, region.start), + min(elmt_n_frame, region.stop), + 1, + ) + assert clamp_region.start != clamp_region.stop + + if isinstance(elmt, h5py.VirtualSource): + frame_dims = elmt.shape[-2], elmt.shape[-1] + n_frames = clamp_region.stop - clamp_region.start + assert n_frames > 0 + shape = (n_frames, frame_dims[0], frame_dims[1]) + vs = h5py.VirtualSource( + path_or_dataset=elmt.path, + name=elmt.name, + shape=shape, + ) + vs.sel = selection.select(elmt.shape, clamp_region) + return vs + else: + if elmt.data_slice() is None: + data_slice = clamp_region + elif isinstance(elmt.data_slice(), slice): + if elmt.data_slice.step not in (1, None): + raise ValueError("DataUrl with step !=1 are not handled") + else: + data_slice = slice( + elmt.data_slice.start + clamp_region.start, + elmt.data_slice.start + clamp_region.stop, + 1, + ) + else: + raise TypeError( + f"data_slice is expected to be None or a slice. Not {type(elmt.data_slice())}" + ) + return DataUrl( + file_path=elmt.file_path(), + data_path=elmt.data_path(), + scheme=elmt.scheme(), + data_slice=data_slice, + ) + + for slice_raw_data, elmt in construct_slices_elmt_list(): + if intersect(section, slice_raw_data): + res.append( + select( + elmt, + slice( + section.start - slice_raw_data.start, + section.stop - slice_raw_data.start, + 1, + ), + ) + ) + return tuple(res) + + def get_invalid_datasets(self) -> dict: + """ + return a dict of invalid dataset compare to the instrument.detector.data dataset. + Key is the location ? path to the invalid dataset. Value is the reason of the failure. + """ + invalid_datasets = {} + n_frames = self._get_n_frames() + + # check rotation_angle + if self.nx_tomo.sample.rotation_angle is not None: + n_rotation_angles = len(self.nx_tomo.sample.rotation_angle) + if n_rotation_angles != n_frames: + invalid_datasets[ + "sample/rotation_angle" + ] = f"{n_rotation_angles} angles found when {n_frames} expected" + + # check image_key_control (force to have the same number as image_key already so only check one) + if self.nx_tomo.instrument.detector.image_key_control is not None: + n_image_key_control = len( + self.nx_tomo.instrument.detector.image_key_control + ) + if n_image_key_control != n_frames: + invalid_datasets[ + "instrument/detector/image_key_control" + ] = f"{n_image_key_control} image_key_control values found when {n_frames} expected" + + # check x_translation + if self.nx_tomo.sample.x_translation.value is not None: + n_x_translation = len(self.nx_tomo.sample.x_translation.value) + if n_x_translation != n_frames: + invalid_datasets[ + "sample/x_translation" + ] = f"{n_x_translation} x translations found when {n_frames} expected" + + # check y_translation + if self.nx_tomo.sample.y_translation.value is not None: + n_y_translation = len(self.nx_tomo.sample.y_translation.value) + if n_y_translation != n_frames: + invalid_datasets[ + "sample/y_translation" + ] = f"{n_y_translation} y translations found when {n_frames} expected" + + # check z_translation + if self.nx_tomo.sample.z_translation.value is not None: + n_z_translation = len(self.nx_tomo.sample.z_translation.value) + if n_z_translation != n_frames: + invalid_datasets[ + "sample/z_translation" + ] = f"{n_z_translation} z translations found when {n_frames} expected" + + return invalid_datasets + + def _get_n_frames(self) -> Optional[int]: + dataset = self.nx_tomo.instrument.detector.data + if dataset is None: + return None + elif isinstance(dataset, numpy.ndarray): + if not dataset.ndim == 3: + raise ValueError( + f"nx_tomo.instrument.detector.data is expected to be 3D and not {dataset.ndim}D." + ) + else: + return dataset.shape[0] + elif isinstance(dataset, (list, tuple)): + n_frames = 0 + for dataset_elmt in dataset: + if isinstance(dataset_elmt, h5py.VirtualSource): + shape = dataset_elmt.shape + if len(shape) == 3: + n_frames += dataset_elmt.shape[0] + elif len(shape) == 2: + n_frames += 1 + else: + raise ValueError( + f"h5py.VirtualSource shape is expected to be 2D (single frame) or 3D. Not {len(shape)}D." + ) + elif isinstance(dataset_elmt, DataUrl): + data = get_data(dataset_elmt) + if not isinstance(data, numpy.ndarray): + raise TypeError( + f"url: {dataset_elmt.path()} is not pointing to an array" + ) + elif data.ndim == 2: + n_frames += 1 + elif data.ndim == 3: + n_frames += data.shape[0] + else: + raise ValueError( + f"url: {dataset_elmt.path()} is expected to be 2D or 3D. Not {dataset_elmt.ndim} D" + ) + else: + raise TypeError( + f"elements of {type(dataset)} must be h5py.VirtualSource) or silx.io.url.DataUrl and not {type(dataset_elmt)}" + ) + return n_frames + else: + raise TypeError( + f"nx_tomo.instrument.detector.data type ({type(dataset)}) is not handled" + ) diff --git a/nxtomomill/utils/test/test_nx_splitter.py b/nxtomomill/utils/test/test_nx_splitter.py new file mode 100644 index 0000000000000000000000000000000000000000..9d2f2f665f73c3205e1cbad90b91bf0ef4029f73 --- /dev/null +++ b/nxtomomill/utils/test/test_nx_splitter.py @@ -0,0 +1,474 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2016-2022 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ + +__authors__ = ["H. Payno"] +__license__ = "MIT" +__date__ = "17/02/2022" + + +from tempfile import TemporaryDirectory +from nxtomomill.nexus.nxtomo import NXtomo +from nxtomomill.utils.nxsplitter import _NXtomoDetectorDataSplitter +from silx.io.url import DataUrl +import numpy +import pytest +import h5py +import os + + +def test_get_invalid_datasets(): + """test the _NXtomoDetectorDataSplitter get_invalid_datasets function""" + nx_tomo = NXtomo("test") + n_frames = 10 + nx_tomo.instrument.detector.data = numpy.random.random(100 * 100 * 10).reshape( + [n_frames, 100, 100] + ) + splitter = _NXtomoDetectorDataSplitter(nx_tomo=nx_tomo) + assert len(splitter.get_invalid_datasets()) == 0 + + # test rotation angle + nx_tomo.sample.rotation_angle = [12, 13] + assert len(splitter.get_invalid_datasets()) == 1 + nx_tomo.sample.rotation_angle = [0] * n_frames + assert len(splitter.get_invalid_datasets()) == 0 + + # test image_key_control + nx_tomo.instrument.detector.image_key_control = [0] + assert len(splitter.get_invalid_datasets()) == 1 + nx_tomo.instrument.detector.image_key_control = [0] * n_frames + assert len(splitter.get_invalid_datasets()) == 0 + + # test x_translation + nx_tomo.sample.x_translation = [0] + assert len(splitter.get_invalid_datasets()) == 1 + nx_tomo.sample.x_translation = [0] * n_frames + assert len(splitter.get_invalid_datasets()) == 0 + + # test y_translation + nx_tomo.sample.y_translation = [0] + assert len(splitter.get_invalid_datasets()) == 1 + nx_tomo.sample.y_translation = [0] * n_frames + assert len(splitter.get_invalid_datasets()) == 0 + + # test z_translation + nx_tomo.sample.z_translation = [0] + assert len(splitter.get_invalid_datasets()) == 1 + nx_tomo.sample.z_translation = [0] * n_frames + assert len(splitter.get_invalid_datasets()) == 0 + + +def test_spliter_raw_data(): + """test the splitter on a simple non virtual h5py Dataset""" + nx_tomo = NXtomo("test_raw_data") + n_frames = 20 + nx_tomo.instrument.detector.data = numpy.random.random( + 100 * 100 * n_frames + ).reshape([n_frames, 100, 100]) + nx_tomo.sample.rotation_angle = [0, 12] + splitter = _NXtomoDetectorDataSplitter(nx_tomo=nx_tomo) + # check incoherent number of rotation + with pytest.raises(ValueError): + splitter.split(data_slice=slice(0, 100, 1), nb_part=2) + nx_tomo.sample.rotation_angle = numpy.linspace(0, 180, num=n_frames, endpoint=False) + # check slice nb_part < 0 + with pytest.raises(ValueError): + splitter.split(data_slice=slice(0, 100, 1), nb_part=-1) + # check slice step != 1 + with pytest.raises(ValueError): + splitter.split(data_slice=slice(0, 100, 2), nb_part=2) + # check incoherent number of frames + with pytest.raises(ValueError): + splitter.split(data_slice=slice(0, 99, 2), nb_part=2) + + # check x translation + nx_tomo.sample.x_translation = [0, 12] + with pytest.raises(ValueError): + splitter.split(data_slice=slice(0, 100, 1), nb_part=2) + nx_tomo.sample.x_translation = numpy.random.random(n_frames) + nx_tomo.sample.y_translation = numpy.random.random(n_frames) + nx_tomo.sample.z_translation = numpy.random.random(n_frames) + # check image key + nx_tomo.instrument.detector.image_key_control = [0, 2] + with pytest.raises(ValueError): + splitter.split(data_slice=slice(0, 100, 1), nb_part=2) + nx_tomo.instrument.detector.image_key_control = [ + numpy.random.randint(low=-1, high=2) for i in range(n_frames) + ] + + assert splitter.split(data_slice=slice(0, 100, 1), nb_part=1) == [ + nx_tomo, + ] + + # check error if request to split a region bigger that the one (100 vs n_frames) + with pytest.raises(ValueError): + splitted_nx_tomo = splitter.split(data_slice=slice(0, 100, 1), nb_part=2) + + splitted_nx_tomo = splitter.split(data_slice=slice(0, 20, 1), nb_part=2) + assert len(splitted_nx_tomo) == 2 + s_nx_tomo_1, s_nx_tomo_2 = splitted_nx_tomo + # chek rotation_angle + numpy.testing.assert_array_equal( + s_nx_tomo_1.sample.rotation_angle, + nx_tomo.sample.rotation_angle[0 : n_frames // 2], + ) + numpy.testing.assert_array_equal( + s_nx_tomo_2.sample.rotation_angle, + nx_tomo.sample.rotation_angle[n_frames // 2 :], + ) + # check image key and image key + numpy.testing.assert_array_equal( + s_nx_tomo_1.instrument.detector.image_key_control, + nx_tomo.instrument.detector.image_key_control[0 : n_frames // 2], + ) + numpy.testing.assert_array_equal( + s_nx_tomo_2.instrument.detector.image_key_control, + nx_tomo.instrument.detector.image_key_control[n_frames // 2 :], + ) + + # chek x translation + numpy.testing.assert_array_equal( + s_nx_tomo_1.sample.x_translation.value, + nx_tomo.sample.x_translation.value[0 : n_frames // 2], + ) + numpy.testing.assert_array_equal( + s_nx_tomo_2.sample.x_translation.value, + nx_tomo.sample.x_translation.value[n_frames // 2 :], + ) + + # chek y translation + numpy.testing.assert_array_equal( + s_nx_tomo_1.sample.y_translation.value, + nx_tomo.sample.y_translation.value[0 : n_frames // 2], + ) + numpy.testing.assert_array_equal( + s_nx_tomo_2.sample.y_translation.value, + nx_tomo.sample.y_translation.value[n_frames // 2 :], + ) + # chek z translation + numpy.testing.assert_array_equal( + s_nx_tomo_1.sample.z_translation.value, + nx_tomo.sample.z_translation.value[0 : n_frames // 2], + ) + numpy.testing.assert_array_equal( + s_nx_tomo_2.sample.z_translation.value, + nx_tomo.sample.z_translation.value[n_frames // 2 :], + ) + # check detector data + numpy.testing.assert_array_equal( + s_nx_tomo_1.instrument.detector.data, + nx_tomo.instrument.detector.data[0 : n_frames // 2], + ) + numpy.testing.assert_array_equal( + s_nx_tomo_2.instrument.detector.data, + nx_tomo.instrument.detector.data[n_frames // 2 :], + ) + + +def test_spliter_virtual_sources_1(): + """ + test the splitter on a simulated h5py virtual dataset composed of two Virtual source. + It must have the same virtual source in the two splitted NXtomo + test the splitter on a h5py virtual dataset + rotation_angle, [W]_translation and image_key datasets are handle are always numpy arrays + not pointing to any external ressources. This is only the case for + detector.data so this is the only dataset to test here + """ + nx_tomo = NXtomo("test_raw_data") + nx_tomo.instrument.detector.data = [ + h5py.VirtualSource("path_to_dataset_1", name="dataset_1", shape=[10, 100, 100]), + h5py.VirtualSource("path_to_dataset_2", name="dataset_2", shape=[10, 100, 100]), + ] + splitter = _NXtomoDetectorDataSplitter(nx_tomo=nx_tomo) + splitted_nx_tomo = splitter.split(data_slice=slice(0, 20, 1), nb_part=2) + assert len(splitted_nx_tomo) == 2 + s_nx_tomo_1, s_nx_tomo_2 = splitted_nx_tomo + det_dataset_1 = s_nx_tomo_1.instrument.detector.data + det_dataset_2 = s_nx_tomo_2.instrument.detector.data + assert len(det_dataset_1) == 1 + assert len(det_dataset_2) == 1 + + det_dataset_vs1 = det_dataset_1[0] + det_dataset_vs2 = det_dataset_2[0] + assert isinstance(det_dataset_vs1, h5py.VirtualSource) + assert det_dataset_vs1.path == "path_to_dataset_1" + assert det_dataset_vs1.shape == (10, 100, 100) + + assert isinstance(det_dataset_vs2, h5py.VirtualSource) + assert det_dataset_vs2.path == "path_to_dataset_2" + assert det_dataset_vs2.shape == (10, 100, 100) + + +def test_spliter_virtual_sources_2(): + """ + test the splitter on a h5py virtual dataset composed of a single Virtual source. + Must split this one into two VirtualSource + rotation_angle, [W]_translation and image_key datasets are handle are always numpy arrays + not pointing to any external ressources. This is only the case for + detector.data so this is the only dataset to test here + """ + nx_tomo = NXtomo("test_raw_data") + nx_tomo.instrument.detector.data = [ + h5py.VirtualSource( + "path_to_dataset", name="path_to_dataset", shape=[20, 100, 100] + ), + ] + splitter = _NXtomoDetectorDataSplitter(nx_tomo=nx_tomo) + splitted_nx_tomo = splitter.split(data_slice=slice(0, 20, 1), nb_part=2) + assert len(splitted_nx_tomo) == 2 + + splitted_nx_tomo = splitter.split(data_slice=slice(0, 20, 1), nb_part=4) + assert len(splitted_nx_tomo) == 4 + + +def test_spliter_virtual_sources_3(): + """ + test the splitter on a concrete h5py virtual dataset + rotation_angle, [W]_translation and image_key datasets are handle are always numpy arrays + not pointing to any external ressources. This is only the case for + detector.data so this is the only dataset to test here + """ + n_file = 5 + n_frame_per_file = 20 + layout = h5py.VirtualLayout( + shape=(n_file * n_frame_per_file, 100, 100), dtype=float + ) + with TemporaryDirectory() as folder: + for i_file in range(n_file): + file_path = os.path.join(folder, f"file{i_file}.hdf5") + data_path = f"path_to_dataset_{i_file}" + with h5py.File(file_path, mode="w") as h5f: + if i_file == 0: + data = numpy.ones([n_frame_per_file, 100, 100]) + elif i_file == n_file - 1: + data = numpy.ones([n_frame_per_file, 100, 100]) * 2 + else: + start = i_file * 1000.0 + stop = i_file * 1000.0 + (n_frame_per_file * 100 * 100) + data = numpy.arange(start, stop).reshape(n_frame_per_file, 100, 100) + h5f[data_path] = data + vs = h5py.VirtualSource(h5f[data_path]) + layout[i_file * n_frame_per_file : (i_file + 1) * n_frame_per_file] = vs + + master_file = os.path.join(folder, "master_file.hdf5") + with h5py.File(master_file, mode="w") as h5f: + h5f.create_virtual_dataset("data", layout) + original_data = h5f["data"][()] + + nx_tomo = NXtomo("entry0000") + with h5py.File(master_file, mode="r") as h5f: + vs_ = [] + for vs_info in h5f["data"].virtual_sources(): + vs_.append( + h5py.VirtualSource( + vs_info.file_name, + vs_info.dset_name, + shape=(n_frame_per_file, 100, 100), + ) + ) + + nx_tomo.instrument.detector.data = vs_ + + splitter = _NXtomoDetectorDataSplitter(nx_tomo=nx_tomo) + data_slice = slice(10, n_frame_per_file * n_file - 10, 1) + splitted_nx_tomo = splitter.split(data_slice=data_slice, nb_part=2) + assert len(splitted_nx_tomo) == 2 + # check the two dataset created + s_nx_tomo_1, s_nx_tomo_2 = splitted_nx_tomo + output_file_1 = os.path.join(folder, "output_file_1.nx") + + # data must contains a common section between the two nxtomo: the first 10 and last 10 frames + # then the rest must be splitted between the two NXtomo + assert len(s_nx_tomo_1.instrument.detector.data) == 5 + assert s_nx_tomo_1.instrument.detector.data[0].shape[0] == 10 + assert s_nx_tomo_1.instrument.detector.data[1].shape[0] == 10 + assert s_nx_tomo_1.instrument.detector.data[2].shape[0] == 20 + assert s_nx_tomo_1.instrument.detector.data[3].shape[0] == 10 + assert s_nx_tomo_1.instrument.detector.data[4].shape[0] == 10 + s_nx_tomo_1.save(output_file_1) + + output_file_2 = os.path.join(folder, "output_file_2.nx") + assert len(s_nx_tomo_2.instrument.detector.data) == 5 + assert s_nx_tomo_2.instrument.detector.data[0].shape[0] == 10 + assert s_nx_tomo_2.instrument.detector.data[1].shape[0] == 10 + assert s_nx_tomo_2.instrument.detector.data[2].shape[0] == 20 + assert s_nx_tomo_2.instrument.detector.data[3].shape[0] == 10 + assert s_nx_tomo_2.instrument.detector.data[4].shape[0] == 10 + s_nx_tomo_2.save(output_file_2) + + # check final datasets are correctly formed + with h5py.File(output_file_1, mode="r") as h5f: + nx_1_data = h5f["entry0000/instrument/detector/data"][()] + assert nx_1_data.shape[0] == 60 + + # check final datasets are correctly formed + with h5py.File(output_file_2, mode="r") as h5f: + nx_2_data = h5f["entry0000/instrument/detector/data"][()] + assert nx_2_data.shape[0] == 60 + + # first 10 frames (common between the three nxtomo) + numpy.testing.assert_array_equal( + nx_1_data[0:10], + nx_2_data[0:10], + ) + numpy.testing.assert_array_equal( + nx_1_data[0:10], + original_data[0:10], + ) + + # last 10 frames (common between the three nxtomo) + numpy.testing.assert_array_equal( + nx_1_data[-10:], + nx_2_data[-10:], + ) + numpy.testing.assert_array_equal( + nx_1_data[-10:], + original_data[-10:], + ) + + # test nx_1_data unique region + numpy.testing.assert_array_equal( + nx_1_data[10:50], + original_data[10:50], + ) + + # test nx_2_data unique region + numpy.testing.assert_array_equal( + nx_2_data[10:50], + original_data[50:90], + ) + + +def test_spliter_data_url(): + """ + test the splitter on a list of DataUrl + rotation_angle, [W]_translation and image_key datasets are handle are always numpy arrays + not pointing to any external ressources. This is only the case for + detector.data so this is the only dataset to test here + """ + urls = [] + n_frame_per_file = 20 + n_file = 5 + original_data = [] + with TemporaryDirectory() as folder: + for i_file in range(n_file): + file_path = os.path.join(folder, f"file{i_file}.hdf5") + data_path = f"path_to_dataset_{i_file}" + with h5py.File(file_path, mode="w") as h5f: + if i_file == 0: + data = numpy.ones([n_frame_per_file, 100, 100]) + elif i_file == n_file - 1: + data = numpy.ones([n_frame_per_file, 100, 100]) * 2 + else: + start = i_file * 1000.0 + stop = i_file * 1000.0 + (n_frame_per_file * 100 * 100) + data = numpy.arange(start, stop).reshape(n_frame_per_file, 100, 100) + h5f[data_path] = data + original_data.append(data) + + urls.append( + DataUrl( + file_path=file_path, + data_path=data_path, + scheme="silx", + ) + ) + + original_data = numpy.concatenate(original_data) + + nx_tomo = NXtomo("entry0000") + nx_tomo.instrument.detector.data = urls + + splitter = _NXtomoDetectorDataSplitter(nx_tomo=nx_tomo) + data_slice = slice(10, n_frame_per_file * n_file - 10, 1) + data_slice = slice(10, n_frame_per_file * n_file - 10, 1) + splitted_nx_tomo = splitter.split(data_slice=data_slice, nb_part=2) + assert len(splitted_nx_tomo) == 2 + # check the two dataset created + s_nx_tomo_1, s_nx_tomo_2 = splitted_nx_tomo + output_file_1 = os.path.join(folder, "output_file_1.nx") + + # data must contains a common section between the two nxtomo: the first 10 and last 10 frames + # then the rest must be splitted between the two NXtomo + def n_elmt(slice_): + return slice_.stop - slice_.start + + assert len(s_nx_tomo_1.instrument.detector.data) == 5 + assert n_elmt(s_nx_tomo_1.instrument.detector.data[0].data_slice()) == 10 + assert n_elmt(s_nx_tomo_1.instrument.detector.data[1].data_slice()) == 10 + assert n_elmt(s_nx_tomo_1.instrument.detector.data[2].data_slice()) == 20 + assert n_elmt(s_nx_tomo_1.instrument.detector.data[3].data_slice()) == 10 + assert n_elmt(s_nx_tomo_1.instrument.detector.data[4].data_slice()) == 10 + s_nx_tomo_1.save(output_file_1) + + output_file_2 = os.path.join(folder, "output_file_2.nx") + assert len(s_nx_tomo_2.instrument.detector.data) == 5 + assert n_elmt(s_nx_tomo_2.instrument.detector.data[0].data_slice()) == 10 + assert n_elmt(s_nx_tomo_2.instrument.detector.data[1].data_slice()) == 10 + assert n_elmt(s_nx_tomo_2.instrument.detector.data[2].data_slice()) == 20 + assert n_elmt(s_nx_tomo_2.instrument.detector.data[3].data_slice()) == 10 + assert n_elmt(s_nx_tomo_2.instrument.detector.data[4].data_slice()) == 10 + s_nx_tomo_2.save(output_file_2) + + # check final datasets are correctly formed + with h5py.File(output_file_1, mode="r") as h5f: + nx_1_data = h5f["entry0000/instrument/detector/data"][()] + assert nx_1_data.shape[0] == 60 + + # check final datasets are correctly formed + with h5py.File(output_file_2, mode="r") as h5f: + nx_2_data = h5f["entry0000/instrument/detector/data"][()] + assert nx_2_data.shape[0] == 60 + + # first 10 frames (common between the three nxtomo) + numpy.testing.assert_array_equal( + nx_1_data[0:10], + nx_2_data[0:10], + ) + numpy.testing.assert_array_equal( + nx_1_data[0:10], + original_data[0:10], + ) + + # last 10 frames (common between the three nxtomo) + numpy.testing.assert_array_equal( + nx_1_data[-10:], + nx_2_data[-10:], + ) + numpy.testing.assert_array_equal( + nx_1_data[-10:], + original_data[-10:], + ) + + # test nx_1_data unique region + numpy.testing.assert_array_equal( + nx_1_data[10:50], + original_data[10:50], + ) + + # test nx_2_data unique region + numpy.testing.assert_array_equal( + nx_2_data[10:50], + original_data[50:90], + ) diff --git a/nxtomomill/utils/utils.py b/nxtomomill/utils/utils.py index ab1884e234cb6558736bd3a6284fdfcb4a737782..35ca4efc9813b1a219ad676497896e02cfd3f869 100644 --- a/nxtomomill/utils/utils.py +++ b/nxtomomill/utils/utils.py @@ -29,7 +29,6 @@ __license__ = "MIT" __date__ = "29/04/2019" -import enum import typing import numpy import os @@ -40,6 +39,7 @@ import logging from silx.utils.deprecation import deprecated from .context import cwd_context from .frameappender import FrameAppender +from silx.utils.enum import Enum as _Enum try: import hdf5plugin # noqa F401 @@ -57,51 +57,6 @@ except ImportError: from tomoscan.esrf.hdf5scan import HDF5TomoScan -class Enum(enum.Enum): - """Enum with additional class methods.""" - - @classmethod - def from_value(cls, value): - """Convert a value to corresponding Enum member - - :param value: The value to compare to Enum members - If it is already a member of Enum, it is returned directly. - :return: The corresponding enum member - :rtype: Enum - :raise ValueError: In case the conversion is not possible - """ - if isinstance(value, cls): - return value - for member in cls: - if value == member.value: - return member - raise ValueError("Cannot convert: %s" % value) - - @classmethod - def members(cls): - """Returns a tuple of all members. - - :rtype: Tuple[Enum] - """ - return tuple(member for member in cls) - - @classmethod - def names(cls): - """Returns a tuple of all member names. - - :rtype: Tuple[str] - """ - return tuple(member.name for member in cls) - - @classmethod - def values(cls): - """Returns a tuple of all member values. - - :rtype: Tuple - """ - return tuple(member.value for member in cls) - - def embed_url(url: DataUrl, output_file: str) -> DataUrl: """ Create a dataset under duplicate_data and with a random name @@ -120,8 +75,7 @@ def embed_url(url: DataUrl, output_file: str) -> DataUrl: return url else: embed_data_path = "/".join(("/duplicate_data", str(uuid.uuid1()))) - with cwd_context(): - os.chdir(os.path.dirname(os.path.abspath(output_file))) + with cwd_context(os.path.dirname(os.path.abspath(output_file))): with HDF5File(output_file, "a") as h5s: h5s[embed_data_path] = get_data(url) h5s[embed_data_path].attrs["original_url"] = url.path() @@ -130,13 +84,13 @@ def embed_url(url: DataUrl, output_file: str) -> DataUrl: ) -class FileExtension(Enum): +class FileExtension(_Enum): H5 = ".h5" HDF5 = ".hdf5" NX = ".nx" -class Format(Enum): +class Format(_Enum): STANDARD = "standard" XRD_CT = "xrd-ct" XRD_3D = "3d-xrd"