diff --git a/radiography/__init__.py b/radiography/__init__.py new file mode 100644 index 0000000..0f2add8 --- /dev/null +++ b/radiography/__init__.py @@ -0,0 +1 @@ +from .operations import fits_to_variable diff --git a/radiography/operations.py b/radiography/operations.py new file mode 100644 index 0000000..81ec720 --- /dev/null +++ b/radiography/operations.py @@ -0,0 +1,46 @@ +import glob +import os +import numpy as np +import scipp as sc +from astropy.io import fits + + +def _load_fits(fit_dir): + if not os.path.isdir(fit_dir): + raise RuntimeError(fit_dir + " is not directory") + stack = [] + path_length = len(fit_dir) + 1 + filenames = sorted(glob.glob(fit_dir + "/*.fits")) + nfiles = len(filenames) + count = 0 + print(f"Loading {nfiles} files from '{fit_dir}'") + for filename in filenames: + count += 1 + print('\r{0}: Image {1}, of {2}'.format(filename[path_length:], count, nfiles), end="") + with fits.open(os.path.join(fit_dir, filename)) as hdu: + data = hdu[0].data + stack.append(np.flipud(data)) + + if len(stack) == 1: + # Fold away a dim + stack = stack[0] + + return np.array(stack, dtype=np.float64) + + +def fits_to_variable(fits_dir, average_data=False): + """ + Loads all fits images from the directory into a scipp Variable. + """ + stack = _load_fits(fits_dir) + if average_data: + stack = stack.mean(axis=0) + + if len(stack.shape) == 3: + return sc.Variable(["slice", "x", "y"], values=stack, variances=stack) + + elif len(stack.shape) == 2: + return sc.Variable(["x", "y"], values=stack, variances=stack) + else: + raise IndexError("Expected 2 or 3 dimensions," + f" found {len(stack.shape)}") diff --git a/scipp-radiography.ipynb b/scipp-radiography.ipynb new file mode 100644 index 0000000..073e345 --- /dev/null +++ b/scipp-radiography.ipynb @@ -0,0 +1,175 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "collapsed": true, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "# How to start\n", + "\n", + "Before starting you must:\n", + "- Ensure that `scipp` and `mantid` are on your `PYTHONPATH`.\n", + "- Generate the `config.py` file using `make_config.py`. Refer to the `README.md` or `python make_config.py --help` for information.\n", + "\n", + "For `scipp` and `mantid` follow instructions at: https://scipp.readthedocs.io/en/latest/getting-started/installation.html." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Some PYTHONPATH sanity checks\n", + "try:\n", + " import scipp\n", + " print(\"scipp found\")\n", + "except ImportError as e:\n", + " print(\"scipp is not available in the PYTHONPATH\")\n", + " raise e\n", + "\n", + "try:\n", + " import mantid\n", + " print(\"mantid found\")\n", + "except ImportError as e:\n", + " print(\"mantid is not available in the PYTHONPATH\")\n", + " raise e\n", + "\n", + "try:\n", + " import scippconfig\n", + " print(\"config found\")\n", + "except ImportError as e:\n", + " print(\"config is not available. Make sure you have generated it with `make_config.py`.\")\n", + " raise e" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "import os\n", + "import radiography\n", + "import scipp as sc\n", + "import numpy as np" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "experiment_dir = scippconfig.script_root\n", + "\n", + "data_dir = os.path.join(experiment_dir, 'radiography', 'IMAT_sample')\n", + "\n", + "# Since we don't have reference data pretend this folder contains our ref\n", + "dark_field_dir = os.path.join(data_dir, \"angle0\")\n", + "flat_field_dir = os.path.join(data_dir, \"Flats\", \"RB1730044\", \"Tomo_test\", \"OpenBeam_aft1\")\n", + "sample_dir = os.path.join(data_dir, \"angle1\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "source": [ + "raw_data = sc.Dataset()\n", + "raw_data[\"sample\"] = radiography.fits_to_variable(sample_dir)\n", + "\n", + "# Set coords based on sample dims, as reference needs to match\n", + "raw_data.coords[\"projection\"] = sc.Variable([\"projection\"], unit=sc.units.dimensionless, values=np.arange(0, raw_data[\"sample\"].shape[0]))\n", + "raw_data.coords[\"x\"] = sc.Variable([\"x\"], unit=sc.units.dimensionless, values=np.arange(0, raw_data[\"sample\"].shape[1]))\n", + "raw_data.coords[\"y\"] = sc.Variable([\"y\"], unit=sc.units.dimensionless, values=np.arange(0, raw_data[\"sample\"].shape[2]))" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + }, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "dark_field = radiography.fits_to_variable(dark_field_dir, average_data=True)\n", + "flat_field = radiography.fits_to_variable(flat_field_dir, average_data=True)\n", + "\n", + "# This blows through 32GB of memory as raw_data[\"sample\"] is already ~17 GB\n", + "raw_data[\"transmission\"] = (raw_data[\"sample\"] - dark_field) / (flat_field - dark_field)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": true + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# TODO WIP below\n", + "roi_x = (100, 200)\n", + "roi_y = (100, 200)\n", + "\n", + "raw_data[\"roi_sample\"] = raw_data[\"sample\"]\n", + "\n", + "raw_data.masks[\"sample\"] = sc.Variable(['x', 'y'], shape=(raw_data.coords[\"x\"].shape[0], raw_data.coords[\"y\"].shape[0]),\n", + " dtype=sc.dtype.bool)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file