Merge pull request #42 from scipp/normalization

Normalization step in the YMIR image workflow.

Author: YooSunYoung

docs/api-reference/index.md

@@ -26,4 +26,8 @@
    :toctree: ../generated/modules
    :template: module-template.rst
    :recursive:
+
+   normalize
+   io
+   workflow
 ```

docs/user-guide/histogram_mode_detector.ipynb

@@ -6,6 +6,77 @@
    "source": [
     "# Histogram Mode Detector\n",
     "\n",
+    "## TL;DR - Get normalized sample images averaged along the rotation angle."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ess.imaging.data import get_ymir_images_path\n",
+    "from ess.imaging.io import FilePath\n",
+    "from ess.imaging.normalize import NormalizedSampleImages\n",
+    "from ess.imaging.workflow import (\n",
+    "    YmirImageNormalizationWorkflow,\n",
+    "    ImageDetectorName,\n",
+    "    RotationMotionSensorName,\n",
+    ")\n",
+    "\n",
+    "wf = YmirImageNormalizationWorkflow()\n",
+    "wf[FilePath] = (\n",
+    "    get_ymir_images_path()\n",
+    ")  # Replace with the path to your images in nexus file.\n",
+    "wf[ImageDetectorName] = 'orca'\n",
+    "wf[RotationMotionSensorName] = 'motion_cabinet_2'\n",
+    "normalized = wf.compute(NormalizedSampleImages)\n",
+    "normalized"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "normalized['time', 0].plot(\n",
+    "    title=f'Normalized Sample Image at Rotation Angle of {normalized.coords[\"rotation_angle\"].values[0]:.3f}'\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can use a slicer to browse images along the rotation angle.\n",
+    "\n",
+    "```python\n",
+    "%matplotlib widget\n",
+    "from plopp import slicer\n",
+    "\n",
+    "slicer(normalized)\n",
+    "```\n",
+    "\n",
+    "If you want to bin the images by linear steps of angles,\n",
+    "you can use ``sc.groupby`` like below.\n",
+    "\n",
+    "```python\n",
+    "%matplotlib widget\n",
+    "import scipp as sc\n",
+    "from plopp import slicer\n",
+    "\n",
+    "slicer(normalized.groupby(\n",
+    "    'rotation_angle',\n",
+    "    bins=sc.linspace(dim='rotation_angle', start=-0.5, stop=4.5, num=5, unit='deg'),\n",
+    ").mean('rotation_angle'))\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "## IO\n",
     "\n",
     "I/O providers can load white beam images and slice them according to the image key.\n",
@@ -14,7 +85,9 @@
     "\n",
     "Instead, we slice the data according to the ``image_key`` and separate image stacks earlier due to performance limitation of 3-dimensional data binning.\n",
     "\n",
-    "> For example, for 400 images of 2_048 x 2_048 pixels couldn't be done in the regular laptop."
+    "> For example, for 400 images of 2_048 x 2_048 pixels couldn't be done in the regular laptop.\n",
+    "\n",
+    "### Load the data as a stack of images"
    ]
   },
   {
@@ -26,19 +99,19 @@
     "from ess.imaging.data import get_ymir_images_path\n",
     "from ess.imaging.io import (\n",
     "    FilePath,\n",
-    "    SampleImageStacks,\n",
+    "    SampleImageStacksWithLogs,\n",
     "    RawSampleImageStacks,\n",
     "    OpenBeamImageStacks,\n",
     "    DarkCurrentImageStacks,\n",
     "    AllImageStacks,\n",
     ")\n",
     "from ess.imaging.workflow import (\n",
-    "    YmirWorkflow,\n",
+    "    YmirImageNormalizationWorkflow,\n",
     "    ImageDetectorName,\n",
     "    RotationMotionSensorName,\n",
     ")\n",
     "\n",
-    "wf = YmirWorkflow()\n",
+    "wf = YmirImageNormalizationWorkflow()\n",
     "wf[FilePath] = get_ymir_images_path()\n",
     "wf[ImageDetectorName] = 'orca'\n",
     "wf[RotationMotionSensorName] = 'motion_cabinet_2'"
@@ -55,7 +128,7 @@
     "    RawSampleImageStacks,\n",
     "    OpenBeamImageStacks,\n",
     "    DarkCurrentImageStacks,\n",
-    "    SampleImageStacks,\n",
+    "    SampleImageStacksWithLogs,\n",
     ")\n",
     "wf.visualize(image_stack_types)"
    ]
