diff --git a/.gitignore b/.gitignore index a87e63d..ddb5015 100644 --- a/.gitignore +++ b/.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 diff --git a/Project.toml b/Project.toml index c291e8d..80d5b8e 100644 --- a/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index d45f6f7..a23e172 100644 --- a/docs/make.jl +++ b/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",] ] ) diff --git a/docs/src/API.md b/docs/src/API.md new file mode 100644 index 0000000..a08d8bf --- /dev/null +++ b/docs/src/API.md @@ -0,0 +1,3 @@ +```@autodocs +Modules = [DoseCalculations] +``` \ No newline at end of file diff --git a/docs/src/ExternalSurfaces.md b/docs/src/ExternalSurfaces.md new file mode 100644 index 0000000..6389687 --- /dev/null +++ b/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. diff --git a/docs/src/Structures.md b/docs/src/Structures.md new file mode 100644 index 0000000..9c74419 --- /dev/null +++ b/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. diff --git a/src/BeamLimitingDevices.jl b/src/BeamLimitingDevices.jl index 1c5b29d..915e00e 100644 --- a/src/BeamLimitingDevices.jl +++ b/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 -------------------------------------------------------------------------------------------------------------- """ diff --git a/src/Fluence.jl b/src/Fluence.jl index b03a698..b588daa 100644 --- a/src/Fluence.jl +++ b/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) diff --git a/src/utils/interpolation.jl b/src/utils/interpolation.jl index d9ff5f1..58f0f27 100644 --- a/src/utils/interpolation.jl +++ b/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 ------------------------------------------------ """ diff --git a/src/utils/misc.jl b/src/utils/misc.jl index 968d7ee..993f1ea 100644 --- a/src/utils/misc.jl +++ b/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