Fluence from Moving MLC Aperture

.gitignore

@@ -2,6 +2,7 @@
 output/
 *.vt?
 examples/sample-data/
+build
 
 # Created by https://www.toptal.com/developers/gitignore/api/julia
 # Edit at https://www.toptal.com/developers/gitignore?templates=julia

Project.toml

@@ -1,7 +1,7 @@
 name = "DoseCalculations"
 uuid = "b0d14063-c48a-4081-83d5-5a5ab48cb495"
 authors = ["lmejn <33491793+lmejn@users.noreply.github.com>"]
-version = "0.1"
+version = "0.2"
 
 [deps]
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"

docs/make.jl

@@ -5,7 +5,9 @@ makedocs(
     modules = [DoseCalculations],
     sitename="DoseCalculations.jl Documentation",
     pages = [
-        "Dose Calculation Algorithms"=>["ScaledIsoplaneKernel.md",
-        ]
+        "Dose Volume"=>["ExternalSurfaces.md",
+                        "Structures.md"],
+        "Dose Calculation Algorithms"=>["ScaledIsoplaneKernel.md",],
+        "API"=>["API.md",]
     ]
     )

docs/src/API.md

@@ -0,0 +1,3 @@
+```@autodocs
+Modules = [DoseCalculations]
+```

docs/src/ExternalSurfaces.md

