From 5f43abd72d07e79c7e1a72fe004f201275ad82d3 Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Thu, 28 Jul 2022 12:51:22 +1000 Subject: [PATCH 1/9] Add docs/build folder to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) 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 From 8f85bab7f098a6c79797bfbad9b8eb3893b8ad5d Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Thu, 28 Jul 2022 12:53:04 +1000 Subject: [PATCH 2/9] Added ExternalSurfaces documentation --- docs/make.jl | 5 ++- docs/src/API.md | 3 ++ docs/src/ExternalSurfaces.md | 73 ++++++++++++++++++++++++++++++++++++ docs/src/Structures.md | 3 ++ 4 files changed, 82 insertions(+), 2 deletions(-) create mode 100644 docs/src/API.md create mode 100644 docs/src/ExternalSurfaces.md create mode 100644 docs/src/Structures.md diff --git a/docs/make.jl b/docs/make.jl index d45f6f7..317a381 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,7 +5,8 @@ makedocs( modules = [DoseCalculations], sitename="DoseCalculations.jl Documentation", pages = [ - "Dose Calculation Algorithms"=>["ScaledIsoplaneKernel.md", - ] + "Dose Volume"=>["ExternalSurfaces.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..8dfcbd3 --- /dev/null +++ b/docs/src/Structures.md @@ -0,0 +1,3 @@ +# Structures + +Test \ No newline at end of file From fcf9f798d66a7bfd2ecc18908e18168d2c868fe9 Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Thu, 28 Jul 2022 15:55:36 +1000 Subject: [PATCH 3/9] Added Structures documentation --- docs/make.jl | 3 ++- docs/src/Structures.md | 27 ++++++++++++++++++++++++++- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 317a381..a23e172 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,7 +5,8 @@ makedocs( modules = [DoseCalculations], sitename="DoseCalculations.jl Documentation", pages = [ - "Dose Volume"=>["ExternalSurfaces.md",], + "Dose Volume"=>["ExternalSurfaces.md", + "Structures.md"], "Dose Calculation Algorithms"=>["ScaledIsoplaneKernel.md",], "API"=>["API.md",] ] diff --git a/docs/src/Structures.md b/docs/src/Structures.md index 8dfcbd3..9c74419 100644 --- a/docs/src/Structures.md +++ b/docs/src/Structures.md @@ -1,3 +1,28 @@ # Structures -Test \ No newline at end of file +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. From 7bdf299e6b1aa9905001686797abb245561ab14c Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Fri, 5 Aug 2022 17:43:58 +1000 Subject: [PATCH 4/9] Added fluence method for moving MLC aperture - Computes time-weighted fluence as the MLC moves from first to second position - Assumes MLC moves from mlcx1 to mlcx2 linearly --- src/Fluence.jl | 55 ++++++++++++++++++++++++++++++++------------------ 1 file changed, 35 insertions(+), 20 deletions(-) diff --git a/src/Fluence.jl b/src/Fluence.jl index b03a698..2fb410f 100644 --- a/src/Fluence.jl +++ b/src/Fluence.jl @@ -366,33 +366,48 @@ function fluence(bixel::AbstractBixel{T}, index::Int, mlcx) where T<:AbstractFlo overlap(position(bixel, 1), width(bixel, 1), mlcx[1, index], mlcx[2, index]) end -#=-- Moving Aperture fluences ------------------------------------------------------------------------------------------ +#--- Moving Aperture fluences ------------------------------------------------------------------------------------------ - Deprecated. - These functions compute the fluence from a leaf moving linearly from one position to another. -=# +function fluence(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator) where T<:AbstractFloat + j = max(1, locate(mlc, bixel[2])) + xL = bixel[1] - 0.5*width(bixel, 1) + xU = bixel[1] + 0.5*width(bixel, 1) -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) - - if isnan(t1) || isnan(t2) - A_triangle = 0. + x1 = mlcx1[1, j] + x2 = mlcx2[1, j] + if(x2 > x1) + ΨB = area_under_path(x1, x2, xL, xU) + elseif(x1 > x2) + ΨB = area_under_path(x2, x1, xL, xU) else - T = minmax(t1, 0., 1.) + minmax(t2, 0., 1.) - D = min(leafx2, xb2) - max(leafx1, xb1) + ΨB = min(1, max(zero(T), xU-x1)) + end - A_triangle = 0.5*T*D + x1 = mlcx1[2, j] + x2 = mlcx2[2, j] + if(x2 > x1) + ΨA = area_under_path(x1, x2, xL, xU) + elseif(x1 > x2) + ΨA = area_under_path(x2, x1, xL, xU) + else + ΨA = min(1, max(zero(T), xU-x1)) end - A_triangle + max(0., xb2 - leafx2) + max(zero(T), ΨB - ΨA) 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) +tfun(x, x1, x2) = min(1, max(0, (x-x1)/(x2-x1))) +function area_under_path(x1, x2, xL, xU) + + xB = max(x1, xL) + xA = min(x2, xU) + + tL = tfun(xB, x1, x2) + tU = tfun(xA, x1, x2) + + A1 = (tU + tL)/2*(xA - xB) + A2 = xU - xA + (A1 + A2)*(xU-xL) + end From 08b1ed663c695fc82518bfee5ddd94d1de2eae40 Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Fri, 5 Aug 2022 17:51:14 +1000 Subject: [PATCH 5/9] Added docstrings for new fluence methods --- src/Fluence.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/Fluence.jl b/src/Fluence.jl index 2fb410f..f3edfd7 100644 --- a/src/Fluence.jl +++ b/src/Fluence.jl @@ -368,7 +368,14 @@ end #--- Moving Aperture fluences ------------------------------------------------------------------------------------------ +""" + fluence(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator) + +From MLC leaf positions which move from `mlcx1` to `mlcx2`. +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(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator) where T<:AbstractFloat j = max(1, locate(mlc, bixel[2])) xL = bixel[1] - 0.5*width(bixel, 1) @@ -397,7 +404,20 @@ function fluence(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator max(zero(T), ΨB - ΨA) end +""" + tfun(x, x1, x2) + +Linear interpolation of `x` between `(x1, 0)` and `(x2, 1)` + +Used in `area_under_path`. +""" tfun(x, x1, x2) = min(1, max(0, (x-x1)/(x2-x1))) + +""" + area_under_path(x1, x2, xL, xU) + +Area of intersection between MLC leaf trajectory and bixel. +""" function area_under_path(x1, x2, xL, xU) xB = max(x1, xL) From e36419421afd20ff85b71388e385a3d677aa43f9 Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Tue, 9 Aug 2022 17:45:53 +1000 Subject: [PATCH 6/9] Tidied moving aperture fluence functions --- src/Fluence.jl | 75 ++++++++++++++++++++++++++------------------------ 1 file changed, 39 insertions(+), 36 deletions(-) diff --git a/src/Fluence.jl b/src/Fluence.jl index f3edfd7..91b9713 100644 --- a/src/Fluence.jl +++ b/src/Fluence.jl @@ -381,53 +381,56 @@ function fluence(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator xL = bixel[1] - 0.5*width(bixel, 1) xU = bixel[1] + 0.5*width(bixel, 1) - x1 = mlcx1[1, j] - x2 = mlcx2[1, j] - if(x2 > x1) - ΨB = area_under_path(x1, x2, xL, xU) - elseif(x1 > x2) - ΨB = area_under_path(x2, x1, xL, xU) - else - ΨB = min(1, max(zero(T), xU-x1)) - end - - x1 = mlcx1[2, j] - x2 = mlcx2[2, j] - if(x2 > x1) - ΨA = area_under_path(x1, x2, xL, xU) - elseif(x1 > x2) - ΨA = area_under_path(x2, x1, xL, xU) - else - ΨA = min(1, max(zero(T), xU-x1)) - end + ΨB = fluence_onesided(mlcx1[1, j], mlcx2[1, j], xL, xU) + ΨA = fluence_onesided(mlcx1[2, j], mlcx2[2, j], xL, xU) - max(zero(T), ΨB - ΨA) + max(zero(T), ΨB - ΨA) # This is needed for cases where xA < xB end """ - tfun(x, x1, x2) + fluence_onesided(xs, xf, xL, xU) -Linear interpolation of `x` between `(x1, 0)` and `(x2, 1)` +Compute fluence for 1 leaf position, assuming the other is infinitely far away -Used in `area_under_path`. +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. """ -tfun(x, x1, x2) = min(1, max(0, (x-x1)/(x2-x1))) +function fluence_onesided(xs, xf, xL, xU) -""" - area_under_path(x1, x2, xL, xU) + # If leaf positions are the same, compute static fluence + xs == xf && return (xU - minmax(xs, xL, xU))/(xU-xL) -Area of intersection between MLC leaf trajectory and bixel. -""" -function area_under_path(x1, x2, xL, xU) - xB = max(x1, xL) - xA = min(x2, xU) + # Otherwise compute moving fluence + + x1, x2 = min(xs, xf), max(xs, xf) # If xf < xs, swap - tL = tfun(xB, x1, x2) - tU = tfun(xA, x1, x2) + 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 = (tU + tL)/2*(xA - xB) - A2 = xU - xA - (A1 + A2)*(xU-xL) + 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 + +""" + minmax(x, l, u) + +Return `x` if `l<=x<=u`, `l` if `x Date: Wed, 10 Aug 2022 09:17:09 +1000 Subject: [PATCH 7/9] Moved minmax function to misc.jl - Also removed minmax function from interpolation.jl --- src/Fluence.jl | 7 ------- src/utils/interpolation.jl | 7 ------- src/utils/misc.jl | 7 +++++++ 3 files changed, 7 insertions(+), 14 deletions(-) diff --git a/src/Fluence.jl b/src/Fluence.jl index 91b9713..3ba21b1 100644 --- a/src/Fluence.jl +++ b/src/Fluence.jl @@ -419,13 +419,6 @@ function fluence_onesided(xs, xf, xL, xU) end -""" - minmax(x, l, u) - -Return `x` if `l<=x<=u`, `l` if `x Date: Wed, 10 Aug 2022 10:25:16 +1000 Subject: [PATCH 8/9] Added fluence methods for moving leaf apertures - Added aliases for MLC Aperture and Sequence --- src/BeamLimitingDevices.jl | 5 +++++ src/Fluence.jl | 40 ++++++++++++++++++++++++++++++-------- 2 files changed, 37 insertions(+), 8 deletions(-) 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 3ba21b1..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,27 +362,52 @@ 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 ------------------------------------------------------------------------------------------ """ - fluence(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator) + 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 + +""" + fluence(bixel::AbstractBixel{T}, index::Int, mlcx::AbstractMLCSequence) + +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 + +""" + fluence_from_moving_aperture(bixel::AbstractBixel{T}, mlcx1, mlcx2) From MLC leaf positions which move from `mlcx1` to `mlcx2`. 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(bixel::AbstractBixel{T}, mlcx1, mlcx2, mlc::MultiLeafCollimator) where T<:AbstractFloat - j = max(1, locate(mlc, bixel[2])) +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, j], mlcx2[1, j], xL, xU) - ΨA = fluence_onesided(mlcx1[2, j], mlcx2[2, j], xL, xU) + Ψ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 @@ -401,7 +426,6 @@ 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 From 5ba2d438a195ee4ecb5a0503557d3dffb71e35c5 Mon Sep 17 00:00:00 2001 From: lmejn <33491793+lmejn@users.noreply.github.com> Date: Wed, 10 Aug 2022 10:42:09 +1000 Subject: [PATCH 9/9] Updated version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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"