initial commit

.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"

.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

.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()'

.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 }}

.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

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.

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"]

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.
+```

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

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)

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

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))))

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