From f69e5fa46f7f1b7c56d9efc753e3f48ff3bebb36 Mon Sep 17 00:00:00 2001
From: "Pedro L. Esquinas" <pedro.esquinasfernandez@bccancer.bc.ca>
Date: Sat, 4 Oct 2025 22:18:49 -0700
Subject: [PATCH 1/4] Implement CSV to voxelS conversion

---
 pytheranostics/MiscTools/Tools.py | 328 +++++++++++++++++++++++++++++-
 1 file changed, 327 insertions(+), 1 deletion(-)

diff --git a/pytheranostics/MiscTools/Tools.py b/pytheranostics/MiscTools/Tools.py
index 735d06f..df42a54 100644
--- a/pytheranostics/MiscTools/Tools.py
+++ b/pytheranostics/MiscTools/Tools.py
@@ -1,9 +1,13 @@
+import math
 from datetime import datetime
-from typing import Dict, Tuple
+from typing import Dict, List, Tuple
 
 import numpy
+import pandas
 from scipy.ndimage import median_filter
 
+MEV_PER_G_TO_GY = 1.602176634e-10  # Gy per (MeV/g)
+
 
 def hu_to_rho(hu: numpy.ndarray) -> numpy.ndarray:
     """Convert a CT array, in HU into a density map in g/cc
@@ -324,3 +328,325 @@ def initialize_biokinetics_from_prior_cycle(
             print("")
 
     return config
+
+
+# Functions to generate voxel-s kernels
+# ----------------------------
+# CSV I/O
+# ----------------------------
+def read_dpk_csv(path: str) -> Tuple[numpy.ndarray, numpy.ndarray]:
+    """Read a DPK CSV (from Graves et al. 2019 https://aapm.onlinelibrary.wiley.com/doi/10.1002/mp.13789)
+    and return (r_mm, K_Gy_per_decay) as 1D arrays.
+
+    Parameters
+    ----------
+    path : str
+        Path to CSV file.
+
+    Returns
+    -------
+    Tuple[numpy.ndarray, numpy.ndarray]
+        Radius in mm, Dose Kernel in Gy/decay.
+
+    Raises
+    ------
+    ValueError
+        Expected CSV format not found.
+    ValueError
+        Expected columns not found.
+    """
+    # Read entire file, skip first line (metadata)
+    with open(path, "r", encoding="utf-8") as f:
+        lines = f.readlines()
+    if len(lines) < 3:
+        raise ValueError(f"File {path} doesn't look like expected CSV (too few lines).")
+
+    # Pandas read_csv from the second line (header is on line index 1)
+    from io import StringIO
+
+    buf = StringIO("".join(lines[1:]))
+    df = pandas.read_csv(buf)
+
+    # Required columns:
+    key_r = "Outer Radius of Bin (cm)"
+    key_d = "Dose per decay (MeV/g)"
+
+    df_columns = list(df.columns)
+    if key_r not in df_columns or key_d not in df_columns:
+        raise ValueError(
+            f"{path}: Required columns not found. Got columns: {df_columns}"
+        )
+
+    # Convert radii to mm
+    r_cm = df[key_r].to_numpy(dtype=float)
+    r_mm = r_cm * 10.0
+
+    # Convert MeV/g to Gy
+    d_mev_per_g = df[key_d].to_numpy(dtype=float)
+    K_Gy = d_mev_per_g * MEV_PER_G_TO_GY
+
+    # Ensure monotonic radii and non-negative K
+    order = numpy.argsort(r_mm)
+    r_mm = r_mm[order]
+    K_Gy = numpy.maximum(K_Gy[order], 0.0)
+
+    return r_mm, K_Gy
+
+
+def merge_multiple_csvs(
+    csvs: List[str], df_mm: float = 0.1, rmax_mm: float = 200.0
+) -> Tuple[numpy.ndarray, numpy.ndarray]:
+    """Read multiple CSVs and sum their K(r). Return common fine-grid r_mm and summed K(r).
+        - Interpolates each K onto a common 0.1 mm grid from r=0 to rmax_mm.
+        - Zero beyond the last provided radius in each file.
+
+    Parameters
+    ----------
+    csvs : List[str]
+        List of pahts to CSV files.
+    df_mm : float, optional
+        grid step in mm, by default 0.1
+    rmax_mm : float, optional
+        maximum radius in mm, by default 200.0
+
+    Returns
+    -------
+    Tuple[numpy.ndarray, numpy.ndarray]
+        Radius in mm, Summed Dose Kernel in Gy/decay.
+    """
+    # Read first to determine radius grid
+    r_mm, K_Gy = read_dpk_csv(csvs[0])
+
+    if len(csvs) == 1:
+        return r_mm, K_Gy
+
+    for p in csvs[1:]:
+        r_mm_, K_Gy_ = read_dpk_csv(p)
+
+        if not (r_mm_ == r_mm).all():
+            raise ValueError("All CSVs must have the same radius grid.")
+        K_Gy += K_Gy_
+
+    return r_mm, K_Gy
+
+
+# ----------------------------
+# Kernel construction
+# ----------------------------
+def build_3d_field_from_radial(
+    K_r_mm: numpy.ndarray,
+    r_mm: numpy.ndarray,
+    df_mm: float = 0.5,
+    Rmax_mm: float = 200.0,
+) -> numpy.ndarray:
+    """Build an isotropic 3-D field on a fine grid (spacing df_mm) by sampling K(r).
+    The field spans [-Rmax_mm, +Rmax_mm] along each axis.
+
+    Parameters
+    ----------
+    K_r_mm : numpy.ndarray
+        Dose Kernel in Gy/decay.
+    r_mm : numpy.ndarray
+        Radius in mm.
+    df_mm : float
+        Grid step in mm.
+    Rmax_mm : float
+        Maximum radius in mm.
+
+    Returns
+    -------
+    numpy.ndarray
+        3-D Dose Kernel field in Gy/decay.
+    """
+    # Determine N (odd) so that extent covers Rmax_mm
+    N = int(numpy.floor(2 * Rmax_mm / df_mm)) + 1
+    if N % 2 == 0:
+        N += 1
+    ax = (numpy.arange(N, dtype=float) - N // 2) * df_mm
+    X, Y, Z = numpy.meshgrid(ax, ax, ax, indexing="ij")
+    R = numpy.sqrt(X * X + Y * Y + Z * Z)
+
+    # Interpolate K(r) for all R (vectorized)
+    K_field = numpy.interp(
+        R.ravel(), r_mm, K_r_mm, left=K_r_mm[0] if r_mm[0] == 0 else 0.0, right=0.0
+    )
+    K_field = K_field.reshape(R.shape)
+    return K_field
+
+
+def box_filter_1d(arr: numpy.ndarray, M: int, axis: int) -> numpy.ndarray:
+    """Separable 1-D box filter (uniform average) along a given axis using cumulative sums.
+     Handles zero-padding at the boundaries.
+
+    Parameters
+    ----------
+    arr : numpy.ndarray
+        3-D dose kernel field.
+    M : int
+        Window length
+    axis : int
+        Axis along which to apply the filter (0, 1, or 2)
+
+    Returns
+    -------
+    numpy.ndarray
+        Filtered array.
+    """
+    if M <= 1:
+        return arr.copy()
+    # Move axis to front
+    arr_swapped = numpy.moveaxis(arr, axis, 0)
+
+    # Zero-pad by floor(M/2) on both sides
+    pad = M // 2
+    pad_before = pad
+    pad_after = M - 1 - pad
+    padded = numpy.pad(
+        arr_swapped,
+        ((pad_before, pad_after),) + tuple((0, 0) for _ in range(arr_swapped.ndim - 1)),
+        mode="constant",
+        constant_values=0.0,
+    )
+    # Cumulative sum along leading axis
+    csum = numpy.cumsum(padded, axis=0, dtype=float)
+    # Windowed sum: s[i] = csum[i+M] - csum[i]
+    s = csum[M:] - csum[:-M]
+    out = s / float(M)
+    # Move axis back
+    out = numpy.moveaxis(out, 0, axis)
+    return out
+
+
+def double_box_average_3d(
+    K_field: numpy.ndarray, L_mm: float, df_mm: float
+) -> Tuple[numpy.ndarray, int]:
+    """Apply two box averages of width L_mm (source & target) to the 3-D field.
+    Implemented as separable 1-D filters along x,y,z.
+
+    Parameters
+    ----------
+    K_field : numpy.ndarray
+        3-D dose kernel field.
+    L_mm : float
+        Box Width in mm. This is the voxel size of the coarse lattice. (i.e., SPECT voxel Size)
+    df_mm : float
+        Grid step in mm.
+
+    Returns
+    -------
+    Tuple[numpy.ndarray, int]
+        Averaged 3-D field, M (box length in voxels).
+    """
+    M = max(1, int(round(L_mm / df_mm)))
+    out = K_field
+    # First box (e.g., target average)
+    for ax in (0, 1, 2):
+        out = box_filter_1d(out, M, axis=ax)
+    # Second box (e.g., source average)
+    for ax in (0, 1, 2):
+        out = box_filter_1d(out, M, axis=ax)
+    return out, M
+
+
+def sample_on_coarse_lattice(
+    K_avg: numpy.ndarray, L_mm: float, df_mm: float, Rmax_mm: float
+) -> Tuple[numpy.ndarray, int]:
+    """Sample the averaged fine field at coarse lattice points (iL, jL, kL).
+
+    Parameters
+    ----------
+    K_avg : numpy.ndarray
+        Averaged 3-D dose kernel field.
+    L_mm : float
+        Box width in mm. This is the voxel size of the coarse lattice. (i.e., SPECT voxel Size)
+    df_mm : float
+        Grid step in mm.
+    Rmax_mm : float
+        Maximum radius in mm.
+
+    Returns
+    -------
+    Tuple[numpy.ndarray, int]
+        h: 3-D kernel (odd-sized cube)
+        Nc: radius in coarse voxels (so size = 2*Nc+1)
+    """
+    # Determine stride in voxels
+    stride = int(round(L_mm / df_mm))
+    N = K_avg.shape[0]
+    center = N // 2
+
+    # Determine Nc so that Nc*L_mm <= Rmax_mm (exclusive on the next)
+    Nc = int(numpy.floor(Rmax_mm / L_mm + 1e-9))
+    # Indices along one axis
+    idx = center + numpy.arange(-Nc, Nc + 1) * stride
+    # Guard within bounds
+    idx = idx[(idx >= 0) & (idx < N)]
+    # Build 3D sub-sampling
+    h = K_avg[numpy.ix_(idx, idx, idx)].copy()
+    # Ensure it's odd
+    assert h.shape[0] == h.shape[1] == h.shape[2], "Kernel shape must be cubic"
+    return h, (len(idx) - 1) // 2
+
+
+# ----------------------------
+# Sanity checks
+# ----------------------------
+def spherical_integral_K(r_mm: numpy.ndarray, K_r: numpy.ndarray) -> float:
+    """
+    Approximate ∫ K(r) dV over 0..∞ using discrete shells on the provided r grid (mm).
+    Returns value in Gy * mm^3 (convert to Gy*m^3 by *1e-9).
+    """
+    # Trapezoidal in r with shell volume 4π r^2 dr
+    r = r_mm
+    K = K_r
+    dr = numpy.diff(r)
+    r_mid = 0.5 * (r[:-1] + r[1:])
+    shell = 4.0 * math.pi * (r_mid**2) * ((K[:-1] + K[1:]) * 0.5) * dr
+    return float(numpy.sum(shell))  # Gy * mm^3
+
+
+def kernel_volume_sum(h: numpy.ndarray, L_mm: float) -> float:
+    """
+    Sum(h) * voxel_volume (Gy * mm^3), comparable to spherical_integral_K.
+    """
+    V_vox_mm3 = L_mm**3
+    return float(numpy.sum(h) * V_vox_mm3)
+
+
+def radius_for_fraction(
+    r_mm: numpy.ndarray, K_r: numpy.ndarray, frac: float = 0.995
+) -> float:
+    """
+    Return the first *tabulated* radius r[j] (mm) at which the cumulative deposited dose
+    (volume integral) exceeds `frac` of the total (default 99.5%).
+
+    This uses trapezoidal integration over W(r)=4π r^2 K(r) *per bin* and returns the
+    outer edge of the first bin whose cumulative integral crosses the threshold.
+    No sub-bin interpolation is performed.
+    """
+    if not (0.0 < frac < 1.0):
+        raise ValueError("frac must be in (0,1).")
+
+    r = numpy.asarray(r_mm, dtype=float)
+    K = numpy.asarray(K_r, dtype=float)
+
+    if r.ndim != 1 or K.ndim != 1 or r.size != K.size or r.size < 2:
+        raise ValueError("r_mm and K_r must be 1D arrays of the same length >= 2.")
+    if numpy.any(numpy.diff(r) <= 0):
+        raise ValueError("r_mm must be strictly increasing.")
+
+    W = 4.0 * numpy.pi * (r**2) * K
+    dr = numpy.diff(r)
+    bin_int = 0.5 * (W[:-1] + W[1:]) * dr
+    total = float(numpy.sum(bin_int))
+    if total <= 0.0:
+        return float(r[0])
+
+    target = frac * total
+    cum = 0.0
+    for i, Ti in enumerate(bin_int):
+        cum += Ti
+        if cum >= target:
+            return float(r[i + 1])  # outer radius of the crossing bin
+
+    return float(r[-1])

From 9df0d8cd276ad5840b2ee5fdea96d8787f957937 Mon Sep 17 00:00:00 2001
From: "Pedro L. Esquinas" <pedro.esquinasfernandez@bccancer.bc.ca>
Date: Thu, 30 Oct 2025 10:05:27 -0700
Subject: [PATCH 2/4] replace numpy.ndarray with NDArray

---
 pytheranostics/MiscTools/Tools.py | 63 ++++++++++++++-----------------
 1 file changed, 28 insertions(+), 35 deletions(-)

diff --git a/pytheranostics/MiscTools/Tools.py b/pytheranostics/MiscTools/Tools.py
index df42a54..19edf7e 100644
--- a/pytheranostics/MiscTools/Tools.py
+++ b/pytheranostics/MiscTools/Tools.py
@@ -4,20 +4,21 @@
 
 import numpy
 import pandas
+from numpy.typing import NDArray
 from scipy.ndimage import median_filter
 
 MEV_PER_G_TO_GY = 1.602176634e-10  # Gy per (MeV/g)
 
 
-def hu_to_rho(hu: numpy.ndarray) -> numpy.ndarray:
+def hu_to_rho(hu: NDArray) -> NDArray:
     """Convert a CT array, in HU into a density map in g/cc
     Conversion based on Schneider et al. 2000 (using GATE's material db example)
 
     Args:
-        hu (numpy.ndarray): _description_
+        hu (NDArray): _description_
 
     Returns:
-        numpy.ndarray: _description_
+        NDArray: _description_
     """
     # Define the bin edges for HU values
     bins = numpy.array(
@@ -334,7 +335,7 @@ def initialize_biokinetics_from_prior_cycle(
 # ----------------------------
 # CSV I/O
 # ----------------------------
-def read_dpk_csv(path: str) -> Tuple[numpy.ndarray, numpy.ndarray]:
+def read_dpk_csv(path: str) -> Tuple[NDArray, NDArray]:
     """Read a DPK CSV (from Graves et al. 2019 https://aapm.onlinelibrary.wiley.com/doi/10.1002/mp.13789)
     and return (r_mm, K_Gy_per_decay) as 1D arrays.
 
@@ -345,7 +346,7 @@ def read_dpk_csv(path: str) -> Tuple[numpy.ndarray, numpy.ndarray]:
 
     Returns
     -------
-    Tuple[numpy.ndarray, numpy.ndarray]
+    Tuple[NDArray, NDArray]
         Radius in mm, Dose Kernel in Gy/decay.
 
     Raises
@@ -393,9 +394,7 @@ def read_dpk_csv(path: str) -> Tuple[numpy.ndarray, numpy.ndarray]:
     return r_mm, K_Gy
 
 
-def merge_multiple_csvs(
-    csvs: List[str], df_mm: float = 0.1, rmax_mm: float = 200.0
-) -> Tuple[numpy.ndarray, numpy.ndarray]:
+def merge_multiple_csvs(csvs: List[str]) -> Tuple[NDArray, NDArray]:
     """Read multiple CSVs and sum their K(r). Return common fine-grid r_mm and summed K(r).
         - Interpolates each K onto a common 0.1 mm grid from r=0 to rmax_mm.
         - Zero beyond the last provided radius in each file.
@@ -404,14 +403,10 @@ def merge_multiple_csvs(
     ----------
     csvs : List[str]
         List of pahts to CSV files.
-    df_mm : float, optional
-        grid step in mm, by default 0.1
-    rmax_mm : float, optional
-        maximum radius in mm, by default 200.0
 
     Returns
     -------
-    Tuple[numpy.ndarray, numpy.ndarray]
+    Tuple[NDArray, NDArray]
         Radius in mm, Summed Dose Kernel in Gy/decay.
     """
     # Read first to determine radius grid
@@ -434,19 +429,19 @@ def merge_multiple_csvs(
 # Kernel construction
 # ----------------------------
 def build_3d_field_from_radial(
-    K_r_mm: numpy.ndarray,
-    r_mm: numpy.ndarray,
+    K_r_mm: NDArray,
+    r_mm: NDArray,
     df_mm: float = 0.5,
     Rmax_mm: float = 200.0,
-) -> numpy.ndarray:
+) -> NDArray:
     """Build an isotropic 3-D field on a fine grid (spacing df_mm) by sampling K(r).
     The field spans [-Rmax_mm, +Rmax_mm] along each axis.
 
     Parameters
     ----------
-    K_r_mm : numpy.ndarray
+    K_r_mm : NDArray
         Dose Kernel in Gy/decay.
-    r_mm : numpy.ndarray
+    r_mm : NDArray
         Radius in mm.
     df_mm : float
         Grid step in mm.
@@ -455,7 +450,7 @@ def build_3d_field_from_radial(
 
     Returns
     -------
-    numpy.ndarray
+    NDArray
         3-D Dose Kernel field in Gy/decay.
     """
     # Determine N (odd) so that extent covers Rmax_mm
@@ -474,13 +469,13 @@ def build_3d_field_from_radial(
     return K_field
 
 
-def box_filter_1d(arr: numpy.ndarray, M: int, axis: int) -> numpy.ndarray:
+def box_filter_1d(arr: NDArray, M: int, axis: int) -> NDArray:
     """Separable 1-D box filter (uniform average) along a given axis using cumulative sums.
      Handles zero-padding at the boundaries.
 
     Parameters
     ----------
-    arr : numpy.ndarray
+    arr : NDArray
         3-D dose kernel field.
     M : int
         Window length
@@ -489,7 +484,7 @@ def box_filter_1d(arr: numpy.ndarray, M: int, axis: int) -> numpy.ndarray:
 
     Returns
     -------
-    numpy.ndarray
+    NDArray
         Filtered array.
     """
     if M <= 1:
@@ -518,14 +513,14 @@ def box_filter_1d(arr: numpy.ndarray, M: int, axis: int) -> numpy.ndarray:
 
 
 def double_box_average_3d(
-    K_field: numpy.ndarray, L_mm: float, df_mm: float
-) -> Tuple[numpy.ndarray, int]:
+    K_field: NDArray, L_mm: float, df_mm: float
+) -> Tuple[NDArray, int]:
     """Apply two box averages of width L_mm (source & target) to the 3-D field.
     Implemented as separable 1-D filters along x,y,z.
 
     Parameters
     ----------
-    K_field : numpy.ndarray
+    K_field : NDArray
         3-D dose kernel field.
     L_mm : float
         Box Width in mm. This is the voxel size of the coarse lattice. (i.e., SPECT voxel Size)
@@ -534,7 +529,7 @@ def double_box_average_3d(
 
     Returns
     -------
-    Tuple[numpy.ndarray, int]
+    Tuple[NDArray, int]
         Averaged 3-D field, M (box length in voxels).
     """
     M = max(1, int(round(L_mm / df_mm)))
@@ -549,13 +544,13 @@ def double_box_average_3d(
 
 
 def sample_on_coarse_lattice(
-    K_avg: numpy.ndarray, L_mm: float, df_mm: float, Rmax_mm: float
-) -> Tuple[numpy.ndarray, int]:
+    K_avg: NDArray, L_mm: float, df_mm: float, Rmax_mm: float
+) -> Tuple[NDArray, int]:
     """Sample the averaged fine field at coarse lattice points (iL, jL, kL).
 
     Parameters
     ----------
-    K_avg : numpy.ndarray
+    K_avg : NDArray
         Averaged 3-D dose kernel field.
     L_mm : float
         Box width in mm. This is the voxel size of the coarse lattice. (i.e., SPECT voxel Size)
@@ -566,7 +561,7 @@ def sample_on_coarse_lattice(
 
     Returns
     -------
-    Tuple[numpy.ndarray, int]
+    Tuple[NDArray, int]
         h: 3-D kernel (odd-sized cube)
         Nc: radius in coarse voxels (so size = 2*Nc+1)
     """
@@ -591,7 +586,7 @@ def sample_on_coarse_lattice(
 # ----------------------------
 # Sanity checks
 # ----------------------------
-def spherical_integral_K(r_mm: numpy.ndarray, K_r: numpy.ndarray) -> float:
+def spherical_integral_K(r_mm: NDArray, K_r: NDArray) -> float:
     """
     Approximate ∫ K(r) dV over 0..∞ using discrete shells on the provided r grid (mm).
     Returns value in Gy * mm^3 (convert to Gy*m^3 by *1e-9).
@@ -605,7 +600,7 @@ def spherical_integral_K(r_mm: numpy.ndarray, K_r: numpy.ndarray) -> float:
     return float(numpy.sum(shell))  # Gy * mm^3
 
 
-def kernel_volume_sum(h: numpy.ndarray, L_mm: float) -> float:
+def kernel_volume_sum(h: NDArray, L_mm: float) -> float:
     """
     Sum(h) * voxel_volume (Gy * mm^3), comparable to spherical_integral_K.
     """
@@ -613,9 +608,7 @@ def kernel_volume_sum(h: numpy.ndarray, L_mm: float) -> float:
     return float(numpy.sum(h) * V_vox_mm3)
 
 
-def radius_for_fraction(
-    r_mm: numpy.ndarray, K_r: numpy.ndarray, frac: float = 0.995
-) -> float:
+def radius_for_fraction(r_mm: NDArray, K_r: NDArray, frac: float = 0.995) -> float:
     """
     Return the first *tabulated* radius r[j] (mm) at which the cumulative deposited dose
     (volume integral) exceeds `frac` of the total (default 99.5%).

From e3fe72ef288acfb75c5f92da8655a2e382c62d0f Mon Sep 17 00:00:00 2001
From: "Pedro L. Esquinas" <pedro.esquinasfernandez@bccancer.bc.ca>
Date: Thu, 30 Oct 2025 17:02:07 -0700
Subject: [PATCH 3/4] Add Lu177 kernels from Graves, data via LFS

---
 .gitattributes                                |   1 +
 pytheranostics/MiscTools/Tools.py             | 348 ++++--------------
 .../Lu177-10-mm-mGyperMBqs-Soft.csv           |   3 +
 .../Lu177-2.33-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-2.40-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-2.46-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-3.00-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-3.59-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-4.80-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-5.00-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-6.00-mm-mGyperMBqs-Soft.csv         |   3 +
 .../Lu177-9.28-mm-mGyperMBqs-Soft.csv         |   3 +
 pytheranostics/dosimetry/dvk.py               |  36 +-
 13 files changed, 118 insertions(+), 297 deletions(-)
 create mode 100644 .gitattributes
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-10-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-2.33-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-2.40-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-2.46-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-3.00-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-3.59-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-4.80-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-5.00-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-6.00-mm-mGyperMBqs-Soft.csv
 create mode 100644 pytheranostics/data/voxel_kernels/Lu177-9.28-mm-mGyperMBqs-Soft.csv

diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..7a697a3
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1 @@
+pytheranostics/data/voxel_kernels/*.csv filter=lfs diff=lfs merge=lfs -text
diff --git a/pytheranostics/MiscTools/Tools.py b/pytheranostics/MiscTools/Tools.py
index 19edf7e..adc92a5 100644
--- a/pytheranostics/MiscTools/Tools.py
+++ b/pytheranostics/MiscTools/Tools.py
@@ -1,6 +1,6 @@
-import math
 from datetime import datetime
-from typing import Dict, List, Tuple
+from pathlib import Path
+from typing import Dict, Tuple
 
 import numpy
 import pandas
@@ -335,311 +335,103 @@ def initialize_biokinetics_from_prior_cycle(
 # ----------------------------
 # CSV I/O
 # ----------------------------
-def read_dpk_csv(path: str) -> Tuple[NDArray, NDArray]:
-    """Read a DPK CSV (from Graves et al. 2019 https://aapm.onlinelibrary.wiley.com/doi/10.1002/mp.13789)
-    and return (r_mm, K_Gy_per_decay) as 1D arrays.
+def load_kernel_from_csv(path: Path) -> NDArray:
+    """Read Voxel Kernels from
+    Graves, S., Tiwari, A., Merrick, M., Hyer, D., Flynn, R.,
+    Kruzer, A., Nelson, A., Dewaraja, Y., Mirando, D.,
+    & Sunderland, J. (2023). Accurate resampling of radial dose point
+    kernels to a Cartesian matrix for voxelwise dose calculation (1.1)
+    [Data set]. Zenodo. https://doi.org/10.5281/zenodo.7596345
 
     Parameters
     ----------
-    path : str
-        Path to CSV file.
+    path : Path
+        Path to .csv file containing voxel-kernel values for positive Octant.
 
     Returns
     -------
-    Tuple[NDArray, NDArray]
-        Radius in mm, Dose Kernel in Gy/decay.
+    NDArray
+        array containing kernel values.
 
     Raises
     ------
     ValueError
-        Expected CSV format not found.
-    ValueError
-        Expected columns not found.
-    """
-    # Read entire file, skip first line (metadata)
-    with open(path, "r", encoding="utf-8") as f:
-        lines = f.readlines()
-    if len(lines) < 3:
-        raise ValueError(f"File {path} doesn't look like expected CSV (too few lines).")
-
-    # Pandas read_csv from the second line (header is on line index 1)
-    from io import StringIO
-
-    buf = StringIO("".join(lines[1:]))
-    df = pandas.read_csv(buf)
-
-    # Required columns:
-    key_r = "Outer Radius of Bin (cm)"
-    key_d = "Dose per decay (MeV/g)"
-
-    df_columns = list(df.columns)
-    if key_r not in df_columns or key_d not in df_columns:
-        raise ValueError(
-            f"{path}: Required columns not found. Got columns: {df_columns}"
-        )
-
-    # Convert radii to mm
-    r_cm = df[key_r].to_numpy(dtype=float)
-    r_mm = r_cm * 10.0
-
-    # Convert MeV/g to Gy
-    d_mev_per_g = df[key_d].to_numpy(dtype=float)
-    K_Gy = d_mev_per_g * MEV_PER_G_TO_GY
-
-    # Ensure monotonic radii and non-negative K
-    order = numpy.argsort(r_mm)
-    r_mm = r_mm[order]
-    K_Gy = numpy.maximum(K_Gy[order], 0.0)
-
-    return r_mm, K_Gy
-
-
-def merge_multiple_csvs(csvs: List[str]) -> Tuple[NDArray, NDArray]:
-    """Read multiple CSVs and sum their K(r). Return common fine-grid r_mm and summed K(r).
-        - Interpolates each K onto a common 0.1 mm grid from r=0 to rmax_mm.
-        - Zero beyond the last provided radius in each file.
-
-    Parameters
-    ----------
-    csvs : List[str]
-        List of pahts to CSV files.
-
-    Returns
-    -------
-    Tuple[NDArray, NDArray]
-        Radius in mm, Summed Dose Kernel in Gy/decay.
+        _description_
     """
-    # Read first to determine radius grid
-    r_mm, K_Gy = read_dpk_csv(csvs[0])
-
-    if len(csvs) == 1:
-        return r_mm, K_Gy
-
-    for p in csvs[1:]:
-        r_mm_, K_Gy_ = read_dpk_csv(p)
-
-        if not (r_mm_ == r_mm).all():
-            raise ValueError("All CSVs must have the same radius grid.")
-        K_Gy += K_Gy_
-
-    return r_mm, K_Gy
-
+    df = pandas.read_csv(
+        path,
+        header=None,
+        skip_blank_lines=False,  # <- important
+    )
 
-# ----------------------------
-# Kernel construction
-# ----------------------------
-def build_3d_field_from_radial(
-    K_r_mm: NDArray,
-    r_mm: NDArray,
-    df_mm: float = 0.5,
-    Rmax_mm: float = 200.0,
-) -> NDArray:
-    """Build an isotropic 3-D field on a fine grid (spacing df_mm) by sampling K(r).
-    The field spans [-Rmax_mm, +Rmax_mm] along each axis.
+    # rows that are completely blank will be all NaN
+    blank_mask = df.isna().all(axis=1)
 
-    Parameters
-    ----------
-    K_r_mm : NDArray
-        Dose Kernel in Gy/decay.
-    r_mm : NDArray
-        Radius in mm.
-    df_mm : float
-        Grid step in mm.
-    Rmax_mm : float
-        Maximum radius in mm.
-
-    Returns
-    -------
-    NDArray
-        3-D Dose Kernel field in Gy/decay.
-    """
-    # Determine N (odd) so that extent covers Rmax_mm
-    N = int(numpy.floor(2 * Rmax_mm / df_mm)) + 1
-    if N % 2 == 0:
-        N += 1
-    ax = (numpy.arange(N, dtype=float) - N // 2) * df_mm
-    X, Y, Z = numpy.meshgrid(ax, ax, ax, indexing="ij")
-    R = numpy.sqrt(X * X + Y * Y + Z * Z)
-
-    # Interpolate K(r) for all R (vectorized)
-    K_field = numpy.interp(
-        R.ravel(), r_mm, K_r_mm, left=K_r_mm[0] if r_mm[0] == 0 else 0.0, right=0.0
-    )
-    K_field = K_field.reshape(R.shape)
-    return K_field
+    # make a group id that increments every time we see a blank row
+    # e.g. rows -> 0..172 (block 0), blank, 174..346 (block 1), ...
+    block_ids = blank_mask.cumsum()
 
+    # drop the blank rows themselves
+    df_data = df[~blank_mask].reset_index(drop=True)
+    block_ids = block_ids[~blank_mask].to_numpy()
 
-def box_filter_1d(arr: NDArray, M: int, axis: int) -> NDArray:
-    """Separable 1-D box filter (uniform average) along a given axis using cumulative sums.
-     Handles zero-padding at the boundaries.
+    # infer N
+    # number of rows in the first block
+    first_block_rows = (block_ids == block_ids[0]).sum()
+    N = first_block_rows
 
-    Parameters
-    ----------
-    arr : NDArray
-        3-D dose kernel field.
-    M : int
-        Window length
-    axis : int
-        Axis along which to apply the filter (0, 1, or 2)
+    # now we have (num_blocks * N) rows, each with N columns
+    # we can groupby block_id and build the 3D array
+    blocks = []
+    for b in numpy.unique(block_ids):
+        block_df = df_data[block_ids == b]
+        arr = block_df.to_numpy(dtype=float)
+        if arr.shape != (N, N):
+            raise ValueError(f"Block {b} has shape {arr.shape}, expected {(N, N)}")
+        blocks.append(arr)
 
-    Returns
-    -------
-    NDArray
-        Filtered array.
-    """
-    if M <= 1:
-        return arr.copy()
-    # Move axis to front
-    arr_swapped = numpy.moveaxis(arr, axis, 0)
-
-    # Zero-pad by floor(M/2) on both sides
-    pad = M // 2
-    pad_before = pad
-    pad_after = M - 1 - pad
-    padded = numpy.pad(
-        arr_swapped,
-        ((pad_before, pad_after),) + tuple((0, 0) for _ in range(arr_swapped.ndim - 1)),
-        mode="constant",
-        constant_values=0.0,
-    )
-    # Cumulative sum along leading axis
-    csum = numpy.cumsum(padded, axis=0, dtype=float)
-    # Windowed sum: s[i] = csum[i+M] - csum[i]
-    s = csum[M:] - csum[:-M]
-    out = s / float(M)
-    # Move axis back
-    out = numpy.moveaxis(out, 0, axis)
-    return out
-
-
-def double_box_average_3d(
-    K_field: NDArray, L_mm: float, df_mm: float
-) -> Tuple[NDArray, int]:
-    """Apply two box averages of width L_mm (source & target) to the 3-D field.
-    Implemented as separable 1-D filters along x,y,z.
+    return expand_octant_to_full(numpy.stack(blocks, axis=0))
 
-    Parameters
-    ----------
-    K_field : NDArray
-        3-D dose kernel field.
-    L_mm : float
-        Box Width in mm. This is the voxel size of the coarse lattice. (i.e., SPECT voxel Size)
-    df_mm : float
-        Grid step in mm.
 
-    Returns
-    -------
-    Tuple[NDArray, int]
-        Averaged 3-D field, M (box length in voxels).
-    """
-    M = max(1, int(round(L_mm / df_mm)))
-    out = K_field
-    # First box (e.g., target average)
-    for ax in (0, 1, 2):
-        out = box_filter_1d(out, M, axis=ax)
-    # Second box (e.g., source average)
-    for ax in (0, 1, 2):
-        out = box_filter_1d(out, M, axis=ax)
-    return out, M
-
-
-def sample_on_coarse_lattice(
-    K_avg: NDArray, L_mm: float, df_mm: float, Rmax_mm: float
-) -> Tuple[NDArray, int]:
-    """Sample the averaged fine field at coarse lattice points (iL, jL, kL).
+def expand_octant_to_full(octant: NDArray) -> NDArray:
+    """Given a 3D array of shape (H, H, H) that represents the center voxel
+    at [0,0,0] and the +x, +y, +z directions (i.e. the positive octant),
+    reconstruct the full symmetric kernel of shape (2H-1, 2H-1, 2H-1).
 
     Parameters
     ----------
-    K_avg : NDArray
-        Averaged 3-D dose kernel field.
-    L_mm : float
-        Box width in mm. This is the voxel size of the coarse lattice. (i.e., SPECT voxel Size)
-    df_mm : float
-        Grid step in mm.
-    Rmax_mm : float
-        Maximum radius in mm.
+    octant : numpy.ndarray
+        Positive octant of the kernel, shape (H, H, H)
 
     Returns
     -------
-    Tuple[NDArray, int]
-        h: 3-D kernel (odd-sized cube)
-        Nc: radius in coarse voxels (so size = 2*Nc+1)
-    """
-    # Determine stride in voxels
-    stride = int(round(L_mm / df_mm))
-    N = K_avg.shape[0]
-    center = N // 2
-
-    # Determine Nc so that Nc*L_mm <= Rmax_mm (exclusive on the next)
-    Nc = int(numpy.floor(Rmax_mm / L_mm + 1e-9))
-    # Indices along one axis
-    idx = center + numpy.arange(-Nc, Nc + 1) * stride
-    # Guard within bounds
-    idx = idx[(idx >= 0) & (idx < N)]
-    # Build 3D sub-sampling
-    h = K_avg[numpy.ix_(idx, idx, idx)].copy()
-    # Ensure it's odd
-    assert h.shape[0] == h.shape[1] == h.shape[2], "Kernel shape must be cubic"
-    return h, (len(idx) - 1) // 2
-
+    numpy.ndarray
+        The full symmetric kernel, shape (2H-1, 2H-1, 2H-1)
 
-# ----------------------------
-# Sanity checks
-# ----------------------------
-def spherical_integral_K(r_mm: NDArray, K_r: NDArray) -> float:
-    """
-    Approximate ∫ K(r) dV over 0..∞ using discrete shells on the provided r grid (mm).
-    Returns value in Gy * mm^3 (convert to Gy*m^3 by *1e-9).
+    Raises
+    ------
+    ValueError
+        If the input octant is not 3D or not cubic.
     """
-    # Trapezoidal in r with shell volume 4π r^2 dr
-    r = r_mm
-    K = K_r
-    dr = numpy.diff(r)
-    r_mid = 0.5 * (r[:-1] + r[1:])
-    shell = 4.0 * math.pi * (r_mid**2) * ((K[:-1] + K[1:]) * 0.5) * dr
-    return float(numpy.sum(shell))  # Gy * mm^3
+    if octant.ndim != 3:
+        raise ValueError("octant must be 3D")
 
+    H = octant.shape[0]
+    if not (octant.shape[1] == H and octant.shape[2] == H):
+        raise ValueError("octant must be cubic (H×H×H)")
 
-def kernel_volume_sum(h: NDArray, L_mm: float) -> float:
-    """
-    Sum(h) * voxel_volume (Gy * mm^3), comparable to spherical_integral_K.
-    """
-    V_vox_mm3 = L_mm**3
-    return float(numpy.sum(h) * V_vox_mm3)
+    # mirror in x (axis=0): [-x | 0..+x]
+    # octant[1:][::-1, :, :] gives slices 1..H-1 reversed → x=-1, -2, ...
+    full_x = numpy.concatenate([octant[1:][::-1, :, :], octant], axis=0)  # (2H-1, H, H)
 
+    # mirror in y (axis=1): [-y | 0..+y]
+    full_xy = numpy.concatenate(
+        [full_x[:, 1:][:, ::-1, :], full_x], axis=1
+    )  # (2H-1, 2H-1, H)
 
-def radius_for_fraction(r_mm: NDArray, K_r: NDArray, frac: float = 0.995) -> float:
-    """
-    Return the first *tabulated* radius r[j] (mm) at which the cumulative deposited dose
-    (volume integral) exceeds `frac` of the total (default 99.5%).
+    # mirror in z (axis=2): [-z | 0..+z]
+    full_xyz = numpy.concatenate(
+        [full_xy[:, :, 1:][:, :, ::-1], full_xy], axis=2
+    )  # (2H-1, 2H-1, 2H-1)
 
-    This uses trapezoidal integration over W(r)=4π r^2 K(r) *per bin* and returns the
-    outer edge of the first bin whose cumulative integral crosses the threshold.
-    No sub-bin interpolation is performed.
-    """
-    if not (0.0 < frac < 1.0):
-        raise ValueError("frac must be in (0,1).")
-
-    r = numpy.asarray(r_mm, dtype=float)
-    K = numpy.asarray(K_r, dtype=float)
-
-    if r.ndim != 1 or K.ndim != 1 or r.size != K.size or r.size < 2:
-        raise ValueError("r_mm and K_r must be 1D arrays of the same length >= 2.")
-    if numpy.any(numpy.diff(r) <= 0):
-        raise ValueError("r_mm must be strictly increasing.")
-
-    W = 4.0 * numpy.pi * (r**2) * K
-    dr = numpy.diff(r)
-    bin_int = 0.5 * (W[:-1] + W[1:]) * dr
-    total = float(numpy.sum(bin_int))
-    if total <= 0.0:
-        return float(r[0])
-
-    target = frac * total
-    cum = 0.0
-    for i, Ti in enumerate(bin_int):
-        cum += Ti
-        if cum >= target:
-            return float(r[i + 1])  # outer radius of the crossing bin
-
-    return float(r[-1])
+    return full_xyz
diff --git a/pytheranostics/data/voxel_kernels/Lu177-10-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-10-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..b4138a1
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-10-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c1169fe687e1a8dd8de9cceddec43cdeeff4d69d000ca6a032a36e83caf58683
+size 810411
diff --git a/pytheranostics/data/voxel_kernels/Lu177-2.33-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-2.33-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..3f684f2
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-2.33-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:23cae7c63bf570bcbe25e6ce1dfffb0e351c2d8e2a6d93a4a4d08d4fae172d74
+size 60440662
diff --git a/pytheranostics/data/voxel_kernels/Lu177-2.40-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-2.40-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..5c05926
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-2.40-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:72f2cbc4ea2a10126558aaf21004f5ee937d3973f2df7e31ba4a06606f8310be
+size 55465783
diff --git a/pytheranostics/data/voxel_kernels/Lu177-2.46-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-2.46-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..c38b47e
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-2.46-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:66f31a3a2306b8091b06e1ecda4c0f994b7a64a9cc048ae72e0e7ef4a73ef1b8
+size 51687432
diff --git a/pytheranostics/data/voxel_kernels/Lu177-3.00-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-3.00-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..4f75b3f
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-3.00-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2668bbf2e9d6533c18412b47d32fb4ba9a6e8e0cc7bc8520684a8b83ee7ab3a2
+size 28655946
diff --git a/pytheranostics/data/voxel_kernels/Lu177-3.59-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-3.59-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..c8dd0a3
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-3.59-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:21dc4b04850c1343a9d1a6e862f9cfbab95997d3a61c7ad0818d7b4cadfd0462
+size 16584954
diff --git a/pytheranostics/data/voxel_kernels/Lu177-4.80-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-4.80-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..95d84c1
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-4.80-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34a97c1053e8c54a4b9109bf7832018384f9479166a0f54e7c42924fe63488ea
+size 7184417
diff --git a/pytheranostics/data/voxel_kernels/Lu177-5.00-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-5.00-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..3caf654
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-5.00-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:426771c730fd36e4781ea230ff69b9c779b9391b597ef4caad34b7b630ecfc0d
+size 6239147
diff --git a/pytheranostics/data/voxel_kernels/Lu177-6.00-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-6.00-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..66b6dab
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-6.00-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2db81acd6c1f28c8aaf565c2b9a359a9c58dbb3129d24249a31fe47810c81d08
+size 3586587
diff --git a/pytheranostics/data/voxel_kernels/Lu177-9.28-mm-mGyperMBqs-Soft.csv b/pytheranostics/data/voxel_kernels/Lu177-9.28-mm-mGyperMBqs-Soft.csv
new file mode 100644
index 0000000..30b49b9
--- /dev/null
+++ b/pytheranostics/data/voxel_kernels/Lu177-9.28-mm-mGyperMBqs-Soft.csv
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:520fc078f912a8f77e1fd8f32993fbea40a1646c3efba00bca98610efaa03187
+size 996357
diff --git a/pytheranostics/dosimetry/dvk.py b/pytheranostics/dosimetry/dvk.py
index 282fbbd..2ed5a52 100644
--- a/pytheranostics/dosimetry/dvk.py
+++ b/pytheranostics/dosimetry/dvk.py
@@ -1,38 +1,36 @@
 import os
+from pathlib import Path
 from typing import Optional
 
 import numpy
 from scipy import signal
 
-from pytheranostics.MiscTools.Tools import hu_to_rho
+from pytheranostics.MiscTools.Tools import hu_to_rho, load_kernel_from_csv
 
 
 class DoseVoxelKernel:
     def __init__(self, isotope: str, voxel_size_mm: float) -> None:
-        """_summary_
+        """Initialize Dose Voxel-Kernel for convolution-based dosimetry.
 
         Args:
-            isotope (str): _description_
-            voxel_size_mm (float): _description_
+            isotope (str): Isotope name, e.g., "Lu177".
+            voxel_size_mm (float): Voxel size in mm, e.g., 4.80
+        Raises:
+            FileNotFoundError: If Voxel-Kernel file is not found.
         """
-        try:
-            self.kernel = numpy.fromfile(
-                os.path.dirname(__file__)
-                + f"/../data/voxel_kernels/{isotope}-{voxel_size_mm:1.2f}-mm-mGyperMBqs-SoftICRP.img",
-                dtype=numpy.float32,
-            )
-        except FileNotFoundError:
-            print(
-                f" >> Voxel Kernel for SPECT voxel size ({voxel_size_mm:2.2f} mm) not found. Using default kernel for 4.8 mm voxels..."
-            )
 
-            self.kernel = numpy.fromfile(
-                os.path.dirname(__file__)
-                + f"/../data/voxel_kernels/{isotope}-4.80-mm-mGyperMBqs-SoftICRP.img",
-                dtype=numpy.float32,
+        # Set file path for Voxel-Kernel.
+        kernel_file = Path(
+            os.path.dirname(__file__)
+            + f"/../data/voxel_kernels/{isotope}-{voxel_size_mm:2.2f}-mm-mGyperMBqs-Soft.csv"
+        )
+
+        if not kernel_file.exists():
+            raise FileNotFoundError(
+                f" >> Voxel Kernel for SPECT voxel size ({voxel_size_mm:2.2f} mm) not found."
             )
 
-        self.kernel = self.kernel.reshape((51, 51, 51)).astype(numpy.float64)
+        self.kernel = load_kernel_from_csv(path=kernel_file)
 
     def tia_to_dose(
         self, tia_mbq_s: numpy.ndarray, ct: Optional[numpy.ndarray] = None

From ddcd04e8867b0f9efc477a5d2b2b5826f6bb184a Mon Sep 17 00:00:00 2001
From: "Pedro L. Esquinas" <pedro.esquinasfernandez@bccancer.bc.ca>
Date: Thu, 30 Oct 2025 22:27:43 -0700
Subject: [PATCH 4/4] handle data with voxel size different than existing
 voxel-s kernels

---
 pytheranostics/dosimetry/dvk.py | 29 +++++++++++++++++++++--------
 1 file changed, 21 insertions(+), 8 deletions(-)

diff --git a/pytheranostics/dosimetry/dvk.py b/pytheranostics/dosimetry/dvk.py
index 2ed5a52..e8d0f30 100644
--- a/pytheranostics/dosimetry/dvk.py
+++ b/pytheranostics/dosimetry/dvk.py
@@ -10,26 +10,39 @@
 
 class DoseVoxelKernel:
     def __init__(self, isotope: str, voxel_size_mm: float) -> None:
-        """Initialize Dose Voxel-Kernel for convolution-based dosimetry.
+        """Initialize Dose Voxel-Kernel for convolution-based dosimetry. If kernel for specified voxel size is not found, the closest
+        available size will be used.
 
         Args:
             isotope (str): Isotope name, e.g., "Lu177".
             voxel_size_mm (float): Voxel size in mm, e.g., 4.80
-        Raises:
-            FileNotFoundError: If Voxel-Kernel file is not found.
         """
 
         # Set file path for Voxel-Kernel.
+        kernel_dir = Path(os.path.dirname(__file__) + "/../data/voxel_kernels/")
+        kernel_files = list(kernel_dir.glob(f"{isotope}-*-mm-mGyperMBqs-Soft.csv"))
+        kernels = {
+            size: file
+            for file in kernel_files
+            for size in [float(file.stem.split("-")[1])]
+        }
+
+        if voxel_size_mm not in kernels:
+            # get closest available kernel size
+            available_sizes = numpy.array(list(kernels.keys()))
+            closest_size = available_sizes[
+                numpy.argmin(numpy.abs(available_sizes - voxel_size_mm))
+            ]
+            print(
+                f" >> Voxel Kernel for SPECT voxel size ({voxel_size_mm:2.2f} mm) not found. Using closest available size: {closest_size:2.2f} mm."
+            )
+            voxel_size_mm = closest_size
+
         kernel_file = Path(
             os.path.dirname(__file__)
             + f"/../data/voxel_kernels/{isotope}-{voxel_size_mm:2.2f}-mm-mGyperMBqs-Soft.csv"
         )
 
-        if not kernel_file.exists():
-            raise FileNotFoundError(
-                f" >> Voxel Kernel for SPECT voxel size ({voxel_size_mm:2.2f} mm) not found."
-            )
-
         self.kernel = load_kernel_from_csv(path=kernel_file)
 
     def tia_to_dose(