@@ -0,0 +1,73 @@
+# External Surfaces
+
+This library provides an interface to provide an external surface, which is used for source-surface distance (SSD) and depth calculations.
+Physically, the external surface denotes the boundary between air and the dose absorbing medium (*e.g.* the patient's skin or the surface of a phantom).
+
+All external surfaces are subtypes of `AbstractExternalSurface`, which exposes two methods:
+
+- [`getSSD(surf, p)`](@ref): Returns distance between the surface and the source for a given position in the dose volume (`p`)
+- [`getdepth(surf, p)`](@ref): Returns the depth of a given position in the dose volume (`p`)
+
+Both of these methods take an external surface (`surf`) and a dose point position (`pᵢ`). The position is usually a location where dose is computed, and is a three-element vector in the IEC Beam-Limiting-Device (BLD) coordinate system.
+
+All surfaces are assumed to be in the IEC BLD coordinate system. This means that they must undergo coordinate transformation in order to account for gantry and/or collimator rotation.
+
+## Types of External Surface
+
+### Plane Surface
+
+A simple plane surface, positioned at a constant source-surface distance along the central beam axis with normal pointed towards the source. It accounts the fact that SSD increases with off-axis position.
+
+```@repl
+using DoseCalculations
+surf = PlaneSurface(800.)
+getSSD(surf, [0., 0., 810.])
+getSSD(surf, [10., 20., 810.])
+```
+
+### Mesh Surface
+
+A general surface as defined by a 3D mesh provided by the user.
+
+This surface is created by suppling a 3D `mesh`,
+```julia
+MeshSurface(mesh)
+```
+which can be read from a `.ply` file, as described in [Structures](@ref)
+
+This is considered a slow method: for each position in the dose volume, it iterates through every face in the mesh (see [`DoseCalculations.intersect_mesh`](@ref)). A better alternative would be to create an [Isoplane Surface](@ref).
+
+
+### Isoplane Surface
+
+An isoplane surface collapses a general surface onto a plane centred at the isocentre with normal towards the source. It stores source-surface distances for points on the isoplane. The position of these points are determined by a pre-defined grid.
+
+This makes use of the fact that the distance between source and surface does not change along ray-lines.
+
+The Isoplane surface can be constructed by specifying either the `x` and `y` axes of the grid,
+```julia
+IsoplaneSurface(x, y, SAD)
+```
+or by specifying dose volume positions (`pos`), and plane grid spacings `Δx` and `Δy`.
+```julia
+IsoplaneSurface(pos, Δx, Δy, SAD)
+```
+In this case, the extent of the grid is computed from dose volume positions.
+
+The SSD values on the plane are computed from a `mesh_surf::DoseCalculations.MeshSurface`,
+```julia
+compute_SSD!(surf, mesh_surf)
+```
+
+### Constant Surface
+
+A constant surface simply returns the provided source-surface distance, regardless of position.
+
+```@repl
+using DoseCalculations
+surf = ConstantSurface(800.)
+getSSD(surf, [0., 0., 810.])
+getSSD(surf, [10., 20., 810.])
+```
+
+While such a surface is technically possible to implement in reality, this surface is largely used for debug purposes.

docs/src/Structures.md

@@ -0,0 +1,28 @@
+# Structures
+
+Dose volumes often contain structures of interest (*e.g.* the PTV or an OAR) where one may want to compute dose, or determine various metrics.
+
+This section describes how to use structures to:
+* Load structures from file.
+* Create dose points within a given structure 
+* Tag dose points based on which structure they are in.
+
+## Creating Structures
+
+One can load a mesh file using [`load_structure_from_ply`](@ref). Currently, only `.ply` files are supported.
+
+```julia
+structure = load_structure_from_ply("mesh.ply")
+```
+This returns a `Meshes.SimpleMesh`, as provided by [Meshes.jl](https://juliageometry.github.io/Meshes.jl/stable/meshes.html).
+
+Currently, loading structures direct from a DICOM RT Structure Set is not supported. Structures stored in DICOM files are not in a "friendly" mesh format, and usually require some modification before they can be used. There are many tools available which convert DICOM structures into a suitable mesh file (*e.g.* [3D Slicer](https://www.slicer.org/)).
+
+## Coordinate Transformations
+
+One can apply general transformations to any structure by calling the [`transform`](@ref) or [`transform!`](@ref) function:
+```julia
+T = fixed_to_bld(0., 0., 1000.)
+transform(structure, T)
+```
+`T` is a general transformation as provided by [CoordinateTransformations.jl](https://github.com/JuliaGeometry/CoordinateTransformations.jl). In this example, the helper function [`fixed_to_bld`](@ref) is used to create a transformation from the IEC Fixed to IEC BLD coordinate system. See [Coordinate Transformations](@ref) for more details.

src/BeamLimitingDevices.jl

@@ -157,6 +157,11 @@ end
 
 extract_subset(mlcx, mlc::MultiLeafCollimator, ylim::AbstractVector) = extract_subset(mlcx, mlc, ylim...)
 
+#--- MultiLeafCollimator Positions -------------------------------------------------------------------------------------
+
+AbstractMLCAperture{T} = AbstractMatrix{T}
+AbstractMLCSequence{T} = AbstractArray{T, 3}
+
 #--- Jaws --------------------------------------------------------------------------------------------------------------
 
 """ 

src/Fluence.jl

@@ -341,7 +341,7 @@ fluence(bixel::Bixel, jaws::Jaws) = fluence_from_rectangle(bixel, jaws.x, jaws.y
 
 From an MLC aperture.
 """
-function fluence(bixel::AbstractBixel{T}, mlcx, mlc::MultiLeafCollimator) where T<:AbstractFloat
+function fluence(bixel::AbstractBixel{T}, mlcx::AbstractMLCAperture, mlc::MultiLeafCollimator) where T<:AbstractFloat
 
     hw = 0.5*width(bixel, 2)
     i1 = max(1, locate(mlc, bixel[2]-hw))
@@ -362,37 +362,92 @@ From an MLC aperture using a given leaf index.
 This method assumes the bixel is entirely within the `i`th leaf track, and does
 not overlap with other leaves. Does not check whether these assumptions are true.
 """
-function fluence(bixel::AbstractBixel{T}, index::Int, mlcx) where T<:AbstractFloat
+function fluence(bixel::AbstractBixel{T}, index::Int, mlcx::AbstractMLCAperture) where T<:AbstractFloat
     overlap(position(bixel, 1), width(bixel, 1), mlcx[1, index], mlcx[2, index])
 end
 
-#=-- Moving Aperture fluences ------------------------------------------------------------------------------------------
+#--- Moving Aperture fluences ------------------------------------------------------------------------------------------
 
-    Deprecated.
+"""
+    fluence(bixel::AbstractBixel{T}, mlcx::AbstractMLCSequence, mlc::MultiLeafCollimator)
+
+From an MLC aperture sequence.
+"""
+function fluence(bixel::AbstractBixel{T}, mlcx::AbstractMLCSequence, mlc::MultiLeafCollimator) where T<:AbstractFloat
+    hw = 0.5*width(bixel, 2)
+    i1 = max(1, locate(mlc, bixel[2]-hw))
+    i2 = min(length(mlc), locate(mlc, bixel[2]-hw))
+
+    Ψ = zero(T)
+    @inbounds for j=i1:i2
+        Ψ += fluence_from_moving_aperture(bixel, (@view mlcx[:, j, 1]), (@view mlcx[:, j, 2]))
+    end
+    Ψ
+end
 
-    These functions compute the fluence from a leaf moving linearly from one position to another.
-=#
+"""
+    fluence(bixel::AbstractBixel{T}, index::Int, mlcx::AbstractMLCSequence)
 
-function fluence_from_moving_leaf(leafx_start, leafx_end, xb1, xb2)
-    leafx1, leafx2 = min(leafx1, leafx2), max(leafx1, leafx2)
-    t1 = (xb1 - leafx1)/(leafx2 - leafx1)
-    t2 = (xb2 - leafx1)/(leafx2 - leafx1)
+From an MLC aperture sequence using a given leaf index.
+"""
+function fluence(bixel::AbstractBixel{T}, index::Int, mlcx::AbstractMLCSequence) where T<:AbstractFloat
+    fluence_from_moving_aperture(bixel, (@view mlcx[:, index, 1]), (@view mlcx[:, index, 2]))
+end
 
-    if isnan(t1) || isnan(t2)
-        A_triangle = 0.
-    else
-        T = minmax(t1, 0., 1.) + minmax(t2, 0., 1.)
-        D = min(leafx2, xb2) - max(leafx1, xb1)
+"""
+    fluence_from_moving_aperture(bixel::AbstractBixel{T}, mlcx1, mlcx2)
 
-        A_triangle = 0.5*T*D
-    end
+From MLC leaf positions which move from `mlcx1` to `mlcx2`.
 
-    A_triangle + max(0., xb2 - leafx2)
+Computes the time-weighted fluence as the MLC moves from position `mlcx1` to `mlcx2`.
+Assumes the MLC leaves move in a straight line.
+"""
+function fluence_from_moving_aperture(bixel::AbstractBixel{T}, mlcx1, mlcx2) where T<:AbstractFloat
+    xL = bixel[1] - 0.5*width(bixel, 1)
+    xU = bixel[1] + 0.5*width(bixel, 1)
+
+    ΨB = fluence_onesided(mlcx1[1], mlcx2[1], xL, xU)
+    ΨA = fluence_onesided(mlcx1[2], mlcx2[2], xL, xU)
+
+    max(zero(T), ΨB - ΨA) # This is needed for cases where xA < xB
 end
 
-function fluence_from_moving_leaf(xb, Δx, mlcx1::Vector, mlcx2::Vector)
-    xb2 = xb + Δx
-    ΨA = fluence_from_moving_leaf(mlcx1[2], mlcx2[2], xb1, xb2)
-    ΨB = fluence_from_moving_leaf(mlcx1[1], mlcx2[1], xb1, xb2)
-    max(0., ΨB - ΨA)
+"""
+    fluence_onesided(xs, xf, xL, xU)
+
+Compute fluence for 1 leaf position, assuming the other is infinitely far away
+
+Computes the fluence for a leaf trajectory from `xs` to `xf`, over a bixel from
+`xL` to `xU`. The aperture is considered open to the right of the leaf position
+(*i.e.* leaf on the B bank). Assumes the bixel is fully in the leaf track.
+"""
+function fluence_onesided(xs, xf, xL, xU)
+
+    # If leaf positions are the same, compute static fluence
+    xs == xf && return (xU - minmax(xs, xL, xU))/(xU-xL)
+
+    # Otherwise compute moving fluence
+
+    x1, x2 = min(xs, xf), max(xs, xf) # If xf < xs, swap
+
+    xl = max(x1, xL) # Get left, centre and right segment positions
+    xc = min(x2, xU)
+    xr = xU
+
+    tl = leaf_trajectory(xl, x1, x2) # Compute the height of each segment pos.
+    tr = leaf_trajectory(xr, x1, x2)
+
+    A1 = 0.5*(tl + tr)*(xc-xl)  # Calculate trapezium segment are
+    A2 = xr-xc  # Calculate the remaining rectangular area
+    (A1 + A2)/(xU-xL)   # Add and scale
+
 end
+
+"""
+    leaf_trajectory(x, x1, x2)
+
+Compute the height of position `x` between `(x1, 0)` and `(x2, 1)`
+
+Used in `intersection_area`.
+"""
+leaf_trajectory(x, x1, x2) = minmax((x-x1)/(x2-x1), 0, 1)

src/utils/interpolation.jl

@@ -13,13 +13,6 @@ Scale the position `xi` within the positions `x1` and `x2`.
 """
 scale_to_cell(x1, x2, xi) = (xi - x1)/(x2 - x1)
 
-"""
-    minmax(xmin, xmax, xi)
-
-Limit the position `xi` between `xmin` and `xmax`.
-"""
-minmax(xmin, xmax, x) = max(xmin, min(xmax, x))
-
 #--- Non-Uniform Grid Location ------------------------------------------------
 
 """

src/utils/misc.jl

@@ -26,3 +26,10 @@ function scale_to_isoplane(pᵢ, z_plane)
     α = z_plane/z
     SVector(α*x, α*y)
 end
+
+"""
+    minmax(x, l, u)
+
+Return `x` if `l<=x<=u`, `l` if `x<l` or `u` if `u<x`
+"""
+minmax(x, l, u) = max(l, min(x, u)) # Ensures l<=x<=u