{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create an NXtomo (from scratch)\n",
    "\n",
    "The goal of this tutorial is to build an [NXtomo](https://manual.nexusformat.org/classes/applications/NXtomo.html) file from scratch, without converting it directly from a BLISS scan.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example description\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create an NXtomo that matches the following sequence:\n",
    "\n",
    "| Frame index | Rotation angle (degrees) | Frame type | Control image key | Image key |\n",
    "| ----------- | ------------------------ | ---------- | ----------------- | --------- |\n",
    "| 0           | 0.0                      | Dark       | 2                 | 2         |\n",
    "| 1           | 0.0                      | Flat       | 1                 | 1         |\n",
    "| 2-201       | 0.0 - 89.9               | Projection | 0                 | 0         |\n",
    "| 202         | 90.0                     | Flat       | 1                 | 1         |\n",
    "| 203-402     | 90.0 - 180.0             | Projection | 0                 | 0         |\n",
    "| 403         | 180.0                    | Flat       | 1                 | 1         |\n",
    "| 404         | 90.0                     | Alignment  | -1                | 0         |\n"
   ]
  },
  {
   "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",
    "            To simplify the setup we create a single dark frame and a single flat frame that we reuse across the tutorial. Remember that these are raw frames.\n",
    "In real acquisitions you usually record several dark and flat frames, and they vary depending on when they are acquired.\n",
    "        </div>\n",
    "    </div>\n",
    "</div>\n"
   ]
  },
  {
   "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 projection frames\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phantom = shepp_logan_phantom()\n",
    "projections = {}\n",
    "proj_rotation_angles = numpy.linspace(0.0, 180.0, 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 dark frames\n"
   ]
  },
  {
   "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 flat frames\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat = numpy.ones((20, max_shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Add noise to radiographs\n"
   ]
  },
  {
   "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 alignment frames\n",
    "To keep things simple we reuse one of the previously generated radiographs.\n"
   ]
  },
  {
   "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": [
    "## Build an NXtomo that matches the sequence\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "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 includes:\n",
    "* detector frames: raw data\n",
    "* image-key (control): the frame type for each image (projections, flats, darks, alignment)\n",
    "* rotation angles: the rotation angle in degrees for every frame\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Detector frames\n",
    "\n",
    "To follow the sequence described earlier we create the series: dark, flat, first half of the projections, flat, second half of the projections, flat, and finally the alignment frame.\n",
    "\n",
    "All frames must be provided as a three-dimensional NumPy array.\n"
   ]
  },
  {
   "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",
    "    [\n",
    "        darks_stack,\n",
    "        flats_stack,\n",
    "        radios[:200],\n",
    "        flats_stack,\n",
    "        radios[200:],\n",
    "        flats_stack,\n",
    "        alignment_stack,\n",
    "    ]\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\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nxtomo.nxobject.nxdetector import ImageKey\n",
    "\n",
    "image_key_control = numpy.concatenate(\n",
    "    [\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",
    "\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 angles\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pint\n",
    "\n",
    "ureg = pint.get_application_registry()\n",
    "\n",
    "rotation_angle = numpy.concatenate(\n",
    "    [\n",
    "        [\n",
    "            0.0,\n",
    "        ],\n",
    "        [\n",
    "            0.0,\n",
    "        ],\n",
    "        proj_rotation_angles[:200],\n",
    "        [\n",
    "            90.0,\n",
    "        ],\n",
    "        proj_rotation_angles[200:],\n",
    "        [\n",
    "            180.0,\n",
    "        ],\n",
    "        [\n",
    "            90.0,\n",
    "        ],\n",
    "    ]\n",
    ")\n",
    "assert len(rotation_angle) == len(data)\n",
    "# register rotation angle to the sample\n",
    "my_nxtomo.sample.rotation_angle = rotation_angle * ureg.degree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Field of view\n",
    "The field of view can be either `Half` or `Full`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_nxtomo.instrument.detector.field_of_view = \"Full\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pixel size\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_nxtomo.instrument.detector.x_pixel_size = (\n",
    "    my_nxtomo.instrument.detector.y_pixel_size\n",
    ") = (10 * ureg.micrometer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a unit is provided it is stored as a dataset attribute. The software that reads the NXtomo file must interpret that unit.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save the NXtomo to disk\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\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 the saved data\n",
    "Use the tomoscan validator to ensure we have enough information for nabu to process the dataset.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import tomoscan  # NOQA\n",
    "except ImportError:\n",
    "    has_tomoscan = False\n",
    "else:\n",
    "    from tomoscan.esrf import NXtomoScan\n",
    "    from tomoscan.validator import ReconstructionValidator\n",
    "\n",
    "    has_tomoscan = True\n",
    "\n",
    "if has_tomoscan:\n",
    "    scan = NXtomoScan(nx_tomo_file_path, entry=\"entry\")\n",
    "    validator = ReconstructionValidator(\n",
    "        scan, check_phase_retrieval=False, check_values=True\n",
    "    )\n",
    "    assert validator.is_valid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also inspect the file layout to confirm it looks correct.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from h5glance import H5Glance\n",
    "\n",
    "H5Glance(nx_tomo_file_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A good practice is to verify the frames, image keys, and rotation angles to ensure the values are consistent.\n"
   ]
  },
  {
   "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 can reconstruct it with [nabu](https://gitlab.esrf.fr/tomotools/nabu).\n",
    "\n",
    "Create a nabu configuration file for contrast tomography, named `nabu-ct.conf`, to reconstruct one slice of the volume.\n"
   ]
  },
  {
   "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 the configuration file make sure to disable `take_logarithm`, because the dataset already contains processed values.\n",
    "        </div>\n",
    "    </div>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `nabu` is installed you can run it:\n"
   ]
  },
  {
   "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",
    "To perform phase tomography you must also record:\n",
    "* incoming beam energy (in keV)\n",
    "* sample-to-detector distance (in metres)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can reuse the existing `my_nxtomo` and add the missing elements.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_nxtomo.energy = 12.5 * ureg.keV  # in keV by default\n",
    "my_nxtomo.instrument.detector.distance = 0.2 * ureg.meter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then you can reconstruct it with phase retrieval by modifying the nabu configuration file accordingly.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Provide more metadata\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also store the sample translations along x, y, and z during the acquisition.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_nxtomo.sample.x_translation = [0, 12] * ureg.millimeter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add useful contextual information such as the sample name, source details, and acquisition start and end times.\n"
   ]
  },
  {
   "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",
    "\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": ".venv",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}