Repository: agencyenterprise/neurotechdevkit (language: python)

Issue:
Add ability to use custom CT scans
## Describe the new feature or enhancement

Add option to load CT scans from NIfTI or DICOM formats to be used within a simulation

## Describe your proposed implementation

This would be similar to what is currently used in scenario 2, where the skull is loaded and used as a (binary) mask for where the skull is, and material properties are homogeneous within the skull; everything inside the skull is brain, outside is water. 


Discussion / hints:
For someone just learning how tweaking parameters affect focused ultrasound, it could also be nice to have an image-interface (e.g., with a bitmap image) so that someone could manipulate the images from one of our examples:

![image](https://github.com/agencyenterprise/neurotechdevkit/assets/3221512/ed9a8f53-96c0-4a53-8646-9081e108c610)

e.g., just in MS Paint, and then re-upload that image (plus a legend) as a scenario.

Fix patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 96102a68..278f9f8d 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -26,7 +26,7 @@ jobs:
       - run:
           name: Install stride
           command: |
-            poetry run pip install git+https://github.com/trustimaging/stride
+            poetry run pip install git+https://github.com/trustimaging/stride@2520c59
 
       - restore_cache:
           name: Restore docs cache
@@ -73,7 +73,7 @@ jobs:
       - run:
           name: Install stride
           command: |
-            poetry run pip install git+https://github.com/trustimaging/stride
+            poetry run pip install git+https://github.com/trustimaging/stride@2520c59
 
       - run:
           name: Build docs
diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml
index bf3cf14f..5fa4fb50 100644
--- a/.github/workflows/lint.yml
+++ b/.github/workflows/lint.yml
@@ -48,7 +48,7 @@ jobs:
         run: poetry install --no-interaction --no-root
 
       - name: Install stride
-        run: poetry run pip install git+https://github.com/trustimaging/stride
+        run: poetry run pip install git+https://github.com/trustimaging/stride@2520c59
 
       - name: Install napari
         run: poetry run pip install "napari[all,pyqt6_experimental]"
diff --git a/Dockerfile b/Dockerfile
index 11f3fca5..bafdf1cd 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -17,7 +17,7 @@ RUN python3 -m venv ./venv && \
 
 ENV POETRY_VIRTUALENVS_CREATE=false
 RUN ./venv/bin/poetry install && \
-    ./venv/bin/pip install git+https://github.com/trustimaging/stride
+    ./venv/bin/pip install git+https://github.com/trustimaging/stride@2520c59
 
 FROM python:3.10.0-slim
 
diff --git a/Dockerfile-gpu b/Dockerfile-gpu
index 5a03a4e6..76ce1598 100644
--- a/Dockerfile-gpu
+++ b/Dockerfile-gpu
@@ -16,7 +16,7 @@ RUN python3 -m venv ./venv && \
 
 ENV POETRY_VIRTUALENVS_CREATE=false
 RUN ./venv/bin/poetry install && \
-    ./venv/bin/pip install git+https://github.com/trustimaging/stride
+    ./venv/bin/pip install git+https://github.com/trustimaging/stride@2520c59
 
 
 FROM nvcr.io/nvidia/nvhpc:23.5-devel-cuda_multi-ubuntu22.04 
diff --git a/Makefile b/Makefile
index dcfb98a7..265f8310 100644
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,7 @@
-.PHONY:help lint lint-check test test-coverage test-unit test-integration docs web web-development
+.PHONY:help lint lint-check test test-coverage test-unit test-integration docs web
 
 help:
-	@echo "Available commands are: \n*lint, lint-check, spellcheck, test, test-unit, test-integration docs web web-development"
+	@echo "Available commands are: \n*lint, lint-check, spellcheck, test, test-unit, test-integration docs web"
 
 lint:
 	poetry run isort src tests docs/examples
@@ -39,8 +39,5 @@ test-integration:
 docs:
 	poetry run mkdocs build
 
-web-development:
-	poetry run ndk-ui
-
 web:
-	poetry run gunicorn -b 0.0.0.0 -w 1 --timeout 99999999 'src.web.app:app'
\ No newline at end of file
+	poetry run ndk-ui
diff --git a/docs/contributing.md b/docs/contributing.md
index 14a84177..37658c6a 100644
--- a/docs/contributing.md
+++ b/docs/contributing.md
@@ -70,7 +70,7 @@ This will resolve and install the dependencies from `poetry.lock` and will insta
 Install stride with:
 
 ```bash
-$ poetry run pip install git+https://github.com/trustimaging/stride
+$ poetry run pip install git+https://github.com/trustimaging/stride@2520c59
 ```
 
 Follow the steps described in [Setting up a compiler](installation.md#setting-up-a-compiler).
diff --git a/docs/installation.md b/docs/installation.md
index 5c6026bc..331d13c8 100644
--- a/docs/installation.md
+++ b/docs/installation.md
@@ -33,12 +33,12 @@ pip install neurotechdevkit
 You also have to install stride, it can be done running:
 
 ```bash
-pip install git+https://github.com/trustimaging/stride
+pip install git+https://github.com/trustimaging/stride@2520c59
 ```
 
 ## Setting up a compiler
 
-NDK uses [devito](https://www.devitoproject.org/devito/) to perform the heavy computational operations. Devito generates, compiles and runs C code to achieve better performance.
+NDK uses [devito](https://www.devitoproject.org/) to perform the heavy computational operations. Devito generates, compiles and runs C code to achieve better performance.
 The compiler used by Devito has to be selected, and paths for the linker might also be added as environment variables.
 
 As a last step **before running NDK**, follow the instructions below depending on your OS.
diff --git a/docs/user_interface.md b/docs/user_interface.md
index 0d0248cf..bf0efafc 100644
--- a/docs/user_interface.md
+++ b/docs/user_interface.md
@@ -33,7 +33,7 @@ There are multiple ways of running the web server, you can choose between:
 1. Open a terminal in the project location and run the following command
 
     ```
-    make web-development
+    make web
     ```
 
 ### Docker
diff --git a/poetry.lock b/poetry.lock
index 59b6f821..123d3620 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1165,26 +1165,6 @@ files = [
 [package.dependencies]
 colorama = ">=0.4"
 
-[[package]]
-name = "gunicorn"
-version = "21.2.0"
-description = "WSGI HTTP Server for UNIX"
-optional = false
-python-versions = ">=3.5"
-files = [
-    {file = "gunicorn-21.2.0-py3-none-any.whl", hash = "sha256:3213aa5e8c24949e792bcacfc176fef362e7aac80b76c56f6b5122bf350722f0"},
-    {file = "gunicorn-21.2.0.tar.gz", hash = "sha256:88ec8bff1d634f98e61b9f65bc4bf3cd918a90806c6f5c48bc5603849ec81033"},
-]
-
-[package.dependencies]
-packaging = "*"
-
-[package.extras]
-eventlet = ["eventlet (>=0.24.1)"]
-gevent = ["gevent (>=1.4.0)"]
-setproctitle = ["setproctitle"]
-tornado = ["tornado (>=0.2)"]
-
 [[package]]
 name = "h5py"
 version = "3.9.0"
@@ -2761,6 +2741,35 @@ files = [
     {file = "nest_asyncio-1.5.8.tar.gz", hash = "sha256:25aa2ca0d2a5b5531956b9e273b45cf664cae2b145101d73b86b199978d48fdb"},
 ]
 
+[[package]]
+name = "nibabel"
+version = "5.1.0"
+description = "Access a multitude of neuroimaging data formats"
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "nibabel-5.1.0-py3-none-any.whl", hash = "sha256:b3deb8130c835b9d26e80880b0d5e443d9e3f30972b3b0302dd2fafa3ca629f8"},
+    {file = "nibabel-5.1.0.tar.gz", hash = "sha256:ce73ca5e957209e7219a223cb71f77235c9df2acf4d3f27f861ba38e9481ac53"},
+]
+
+[package.dependencies]
+numpy = ">=1.19"
+packaging = ">=17"
+
+[package.extras]
+all = ["nibabel[dev,dicomfs,doc,minc2,spm,style,test,zstd]"]
+dev = ["gitpython", "nibabel[style]", "twine"]
+dicom = ["pydicom (>=1.0.0)"]
+dicomfs = ["nibabel[dicom]", "pillow"]
+doc = ["matplotlib (>=1.5.3)", "numpydoc", "sphinx (>=5.3,<6.0)", "texext", "tomli"]
+doctest = ["nibabel[doc,test]"]
+minc2 = ["h5py"]
+spm = ["scipy"]
+style = ["blue", "flake8", "isort"]
+test = ["coverage", "pytest (!=5.3.4)", "pytest-cov", "pytest-doctestplus", "pytest-httpserver", "pytest-xdist"]
+typing = ["importlib-resources", "mypy", "pydicom", "pytest", "pyzstd", "types-pillow", "types-setuptools"]
+zstd = ["pyzstd (>=0.14.3)"]
+
 [[package]]
 name = "nodeenv"
 version = "1.8.0"
@@ -3332,6 +3341,20 @@ typing-extensions = ">=4.2.0"
 dotenv = ["python-dotenv (>=0.10.4)"]
 email = ["email-validator (>=1.0.3)"]
 
+[[package]]
+name = "pydicom"
+version = "2.4.3"
+description = "A pure Python package for reading and writing DICOM data"
+optional = false
+python-versions = ">=3.7"
+files = [
+    {file = "pydicom-2.4.3-py3-none-any.whl", hash = "sha256:797e84f7b22e5f8bce403da505935b0787dca33550891f06495d14b3f6c70504"},
+    {file = "pydicom-2.4.3.tar.gz", hash = "sha256:51906e0b9fb6e184a0f56298cb43ed716b7cf7edc00f6b71d5c769bc1f982402"},
+]
+
+[package.extras]
+docs = ["matplotlib", "numpy", "numpydoc", "pillow", "sphinx", "sphinx-copybutton", "sphinx-gallery", "sphinx_rtd_theme", "sphinxcontrib-napoleon"]
+
 [[package]]
 name = "pydocstyle"
 version = "6.3.0"
@@ -4867,4 +4890,4 @@ testing = ["big-O", "jaraco.functools", "jaraco.itertools", "more-itertools", "p
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.9,<3.12"
-content-hash = "15b42644a1569d3f93cca06af06265bdbf9148c506cc1f10d37ab53ca9947d76"
+content-hash = "ddf993bdb9dda1bfd6ef8dcab20f4287e611c09781de9c26930f124a6c641f58"
diff --git a/pyproject.toml b/pyproject.toml
index 38129113..82b136e0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "neurotechdevkit"
-version = "v0.5.3"
+version = "v0.6.0"
 description = "Neurotech Development Kit: an open-source software library designed to enhance accessibility to cutting-edge neurotechnology"
 authors = ["AE Studio <bci@ae.studio>"]
 maintainers = ["AE Studio <bci@ae.studio>"]
@@ -43,8 +43,9 @@ mkdocs-gallery = "0.7.6"
 xarray = ">=2022.3.0"
 flask = {extras = ["async"], version = "^2.3.3"}
 pydantic = "<2"
-gunicorn = "^21.2.0"
 pyrevolve = ">=2.2.3"
+pydicom = "^2.4.3"
+nibabel = "^5.1.0"
 
 [tool.poetry.group.dev.dependencies]
 black = {extras = ["jupyter"], version = "^23.9.1"}
diff --git a/src/neurotechdevkit/grid.py b/src/neurotechdevkit/grid.py
index 15657c2e..e1f92b38 100644
--- a/src/neurotechdevkit/grid.py
+++ b/src/neurotechdevkit/grid.py
@@ -34,20 +34,45 @@ class Grid(stride.Grid):
 
     @staticmethod
     def make_shaped_grid(
-        shape: Tuple[int, int],
-        spacing: float,
+        shape: Tuple[int, ...],
+        spacing: Union[float, Tuple[float, ...]],
+        extra: Union[int, Iterable[int]] = 50,
+        absorbing: Union[int, Iterable[int]] = 40,
     ) -> "Grid":
-        """Create an NDK Grid with the given shape.
+        """Create an NDK Grid with a given shape.
+
+        Note that the time component of the grid is not defined here. That is
+        created at simulation time because it depends on simulation parameters.
 
         Args:
-            shape: a tuple of two integers representing the shape of the grid.
+            shape: a tuple of integers representing the shape of the grid.
             spacing: a float representing the axis-wise spacing of the grid, in
                 meters.
+            extra: an integer or an iterable of integers representing the number
+                of gridpoints to add as boundary layers on each side of the grid.
+                Extras are added both before and after the grid on each axis.
+                Default is 50.
+            absorbing: an integer or an iterable of integers representing the
+                number of gridpoints within the boundary layers that are
+                absorbing. Default is 40.
 
         Returns:
             Grid: the Grid object.
         """
-        space = stride.Space(shape=shape, spacing=spacing)
+        n_dims = len(shape)
+
+        if isinstance(extra, int):
+            extra = (extra,) * n_dims
+        else:
+            extra = tuple(extra)
+        if isinstance(absorbing, int):
+            absorbing = (absorbing,) * n_dims
+        else:
+            absorbing = tuple(absorbing)
+
+        space = stride.Space(
+            shape=shape, extra=extra, absorbing=absorbing, spacing=spacing
+        )
         return Grid(space=space, time=None)
 
     @staticmethod
@@ -84,18 +109,7 @@ def make_grid(
             Grid: the Grid object.
         """
         _extent = np.array(extent, dtype=float)
-        n_dims = len(_extent)
         dx = speed_water / center_frequency / ppw  # m
         shape = _compute_grid_shape(_extent, dx)
 
-        if isinstance(extra, int):
-            extra = (extra,) * n_dims
-        else:
-            extra = tuple(extra)
-        if isinstance(absorbing, int):
-            absorbing = (absorbing,) * n_dims
-        else:
-            absorbing = tuple(absorbing)
-
-        space = stride.Space(shape=shape, extra=extra, absorbing=absorbing, spacing=dx)
-        return Grid(space=space, time=None)
+        return Grid.make_shaped_grid(shape, dx, extra, absorbing)
diff --git a/src/neurotechdevkit/results/_metrics.py b/src/neurotechdevkit/results/_metrics.py
index bbdb16cf..b7f2ac38 100644
--- a/src/neurotechdevkit/results/_metrics.py
+++ b/src/neurotechdevkit/results/_metrics.py
@@ -1,14 +1,17 @@
 from __future__ import annotations
 
+from typing import List, Tuple
+
 import numpy as np
 import numpy.typing as npt
+import scipy.ndimage
 
 from . import _results as results
 
 
 def calculate_all_metrics(
     result: results.SteadyStateResult,
-) -> dict[str, dict[str, float | str]]:
+) -> dict[str, dict[str, float | int | str]]:
     """Calculate all metrics for the steady-state result and return as a dictionary.
 
     The keys for the dictionary are the names of the metrics. The value for each metric
@@ -25,6 +28,22 @@ def calculate_all_metrics(
         The dictionary structure containing metrics and their descriptions.
     """
     return {
+        "focal_pressure": {
+            "value": calculate_focal_pressure(result, layer="brain"),
+            "unit-of-measurement": "Pa",
+            "description": ("The peak pressure amplitude within the brain."),
+        },
+        "focal_volume": {
+            "value": calculate_focal_volume(result, layer="brain"),
+            "unit-of-measurement": "voxels",
+            "description": (
+                "The region of the brain where the pressure amplitude is above"
+                " 50%% of the maximum pressure amplitude. If more than one"
+                " region is above 50%% of the maximum pressure amplitude, the"
+                " volume of largest connected region is returned. Also called"
+                " the -6 dB focal volume."
+            ),
+        },
         "focal_gain": {
             "value": calculate_focal_gain(result),
             "unit-of-measurement": "dB",
@@ -35,6 +54,19 @@ def calculate_all_metrics(
                 " amplitude inside the brain but excluding the target."
             ),
         },
+        **{
+            f"FWHM_{axis}": {
+                "value": calculate_focal_fwhm(result, axis=axis, layer="brain"),
+                "unit-of-measurement": "grid steps",
+                "description": (
+                    "The full-width at half-maximum-pressure (FWHM) along the"
+                    f" {axis}-axis. The FWHM is calculated by taking a 1-D"
+                    " profile through the focal point. Also called the -6 dB"
+                    " focal width."
+                ),
+            }
+            for axis in ("x", "y", "z")[: result.get_steady_state().ndim]
+        },
         "I_ta_target": {
             "value": Conversions.convert(
                 from_uom="W/m²", to_uom="mW/cm²", value=calculate_i_ta_target(result)
@@ -118,6 +150,50 @@ def calculate_all_metrics(
     }
 
 
+def calculate_focal_pressure(
+    result: results.SteadyStateResult, layer: str | None = None
+) -> float:
+    """Calculate the focal pressure of the simulation result for a particular layer.
+
+    Focal pressure is also known as peak pressure.
+
+    Args:
+        result: the Result object containing the simulation results.
+        layer: the layer within which to calculate the focal pressure. If None, the
+            default, the focal pressure is calculated over the entire simulation
+            space.
+
+    Returns:
+        The focal pressure (in Pa)
+    """
+    ss_amp_masked = _get_steady_state_in_layer(result, layer=layer)
+    focal_pressure = np.max(ss_amp_masked)
+    return focal_pressure
+
+
+def calculate_focal_position(
+    result: results.SteadyStateResult, layer: str | None = None
+) -> Tuple[np.int_, ...]:
+    """Calculate the focal position of the simulation result for a particular layer.
+
+    The focal position is the position of peak pressure.
+
+    Args:
+        result: the Result object containing the simulation results.
+        layer: the layer within which to calculate the focal position. If None, the
+            default, the focal position is calculated over the entire simulation
+            space.
+
+    Returns:
+        The focal position (as grid index tuple)
+    """
+    ss_amp_masked = _get_steady_state_in_layer(result, layer=layer)
+    focal_position = np.unravel_index(
+        np.argmax(ss_amp_masked, axis=None), ss_amp_masked.shape
+    )
+    return focal_position
+
+
 def calculate_mechanical_index(
     result: results.SteadyStateResult, layer: str | None = None
 ) -> float:
@@ -135,14 +211,8 @@ def calculate_mechanical_index(
     Returns:
         The mechanical index (in Pa √s̅)
     """
-    if layer is None:
-        mask = np.ones_like(result.get_steady_state(), dtype=bool)
-    else:
-        mask = result.scenario.material_masks[layer]
-    ss_amp_in_brain: npt.NDArray[np.float_] = np.ma.masked_array(
-        result.get_steady_state(), mask=~mask
-    )
-    peak_neg_pressure = np.max(ss_amp_in_brain)
+    ss_amp_masked = _get_steady_state_in_layer(result, layer=layer)
+    peak_neg_pressure = np.max(ss_amp_masked)
     return peak_neg_pressure / np.sqrt(result.center_frequency)
 
 
@@ -289,6 +359,87 @@ def calculate_i_pa_off_target(result: results.SteadyStateResult) -> float:
     return calculate_i_ta_off_target(result)
 
 
+def calculate_focal_volume(
+    result: results.SteadyStateResult,
+    layer: str | None = "brain",
+) -> float:
+    """Calculate the focal volume of the simulation result.
+
+    From https://doi.org/10.1121/10.0013426:
+    > The focal volume is calculated by thresholding the pressure field inside
+    > the brain to 50% of the maximum value and then counting the voxels in the
+    > largest connected component.
+
+    Args:
+        result: the Result object containing the simulation results.
+        layer: the layer within which to calculate the focal volume.
+            If None, the focal volume is calculated over the entire grid.
+
+    Returns:
+        The focal volume (in voxels).
+    """
+    ss_amp_brain_masked = _get_steady_state_in_layer(result, layer=layer)
+    above_threshold = ss_amp_brain_masked >= (0.5 * ss_amp_brain_masked.max())
+    # `scipy.ndimage.label` expects a normal numpy array, so let's fill the masked array
+    if isinstance(above_threshold, np.ma.MaskedArray):
+        above_threshold = above_threshold.filled(False)
+    # Get contiguous regions
+    labeled_mask, _ = scipy.ndimage.label(above_threshold)
+    # Count the occurrences of each label in the labeled array
+    unique_labels, label_counts = np.unique(labeled_mask, return_counts=True)
+    # Exclude background, which `scipy.ndimage.label` labels as 0
+    label_counts = label_counts[unique_labels != 0]
+    assert label_counts.sum(), "Expected >=1 connected regions above 50% threshold"
+    # Get size of largest connected component
+    return label_counts.max()
+
+
+def calculate_focal_fwhm(
+    result: results.SteadyStateResult,
+    *,
+    axis: str | int,
+    layer: str | None = "brain",
+) -> float:
+    """Calculate the full-width at half-maximum-pressure (FWHM) along an axis.
+
+    Uses the 1-D profile through the focal position.
+
+    Args:
+        result: the Result object containing the simulation results.
+        axis: the axis along which to calculate the FWHM. Can be specified as a string
+            ("x", "y", or "z") or as an integer (0, 1, or 2).
+        layer: the layer within which to calculate the FWHM
+            If None, the FWHM is calculated over the entire grid.
+
+    Returns:
+        The full-width at half-maximum (in integer grid steps).
+    """
+    if isinstance(axis, str):
+        # Convert axis to index
+        axis = ["x", "y", "z"].index(axis.lower())
+
+    # Extract 1-D slice going through the focal position along `axis`
+    focal_position_idx = calculate_focal_position(result, layer=layer)
+    ss_amp_masked = _get_steady_state_in_layer(result, layer=layer)
+    slicer: List[np.int_ | slice] = list(focal_position_idx)
+    slicer[axis] = slice(None)
+    focal_slice = ss_amp_masked[tuple(slicer)]
+
+    # Find which indices are above the half-maximum-pressure
+    focal_position_axis_idx = focal_position_idx[axis]
+    above_half_max = focal_slice >= (0.5 * focal_slice[focal_position_axis_idx])
+
+    # `scipy.ndimage.label` expects a normal numpy array, so let's fill the masked array
+    if isinstance(above_half_max, np.ma.MaskedArray):
+        above_half_max = above_half_max.filled(False)
+    labeled_mask, _ = scipy.ndimage.label(above_half_max)
+    # We only care about the connected region containing the focal position
+    fwhm_group_label = labeled_mask[focal_position_axis_idx]
+    fwhm = np.sum(labeled_mask == fwhm_group_label)
+
+    return fwhm
+
+
 class Conversions:
     """A set of unit-of-measurement conversion tools."""
 
@@ -327,3 +478,26 @@ def convert(from_uom, to_uom, value):
             return value * 1e-4
 
         raise ValueError(f"Unsupported conversion: from '{from_uom}' to '{to_uom}'")
+
+
+def _get_steady_state_in_layer(
+    result: results.SteadyStateResult,
+    layer: str | None = None,
+) -> npt.NDArray | np.ma.MaskedArray:
+    """Get the steady-state pressure amplitude within a particular layer.
+
+    Args:
+        result: the Result object containing the simulation results.
+        layer: the layer to extract
+
+    Returns:
+        The steady-state pressure amplitude, with non-layer values masked.
+    """
+    steady_state = result.get_steady_state()
+    if layer is None:
+        return steady_state
+    else:
+        return np.ma.masked_array(
+            steady_state,
+            mask=~result.scenario.material_masks[layer],
+        )
diff --git a/src/neurotechdevkit/scenarios/built_in/_scenario_2.py b/src/neurotechdevkit/scenarios/built_in/_scenario_2.py
index e880fa56..8f989d2a 100644
--- a/src/neurotechdevkit/scenarios/built_in/_scenario_2.py
+++ b/src/neurotechdevkit/scenarios/built_in/_scenario_2.py
@@ -365,7 +365,7 @@ def _create_scenario_2_mask(
         brain_mask = brain_mask[:, :, z_center_idx]
 
     if skull_mask.shape != output_shape:
-        # resample the mask to match the grid
+        # re-sample the mask to match the grid
         scale_factor = np.array(output_shape) / np.array(skull_mask.shape)
         skull_mask = scipy.ndimage.zoom(skull_mask, scale_factor, order=0)
         brain_mask = scipy.ndimage.zoom(brain_mask, scale_factor, order=0)
diff --git a/src/web/app.py b/src/web/app.py
index cd85c082..daf4ef8e 100644
--- a/src/web/app.py
+++ b/src/web/app.py
@@ -1,9 +1,14 @@
 """Flask app factory."""
+import os
+import pathlib
+
 import matplotlib
 from flask import Flask
 from matplotlib import rc
 from web.simulation_runner import SimulationRunner
 
+CT_FOLDER = pathlib.Path("./ct-scans/")
+
 
 def create_app():
     """Create the Flask app."""
@@ -23,11 +28,19 @@ def create_app():
     return app
 
 
+def assign_CT_FOLDER(app):
+    """Assign the CT_FOLDER to the app."""
+    if not os.path.exists(CT_FOLDER):
+        os.makedirs(CT_FOLDER)
+    app.config["CT_FOLDER"] = CT_FOLDER
+
+
 app = create_app()
 
 
 def run():
     """Run the web app."""
+    assign_CT_FOLDER(app)
     app.run(host="0.0.0.0", port=5000, debug=True)
 
 
diff --git a/src/web/computed_tomography.py b/src/web/computed_tomography.py
new file mode 100644
index 00000000..d7ef2208
--- /dev/null
+++ b/src/web/computed_tomography.py
@@ -0,0 +1,224 @@
+"""Contains functions for loading Computed Tomography images and their masks."""
+import json
+import pathlib
+import tempfile
+import zipfile
+from dataclasses import dataclass
+from typing import Dict, List, Optional, Tuple, Union
+
+import nibabel as nib
+import numpy as np
+import pydicom
+from nibabel.nifti1 import Nifti1Image
+from nibabel.nifti2 import Nifti2Image
+from web.messages.settings import Axis
+
+from neurotechdevkit.materials import DEFAULT_MATERIALS
+
+
+@dataclass
+class ComputedTomographyImage:
+    """A Computed Tomography image with its voxel spacing and material masks."""
+
+    data: np.ndarray
+    spacing: Tuple[float, float, float]  # voxel sizes in millimeter
+    material_masks: Dict[str, np.ndarray]
+
+    @property
+    def spacing_in_meters(self) -> Tuple[float, ...]:
+        """The voxel spacing in meters."""
+        return tuple(s / 1000 for s in self.spacing)
+
+
+@dataclass
+class ComputedTomographyInfo:
+    """Information about a Computed Tomography."""
+
+    filename: str
+    shape: Tuple[int, ...]
+    spacing: str  # We need to convert the tuple to a string to be able to serialize it.
+
+
+def _assign_masks(mapping: dict[str, int], data: np.ndarray) -> Dict[str, np.ndarray]:
+    """Assign the masks to the materials.
+
+    Args:
+        mapping: A dictionary mapping material names to label values.
+        data: The Computed Tomography image data.
+
+    Returns:
+        A dictionary of material masks with material names as keys.
+    """
+    mask_materials: Dict[str, np.ndarray] = {}
+    for material, label_value in mapping.items():
+        label_mask = data == label_value
+        mask_materials[material] = label_mask.astype(np.bool_)
+
+    return mask_materials
+
+
+def _load_nii(filepath: pathlib.Path) -> Tuple[np.ndarray, Tuple[float, float, float]]:
+    """Load a NIfTI file.
+
+    Args:
+        filepath: Path to the NIfTI file.
+
+    Returns:
+        A tuple containing the image data and the voxel spacing.
+    """
+    image: Union[Nifti1Image, Nifti2Image] = nib.load(filepath)  # type: ignore
+    return image.get_fdata(), image.header.get_zooms()
+
+
+def _load_compacted_dicom(
+    filepath: pathlib.Path,
+) -> Tuple[np.ndarray, Tuple[float, float, float]]:
+    """Extract and load the dicom files from a zip file.
+
+    Args:
+        filepath: Path to the zip file.
+
+    Returns:
+        A tuple containing the image data and the voxel spacing.
+    """
+    with tempfile.TemporaryDirectory() as tmpdir:
+        with zipfile.ZipFile(filepath, "r") as zip_ref:
+            zip_ref.extractall(tmpdir)
+        dicom_files = list(pathlib.Path(tmpdir).glob("*.dcm"))
+        dicom_files.sort()
+        dicom_data = []
+        for dicom_file in dicom_files:
+            data = pydicom.dcmread(dicom_file, force=True)
+            dicom_data.append(data)
+
+        image = np.stack([dicom.pixel_array for dicom in dicom_data])
+        spacing = (
+            float(dicom_data[0].PixelSpacing[0]),
+            float(dicom_data[0].PixelSpacing[1]),
+            float(dicom_data[0].SliceThickness),
+        )
+        return image, spacing
+
+
+def _load_ct_file(
+    ct_path: pathlib.Path,
+) -> Tuple[np.ndarray, Tuple[float, float, float]]:
+    """Load a CT image file.
+
+    Args:
+        ct_path: Path to the CT image file.
+
+    Returns:
+        A tuple containing the image data and the voxel spacing.
+    """
+    if ct_path.suffix == ".nii":
+        return _load_nii(ct_path)
+    elif ct_path.suffix == ".zip":
+        return _load_compacted_dicom(ct_path)
+    else:
+        raise ValueError(f"Unsupported CT file type: {ct_path.suffix}")
+
+
+def _load_mapping_file(mapping_path: pathlib.Path) -> dict[str, int]:
+    """Load the mapping file.
+
+    Args:
+        mapping_path: Path to the mapping file.
+
+    Returns:
+        A dictionary mapping material names to label values.
+    """
+    with open(mapping_path, "r") as f:
+        try:
+            return json.load(f)
+        except json.JSONDecodeError:
+            raise ValueError(
+                f"Could not load masks file {mapping_path.name} as JSON."
+            ) from None
+
+
+def get_ct_image(
+    ct_path: pathlib.Path,
+    slice_axis: Optional[Axis],
+    slice_position: Optional[float],
+) -> ComputedTomographyImage:
+    """Load a CT image and its masks.
+
+    Args:
+        ct_path: Path to the CT image file.
+        slice_axis: The axis along which to slice the image.
+        slice_position: The position along the slice axis.
+
+    Returns:
+        A ComputedTomographyImage object.
+    """
+    array, spacing = _load_ct_file(ct_path)
+    mapping = _load_mapping_file(ct_path.with_suffix(".json"))
+
+    if slice_axis is not None and slice_position is not None:
+        if slice_axis == Axis.x:
+            position = int(slice_position / spacing[0] * 1000)
+            array = array[position, :, :]
+        elif slice_axis == Axis.y:
+            position = int(slice_position / spacing[1] * 1000)
+            array = array[:, position, :]
+        elif slice_axis == Axis.z:
+            position = int(slice_position / spacing[2] * 1000)
+            array = array[:, :, position]
+
+    material_masks = _assign_masks(mapping, array)
+    return ComputedTomographyImage(
+        data=array, spacing=spacing, material_masks=material_masks
+    )
+
+
+def get_available_cts(cts_folder: pathlib.Path) -> List[ComputedTomographyInfo]:
+    """Get the list of available Computed Tomography's in the given directory.
+
+    Args:
+        cts_folder: Path to the folder containing the Computed Tomography's.
+
+    Returns:
+        The list of available Computed Tomography's.
+    """
+    files: List[ComputedTomographyInfo] = []
+    for file in cts_folder.glob("**/*"):
+        if file.is_file() and file.suffix in [".nii", ".zip"]:
+            ct, spacing = _load_ct_file(file)
+            ct_info = ComputedTomographyInfo(
+                filename=file.name, shape=ct.shape, spacing=str(list(spacing))
+            )
+            files.append(ct_info)
+    return files
+
+
+def validate_ct(directory: pathlib.Path, files: list[str]) -> ComputedTomographyInfo:
+    """Validate the CT image and its masks.
+
+    Args:
+        directory: Path to the directory containing the CT image and its masks.
+        files: List of files in the directory to be checked.
+
+    Returns:
+        The ComputedTomographyInfo object.
+    """
+    assert (
+        len(files) == 2
+    ), f"Expected exactly two files but found {len(files)}: {files}"
+    ct_file, mapping_file = (
+        (files[0], files[1]) if files[1].endswith(".json") else (files[1], files[0])
+    )
+    data, spacing = _load_ct_file(directory / ct_file)
+    mapping = _load_mapping_file(directory / mapping_file)
+    data_layers = set(np.unique(data).astype(np.int_))
+    assert data_layers == set(mapping.values()), (
+        "The labels in the CT image do not match the labels in the masks file."
+        f"Expected labels: {data_layers}"
+    )
+
+    for material in mapping.keys():
+        if material not in DEFAULT_MATERIALS:
+            raise ValueError(f"Material {material} is not supported.")
+    return ComputedTomographyInfo(
+        filename=ct_file, shape=data.shape, spacing=str(list(spacing))
+    )
diff --git a/src/web/controller.py b/src/web/controller.py
index 650ac965..a3ac9780 100644
--- a/src/web/controller.py
+++ b/src/web/controller.py
@@ -5,12 +5,18 @@
 from typing import Dict, Optional, Tuple, Union
 
 import matplotlib.pyplot as plt
+from flask import current_app
+from web.computed_tomography import get_ct_image
+from web.messages.material_properties import Material, MaterialName, MaterialProperties
 from web.messages.requests import (
     IndexBuiltInScenario,
     RenderLayoutRequest,
     SimulateRequest,
 )
+from web.messages.settings import Axis
 
+from neurotechdevkit.grid import Grid
+from neurotechdevkit.materials import DEFAULT_MATERIALS, get_material
 from neurotechdevkit.results import SteadyStateResult2D, SteadyStateResult3D
 from neurotechdevkit.scenarios import Scenario2D, Scenario3D
 from neurotechdevkit.scenarios.built_in import BUILT_IN_SCENARIOS
@@ -84,6 +90,26 @@ def get_scenario_layout(config: RenderLayoutRequest) -> str:
         config.scenarioSettings.scenario_id,
         config.simulationSettings.centerFrequency,
     )
+    if not config.scenarioSettings.isPreBuilt:
+        assert config.scenarioSettings.ctFile is not None
+        ct_image = get_ct_image(
+            ct_path=current_app.config["CT_FOLDER"] / config.scenarioSettings.ctFile,
+            slice_axis=config.scenarioSettings.ctSliceAxis,
+            slice_position=config.scenarioSettings.ctSlicePosition,
+        )
+        spacing = ct_image.spacing_in_meters
+        if config.is2d:
+            assert config.scenarioSettings.ctSliceAxis is not None
+            if config.scenarioSettings.ctSliceAxis == Axis.x:
+                spacing = (spacing[1], spacing[2])
+            elif config.scenarioSettings.ctSliceAxis == Axis.y:
+                spacing = (spacing[0], spacing[2])
+            else:
+                spacing = spacing = (spacing[0], spacing[1])
+        scenario.grid = Grid.make_shaped_grid(
+            shape=ct_image.data.shape, spacing=spacing
+        )
+        scenario.material_masks = ct_image.material_masks
     _configure_scenario(scenario, config)
     fig = scenario.render_layout(show_sources=len(scenario.sources) > 0)
     buf = io.BytesIO()
@@ -137,16 +163,25 @@ def _instantiate_scenario(
             grid=builtin_scenario.grid,
         )
     else:
-        raise NotImplementedError
+        scenario_class = Scenario2D if is_2d else Scenario3D
+        scenario = scenario_class(
+            sources=[],
+            origin=[0, 0] if is_2d else [0, 0, 0],
+            material_outline_upsample_factor=16,  # TODO: make this configurable
+        )
+
     return scenario
 
 
-async def get_simulation_image(config: SimulateRequest) -> Tuple[str, str]:
+async def get_simulation_image(
+    config: SimulateRequest, app_config: dict
+) -> Tuple[str, str]:
     """
     Simulate a scenario and return the result as a base64 GIF or PNG.
 
     Args:
         config (SimulateRequest): The configuration for the scenario.
+        app_config(dict): The configuration of the app.
 
     Returns:
         Tuple[str, str]: The base64 encoded image and the image format.
@@ -157,6 +192,17 @@ async def get_simulation_image(config: SimulateRequest) -> Tuple[str, str]:
         config.scenarioSettings.scenario_id,
         config.simulationSettings.centerFrequency,
     )
+    if not config.scenarioSettings.isPreBuilt:
+        assert config.scenarioSettings.ctFile is not None
+        ct_image = get_ct_image(
+            ct_path=app_config["CT_FOLDER"] / config.scenarioSettings.ctFile,
+            slice_axis=config.scenarioSettings.ctSliceAxis,
+            slice_position=config.scenarioSettings.ctSlicePosition,
+        )
+        scenario.grid = Grid.make_shaped_grid(
+            shape=ct_image.data.shape, spacing=ct_image.spacing[0] / 1000
+        )
+        scenario.material_masks = ct_image.material_masks
     _configure_scenario(scenario, config)
     scenario.compile_problem()
 
@@ -182,11 +228,36 @@ async def get_simulation_image(config: SimulateRequest) -> Tuple[str, str]:
             return data, "gif"
 
 
+def get_default_material_properties(center_frequency: float) -> MaterialProperties:
+    """Get the default material properties with center frequency.
+
+    Args:
+        center_frequency: The center frequency of the scenario.
+
+    Returns:
+        The default material properties.
+    """
+    material_properties = {}
+    for ndk_material_name in DEFAULT_MATERIALS:
+        ndk_material = get_material(ndk_material_name, center_frequency)
+        material = Material.from_ndk_material(ndk_material)
+        material_name = MaterialName.get_material_name(ndk_material_name)
+        material_properties[material_name] = material
+
+    return MaterialProperties(**material_properties)
+
+
 def _configure_scenario(
     scenario: Union[Scenario2D, Scenario3D],
     config: Union[RenderLayoutRequest, SimulateRequest],
 ):
-    """Configure a scenario based on the given configuration."""
+    """Configure a scenario based on the given configuration.
+
+    Args:
+        scenario: The scenario to configure.
+        config: The configuration.
+    """
+    scenario.center_frequency = config.simulationSettings.centerFrequency
     if config_target := config.target:
         scenario.target = config_target.to_ndk_target()
 
diff --git a/src/web/messages/settings.py b/src/web/messages/settings.py
index fe60a0f6..447bea4c 100644
--- a/src/web/messages/settings.py
+++ b/src/web/messages/settings.py
@@ -11,7 +11,7 @@
 
 
 class Axis(str, Enum):
-    """Axis enum for the axis of CT scans."""
+    """Axis enum for the axis of Computed Tomography scans."""
 
     x = "x"
     y = "y"
@@ -82,6 +82,9 @@ class ScenarioSettings(BaseModel):
 
     isPreBuilt: bool
     scenarioId: Optional[str]
+    ctFile: Optional[str]
+    ctSliceAxis: Optional[Axis]
+    ctSlicePosition: Optional[float]
 
     @validator("scenarioId", always=True)
     def _validate(cls, value, values):
@@ -90,9 +93,8 @@ def _validate(cls, value, values):
         return value
 
     @property
-    def scenario_id(self) -> str:
+    def scenario_id(self) -> Optional[str]:
         """Get the scenario id."""
-        assert self.scenarioId
         return self.scenarioId
 
 
diff --git a/src/web/templates/components/parameters.html b/src/web/templates/components/parameters.html
index 12a46f7b..63cc3838 100644
--- a/src/web/templates/components/parameters.html
+++ b/src/web/templates/components/parameters.html
@@ -74,23 +74,42 @@
     });
   }
 
+  function makeFieldsRequired(querySelector, required){
+    document.querySelectorAll(querySelector).forEach((element) => {
+      if (required){
+        element.setAttribute("required", true)
+      } else {
+        element.removeAttribute("required")
+      }
+    })
+  }
+  function hideFields(querySelector, shouldHide){
+    document.querySelectorAll(querySelector).forEach((element) => {
+      if (shouldHide){
+        element.setAttribute("hidden", true)
+      } else {
+        element.removeAttribute("hidden")
+      }
+    })
+  }
+
   function changeSimulationType(is2d) {
     if (is2d) {
-      document.querySelectorAll('[data-3d-only="true"]').forEach((element) => {
-        element.setAttribute("hidden", true)
-      })
-      document.querySelectorAll('[data-3d-input-only="true"]').forEach((element) => {
-        element.removeAttribute("required")
-      })
+      makeFieldsRequired('[data-3d-input-only="true"]', false)
+      makeFieldsRequired('[data-2d-input-only="true"]', true)
+      hideFields('[data-3d-only="true"]', true)
+      hideFields('[data-2d-only="true"]', false)
     } else {
-      document.querySelectorAll('[data-3d-only="true"]').forEach((element) => {
-        element.removeAttribute("hidden")
-      })
-      document.querySelectorAll('[data-3d-input-only="true"]').forEach((element) => {
-        element.setAttribute("required", true)
-      })
+      makeFieldsRequired('[data-3d-input-only="true"]', true)
+      makeFieldsRequired('[data-2d-input-only="true"]', false)
+      hideFields('[data-2d-only="true"]', true)
+      hideFields('[data-3d-only="true"]', false)
     }
+    resetAll(is2d)
     collapseParameters()
+  }
+
+  function resetAll(is2d){
     clearPlotArea()
     disallowRenderLayout()
     disallowRunSimulation()
@@ -124,6 +143,14 @@
     const isRunningSimulation = {{is_running_simulation | tojson}}
     const currentConfiguration = {{configuration | tojson}}
 
+    const defaultMaterialProperties = {{ all_material_properties | tojson }}
+    setMaterialProperties(defaultMaterialProperties)
+
+    const availableCTs = {{ available_cts | tojson }}
+    const selectedCT = null
+    loadCTs(availableCTs, selectedCT)
+
+
     setScenariosOptions((is2d = true)) // default is 2D
     if (currentConfiguration){
       setScenariosOptions(currentConfiguration.is2d)
@@ -132,8 +159,7 @@
       } else {
         document.getElementById("is3d").checked = true
       }
-      const scenarioSelect = document.getElementById("scenario")
-      scenarioSelect.value = currentConfiguration.scenarioSettings.scenarioId
+      
       setScenarioValues(currentConfiguration)
     }
 
@@ -150,7 +176,6 @@
     if(isRunningSimulation){
       showCancelSimulation()
     }
-
   }
 </script>
 <h1>
@@ -166,8 +191,7 @@ <h1>
         id="is2d"
         autocomplete="off"
         onchange="changeSimulationType(is2d=true)"
-        checked
-      />
+        checked />
       <label class="btn btn-primary" for="is2d">2D</label>
 
       <input
@@ -176,12 +200,11 @@ <h1>
         name="options"
         id="is3d"
         autocomplete="off"
-        onchange="changeSimulationType(is2d=false)"
-      />
+        onchange="changeSimulationType(is2d=false)" />
       <label class="btn btn-primary" for="is3d">3D</label>
     </div>
   </fieldset>
-  <form id="parameters-form">
+  <form id="parametersForm">
     <div class="accordion accordion-flush" id="parameters-accordion">
       <div class="accordion-item">
         <h2 class="accordion-header">
@@ -191,16 +214,14 @@ <h2 class="accordion-header">
             data-bs-toggle="collapse"
             data-bs-target="#scenario-accordion-item"
             aria-expanded="false"
-            aria-controls="scenario-accordion-item"
-          >
+            aria-controls="scenario-accordion-item">
             Scenario
           </button>
         </h2>
         <div
           id="scenario-accordion-item"
           class="accordion-collapse collapse show"
-          data-bs-parent="#parameters-accordion"
-        >
+          data-bs-parent="#parameters-accordion">
           <div class="accordion-body">
             {% include 'components/parameters/scenario.html' %}
           </div>
@@ -215,16 +236,14 @@ <h2 class="accordion-header">
             data-bs-target="#display-accordion-item"
             aria-expanded="false"
             aria-controls="display-accordion-item"
-            disabled
-          >
+            disabled>
             Display
           </button>
         </h2>
         <div
           id="display-accordion-item"
           class="accordion-collapse collapse"
-          data-bs-parent="#parameters-accordion"
-        >
+          data-bs-parent="#parameters-accordion">
           <div class="accordion-body">
             {% include 'components/parameters/display.html' %}
           </div>
@@ -240,16 +259,14 @@ <h2 class="accordion-header">
             data-bs-target="#transducers-accordion-item"
             aria-expanded="false"
             aria-controls="transducers-accordion-item"
-            disabled
-          >
+            disabled>
             Transducers
           </button>
         </h2>
         <div
           id="transducers-accordion-item"
           class="accordion-collapse collapse"
-          data-bs-parent="#parameters-accordion"
-        >
+          data-bs-parent="#parameters-accordion">
           <div class="accordion-body">
             {% include 'components/parameters/transducers.html' %}
           </div>
@@ -264,16 +281,14 @@ <h2 class="accordion-header">
             data-bs-target="#target-accordion-item"
             aria-expanded="false"
             aria-controls="target-accordion-item"
-            disabled
-          >
+            disabled>
             Target
           </button>
         </h2>
         <div
           id="target-accordion-item"
           class="accordion-collapse collapse"
-          data-bs-parent="#parameters-accordion"
-        >
+          data-bs-parent="#parameters-accordion">
           <div class="accordion-body">
             {% include 'components/parameters/target.html' %}
           </div>
@@ -288,16 +303,14 @@ <h2 class="accordion-header">
             data-bs-target="#simulation-settings-accordion-item"
             aria-expanded="false"
             aria-controls="simulation-settings-accordion-item"
-            disabled
-          >
+            disabled>
             Simulation Settings
           </button>
         </h2>
         <div
           id="simulation-settings-accordion-item"
           class="accordion-collapse collapse"
-          data-bs-parent="#parameters-accordion"
-        >
+          data-bs-parent="#parameters-accordion">
           <div class="accordion-body">
             {% include 'components/parameters/simulation_settings.html' %}
           </div>
diff --git a/src/web/templates/components/parameters/materials.html b/src/web/templates/components/parameters/materials.html
index 3b821222..8a9f2413 100644
--- a/src/web/templates/components/parameters/materials.html
+++ b/src/web/templates/components/parameters/materials.html
@@ -1,5 +1,5 @@
 <script>
-  function setMaterialPropertiesPreBuiltValues(materialProperties) {
+  function setMaterialProperties(materialProperties) {
     for (let key in materialProperties) {
       if (!materialProperties[key]){
         // hiding the material as it is not part of the scenario
diff --git a/src/web/templates/components/parameters/scenario.html b/src/web/templates/components/parameters/scenario.html
index 5d19e599..2449e846 100644
--- a/src/web/templates/components/parameters/scenario.html
+++ b/src/web/templates/components/parameters/scenario.html
@@ -1,28 +1,34 @@
 <script>
-  let scenarioInfos = {{built_in_scenarios | tojson}};
+  let scenarioInfos = {{ built_in_scenarios | tojson}};
+  let ctDimensionPerAxis = {};
+
   function togglePreBuiltDiv() {
     document.getElementById("preBuiltDiv").hidden = false
     document.getElementById("ctFileDiv").hidden = true
+    resetAll(document.getElementById("is2d").checked)
   }
 
   function toggleCTFileDiv() {
     document.getElementById("preBuiltDiv").hidden = true
     document.getElementById("ctFileDiv").hidden = false
+    resetAll(document.getElementById("is2d").checked)
   }
 
   function valueChanged() {
     let hasError = false
-    if (!validateTargetValues()){
+    if (!validateCTValues()) {
+      hasError = true
+    } else if (!validateTargetValues()) {
       hasError = true
-    } else if (!validateDisplayValues()){
+    } else if (!validateDisplayValues()) {
       hasError = true
-    } else if(!validateSimulationSettingsValues()) {
+    } else if (!validateSimulationSettingsValues()) {
       hasError = true
-    } else if (!validateMaterialValues()){
+    } else if (!validateMaterialValues()) {
       hasError = true
     }
 
-    if (hasError){
+    if (hasError) {
       disallowRunSimulation()
       clearPlotArea()
     } else {
@@ -30,15 +36,25 @@
     }
   }
 
-  function scenarioSelected(scenarioInput){
+  function validateCTValues() {
+    let form = document.getElementById("ctForm")
+    if (!form.checkValidity()) {
+      $("#scenario-accordion-item").collapse("show")
+      form.reportValidity()
+      return false
+    }
+    return true
+  }
+
+  function scenarioSelected(scenarioInput) {
     clearPlotArea()
 
     const selection = scenarioInput.value;
-    if (selection===""){
+    if (selection === "") {
       disallowRunSimulation()
       disallowRenderLayout()
       clearAllTransducers()
-      document.getElementById("parameters-form").reset()
+      document.getElementById("parametersForm").reset()
       return;
     }
     const selectedScenario = scenarioInfos[selection]
@@ -47,32 +63,58 @@
     allowViewScenarioParameters()
   }
 
-  function setScenarioValues(scenario){
-    setTargetValues(scenario.target)
-    setSimulationSettingsValues(scenario.simulationSettings)
-    setTransducersValues(scenario.transducers)
-    if (!scenario.is2d){
-      setDisplayValues(scenario.displaySettings)
+  function setScenarioValues(configuration) {
+    if (configuration.scenarioSettings.scenarioId){
+        // running a built-in scenario
+        const scenarioSelect = document.getElementById("scenario")
+        scenarioSelect.value = configuration.scenarioSettings.scenarioId
+    } else if (configuration.scenarioSettings.ctFile){
+        // running a ct scan
+        document.getElementById("isCtScan").checked = true
+        toggleCTFileDiv()
+        const loadedCTsElement = document.getElementById("loadedCTs")
+        loadedCTsElement.value = configuration.scenarioSettings.ctFile
+        if (configuration.is2d){
+          const ctSliceAxis = document.getElementById("ctSliceAxis")
+          ctSliceAxis.value = configuration.scenarioSettings.ctSliceAxis
+          const ctSlicePosition = document.getElementById("ctSlicePosition")
+          ctSlicePosition.value = configuration.scenarioSettings.ctSlicePosition
+        }
+    }
+    setTargetValues(configuration.target)
+    setSimulationSettingsValues(configuration.simulationSettings)
+    setTransducersValues(configuration.transducers)
+    if (!configuration.is2d) {
+      setDisplayValues(configuration.displaySettings)
     }
   }
 
   function getScenarioParams() {
     const isPreBuilt = document.getElementById("isPreBuilt").checked
+    const is2d = document.getElementById("is2d").checked
     let params = {
       isPreBuilt: isPreBuilt,
     }
     if (isPreBuilt) {
       params["scenarioId"] = document.getElementById("scenario").value
     } else {
-      params["ctFile"] = document.getElementById("ctFile").value
+      // ct file was loaded
+      params["ctFile"] = document.getElementById("loadedCTs").value
+      if (is2d){
+        // these attributes only exist for 2D
+        params["ctSliceAxis"] = document.getElementById("ctSliceAxis").value
+        params["ctSlicePosition"] = document.getElementById("ctSlicePosition").value
+      }
     }
     return params
   }
 
-  function setScenariosOptions(is2d){
+  function setScenariosOptions(is2d) {
     const scenarioSelect = document.getElementById("scenario")
     scenarioSelect.innerHTML = ""
 
+    document.getElementById("ctForm").reset()
+
     const option = document.createElement("option")
     option.value = ""
     option.text = "Select a scenario"
@@ -87,63 +129,214 @@
       }
     }
   }
+
+  function filesChosen(fileInput) {
+    clearPlotArea()
+    if (fileInput.files.length == 0){
+      document.getElementById("fileUploadButton").disabled = true
+      return
+    }
+    if (fileInput.files.length !== 2) {
+      document.getElementById("plotArea").innerHTML = "Please select two files"
+      document.getElementById("fileUploadButton").disabled = true
+      return
+    }
+    const firstFileName = fileInput.files[0].name
+    const secondFileName = fileInput.files[1].name
+    if (firstFileName.split('.')[0] !== secondFileName.split('.')[0]) {
+      document.getElementById("plotArea").innerHTML = 'Files must have the same name'
+      document.getElementById("fileUploadButton").disabled = true
+      return
+    }
+    document.getElementById("fileUploadButton").disabled = false
+  }
+
+  function fillDefaultValues() {
+    setSimulationSettingsValues({
+      "simulationPrecision": 6,
+      "centerFrequency": 5e5,
+      "isSteadySimulation": true,
+    })
+    const is3d = document.getElementById("is3d").checked
+
+    setTargetValues({
+      "centerX": ctDimensionPerAxis["x"] / 2,
+      "centerY": ctDimensionPerAxis["y"] / 2,
+      "centerZ": is3d ? ctDimensionPerAxis["z"] / 2 : '',
+      "radius": 0.01,
+    })
+  }
+
+  function fillCTDimensions(shape, spacing){
+    // Calculate the max value for each axis
+    ctDimensionPerAxis["x"] = (shape[0] * spacing[0]) / 1000
+    ctDimensionPerAxis["y"] = (shape[1] * spacing[1]) / 1000
+    ctDimensionPerAxis["z"] = (shape[2] * spacing[2]) / 1000
+  }
+
+  function ctSliceAxisChanged(input){
+    let ctSliceAxisValue = input.value
+    let ctSlicePositionInput = document.getElementById("ctSlicePosition")
+    ctSlicePositionInput.max = ctDimensionPerAxis[ctSliceAxisValue]
+    valueChanged()
+  }
+
+  function ctSelected(ctInput) {
+    const selectedCT = ctInput.value
+    if (selectedCT === '') {
+      return
+    }
+    document.getElementById("ctSlicePosition").value = ""
+    document.getElementById("ctSlicePosition").removeAttribute("max")
+    document.getElementById("ctSliceAxis").value = ""
+    const shape = JSON.parse(ctInput.selectedOptions[0].dataset.shape)
+    const spacing = JSON.parse(ctInput.selectedOptions[0].dataset.spacing)
+    fillCTDimensions(shape, spacing)
+    fillDefaultValues()
+    allowViewScenarioParameters()
+    valueChanged()
+  }
+
+  function fileUploadClicked(button) {
+    BtnLoading(button)
+    var input = document.getElementById("ctFiles")
+
+    var data = new FormData()
+    data.append('file_0', input.files[0], input.files[0].name);
+    data.append('file_1', input.files[1], input.files[1].name);
+
+    const url = "{{ url_for('main.ct_scan') }}"
+
+    fetch(url, {
+      method: 'POST',
+      body: data
+    }).then(response => {
+        if (!response.ok) {
+            return Promise.reject(response);
+        }
+        return response.json();
+    })
+      .then(message => {
+        const selectedCT = message.selected_ct
+        const availableCTs = message.available_cts
+        loadCTs(availableCTs, selectedCT.filename)
+        const shape = selectedCT.shape
+        const spacing = JSON.parse(selectedCT.spacing)
+
+        document.getElementById("ctFiles").value = ""
+        fillCTDimensions(shape, spacing)
+        fillDefaultValues()
+        allowViewScenarioParameters()
+        valueChanged()
+        BtnReset(button)
+        document.getElementById("fileUploadButton").disabled = true
+      })
+      .catch(response => {
+        response.json().then((error) => {
+          clearPlotArea()
+          document.getElementById("plotArea").innerHTML = error.error
+        })
+        BtnReset(button)
+        document.getElementById("fileUploadButton").disabled = true
+      });
+  }
+
+  function loadCTs(availableCTs, selectedCT) {
+    const loadedCTsElement = document.getElementById("loadedCTs")
+    loadedCTsElement.innerHTML = ""
+    for (const ct of availableCTs) {
+      if (ct.filename === selectedCT) {
+        loadedCTsElement.innerHTML += `<option value="${ct.filename}" data-shape="[${ct.shape}]" data-spacing="${ct.spacing}" selected>${ct.filename}</option>`
+      } else {
+        loadedCTsElement.innerHTML += `<option value="${ct.filename}" data-shape="[${ct.shape}]" data-spacing="${ct.spacing}">${ct.filename}</option>`
+      }
+    }
+    if (loadedCTs.length === 0) {
+      loadedCTsElement.innerHTML += `<option value="">No CTs loaded</option>`
+      loadedCTsElement.setAttribute("disabled", "disabled")
+      document.getElementById("ctSliceAxis").setAttribute("disabled", "disabled")
+      document.getElementById("ctSlicePosition").setAttribute("disabled", "disabled")
+    } else {
+      loadedCTsElement.removeAttribute("disabled")
+      document.getElementById("ctSliceAxis").removeAttribute("disabled")
+      document.getElementById("ctSlicePosition").removeAttribute("disabled")
+    }
+  }
+
 </script>
-<form>
-  <fieldset>
-    <div class="mb-3">
-      <div class="btn-group">
-        <input
-          type="radio"
-          class="btn-check"
-          name="options"
-          id="isPreBuilt"
-          autocomplete="off"
-          onchange="togglePreBuiltDiv()"
-          checked
-        />
-        <label
-          class="btn btn-primary"
-          for="isPreBuilt"
-          data-bs-toggle="tooltip"
-          data-bs-placement="right"
-          title="Use a prebuilt scenario with a target and transducers"
-          >Prebuilt</label
-        >
-
-        <input
-          type="radio"
-          class="btn-check"
-          name="options"
-          id="isCtScan"
-          autocomplete="off"
-          onchange="toggleCTFileDiv()"
-          disabled
-        />
-        <label
-          class="btn btn-primary"
-          for="isCtScan"
-          data-bs-toggle="tooltip"
-          data-bs-placement="right"
-          title="Load a CT scan from a file"
-          >CT Scan</label
-        >
-      </div>
+<form></form>
+<fieldset>
+  <div class="mb-3">
+    <div class="btn-group">
+      <input type="radio" class="btn-check" name="options" id="isPreBuilt"
+        autocomplete="off"
+        onchange="togglePreBuiltDiv()" checked />
+      <label class="btn btn-primary" for="isPreBuilt" data-bs-toggle="tooltip"
+        data-bs-placement="right"
+        title="Use a prebuilt scenario with a target and transducers">Prebuilt</label>
+
+      <input type="radio" class="btn-check" name="options" id="isCtScan"
+        autocomplete="off"
+        onchange="toggleCTFileDiv()" />
+      <label class="btn btn-primary" for="isCtScan" data-bs-toggle="tooltip"
+        data-bs-placement="right"
+        title="Load a CT scan from a file">CT Scan</label>
     </div>
+  </div>
 
-    <div class="mb-3">
-      <div id="preBuiltDiv">
-        <label for="scenario" class="form-label">Scenarios:</label>
-        <select
-          name="scenario"
-          id="scenario"
-          class="form-select"
-          onchange="scenarioSelected(this)"
-        ></select>
+  <div class="mb-3">
+    <div id="preBuiltDiv">
+      <label for="scenario" class="form-label">Scenarios:</label>
+      <select name="scenario" id="scenario" class="form-select"
+        onchange="scenarioSelected(this)"></select>
+    </div>
+    <div id="ctFileDiv" hidden>
+      <form id="ctForm">
+        <div class="mb-3">
+          <label class="form-label" for="loadedCTs">Available CTs</label>
+          <select class="form-select" size="3" name="loadedCTs" id="loadedCTs"
+            onchange="ctSelected(this)" required>
+          </select>
+        </div>
+        <div class="mb-3" data-2d-only="true">
+          <label for="ctSliceAxis" class="form-label" data-bs-toggle="tooltip"
+            data-bs-placement="right"
+            title="Axis along which to slice the 3D field to be recorded">Axis</label>
+          <select class="form-select" aria-label="Axis" id="ctSliceAxis"
+            data-2d-input-only="true"
+            name="ctSliceAxis" required onchange="ctSliceAxisChanged(this)">
+            <option value selected>Select an axis</option>
+            <option value="x">X</option>
+            <option value="y">Y</option>
+            <option value="z">Z</option>
+          </select>
+        </div>
+        <div class="mb-3" data-2d-only="true">
+          <label for="ctSlicePosition" class="form-label"
+            data-bs-toggle="tooltip" data-bs-placement="right"
+            title="The position (in meters) along the slice axis at which the slice of the 3D field should be made">Distance
+            from Origin (m)</label>
+          <input type="number" step="any" class="form-control"
+            id="ctSlicePosition" name="ctSlicePosition"
+            data-2d-input-only="true"
+            placeholder="0.0" required
+            onchange="valueChanged()" />
+        </div>
+      </form>
+      <div class="mb-3">
+        <label for="ctFiles" class="form-label"
+          title="Select the CT file and the file containing the mapping between layers and masks">CT
+          and mapping
+          Files</label>
+        <input class="form-control" type="file" id="ctFiles"
+          onchange="filesChosen(this)" multiple />
       </div>
-      <div id="ctFileDiv" hidden>
-        <label for="ctFile" class="form-label" for="ctFile">CT File</label>
-        <input class="form-control" type="file" id="ctFile" />
+      <div class="mb-3">
+        <button id="fileUploadButton" class="btn btn-primary" type="button"
+          onclick="fileUploadClicked(this)" disabled>
+          Upload new CT
+        </button>
       </div>
     </div>
-  </fieldset>
-</form>
+  </div>
+</fieldset>
diff --git a/src/web/templates/components/parameters/simulation_settings.html b/src/web/templates/components/parameters/simulation_settings.html
index e3342fb1..5ce39ff2 100644
--- a/src/web/templates/components/parameters/simulation_settings.html
+++ b/src/web/templates/components/parameters/simulation_settings.html
@@ -8,7 +8,7 @@
       !scenarioSettings.isSteadySimulation
     document.getElementById("simulationTypeSteady").checked =
       scenarioSettings.isSteadySimulation
-    setMaterialPropertiesPreBuiltValues(scenarioSettings.materialProperties)
+    setMaterialProperties(scenarioSettings.materialProperties)
   }
 
   function validateSimulationSettingsValues() {
diff --git a/src/web/views.py b/src/web/views.py
index c46aa8e5..826c53b2 100644
--- a/src/web/views.py
+++ b/src/web/views.py
@@ -1,8 +1,14 @@
 """Views for the web app."""
-from flask import Blueprint, jsonify, render_template, request
+import shutil
+import tempfile
+from pathlib import Path
+
+from flask import Blueprint, current_app, jsonify, render_template, request
 from pydantic import ValidationError
+from web.computed_tomography import get_available_cts, validate_ct
 from web.controller import (
     get_built_in_scenarios,
+    get_default_material_properties,
     get_scenario_layout,
     get_simulation_image,
 )
@@ -10,23 +16,29 @@
 from web.messages.requests import RenderLayoutRequest, SimulateRequest
 from web.messages.transducers import TransducerType
 from web.simulation_runner import SimulationRunner
+from werkzeug.utils import secure_filename
 
 bp = Blueprint("main", __name__, url_prefix="/")
 
+DEFAULT_CENTER_FREQUENCY = 5e5
+
 
 @bp.route("/")
 async def index():
     """Render the index page, listing all the built-in scenarios."""
-    title = "Neurotech Web App"
     return render_template(
         "index.html",
-        title=title,
+        title="Neurotech Web App",
         has_simulation=SimulationRunner().has_last_result,
         is_running_simulation=SimulationRunner().is_running,
         configuration=SimulationRunner().configuration,
         built_in_scenarios=get_built_in_scenarios(),
         all_materials=MaterialName.get_material_titles(),
+        all_material_properties=get_default_material_properties(
+            DEFAULT_CENTER_FREQUENCY
+        ).dict(),
         all_transducer_types=TransducerType.get_transducer_titles(),
+        available_cts=get_available_cts(current_app.config["CT_FOLDER"]),
     )
 
 
@@ -39,7 +51,10 @@ def simulate():
         # If the JSON data doesn't match the Pydantic model,
         # return a 400 Bad Request response
         return jsonify({"error": str(e)}), 400
-    result = SimulationRunner().run(get_simulation_image(config), config.dict())
+    result = SimulationRunner().run(
+        get_simulation_image(config=config, app_config=current_app.config),
+        config.dict(),
+    )
     if result.type == "simulation":
         data, image_format = result.data
         return f"<img src='data:image/{image_format};base64,{data}'/>"
@@ -78,3 +93,36 @@ async def render_layout():
     data = get_scenario_layout(config)
     image_format = "png"
     return f"<img src='data:image/{image_format};base64,{data}'/>"
+
+
+@bp.route("/ct_scan", methods=["POST"])
+def ct_scan():
+    """Upload a Computed Tomography scan and return the available scans."""
+    files = list(request.files.values())
+    temp_dir = Path(tempfile.mkdtemp())
+    saved_files = []
+    for file in files:
+        filename = secure_filename(file.filename)
+        file.save(temp_dir / filename)
+        saved_files.append(filename)
+
+    try:
+        ct_info = validate_ct(temp_dir, saved_files)
+        for filename in saved_files:
+            # moving files from the temporary directory to the CT_FOLDER
+            (temp_dir / filename).rename(current_app.config["CT_FOLDER"] / filename)
+    except Exception as e:
+        return jsonify({"error": str(e)}), 400
+    finally:
+        shutil.rmtree(temp_dir)
+
+    return jsonify(
+        {
+            "available_cts": get_available_cts(current_app.config["CT_FOLDER"]),
+            "selected_ct": {
+                "filename": ct_info.filename,
+                "shape": ct_info.shape,
+                "spacing": ct_info.spacing,
+            },
+        }
+    )
diff --git a/whitelist.txt b/whitelist.txt
index 5b7383e7..c6ca15c1 100644
--- a/whitelist.txt
+++ b/whitelist.txt
@@ -25,6 +25,7 @@ coroutine
 dataset
 der
 devito
+dicom
 dimensionality
 directivity
 docstring
@@ -94,6 +95,7 @@ str
 struct
 structs
 suptitle
+thresholding
 timesteps
 transcranial
 undersampled


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 431509c4..3c33e515 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -66,7 +66,7 @@ jobs:
         run: poetry install --no-interaction --no-root
 
       - name: Install stride
-        run: poetry run pip install git+https://github.com/trustimaging/stride
+        run: poetry run pip install git+https://github.com/trustimaging/stride@2520c59
 
       - name: Install neurotechdevkit
         run: poetry install --only-root
diff --git a/tests/neurotechdevkit/scenarios/test_data/PH1-BM4-SC1_STRIDE.mat b/tests/neurotechdevkit/scenarios/test_data/PH1-BM4-SC1_STRIDE.mat
new file mode 100644
index 00000000..11049a1d
Binary files /dev/null and b/tests/neurotechdevkit/scenarios/test_data/PH1-BM4-SC1_STRIDE.mat differ
diff --git a/tests/neurotechdevkit/scenarios/test_data/metrics-BM4-SC1_STRIDE.mat b/tests/neurotechdevkit/scenarios/test_data/metrics-BM4-SC1_STRIDE.mat
new file mode 100644
index 00000000..4bc64a9a
Binary files /dev/null and b/tests/neurotechdevkit/scenarios/test_data/metrics-BM4-SC1_STRIDE.mat differ
diff --git a/tests/neurotechdevkit/scenarios/test_metrics.py b/tests/neurotechdevkit/scenarios/test_metrics.py
index 2279cd43..e0bc0605 100644
--- a/tests/neurotechdevkit/scenarios/test_metrics.py
+++ b/tests/neurotechdevkit/scenarios/test_metrics.py
@@ -1,12 +1,18 @@
+from pathlib import Path
 from types import SimpleNamespace
+from typing import Optional, Union
 
+import hdf5storage
 import numpy as np
 import pytest
 
 from neurotechdevkit.results._metrics import (
     Conversions,
     calculate_all_metrics,
+    calculate_focal_fwhm,
     calculate_focal_gain,
+    calculate_focal_pressure,
+    calculate_focal_volume,
     calculate_i_pa_off_target,
     calculate_i_pa_target,
     calculate_i_ta_off_target,
@@ -14,10 +20,13 @@
     calculate_mechanical_index,
 )
 from neurotechdevkit.results._results import SteadyStateResult2D
+from neurotechdevkit.scenarios.built_in import Scenario1_2D
 
+GRID_SHAPE = (21, 31)
 CENTER_FREQUENCY = 1.5e6
 AMBIENT_PRESSURE = 2.5e6
 PEAK_PRESSURE = 6.3e6
+PEAK_PRESSURE_SHAPE = (2, 2)
 SPEED_OF_SOUND = 1800.0
 DENSITY = 1200.0
 
@@ -30,7 +39,12 @@ def patched_metric_fns(monkeypatch):
     monkeypatch.setattr(
         metrics, "calculate_mechanical_index", lambda result, layer: 0.1234
     )
+    monkeypatch.setattr(
+        metrics, "calculate_focal_pressure", lambda result, layer: 8.2e5
+    )
+    monkeypatch.setattr(metrics, "calculate_focal_volume", lambda result, layer: 123)
     monkeypatch.setattr(metrics, "calculate_focal_gain", lambda result: 2.4)
+    monkeypatch.setattr(metrics, "calculate_focal_fwhm", lambda result, axis, layer: 6)
     monkeypatch.setattr(metrics, "calculate_i_ta_target", lambda result: 2.34e-3)
     monkeypatch.setattr(metrics, "calculate_i_ta_off_target", lambda result: 3.45e-3)
     monkeypatch.setattr(metrics, "calculate_i_pa_target", lambda result: 45.6)
@@ -44,8 +58,10 @@ def fake_steady_state_matrix():
     The array has a shape of (21, 31) with a constant value of 2.5e6 everywhere except
     for four pixels which have a value of 6.3e6 in the lower half.
     """
-    steady_state = AMBIENT_PRESSURE * np.ones((21, 31))
-    steady_state[15:17, 24:26] = PEAK_PRESSURE
+    steady_state = AMBIENT_PRESSURE * np.ones(GRID_SHAPE)
+    steady_state[
+        15 : (15 + PEAK_PRESSURE_SHAPE[0]), 24 : (24 + PEAK_PRESSURE_SHAPE[1])
+    ] = PEAK_PRESSURE
     return steady_state
 
 
@@ -113,13 +129,17 @@ def test_calculate_all_metrics_has_correct_structure(fake_result, patched_metric
     metrics = calculate_all_metrics(fake_result)
     for _, data in metrics.items():
         assert set(data.keys()) == {"value", "unit-of-measurement", "description"}
-        assert isinstance(data["value"], float)
+        assert isinstance(data["value"], (float, int))
 
 
 def test_calculate_all_metrics_has_expected_metrics(fake_result, patched_metric_fns):
     """Verify that the metrics data contains the expected set of metrics"""
     expected_metrics = [
+        "focal_pressure",
+        "focal_volume",
         "focal_gain",
+        "FWHM_x",
+        "FWHM_y",
         "mechanical_index_all",
         "mechanical_index_brain",
         "mechanical_index_water",
@@ -137,7 +157,11 @@ def test_calculate_all_metrics_conversions(fake_result, patched_metric_fns):
     correct.
     """
     metrics = calculate_all_metrics(fake_result)
+    np.testing.assert_allclose(metrics["focal_pressure"]["value"], 8.2e5)
+    assert metrics["focal_volume"]["value"] == 123
     np.testing.assert_allclose(metrics["focal_gain"]["value"], 2.4)
+    assert metrics["FWHM_x"]["value"] == 6
+    assert metrics["FWHM_y"]["value"] == 6
     np.testing.assert_allclose(metrics["mechanical_index_all"]["value"], 0.1234)
     np.testing.assert_allclose(metrics["mechanical_index_brain"]["value"], 0.1234)
     np.testing.assert_allclose(metrics["mechanical_index_water"]["value"], 0.1234)
@@ -172,7 +196,33 @@ def test_calculate_mechanical_index_with_bottom_layer(fake_result):
     np.testing.assert_allclose(metric_value, expected)
 
 
-def test_calculate_focal_gain(fake_result):
+@pytest.mark.parametrize(
+    "layer, expected_focal_pressure",
+    [(None, PEAK_PRESSURE), ("water", AMBIENT_PRESSURE), ("brain", PEAK_PRESSURE)],
+)
+def test_calculate_focal_pressure(
+    fake_result: SteadyStateResult2D,
+    layer: Optional[str],
+    expected_focal_pressure: float,
+):
+    """Verify that the focal pressure is calculated as expected"""
+    metric_value = calculate_focal_pressure(fake_result, layer=layer)
+    np.testing.assert_allclose(metric_value, expected_focal_pressure)
+
+
+@pytest.mark.parametrize(
+    "layer, expected_focal_volume",
+    [(None, np.prod(PEAK_PRESSURE_SHAPE)), ("brain", np.prod(PEAK_PRESSURE_SHAPE))],
+)
+def test_calculate_focal_volume(
+    fake_result: SteadyStateResult2D, layer: Optional[str], expected_focal_volume: int
+):
+    """Verify that the focal pressure is calculated as expected"""
+    metric_value = calculate_focal_volume(fake_result, layer=layer)
+    np.testing.assert_allclose(metric_value, expected_focal_volume)
+
+
+def test_calculate_focal_gain(fake_result: SteadyStateResult2D):
     """Verify that the focal gain is calculated as expected"""
     metric_value = calculate_focal_gain(fake_result)
     # 16 pixels in the target region, 4 with high pressure
@@ -182,6 +232,30 @@ def test_calculate_focal_gain(fake_result):
     np.testing.assert_allclose(metric_value, expected)
 
 
+@pytest.mark.parametrize(
+    "axis,layer,expected_fwhm",
+    [
+        (0, None, PEAK_PRESSURE_SHAPE[0]),
+        (1, None, PEAK_PRESSURE_SHAPE[1]),
+        ("x", "brain", PEAK_PRESSURE_SHAPE[0]),
+        ("y", "brain", PEAK_PRESSURE_SHAPE[1]),
+    ],
+)
+def test_calculate_focal_fwhm(
+    fake_result: SteadyStateResult2D,
+    axis: Union[str, int],
+    layer: Optional[str],
+    expected_fwhm: int,
+):
+    """Verify that the full-width at half-maximum is calculated as expected"""
+    assert PEAK_PRESSURE > (AMBIENT_PRESSURE * 2), (
+        "This test's expected values are invalid if the peak pressure is not"
+        " 2x larger than the ambient pressure"
+    )
+    metric_value = calculate_focal_fwhm(fake_result, axis=axis, layer=layer)
+    np.testing.assert_allclose(metric_value, expected_fwhm)
+
+
 def test_calculate_i_ta_off_target(fake_result):
     """Verify that the time-averaged intensity out of the target is correct"""
     metric_value = calculate_i_ta_off_target(fake_result)
@@ -239,3 +313,65 @@ def test_unknown_conversion():
     """Unknown conversions should raise a ValueError"""
     with pytest.raises(ValueError):
         Conversions.convert("foo", "bar", 2.3)
+
+
+def test_compare_metrics_to_aubry2022():
+    """Tests FWHM calculation against Aubry et al. (2022) code:
+
+    To download the test data, see:
+        https://zenodo.org/record/6020543
+
+    To calculate the reference metrics, see:
+        https://github.com/agencyenterprise/transcranial-ultrasound-benchmarks/blob/master/strideMetrics.m
+    which is forked from:
+        https://github.com/ucl-bug/transcranial-ultrasound-benchmarks
+    """
+    steady_state_file = Path(__file__).parent.joinpath(
+        "test_data", "PH1-BM4-SC1_STRIDE.mat"
+    )
+    steady_state = hdf5storage.loadmat(str(steady_state_file))["p_amp"]
+
+    # _metrics.py expects a Results object
+    # Benchmark 4 corresponds to scenario 1
+    # The simulation resolution was slightly different, but
+    # let's use the object to wrap the data anyway
+    scenario = Scenario1_2D()
+    scenario.make_grid()
+    scenario.compile_problem()
+    result = SteadyStateResult2D(
+        scenario=scenario,
+        steady_state=steady_state,
+        center_frequency=scenario.center_frequency,
+        effective_dt=None,
+        pde=None,
+        shot=None,
+        wavefield=None,
+        traces=None,
+    )
+
+    # Calculate metrics on Aubry et al. 2022 results
+    metrics = calculate_all_metrics(result)
+
+    # Load in metrics calculated with published MATLAB code
+    expected_metrics_file = Path(__file__).parent.joinpath(
+        "test_data", "metrics-BM4-SC1_STRIDE.mat"
+    )
+    expected_metrics = hdf5storage.loadmat(str(expected_metrics_file))
+
+    # Check expected metrics
+    np.testing.assert_approx_equal(
+        metrics["focal_pressure"]["value"],
+        expected_metrics["max_amp_field1"].squeeze(),
+    )
+    assert (
+        metrics["focal_volume"]["value"]
+        == expected_metrics["focal_volume_num_vox_field1"].squeeze()
+    ), "Focal volume does not match expected"
+    for dim in ("x", "y"):
+        # FWHM in NDK is calculated as a whole number of pixels,
+        # while the MATLAB code estimates the FWHM with interpolation
+        np.testing.assert_almost_equal(
+            metrics[f"FWHM_{dim}"]["value"],
+            expected_metrics[f"FWHM_{dim}_field1"].squeeze(),
+            decimal=0,
+        )
diff --git a/tests/web/test_views.py b/tests/web/test_views.py
index c1ee481a..19458434 100644
--- a/tests/web/test_views.py
+++ b/tests/web/test_views.py
@@ -5,9 +5,10 @@
 from web.app import app
 
 
-@pytest.fixture(scope="module")
-def test_client():
+@pytest.fixture
+def test_client(tmp_path):
     """Create a test client for the flask app."""
+    app.config["CT_FOLDER"] = tmp_path
     with app.test_client() as testing_client:
         yield testing_client