{
 "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": [
    "<div style=\"width:image width px; text-align:center; vertical-align:middle; border:solid; background-color:#5cadad\">\n",
    "    <div class=\"two-columns-grid\">\n",
    "        <img src=\"img/note.svg\"\n",
    "             style=\"width: 20px\"\n",
    "             align=\"left\"\n",
    "            />  \n",
    "        <div style=\"text-align:left; vertical-align:middle\">\n",
    "            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 it 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.\n",
    "        </div>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "source": [
    "## create an NXtomo that fits the sequence we want"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import nxtomo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from 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.instrument.detector.data = data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### image key control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nxtomo.nxobject.nxdetector 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 indexes are\", 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": 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.instrument.detector.field_of_view = \"Full\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### pixel size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.makedirs(\"output\", exist_ok=True)\n",
    "nx_tomo_file_path = os.path.join(\"output\", \"nxtomo.nx\")\n",
    "my_nxtomo.save(file_path=nx_tomo_file_path, data_path=\"entry\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### check data saved\n",
    "We can use some validator from tomoscan to insure we have enought data to be treated by nabu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import tomoscan\n",
    "except ImportError:\n",
    "    has_tomoscan = False\n",
    "else:\n",
    "    from tomoscan.esrf import NXtomoScan\n",
    "    from tomoscan.validator import ReconstructionValidator\n",
    "    has_tomoscan = True\n",
    "\n",
    "if has_tomoscan:\n",
    "    scan = NXtomoScan(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 output/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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"width:image width px; text-align:center; vertical-align:middle; border:solid; background-color:#5cadad\">\n",
    "    <div class=\"two-columns-grid\">\n",
    "        <img src=\"img/note.svg\"\n",
    "             style=\"width: 20px\"\n",
    "             align=\"left\"\n",
    "            />  \n",
    "        <div style=\"text-align:left; vertical-align:middle\">\n",
    "            on the configuration you must disable take_logarithm due to the dataset.\n",
    "        </div>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_nxtomo.energy = 12.5  # in keV by default\n",
    "my_nxtomo.instrument.detector.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.sample.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.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.end_time = datetime(2022, 2, 27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_nxtomo.save(\n",
    "    file_path=nx_tomo_file_path,\n",
    "    data_path=\"entry\",\n",
    "    overwrite=True,\n",
    ")"
   ]
  }
 ],
 "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}