@@ -67,14 +140,14 @@
    "outputs": [],
    "source": [
     "results = wf.compute(image_stack_types)\n",
-    "results[SampleImageStacks]"
+    "results[SampleImageStacksWithLogs]"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Save Individual Frame"
+    "### Save Individual Frame"
    ]
   },
   {
@@ -117,7 +190,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Save All Frames in One File"
+    "### Save All Frames in One File"
    ]
   },
   {
@@ -167,11 +240,43 @@
    "source": [
     "!rm -rf images"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Normalization\n",
+    "\n",
+    "See :attr:ess.imaging.normalize.normalize_sample_images for more details."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ess.imaging.normalize import NormalizedSampleImages\n",
+    "\n",
+    "wf[OpenBeamImageStacks] = results[OpenBeamImageStacks]\n",
+    "wf[DarkCurrentImageStacks] = results[DarkCurrentImageStacks]\n",
+    "wf[SampleImageStacksWithLogs] = results[SampleImageStacksWithLogs]\n",
+    "wf.visualize(NormalizedSampleImages)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wf.compute(NormalizedSampleImages)"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "img-dev-310",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -189,5 +294,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

src/ess/imaging/io.py

@@ -22,14 +22,18 @@
     RotationMotionSensorName,
 )
 
+FileLock = NewType("FileLock", bool)
+"""File lock mode for reading nexus file."""
+DEFAULT_FILE_LOCK = FileLock(True)
+
 HistogramModeDetector = NewType("HistogramModeDetector", sc.DataGroup)
 """Histogram mode detector data group."""
 HistogramModeDetectorData = NewType("HistogramModeDetectorData", sc.DataArray)
 """Histogram mode detector data."""
 ImageKeyCoord = NewType("ImageKeyCoord", sc.Variable)
 """Image key coordinate."""
-SampleImageStacks = NewType("SampleImageStacks", sc.DataArray)
-"""Image stacks separated by ImageKey values via timestamp."""
+SampleImageStacksWithLogs = NewType("SampleImageStacksWithLogs", sc.DataArray)
+"""Raw image stacks separated by ImageKey values via timestamp."""
 RotationAngleCoord = NewType("RotationAngleCoord", sc.Variable)
 """Rotation angle coordinate."""
 
@@ -79,10 +83,20 @@ def load_nexus_histogram_mode_detector(
     file_path: FilePath,
     image_detector_name: ImageDetectorName,
     histogram_mode_detectors_path: HistogramModeDetectorsPath = DEFAULT_HISTOGRAM_PATH,
+    locking: FileLock = DEFAULT_FILE_LOCK,
 ) -> HistogramModeDetector:
-    with snx.File(file_path, mode="r") as f:
-        img_path = f"{histogram_mode_detectors_path}/{image_detector_name}"
-        dg: sc.DataGroup = f[img_path][()]
+    try:
+        with snx.File(file_path, mode="r", locking=locking) as f:
+            img_path = f"{histogram_mode_detectors_path}/{image_detector_name}"
+            dg: sc.DataGroup = f[img_path][()]
+    except PermissionError as e:
+        raise PermissionError(
+            f"Permission denied to read the nexus file [{file_path}]. "
+            "Please check the permission of the file or the directory. "
+            "Consider using the `file_lock` parameter to avoid file locking "
+            "if the file system is mounted on a network file system. "
+            "and it is safe to read the file without locking."
+        ) from e
 
     # Manually assign unit to the histogram detector mode data
     img: sc.DataArray = dg['data']
@@ -134,9 +148,10 @@ def separate_image_key_logs(*, dg: HistogramModeDetector) -> ImageKeyLogs:
 def load_nexus_rotation_logs(
     file_path: FilePath,
     motion_sensor_name: RotationMotionSensorName,
+    locking: FileLock = DEFAULT_FILE_LOCK,
 ) -> RotationLogs:
     log_path = f"entry/instrument/{motion_sensor_name}/rotation_stage_readback"
-    with snx.File(file_path, mode="r") as f:
+    with snx.File(file_path, mode="r", locking=locking) as f:
         return RotationLogs(f[log_path][()]['value'])
 
 
@@ -230,7 +245,7 @@ def retrieve_sample_images(
 
 def apply_logs_as_coords(
     samples: RawSampleImageStacks, rotation_angles: RotationLogs
-) -> SampleImageStacks:
+) -> SampleImageStacksWithLogs:
     # Make sure the data has the same range as the rotation angle coordinate
     min_log_time = rotation_angles.coords[TIME_COORD_NAME].min(TIME_COORD_NAME)
     sliced = samples[TIME_COORD_NAME, min_log_time:].copy(deep=False)
@@ -241,7 +256,7 @@ def apply_logs_as_coords(
         )
     rotation_angle_coord = derive_log_coord_by_range(samples, rotation_angles)
     sliced.coords['rotation_angle'] = rotation_angle_coord
-    return SampleImageStacks(sliced)
+    return SampleImageStacksWithLogs(sliced)
 
 
 DEFAULT_IMAGE_NAME_PREFIX_MAP = {
@@ -256,7 +271,7 @@ def dummy_progress_wrapper(core_iterator: Iterable) -> Iterable:
 
 
 def _save_merged_images(
-    *, image_stacks: SampleImageStacks, image_prefix: str, output_dir: Path
+    *, image_stacks: SampleImageStacksWithLogs, image_prefix: str, output_dir: Path
 ) -> None:
     image_path = output_dir / Path(
         f"{image_prefix}_0000_{image_stacks.sizes['time']:04d}.tiff"
@@ -266,7 +281,7 @@ def _save_merged_images(
 
 def _save_individual_images(
     *,
-    image_stacks: SampleImageStacks,
+    image_stacks: SampleImageStacksWithLogs,
     image_prefix: str,
     output_dir: Path,
     progress_wrapper: Callable[[Iterable], Iterable] = dummy_progress_wrapper,
@@ -332,13 +347,13 @@ def export_image_stacks_as_tiff(
     for image_key, cur_images in progress_wrapper(image_stacks.items()):
         if merge_image_by_key:
             _save_merged_images(
-                image_stacks=SampleImageStacks(cur_images),
+                image_stacks=SampleImageStacksWithLogs(cur_images),
                 image_prefix=image_prefix_map[image_key],
                 output_dir=output_path,
             )
         else:
             _save_individual_images(
-                image_stacks=SampleImageStacks(cur_images),
+                image_stacks=SampleImageStacksWithLogs(cur_images),
                 image_prefix=image_prefix_map[image_key],
                 output_dir=output_path,
                 progress_wrapper=progress_wrapper,