{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Manipulating Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "## Notebook bootstrap\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "shared_dir = next(path / \"shared\" for path in [Path.cwd(), *Path.cwd().parents] if (path / \"shared\").is_dir())\n",
    "if str(shared_dir) not in sys.path:\n",
    "    sys.path.append(str(shared_dir))\n",
    "\n",
    "from demo_utils import unzip_folder, mpr_plot\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "## Import Modules\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imfusion\n",
    "import numpy as np\n",
    "from copy import deepcopy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Image geometry\n",
    "A `SharedImage` defines the orientation and size of an image using three coordinate systems:\n",
    "- **Pixel Space**: The discrete voxel grid, where voxels are indexed by: (`i` in `[0, Width]`, `j` in `[0, Height]`, `k` in `[0, Slices]`).\n",
    "- **Image Space**: The voxel grid expressed in metric units (mm), with the image center located at the origin: `O = (x, y, z) = (0, 0, 0)`.\n",
    "The voxel size is defined by the image `spacing`. Image coordinates range from `(-SizeX/2, -SizeY/2, -SizeZ/2)` to `(SizeX/2, SizeY/2, SizeZ/2)`.\n",
    "- **World Space**: The physical coordinate system. Each `SharedImage` stores a transformation matrix that defines the image's orientation and position by mapping image-space coordinates to world-space coordinates.\n",
    "\n",
    "![Coordinate Systems](images/CoordinateSystems.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-05T14:10:19.507578Z",
     "start_time": "2025-03-05T14:10:19.143800Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image matrix:\n",
      " [[  1.          -0.           0.           0.00001776]\n",
      " [  0.           1.          -0.           0.00001776]\n",
      " [ -0.           0.           1.         414.51000977]\n",
      " [  0.          -0.          -0.           1.        ]]\n",
      "Image spacing: Spacing(x=3.6458332538605, y=3.6458332538605, z=3.270000086421467)\n",
      "Image dimensions: [192 192 227]\n",
      "Image extent: [699.99998474 699.99998474 742.29001962]\n",
      "Image coordinate of Voxel (0, 0, 0) [-348.17707574 -348.17707574 -369.51000977]\n",
      "Voxel coordinate of Image Center (0., 0., 0.) [ 95.5  95.5 113. ]\n",
      "World coordinate (aka Position) of Voxel (0, 0, 0): [-348.17709351 -348.17709351 -784.02001953]\n",
      "World coordinate of Image Center (0., 0., 0.): [  -0.00001776   -0.00001776 -414.51000977    1.        ]\n"
     ]
    }
   ],
   "source": [
    "unzipped_folder = unzip_folder(\"../shared/data/pet-ct-rtstruct.zip\")\n",
    "\n",
    "pet_sis, *_ = imfusion.io.load('../shared/data/pet-ct-rtstruct/pet.imf')\n",
    "pet_image = pet_sis[0]\n",
    "\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "print('Image matrix:\\n', pet_image.world_to_image_matrix)\n",
    "print('Image spacing:', pet_image.descriptor.spacing)\n",
    "print('Image dimensions:', pet_image.descriptor.dimensions)\n",
    "print('Image extent:', pet_image.descriptor.extent)\n",
    "print('Image coordinate of Voxel (0, 0, 0)', pet_image.descriptor.pixel_to_image([0, 0, 0]))\n",
    "print('Voxel coordinate of Image Center (0., 0., 0.)', pet_image.descriptor.image_to_pixel([0., 0., 0.]))\n",
    "\n",
    "# get world coordinates of Voxel (0, 0, 0):\n",
    "print(f'World coordinate (aka Position) of Voxel (0, 0, 0): {pet_image.descriptor_world.pixel_to_world([0, 0, 0])}')\n",
    "print(f'World coordinate of Image Center (0., 0., 0.): {np.matmul(pet_image.descriptor_world.matrix_to_world, np.array([0., 0., 0., 1.]))}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Basic image processing\n",
    "\n",
    "Many ImFusion SDK algorithms are available as free functions in `imfusion.processing`. This tutorial uses image resampling as an example, but you can explore `imfusion.processing` to discover the full range of available Python algorithms.\n",
    "\n",
    "> **NOTE**: Algorithms without a dedicated Python function can still be executed through the generic `imfusion.algorithm.execute(...)` or `imfusion.algorithm.create(...)` interface.\n",
    "\n",
    "For example, the following code resamples the image to a spacing of 2 mm in each dimension. The corresponding geometry information stored in the SharedImage is updated automatically. The optional arguments are shown with their default values and control how the resampling is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-05T14:10:21.420649Z",
     "start_time": "2025-03-05T14:10:19.711429Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "before:\t imfusion.BoundImageDescriptor(dimensions: [w:192, h:192, s:227, c:1] | spacing:[3.64583, 3.64583, 3.27] | isMetric:true | type:FLOAT | memory: 31.92 MB)\n",
      "after:\t imfusion.BoundImageDescriptor(dimensions: [w:350, h:350, s:371, c:1] | spacing:[2, 2, 2] | isMetric:true | type:FLOAT | memory: 173.4 MB)\n"
     ]
    }
   ],
   "source": [
    "print(\"before:\\t\", pet_sis[0].descriptor)\n",
    "\n",
    "resampled_pet = imfusion.processing.resample_to_spacing(\n",
    "    pet_sis,\n",
    "    target_spacing=[2., 2., 2.],  # mm\n",
    "    # reduction_mode=imfusion.ReductionMode.AVERAGE,  # Defines how contributing input voxels are combined.\n",
    "    # interpolation_mode=imfusion.InterpolationMode.LINEAR,  # Defines how values are sampled between voxel centers.\n",
    "    # keep_zero_values=False,  # Controls whether zero-valued input voxels remain zero after resampling.\n",
    "    # preserve_extent=False,  # Controls whether output spacing is adjusted to preserve the input extent.\n",
    "    # device_policy=imfusion.processing.DevicePolicy.PREFER_GPU,  # Selects the preferred processing device.\n",
    "    # mode=imfusion.processing.Mode.CREATE_NEW,  # Selects whether processing modifies the input or creates a new image.\n",
    ")\n",
    "\n",
    "imfusion.io.save([resampled_pet], 'output/resampled_pet.imf')\n",
    "\n",
    "print(\"after:\\t\", resampled_pet[0].descriptor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "Since `SharedImage`s store orientation and spacing information, they can also be resampled to the geometry of another reference image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-05T14:10:22.153817Z",
     "start_time": "2025-03-05T14:10:21.430536Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "original CT: imfusion.BoundImageDescriptor(dimensions: [w:167, h:167, s:247, c:1] | spacing:[3, 3, 3] | type:INT | memory: 26.28 MB)\n",
      "\n",
      "CT: imfusion.BoundImageDescriptor(dimensions: [w:350, h:350, s:371, c:1] | spacing:[2, 2, 2] | isMetric:true | type:INT | memory: 173.4 MB)\n",
      "PET: imfusion.BoundImageDescriptor(dimensions: [w:350, h:350, s:371, c:1] | spacing:[2, 2, 2] | isMetric:true | type:FLOAT | memory: 173.4 MB)\n"
     ]
    }
   ],
   "source": [
    "(ct_sis,) = imfusion.io.load('../shared/data/pet-ct-rtstruct/ct_base.imf')\n",
    "\n",
    "print(f'original CT: {ct_sis[0].descriptor}\\n')\n",
    "\n",
    "resampled_ct = imfusion.processing.resample_to_ref_image(ct_sis, resampled_pet)\n",
    "\n",
    "imfusion.io.save([resampled_ct], 'output/resampled_ct.imf')\n",
    "\n",
    "\n",
    "print(f'CT: {resampled_ct[0].descriptor}')\n",
    "print(f'PET: {resampled_pet[0].descriptor}')"
   ]
  }
 ],
 "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}