diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..700707c --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" # Location of package manifests + schedule: + interval: "weekly" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..65b8b7d --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,41 @@ +name: CI +on: + push: + branches: + - master + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read + strategy: + fail-fast: false + matrix: + version: + - '1.10' + - '1.6' + - 'pre' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..cba9134 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,16 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..0cd3114 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,31 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: + inputs: + lookback: + default: "3" +permissions: + actions: read + checks: read + contents: write + deployments: read + issues: read + discussions: read + packages: read + pages: read + pull-requests: read + repository-projects: read + security-events: read + statuses: read +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..21caefb --- /dev/null +++ b/.gitignore @@ -0,0 +1,24 @@ +# Files generated by invoking Julia with --code-coverage +*.jl.cov +*.jl.*.cov + +# Files generated by invoking Julia with --track-allocation +*.jl.mem + +# System-specific files and directories generated by the BinaryProvider and BinDeps packages +# They contain absolute paths specific to the host computer, and so should not be committed +deps/deps.jl +deps/build.log +deps/downloads/ +deps/usr/ +deps/src/ + +# Build artifacts for creating documentation generated by the Documenter package +docs/build/ +docs/site/ + +# File generated by Pkg, the package manager, based on a corresponding Project.toml +# It records a fixed state of all packages used by the project. As such, it should not be +# committed for packages, but should be committed for applications that require a static +# environment. +Manifest.toml \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..dd3f580 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 Roy Sela Arazi + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..d7e1f74 --- /dev/null +++ b/Project.toml @@ -0,0 +1,22 @@ +name = "PTMC" +uuid = "3e88bec1-2450-477a-bc2d-163824aebffa" +authors = ["Roy Sela Arazi"] +version = "0.1.0" + +[deps] +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MAT = "23992714-dd62-5051-b70f-ba57cb901cac" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[extras] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Combinatorics", "BenchmarkTools", "Test"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..3311244 --- /dev/null +++ b/README.md @@ -0,0 +1,12 @@ +# PTMC.jl + +[![Build Status](https://github.com/ryarazi/PTMC.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/ryarazi/PTMC.jl/actions/workflows/CI.yml?query=branch%3Amaster) + +PTMC.jl (Photon Transport Monte Carlo) a Julia code based on the Monte Carlo modeling of photon transport in multi-layered tissues (MCML) algorithm. The original MCML algorithm was developed by Wang et al. in 1995 with a C code implementation available at [https://omlc.org/software/mc/](https://omlc.org/software/mc/). + +The MCML algorithm is a tool for simulating light transport in biological tissues. It is based on the Monte Carlo method, which simulates the path of photons as they travel through the tissue, and is used to calculate the spatial distribution of light fluence, absorption, and scattering in the tissue. + +## License + +This package is released under the MIT License. Please see the [LICENSE](LICENSE) file for more information. +``` \ No newline at end of file diff --git a/src/PTMC.jl b/src/PTMC.jl new file mode 100644 index 0000000..e22218b --- /dev/null +++ b/src/PTMC.jl @@ -0,0 +1,37 @@ +module PTMC + using LinearAlgebra + using Random + using Rotations + using StaticArrays + using MAT + using Interpolations + using Plots + + import Plots: plot + import Base: maximum, minimum, extrema, isapprox, ∈, ∩ + + export Laser, Box, Solid, Turbid, Refractive, Scene, HistoryAggregator + export simulate_particle + + include("linear.jl") + + #geometry + include("geometry/types.jl") + include("geometry/ray.jl") + include("geometry/plane.jl") + include("geometry/sphere.jl") + include("geometry/box.jl") + include("geometry/intersection.jl") + + include("particle.jl") + include("aggregator.jl") + include("laser.jl") + + include("material.jl") + include("solid.jl") + include("scene.jl") + + include("driver.jl") + + include("plots.jl") +end \ No newline at end of file diff --git a/src/aggregator.jl b/src/aggregator.jl new file mode 100644 index 0000000..4660e3c --- /dev/null +++ b/src/aggregator.jl @@ -0,0 +1,21 @@ +""" + abstract type Aggregator + +Abstract type for aggregators of PTMC simulation results. +Concrete types must dispatch the `create_aggregation` and `update_aggregation` functions. +See `HistoryAggregator` for an example. +""" +abstract type Aggregator end + +""" + struct HistoryAggregator <: Aggregator + +An aggregator that stores the history of each photon in the simulation. +""" +struct HistoryAggregator <: Aggregator + history::Vector{Vector{Photon}} + HistoryAggregator() = new([]) +end + +create_aggregation(agg::HistoryAggregator, pstart::Photon) = push!(agg.history, [pstart]) +update_aggregation(agg::HistoryAggregator, pstart::Photon, pend::Photon) = push!(agg.history[end], pend) \ No newline at end of file diff --git a/src/driver.jl b/src/driver.jl new file mode 100644 index 0000000..1d8c824 --- /dev/null +++ b/src/driver.jl @@ -0,0 +1,64 @@ + +""" + simulate_particle(source::Source, scene::Scene, aggregators = nothing; roulette_threshold = 0.001, roulette_chance = 0.1, δs = 1.e-8) + +Simulate a single particle through the scene. The particle is created by sampling the source and then +propagated through the scene until it exits the system or its weight is zero. The particle is then returned. + +# Arguments +- `source::Source`: The source from which the particle is sampled. +- `scene::Scene`: The scene through which the particle is propagated. +- `aggregators`: A collection of aggregators which collect data during the simulation. +- `roulette_threshold::Float64`: The threshold at which the particle weight is considered small enough to be terminated. +- `roulette_chance::Float64`: The probability that a particle is terminated when its weight is below the threshold. +- `δs::Float64`: The distance to move the particle past the boundary before it is reflected or refracted. + +# Returns +- `Photon`: The particle after it has been propagated through the scene. +""" +function simulate_particle(source::Source, scene::Scene, aggregators = nothing; roulette_threshold = 0.001, roulette_chance = 0.1, δs = 1.e-8) + p = sample_source(source) + solid = find_container(p, scene) + + if !isnothing(aggregators) for agg in aggregators create_aggregation(agg, p) end end + + while p.weight != 0. && solid != Outside() + s = randexp() # HOP + + #check boundaries and move between solids + while s != 0. + pstart = p + if solid == Outside() break end #photon is out of the systems + + boundary_intersection, normal = exit_solid(p, solid) + d_boundary = sqrt.(sum((p.r - boundary_intersection).^2)) #calculate norm, faster then norm function + s_boundary = d_boundary * solid.material.μₛ + + if s_boundary > s + d = s/solid.material.μₛ + p = advance(p, d) + p = spin(p, solid.material.g) # SPIN + p = drop(p, d, solid.material.μₐ) # DROP + else + p = Photon(boundary_intersection, p.u, p.weight) #move photon to the surface boundary + solid_next = find_container(advance(p, δs), scene) + isreflect = rand() < unpolarized_reflectance(p, normal, solid.material.n, solid_next.material.n) #decide if reflect or refract + p = isreflect ? reflect(p, normal) : refract(p, normal, solid.material.n, solid_next.material.n) + p = advance(p, δs) #make sure the particle is just past the boundary + p = drop(p, d_boundary, solid.material.μₐ) # DROP + solid = isreflect ? solid : solid_next #decide on next solid + # @assert p∈solid.geometry + end + s = max(s - s_boundary, 0.) + + if !isnothing(aggregators) for agg in aggregators update_aggregation(agg, pstart, p) end end + + #Roulette Method (see https://fa.bianp.net/blog/2022/russian-roulette/) + if p.weight < roulette_threshold + p = roulette(p, roulette_chance) + end + end + end + + return p +end \ No newline at end of file diff --git a/src/geometry/box.jl b/src/geometry/box.jl new file mode 100644 index 0000000..07f7b62 --- /dev/null +++ b/src/geometry/box.jl @@ -0,0 +1,23 @@ +""" + struct Box{D,T} <: Geometry{D,T} + +A box geometry object in the scene. +""" +struct Box{D,T} <: Geometry{D,T} + vmin::SVector{D,T} + vmax::SVector{D,T} + + function Box{D,T}(vmin, vmax) where {D,T} + @assert all(vmin .<= vmax) + new(vmin, vmax) + end +end + +Box(vmin::NTuple{D,T}, vmax::NTuple{D,T}) where {D,T} = Box{D,T}(vmin, vmax) +Box(vmin::Union{AbstractVector,Tuple}, vmax::Union{AbstractVector,Tuple}) = Box(Tuple(vmin), Tuple(vmax)) + +#overload Base methods +minimum(b::Box) = b.vmin +maximum(b::Box) = b.vmax +extrema(b::Box) = b.vmin, b.vmax +∈(x::AbstractVector, b::Box) = all(@. ((b.vmin < x) || isapprox(b.vmin, x; atol=100*eps(Float64))) && ((x < b.vmax) || isapprox(b.vmax, x; atol=100*eps(Float64)))) \ No newline at end of file diff --git a/src/geometry/intersection.jl b/src/geometry/intersection.jl new file mode 100644 index 0000000..1f4897c --- /dev/null +++ b/src/geometry/intersection.jl @@ -0,0 +1,58 @@ +""" + ∩ + +Compute the intersection between a ray and a geometry. The intersection +can br `nothing` if the ray does not intersect the geometry, a single +point or a tuple of two points. +""" +∩ + +function ∩(ray::Ray, plane::Plane) + if ray.direction ⟂ plane.normal + return ray.origin∈plane ? ray : nothing + end + + # solve r = origin + α * direction & (r - p₀) * normal = 0 + α = dot(plane.p₀ - ray.origin, plane.normal) / dot(ray.direction, plane.normal) + if α<0. return nothing end #ray defined only for α≥0 + return ray(α) +end + +function ∩(ray::Ray, sphere::Sphere) + # solve quadratic equation r = origin + α * direction & (r - p₀)^2 = R^2 + a, b, c = dot(ray.direction, ray.direction), 2*dot(ray.direction, ray.origin-sphere.p₀), dot(ray.origin-sphere.p₀,ray.origin-sphere.p₀)-sphere.R^2 + Δ² = b^2 - 4*a*c #squared descriminant + if Δ² < 0 + return nothing + elseif abs(Δ²) < 1000*eps(eltype(Δ²)) + α = -b/2a + return α>=0. ? ray(α) : nothing + else + Δ = sqrt(Δ²) + α₁, α₂ = (-b + Δ)/2a, (-b - Δ)/2a + if α₂>=0 + return ray(α₂), ray(α₁) + elseif α₁>=0 + return ray(α₁) + else + return nothing + end + end +end + +function ∩(ray::Ray, b::Box) + #solve r = origin + α * direction & vmin ≤ r ≤ vmax + α₁, α₂ = (b.vmin - ray.origin) ./ ray.direction, (b.vmax - ray.origin) ./ ray.direction + + #see https://www.cs.cornell.edu/courses/cs4620/2013fa/lectures/03raytracing1.pdf + α_enter = min.(α₁, α₂) + α_exit = max.(α₁, α₂) + + α_interval = maximum(α_enter), minimum(α_exit) + + if α_interval[1] > α_interval[2] return nothing + elseif α_interval[1] < 0. + if α_interval[2] < 0. return nothing + else return ray(α_interval[2]) end + else return ray(α_interval[1]), ray(α_interval[2]) end +end \ No newline at end of file diff --git a/src/geometry/plane.jl b/src/geometry/plane.jl new file mode 100644 index 0000000..3fe76eb --- /dev/null +++ b/src/geometry/plane.jl @@ -0,0 +1,15 @@ +""" + struct Plane{D,T} <: Geometry{D,T} + +A plane geometry object in the scene. +""" +struct Plane{D,T} <: Geometry{D,T} + "A point on the plane" + p₀::SVector{D,T} + "The normal to the plane" + normal::SVector{D,T} +end + +Plane(p₀::NTuple{D,T}, normal::NTuple{D,T}) where {D,T} = Plane{D,T}(p₀, normal) + +∈(x::AbstractVector, plane::Plane) = (x⋅plane.normal ≈ plane.p₀⋅plane.normal) \ No newline at end of file diff --git a/src/geometry/ray.jl b/src/geometry/ray.jl new file mode 100644 index 0000000..7b69553 --- /dev/null +++ b/src/geometry/ray.jl @@ -0,0 +1,19 @@ +""" + struct Ray{D,T} <: Geometry{D,T} + +A ray geometry object in the scene. +""" +struct Ray{D,T} <: Geometry{D,T} + "The origin of the ray" + origin::SVector{D,T} + "The direction vector of the ray" + direction::SVector{D,T} +end + +Ray(origin::NTuple{D,T}, direction::NTuple{D,T}) where {D,T} = Ray{D,T}(origin, direction) +Ray(origin::Union{AbstractVector,Tuple}, direction::Union{AbstractVector,Tuple}) = Ray(Tuple(origin), Tuple(direction)) + +@inline function (r::Ray)(α) + # @assert α>=0. + r.origin + α*r.direction +end \ No newline at end of file diff --git a/src/geometry/sphere.jl b/src/geometry/sphere.jl new file mode 100644 index 0000000..6d75a93 --- /dev/null +++ b/src/geometry/sphere.jl @@ -0,0 +1,16 @@ +""" + struct Sphere{D,T} <: Geometry{D,T} + +A sphere geometry object in the scene. +""" +struct Sphere{D,T} <: Geometry{D,T} + "The center of the sphere" + p₀::SVector{D,T} + "The radius of the sphere" + R::T +end + +Sphere(p₀::NTuple{D,T}, R::T) where {D,T} = Sphere{D,T}(p₀, R) +Sphere(p₀::Union{AbstractVector,Tuple}, R) = Sphere(Tuple(p₀), R) + +∈(x::AbstractVector, sphere::Sphere) = dot(x-sphere.p₀,x-sphere.p₀) < sphere.R^2 || isapprox(dot(x-sphere.p₀,x-sphere.p₀), sphere.R^2; atol=100*eps(Float64)) diff --git a/src/geometry/types.jl b/src/geometry/types.jl new file mode 100644 index 0000000..bc836fa --- /dev/null +++ b/src/geometry/types.jl @@ -0,0 +1,6 @@ +""" + Geometry{D,T} + +Abstract type for all geometry objects in the scene. +""" +abstract type Geometry{D,T} end \ No newline at end of file diff --git a/src/laser.jl b/src/laser.jl new file mode 100644 index 0000000..3e387bf --- /dev/null +++ b/src/laser.jl @@ -0,0 +1,64 @@ +""" + abstract type Source + +An abstract type representing a generic source of particles. +""" +abstract type Source end + +""" + struct Laser <: Source + +A type representing a laser source. + +# Fields +- `σ::Float64`: The standard deviation of the laser beam. +- `direction::SVector{3,Float64}`: The direction vector of the laser beam. +- `center::SVector{3,Float64}`: The center position of the laser beam. + +# Constructor +- `Laser(σ, direction::SVector{3,Float64}, center::SVector{3,Float64})`: Creates a new `Laser` instance with the given standard deviation, direction, and center. The direction vector is normalized. +""" +struct Laser <: Source + σ::Float64 + direction::SVector{3,Float64} + center::SVector{3,Float64} + + Laser(σ, direction::SVector{3,Float64}, center::SVector{3,Float64}) = new(σ, normalize(direction), center) +end + +""" + Laser(σ, origin::Union{AbstractVector,Tuple}, direction::Union{AbstractVector,Tuple}) -> Laser + +Alternative constructor for `Laser` that accepts the origin and direction as vectors or tuples. + +# Arguments +- `σ::Float64`: The standard deviation of the laser beam. +- `origin::Union{AbstractVector,Tuple}`: The origin position of the laser beam. +- `direction::Union{AbstractVector,Tuple}`: The direction vector of the laser beam. + +# Returns +A new `Laser` instance with the given standard deviation, origin, and direction. +""" +Laser(σ, origin::Union{AbstractVector,Tuple}, direction::Union{AbstractVector,Tuple}) = Laser(σ, SVector(Tuple(origin)), SVector(Tuple(direction))) + +""" + sample_source(laser::Laser) + +Samples a particle from the given laser source. + +# Arguments +- `laser::Laser`: The laser source from which to sample a particle. + +# Returns +A particle sampled from the laser source. +""" +function sample_source(laser::Laser) + radii = √2 * laser.σ * sqrt(randexp()) + θ = 2pi*rand() + r = @SVector [radii*cos(θ), radii*sin(θ), 0.] + + z_axis = @SVector [0., 0., 1.] + rot = rotation_a2b(z_axis, laser.direction) + + return Photon(rot*r .+ laser.center, laser.direction) +end \ No newline at end of file diff --git a/src/linear.jl b/src/linear.jl new file mode 100644 index 0000000..9b6eecb --- /dev/null +++ b/src/linear.jl @@ -0,0 +1,50 @@ + +""" + ⟂(x, y) + +Check if two vectors are orthogonal. +""" +⟂(x, y) = abs(x⋅y)<1000*eps(eltype(x)) + +""" + rotate_vector(u::SVector{3, Float64}, cosθ, ψ) + +Rotate a vector `u` by rotaion angles `cosθ` and `ψ`. +""" +function rotate_vector(u::SVector{3, Float64}, cosθ, ψ) + #@assert norm(u) ≈ 1. + sinθ = sqrt(1-cosθ^2) + sinψ, cosψ = sincos(ψ) + + uˣ,uʸ,uᶻ = u + + if abs(abs(uᶻ)-1) < 1000*eps(eltype(uᶻ)) + u_rotated = @SVector [sinθ*cosψ, sinθ*sinψ, sign(uᶻ)*cosθ] + else + u_rotated = @SVector [sinθ * (uˣ*uᶻ*cosψ - uʸ*sinψ)/sqrt(1-uᶻ^2) + uˣ*cosθ, + sinθ * (uʸ*uᶻ*cosψ + uˣ*sinψ)/sqrt(1-uᶻ^2) + uʸ*cosθ, + -sinθ * cosψ * sqrt(1-uᶻ^2) + uᶻ*cosθ] + end + return u_rotated +end + +""" + rotation_a2b(a::SVector{3, Float64}, b::SVector{3, Float64}) + +Calculate the rotation matrix that rotates vector `a` to vector `b`. +""" +function rotation_a2b(a::SVector{3, Float64}, b::SVector{3, Float64}) + anorm, bnorm = norm(a), norm(b) + k = a × b + cosθ = (a⋅b) / anorm / bnorm + if abs(cosθ) ≉ 1 + v_cross = @SMatrix [0. -k[3] k[2] + k[3] 0. -k[1] + -k[2] k[1] 0.] + + return SMatrix{3,3,Float64}(I) + v_cross + (v_cross*v_cross) ./ (1+cosθ) #https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d + else # a is parallel to b + if dot(a,b)>0 return SMatrix{3,3,Float64}(I) + else return -SMatrix{3,3,Float64}(I) end + end +end \ No newline at end of file diff --git a/src/material.jl b/src/material.jl new file mode 100644 index 0000000..40debee --- /dev/null +++ b/src/material.jl @@ -0,0 +1,92 @@ +""" + struct Turbid + +A struct to represent a turbid material +""" +struct Turbid + "material scattering coefficient" + μₛ::Float64 + + "material absorption coefficient" + μₐ::Float64 + + "refractive index" + n::Float64 + + "Henyey-Greenstein coefficient" + g::Float64 + + "total interaction coefficient" + μₜ::Float64 + + Turbid(;μₛ=0.,μₐ=0.,n=1.,g=0.) = new(μₛ, μₐ, n, g, μₛ+μₐ) +end + +""" + mean_free_path(mat::Turbid) + +Return the mean free path of a turbid material +""" +mean_free_path(mat::Turbid) = mat.μₜ==0. ? Inf : 1/mat.μₜ +Base.show(io::IO, t::Turbid) = print(io, "Turbid(μₛ=$(t.μₛ), μₐ=$(t.μₐ), n=$(t.n), g=$(t.g))") + +""" + Refractive(;n=1.) + +Return a refractive material with refractive index `n` +""" +Refractive(;n=1.) = Turbid(n=n) + +""" + Tissue(λₙₘ, B, S, W, M, F, μₛ_reduced_500nm, f, b_mie, g) + +Return a turbid material representing a tissue. The model is based on +"Jacques, S. L. (2013). Optical properties of biological tissues: a review". The +mathematical model is described shortly in https://omlc.org/news/feb15/generic_optics/index.html. +The code is based on `makeTissueList.m` taken from https://omlc.org/software/mc/mcxyz/ + +# Arguments +- `λₙₘ::Float64`: wavelength in nm +- `B::Float64`: blood volume fraction +- `S::Float64`: blood oxygen saturation +- `W::Float64`: water volume fraction +- `M::Float64`: melanosome volume fraction +- `F::Float64`: fat volume fraction +- `μₛ_reduced_500nm::Float64`: reduced scattering coefficient at 500nm +- `f::Float64`: scattering power law +- `b_mie::Float64`: scattering power law +- `g::Float64`: Henyey-Greenstein coefficient + +# Returns +- `Turbid`: a turbid material representing a tissue +""" + +function Tissue(λₙₘ, B, S, W, M, F, μₛ_reduced_500nm, f, b_mie, g) + data = Dict(label=>vec(mat) for (label,mat) in matread("data/spectralLIB.mat")) + @assert minimum(data["nmLIB"]) <= λₙₘ <= maximum(data["nmLIB"]) + + μₐ_HGb_oxygen = linear_interpolation(data["nmLIB"], data["muaoxy"])(λₙₘ) + μₐ_HGb_deoxygen = linear_interpolation(data["nmLIB"], data["muadeoxy"])(λₙₘ) + μₐ_water = linear_interpolation(data["nmLIB"], data["muawater"])(λₙₘ) + μₐ_melanosome = linear_interpolation(data["nmLIB"], data["muamel"])(λₙₘ) + μₐ_fat = linear_interpolation(data["nmLIB"], data["muafat"])(λₙₘ) + + μₐ = B*S*μₐ_HGb_oxygen + B*(1-S)*μₐ_HGb_deoxygen + W*μₐ_water + M*μₐ_melanosome + F*μₐ_fat + μₛ_reduced = μₛ_reduced_500nm * (f*(λₙₘ/500)^(-4) + (1-f)*(λₙₘ/500)^(-b_mie)) + μₛ = μₛ_reduced / (1-g) + + ndry = 1.514 #see "Jacques, S. L. (2013). Optical properties of biological tissues: a review" + nwater = 1.33 + n = ndry - W*(ndry - nwater) + + return Turbid(μₛ=μₛ,μₐ=μₐ,n=n,g=g) +end + +#based on makeTissueList.m taken from https://omlc.org/software/mc/mcxyz/ +#(see https://omlc.org/software/mc/mcxyz/makeTissueList.m) +Blood(λₙₘ) = Tissue(λₙₘ, 1., 0.75, 0.95, 0., 0., 10., 0., 1.0, 0.9) +Dermis(λₙₘ) = Tissue(λₙₘ, 0.002, 0.67, 0.65, 0., 0., 42.4, 0.62, 1.0, 0.9) +Epidermis(λₙₘ) = Tissue(λₙₘ, 0., 0.75, 0.75, 0.03, 0., 40., 0., 1.0, 0.9) +Skull(λₙₘ) = Tissue(λₙₘ, 0.0005, 0.75, 0.35, 0., 0., 30., 0., 1.0, 0.9) +GrayMatter(λₙₘ) = Tissue(λₙₘ, 0.01, 0.75, 0.75, 0., 0., 20., 0.2, 1.0, 0.9) +WhiteMatter(λₙₘ) = Tissue(λₙₘ, 0.01, 0.75, 0.75, 0., 0., 20., 0.2, 1.0, 0.9) \ No newline at end of file diff --git a/src/particle.jl b/src/particle.jl new file mode 100644 index 0000000..298bf18 --- /dev/null +++ b/src/particle.jl @@ -0,0 +1,162 @@ +""" + struct Photon + +A representation of a photon in a medium with a position `r`, a direction `u` and a weight `weight`. +""" +struct Photon + "photon position" + r::SVector{3,Float64} + "photon direction" + u::SVector{3,Float64} + "photon weight" + weight::Float64 + + Photon(r::SVector{3,Float64}, u::SVector{3,Float64}, weight=1.) = new(r, u, weight) +end + +Photon(r::Union{AbstractVector,Tuple}, u::Union{AbstractVector,Tuple}, weight=1.) = Photon(SVector(Tuple(r)), SVector(Tuple(u)), weight) + +""" + advance(p::Photon, distance::Float64) + +Advance a photon `p` by a distance `distance` in the direction of the photon `p.u`. + +# Arguments +- `p::Photon`: the photon to advance +- `distance::Float64`: the distance to advance the photon + +# Returns +- a new photon advanced by `distance` in the direction of `p.u`. +""" +@inline function advance(p::Photon, distance::Float64) + Photon(p.r + distance * p.u, p.u, p.weight) +end + +""" + reflect(p::Photon, normal) + +Reflect a photon `p` on a surface with normal `normal` using the law of reflection. + +# Arguments +- `p::Photon`: the photon to reflect +- `normal::SVector{3,Float64}`: the normal of the surface + +# Returns +- a new photon reflected on the surface +""" +@inline function reflect(p::Photon, normal) + cosθ1 = -dot(normal, p.u) + u_reflect = p.u + 2cosθ1*normal + Photon(p.r, u_reflect, p.weight) +end + +@inline function snell_angles(p::Photon, normal, n1::Float64, n2::Float64) + cosθ1 = -dot(normal, p.u) + sqrt_body = 1 - (n1/n2)^2 * (1-cosθ1^2) + if sqrt_body < 0. #pass critical angle + return cosθ1, nothing + else + cosθ2 = sqrt(sqrt_body) + return cosθ1, cosθ2 + end +end + +""" + refract(p::Photon, normal, n1::Float64, n2::Float64) + +Refract a photon `p` on a surface with normal `normal` where the photon +is coming from a medium with refractive index `n1` to +a medium with refractive index `n2` using Snell's law. + +# Arguments +- `p::Photon`: the photon to refract +- `normal::SVector{3,Float64}`: the normal of the surface +- `n1::Float64`: the refractive index of the medium where the photon is coming from +- `n2::Float64`: the refractive index of the medium where the photon is going to + +# Returns +- a new photon refracted on the surface +""" +@inline function refract(p::Photon, normal, n1::Float64, n2::Float64) + cosθ1, cosθ2 = snell_angles(p, normal, n1, n2) + if isnothing(cosθ2) return nothing end #pass critical angle + u_refract = (n1/n2)*p.u + ((n1/n2)cosθ1 - cosθ2)*normal + Photon(p.r, u_refract, p.weight) +end + +""" + unpolarized_reflectance(p::Photon, normal, n1::Float64, n2::Float64) + +Compute the unpolarized reflectance (Fresnel reflectance) of a photon `p` on a surface with normal `normal` where the photon +is coming from a medium with refractive index `n1` to a medium with refractive index `n2`. + +# Arguments +- `p::Photon`: the photon to reflect +- `normal::SVector{3,Float64}`: the normal of the surface +- `n1::Float64`: the refractive index of the medium where the photon is coming from +- `n2::Float64`: the refractive index of the medium where the photon is going to + +# Returns +- the Fresnel reflectance of the photon +""" +@inline function unpolarized_reflectance(p::Photon, normal, n1::Float64, n2::Float64) + cosθ1, cosθ2 = snell_angles(p, normal, n1, n2) + if isnothing(cosθ2) return 1. end #pass critical angle = always reflect + R = 0.5 * ( ((n1*cosθ1-n2*cosθ2)/(n1*cosθ1+n2*cosθ2))^2 + ((n2*cosθ1-n1*cosθ2)/(n2*cosθ1+n1*cosθ2))^2 ) + R +end + +""" + spin(p::Photon, g) + +Spin a photon `p` with a Henyey-Greenstein coefficient `g` in a new direction. + +# Arguments +- `p::Photon`: the photon to spin +- `g::Float64`: the Henyey-Greenstein coefficient + +# Returns +- a new photon with a new random direction +""" +function spin(p::Photon, g) + if abs(g) < 1000*eps(eltype(g)) + cosθ = 2*rand() - 1 + else + cosθ = (1 + g^2 - (1-g^2)^2/(1-g+2*g*rand())^2) / 2g + end + ψ = 2pi*rand() + u_rotated = rotate_vector(p.u, cosθ, ψ) + return Photon(p.r, u_rotated, p.weight) +end + +""" + drop(p::Photon, distnace, μₐ) + +Decrease the weight of a photon `p` by a distance `distance` as it is +partially absorbed in the medium with absorption coefficient `μₐ`. + +# Arguments +- `p::Photon`: the moving photon +- `distance::Float64`: the distance the photon has traveled +- `μₐ::Float64`: the absorption coefficient of the medium + +# Returns +- a new photon with a decreased weight by a factor of `exp(-distance*μₐ)` +""" +@inline function drop(p::Photon, distnace, μₐ) + return Photon(p.r, p.u, p.weight * exp(-distnace*μₐ)) +end + +""" + roulette(p::Photon, chance) + +Apply Russian roulette to a photon `p` with a chance `chance` to kill the photon. + +# Arguments +- `p::Photon`: the photon to apply Russian roulette +- `chance::Float64`: the chance to kill the photon + +# Returns +- a new photon with a increased weight by a factor of `1/chance` if the photon survives, otherwise a photon with weight 0. +""" +@inline roulette(p::Photon, chance) = rand() "x" + yguide --> "y" + zguide --> "z" + + @series begin + label := :none + linewidth := 0.7 + linecolor := :black + xe,ye,ze + end + + @series begin + fa := 0.3 + linewidth := 0. + seriestype := :mesh3d + connections := connections + xp,yp,zp + end +end + +@recipe function f(s::Solid) + label --> s.name + s.geometry +end + +@recipe function f(scene::Scene) + for i in length(scene.solids)-1:-1:1 #first solid is Outside + @series begin scene.solids[i] end + end +end + +@recipe function f(history::Vector{Photon}; slice = nothing) + marker --> :circle + markersize --> 1. + arrow --> (:closed, 20.0) + label --> :none + + @assert slice∈(nothing, "xy", "yx", "yz", "zy", "xz", "zx") + + r = stack(p.r for p in history) + if isnothing(slice) + return r[1,:], r[2,:], r[3,:] + elseif slice=="xy" || slice=="yx" + return r[1,:], r[2,:] + elseif slice=="yz" || slice=="zy" + return r[2,:], r[3,:] + elseif slice=="zx" || slice=="xz" + return r[1,:], r[3,:] + end +end \ No newline at end of file diff --git a/src/scene.jl b/src/scene.jl new file mode 100644 index 0000000..af0b69b --- /dev/null +++ b/src/scene.jl @@ -0,0 +1,46 @@ +""" + struct Scene + +A collection of solids in the scene. The order of the solids is important, +as the earliest solids "hide" the later ones. +""" +struct Scene + "The ordered collections of solids in the scene, ordered by priority" + solids::Vector{Solid} + Scene() = new([Outside()]) +end + +∈(p::Photon, s::Scene) = any(p.r∈sol.geometry for sol in s.solids) + +""" + add_solid!(s::Scene, solid::Solid) + +Add a solid to the scene. +""" +add_solid!(s::Scene, solid::Solid) = pushfirst!(s.solids, solid) + +""" + add_solid!(s::Scene, geometry, material) + +Add a solid to the scene with the given geometry and material. +""" +add_solid!(s::Scene, geometry, material) = add_solid!(s, Solid(geometry, material)) + +""" + Scene(solid_collection...) + +Create a scene from a collection of solids. +""" +function Scene(solid_collection) + scene = Scene() + pushfirst!(scene.solids, solid_collection...) + scene +end + +""" + find_container(p::Photon, s::Scene) + +Find the solid that contains the photon. The solid is the first solid in the +scene's list that contains the photon. +""" +find_container(p::Photon, s::Scene) = first(sol for sol in s.solids if p.r∈sol.geometry) \ No newline at end of file diff --git a/src/solid.jl b/src/solid.jl new file mode 100644 index 0000000..a2e01cd --- /dev/null +++ b/src/solid.jl @@ -0,0 +1,44 @@ +""" + Solid + +A solid is a geometry with a material. It is the basic building block of a scene. +""" +struct Solid + geometry::Geometry{3,Float64} + material::Turbid + name::String + Solid(geometry, material; name="") = new(geometry, material, name) +end + +""" + Outside() + +Return a solid that represents the outside of the scene. +""" +Outside() = Solid(Box((-Inf,-Inf,-Inf), (Inf,Inf,Inf)), Turbid(), name="Outside") + +∈(p::Photon, s::Geometry{3,Float64}) = p.r∈s +∈(p::Photon, s::Solid) = p∈s.geometry + +""" + exit_solid(p::Photon, s::Solid) + +Return the point where the photon exits the solid and the normal to the surface at that point. + +# Arguments +- `p::Photon`: the photon moving through the solid +- `s::Solid`: the solid the photon is moving through + +# Returns +- `(boundary_intersection, normal)`: the point where the photon exits the solid and the normal to the surface at that point +""" +exit_solid(p::Photon, s::Solid) = exit_solid(p, s.geometry) + +function exit_solid(p::Photon, b::Box{3,Float64}) + # @assert p∈b + ray = Ray(p.r, p.u) + boundary_intersection::SVector{3, Float64} = ray∩b + normal = (abs.(boundary_intersection-b.vmin) .< 1e5*eps(Float64)) - (abs.(boundary_intersection-b.vmax) .< 1e5*eps(Float64)) + # @assert sum(abs.(normal))==1 + return boundary_intersection, normal +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..01b9428 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,226 @@ +using PTMC +import PTMC: Ray, Box, Sphere, Photon +import PTMC: ⟂, sample_source, reflect, unpolarized_reflectance, refract, spin, roulette, rotate_vector, advance, rotation_a2b, find_container, exit_solid +using Test +using StaticArrays +using Statistics +using Random +using LinearAlgebra + +@testset "PTMC.jl" begin + @testset verbose = true "Unit Tests" begin + @testset "linear rotations" begin + u = normalize(@SVector randn(3)) + cosθ = 2*rand() - 1 + ψ = 2pi*rand() + + @test norm(u) ≈ norm(rotate_vector(u, cosθ, ψ)) + @test u ≈ -rotate_vector(u, -1, ψ) + @test dot(rotate_vector(u, cosθ, ψ), u) ≈ cosθ + + v = normalize(@SVector randn(3)) + rot = rotation_a2b(u,v) + @test all(rot*u .≈ v) + @test all(rot*cross(u,v) .≈ cross(u,v)) #axis of rotation stays in the same place + end; + + @testset "geometric calculations" begin + ray = Ray(Tuple(rand(3)), Tuple(normalize(randn(3)))) + + b = Box((0.,0.,0.), (1.,1.,1.)) + + @test rand(3)∈b + @test (rand(3).+1∉b) && (rand(3).-1∉b) + + p = @SVector rand(3) + direction = zeros(3) + direction[rand((1,2,3))] = rand((-1,1)) + ray = Ray(p, direction) + + inter = ray ∩ b + @test !(typeof(inter)<:Tuple) + + normal = isapprox.(inter, b.vmin; atol=eps(Float64)) - isapprox.(inter, b.vmax; atol=eps(Float64)) + @test all(normal .≈ -direction) + + R = 2. + 3*rand() + sphere = Sphere(zeros(3), R) + ray = Ray((@SVector rand(3)), (@SVector randn(3))) + @test ray.origin∈sphere + @test !isnothing(ray∩sphere) && typeof(ray∩sphere)==SVector{3,Float64} && norm(ray∩sphere) ≈ R + + direction = zeros(3) + direction[rand((1,2,3))] = rand((-1,1)) + ray_axis = Ray(zeros(3), direction) + @test ray_axis∩sphere ≈ R .* direction + + p_far = 10 .+ (@SVector rand(3)) + ray = Ray(p_far, -p_far) + @test !isnothing(ray∩sphere) && length(ray∩sphere)==2 + ray_neg = Ray(p_far, p_far) + @test isnothing(ray_neg∩sphere) + end; + + @testset "laser sampling" begin + σ = rand() + 0.1 + direction = normalize(@SVector randn(3)) + center = @SVector rand(3) + las = Laser(σ, direction, center) #laser is pointed to xy plane + photons = [sample_source(las) for _ in 1:100000] + + photons_r = stack((p.r for p in photons)) + + @test isapprox(vec(mean(photons_r, dims=2)), center, atol=5.e-2) #center is ok + @test isapprox(mean(sum((photons_r .- center).^2, dims=1)), (√2σ)^2, atol=5.e-2) #σ is ok + + p1, p2 = rand(photons, 2) + @test (p1.r - p2.r) ⟂ direction + end + + @testset "scene creation" begin + mat1 = Turbid(μₛ=0.5, μₐ=0., n=1., g=0.) + + start_point = rand(3) + + s1 = Solid(Box(start_point, start_point + ones(3)), mat1) + + direction = 2*bitrand(3) .- 1. + mat2 = Turbid(μₛ=1.5, μₐ=0., n=1., g=0.) + s2 = Solid(Box(start_point+direction, start_point+direction+ones(3)), mat2) + + s3 = Solid(Box(min.(minimum(s1.geometry), minimum(s2.geometry)) + ,max.(maximum(s1.geometry), maximum(s2.geometry))), mat1) + + scene = Scene([s1, s2, s3]) + + p1 = Photon((start_point + 0.5*ones(3)), direction+1e-6*rand(3)) + + @test p1∈scene + @test p1∈s1 + @test p1∈s3 + + p2 = advance(p1, 1.) + + @test p2∈s2 + + @test find_container(p1, scene) == s1 + + rbound, normal = exit_solid(p1, find_container(p1, scene)) + @test rbound≈(p1.r+p2.r)/2 atol=1.e-3 + end; + + @testset "particle reflect" begin + #wikipedia example https://en.wikipedia.org/wiki/Snell%27s_law#Vector_form + p = Photon([0, 0, 0.], [1/sqrt(2), -1/sqrt(2), 0.]) + n1 = 1. + n2 = 1. / 0.9 + normal = [0, 1., 0.] + reflect(p, normal) + + @test all(reflect(p, normal).u .≈ [1/sqrt(2), 1/sqrt(2), 0.0]) + @test norm(refract(p, normal, n1, n2).u - [0.636396, -0.771362, 0.0]) < 1.e-5 + end; + + @testset "particle fresnel reflectance" begin + p = Photon(zeros(3), normalize(randn(3))) + normal = normalize(@SVector randn(3)) + normal = dot(normal, p.u)<0 ? normal : -normal + + @test unpolarized_reflectance(p, normal, 1., 1.) ≈ 0. atol=1.e-10 + @test unpolarized_reflectance(p, normal, 1., 1.e8) ≈ 1. atol=1.e-5 + @test unpolarized_reflectance(p, normal, 1.e8, 1.) ≈ 1. atol=1.e-5 + + n1 = rand()*3 + 1. + n2 = rand()*3 + 1. + @test unpolarized_reflectance(p, -p.u, n1, n2) ≈ (n1-n2)^2/(n1+n2)^2 #https://en.wikipedia.org/wiki/Fresnel_equations#Normal_incidence + end; + + @testset "particle snell law" begin + p = Photon(zeros(3), normalize(randn(3))) + + normal = normalize(@SVector randn(3)) + normal = dot(normal, p.u)<0 ? normal : -normal + + n1 = rand()*3 + 1. + n2 = rand()*3 + 1. + + reflected = reflect(p, normal) + refracted = refract(p, normal, n1, n2) + R = unpolarized_reflectance(p, normal, n1, n2) + + @test dot(p.u, normal) ≈ -dot(reflected.u, normal) # incident angle = reflection angle + + hit = [0., 0., 0.] + start = hit - p.u + hit2reflect = hit + reflected.u + if !isnothing(refracted) hit2refract = hit + refracted.u end + + normal_vec = hit + normal + end + + @testset "particle spin" begin + p = Photon(zeros(3), normalize(randn(3))) + g = rand()*0.2+0.7 + cosθ_samp = [dot(spin(p, g).u, p.u) for i in 1:100000] + @test g ≈ mean(cosθ_samp) atol=1.e-2 + + g = 0. + cosθ_samp = [dot(spin(p, g).u, p.u) for i in 1:100000] + @test g ≈ mean(cosθ_samp) atol=1.e-2 + end; + + @testset "russian roulette" begin + roulette_chance = 0.1 + + photons = vcat([Photon(zeros(3), normalize(randn(3)), rand()) for i in 1:1000000], + [Photon(zeros(3), normalize(randn(3)), 2*roulette_chance*rand()) for i in 1:1000000]) + + energy = sum(x->x.weight, photons) + + roulette_photons = [roulette(p,roulette_chance) for p in photons] + energy_roulette = sum(x->x.weight, roulette_photons) + + @test energy≈energy_roulette rtol=1.e-2 #energy conversion + end + end + + Nₚ = 100000 #number of particles + @testset verbose = true "Integral Tests" begin + @testset "PTMC benchmark - finite slab" begin + las = Laser(0.04, [0.,0.,-1.], [0.,0.,0.02]) #laser is pointed in the -̂z direction + + air = Solid(Box((-5000.,-5000.,0.), (5000.,5000.,0.04)), Turbid(); name="air") + scatterer = Solid(Box((-5000.,-5000.,-0.02), (5000.,5000.,0.)), Turbid(μₛ=90, μₐ=10, n=1., g=0.75); name="scatterer") + air2 = Solid(Box((-5000.,-5000.,-0.04), (5000.,5000.,-0.02)), Turbid(); name="air2") + scene = Scene([air, scatterer, air2]) + end_particles = [simulate_particle(las, scene) for i in 1:Nₚ] + + R = 0. + T = 0. + for p in end_particles + if p.r[3] > 0.01 R += p.weight/Nₚ end + if p.r[3] < -0.03 T += p.weight/Nₚ end + end + + #see chapter 6 in https://omlc.org/classroom/ece532/class4/MCMan.pdf + @test R ≈ 0.09739 rtol=5.e-1 + @test T ≈ 0.66096 rtol=5.e-1 + end; + + @testset "PTMC benchmark - infinite half-slab" begin + las = Laser(0.01, [0.,0.,-1.], [0.,0.,0.02]) #laser is pointed in the -̂z direction + + air = Solid(Box((-5000.,-5000.,0.), (5000.,5000.,0.04)), Turbid(); name="air") + scatterer = Solid(Box((-5000.,-5000.,-1000.), (5000.,5000.,0.)), Turbid(μₛ=90, μₐ=10, n=1.5, g=0.); name="scatterer") + scene = Scene([air, scatterer]) + end_particles = [simulate_particle(las, scene) for i in 1:Nₚ] + + R = 0. + for p in end_particles + if p.r[3] > 0.01 R += p.weight/Nₚ end + end + + @test R ≈ 0.26 rtol=5.e-1 #see chapter 6 in https://omlc.org/classroom/ece532/class4/MCMan.pdf + end + end +end;