PS-OCT-sim

Polarization-sensitive OCT simulation framework with radiative transfer equation (RTE/VRTE) forward chain and conservative scattering (ω = 1) operator-theoretic analysis.

MATLAB implementation with package namespaces. Loosely coupled, independently testable modules.

Chinese documentation: README_CN.md.

---

Table of Contents

• Project Structure
• Quick Start
• Models
• ω → 1 vs ω = 1 Analysis
• Output Figures
• Notes on PST θ Inversion
• Dependencies & References

---

Project Structure

PS-OCT-sim/
├── init.m                             Parameters
├── run_demo.m                         Main demo (6 parts)
├── README.md                          This file
├── README_CN.md                       Chinese README
├── .gitignore

├── +rte/                              Layer 1 — Radiative transfer
│   ├── GaussQuadrature.m             Gauss-Legendre quadrature
│   ├── PhaseFunctionHG.m             HG phase function (Legendre expansion)
│   ├── ScalarRTE.m                    Scalar RTE discrete ordinates solver
│   ├── VectorRTE.m                    Vector RTE (Stokes) solver
│   ├── VRTE_OperatorExact.m           Exact ω=1 VRTE operator (SVD kernel)
│   └── BeerLambert.m                  Beer-Lambert law baseline

├── +polarization/                     Layer 2 — Polarization optics
│   ├── JonesMatrix.m                  Uniaxial birefringent Jones matrix
│   ├── MuellerConvert.m              Jones → Mueller (Pauli basis)
│   ├── StokesPropagate.m              Depth-resolved Stokes evolution
│   └── PSTInverse.m                   PST inversion of Δn and θ

├── +oct/                              Layer 3 — OCT signal chain
│   ├── SourceSpectrum.m              Gaussian source (k-space)
│   ├── Interferogram.m               Michelson interferometer + noise
│   ├── DepthReconstruct.m            FFT reconstruction
│   └── PS_OCT_Bscan.m                B-scan simulation

├── +tissue/                           Layer 4 — Tissue models
│   ├── BrainCortexParams.m            Mouse cortex L1–L6 parameters
│   └── BirefringentSample.m           Δn/θ profile generator (5 modes)

├── +analysis/                          Analysis layer
│   ├── ConservativeScatteringTheory.m Operator-level ω=1 analysis
│   └── OmegaToOne.m                   Integration: theory + OCT diagnostics

└── output/                             Auto-generated figures
    ├── fig1_rte_vs_bl.png             RTE vs Beer-Lambert
    ├── fig2_omega_analysis.png        ω→1 analysis
    ├── fig3_stokes_evolution.png      Stokes + DOP
    ├── fig4_pst_inversion.png         PST Δn/θ inversion
    ├── fig5_ps_oct_bscan.png          B-scan H/V/DOP
    └── fig6_vrte_dop.png              VRTE DOP profiles

---

Quick Start

git clone https://github.com/<user>/PS-OCT-sim.git
cd('PS-OCT-sim');
run_demo

Requirements: MATLAB R2020b or later, no toolboxes needed. Paths are auto-added by init().

run_demo runs 6 parts sequentially:

Part	Content	Variables	Figure	Time
1	RTE vs Beer-Lambert	rte_bio, bl_bio	fig1	<1 s
2	ω→1 well-posedness	report	fig2	~5-15 s
3	Stokes depth evolution	S_z, J_cum	fig3	<1 s
4	PST inversion	delta_n_est, theta_est	fig4	<1 s
5	B-scan (128 A-lines)	bscan_result	fig5	~10-30 s
6	VRTE: polarization + scattering	vrte_results	fig6	~10-30 s

Key quantities in the report structure:

report.theory.ker_dim_omega1        % Scalar RTE ω=1 kernel dimension
report.theory.vrte_ker_dim_exact    % VRTE ω=1 kernel dimension
report.theory.scaling_power         % λ_min ~ C(1-ω)^p exponent
report.theory.T_ratio_ref           % RTE/BL transmittance ratio at ω≈0.99
report.comparison.rows              % ω→1 vs ω=1 comparison table
report.oct_diag.key_numbers         % OCT diagnostics (l*, kernel, etc.)

Standalone Module Usage

params = init();

% Scalar RTE
[mu, w] = rte.GaussQuadrature(16);
sol = rte.ScalarRTE(15e2, 0.1e2, 0.9, 1e-3, mu, w);

% Polarization propagation
[~, dn, th] = tissue.BirefringentSample(200, 10e-6, 'Profile', 'layered_random');
z = (0:200)' * 10e-6;
S = polarization.StokesPropagate(z, dn, th, 1310e-9, [1;0;0;1]);

% Interferogram
[I_det, k, spec] = oct.Interferogram(...
    [0; 500e-6; 1000e-6], [0.1; 0.05; 0.02], 1.36, 1310e-9, 100e-9, 2048);
[A_lin, z, A_dB] = oct.DepthReconstruct(I_det, k, 1.36);

---

Models

Layer 1: RTE

The steady-state scalar RTE in plane-parallel geometry:

μ ∂I/∂τ + I = (ω/2) ∫ p(μ,μ') I(τ,μ') dμ'

• τ = μ_t·z (optical depth), μ = cos θ, ω = μ_s/μ_t (albedo, ω = 1 = purely scattering)
• HG phase function via Legendre expansion, truncation L ≈ ⌈1/(1-|g|)⌉
• Numerical method: Gauss-Legendre discrete ordinates, eigenvalue decomposition of the 2N×2N system matrix, boundary-condition solve
• Normalization: DISORT standard ∫p dμ = 2, internally consistent with ω/2 prefactor. Verified: Σ_j w_j P_ij = 2 to within 2.3×10⁻¹⁴

VRTE: Stokes 4-vector + depolarizing 4×4 phase matrix (parameter δ), system matrix 8N×8N.

Layer 2: Polarization Optics

Depth-stratified uniaxial birefringence. Layer k: fast axis θ_k [rad], phase retardance δ_k = 2π Δn_k d_k / λ₀. Unitary Jones matrix, cumulative product J(z) = Π J_k. Pauli-basis Mueller conversion.

PST inversion: Rotation axis and angle of the Stokes trajectory on the Poincaré sphere between adjacent depths directly encode local Δn and θ. Steps: normalize S→[Q,U,V]/I → moving-average smooth → skip DOP < 0.1 → rotation angle δ = 2·arcsin(|Δs|/2|s|) → axis a = s_i × s_{i+1} → Δn = δ·λ₀/(2π·Δz) → θ = ½·atan2(a_U, a_Q), wrapped to [0, π).

Layer 3 & 4

Gaussian source S(k), Michelson interferogram + shot/phase noise, FFT reconstruction (DC removal → windowing → zero-padding → FFT). Mouse cortex L1–L6 layered parameters. 5 birefringence profile modes including random walk (Ornstein-Uhlenbeck).

---

ω → 1 vs ω = 1 Analysis

Tissue ω ≈ 0.99. ω = 1 is a structurally singular point, not merely the endpoint of the ω → 1 limit.

The diagonal cancellation: At ω = 1, ∫p dμ'·ω/2μ = 1/2μ = the extinction coefficient. The diagonal blocks of the system matrix A cancel → nontrivial kernel emerges. The kernel dimension has a step discontinuity at ω = 1 — zero for all ω < 1, ≥1 at ω = 1. This cannot be predicted from any ω → 1 eigenvalue sequence; it requires direct construction of the ω = 1 operator.

Implementation: exact ω=1 operators for scalar RTE (via compute_RTE_operator) and VRTE (via rte.VRTE_OperatorExact). Kernel detected via SVD (singular values < 10⁻⁷). The ω sweep (0.5→0.999) independently validates the discontinuity.

Comparison Table

Property	ω → 1	ω = 1
A diagonal	T + Σ - S, Σ>0	Σ ≡ S cancels on average
ker(A)	{0}	≥1 (scalar); 1 (VRTE)
λ_min	~C(1-ω)^p, finite→0⁺	≡0 (exact degeneracy)
cond(B)	bounded→∞	∞ (singular)
Solution	unique	non-unique (up to kernel)

PS-OCT Implications

1. Depth attenuation deviates from pure exponential beyond l* ≈ 670 µm
2. H/V channel coupling: f(Δn)·g(ω)·h(z/l*), with g(ω) diverging as ω→1
3. Conjecture: DOP decay under conservative scattering is multiple-scattering-limited → PS-OCT imaging depth may be DOP-limited rather than SNR-limited (not experimentally verified)

---

Output Figures

Figure	Content	Expected behavior
fig1	RTE forward/backward flux + BL comparison	RTE flux » BL (multi-scatter); BL is exponential
fig2	λ_min vs ω, cond(B) scaling, VRTE cond	λ_min→0 with ω→1; cond(B) ~ 1/(1-ω)^p
fig3	Stokes 4-vector + DOP vs depth	DOP≡1 (pure Jones, no depolarization); Q/U oscillate
fig4	PST Δn and θ: true vs estimated	Δn RMSE ~10⁻⁵; θ: see known issues
fig5	H/V B-scan + DOP map	H/V envelops similar; DOP stays high
fig6	VRTE DOP profiles + transmittance	DOP decays faster with higher ω; T rises

---

Notes on PST θ Inversion

θ Panel (fig4, right subplot)

The θ true-vs-estimated traces may exhibit local discontinuities (spikes/folds). This is an inherent consequence of the Jones matrix symmetry, not a code defect.

Physical origin: The Jones matrix depends on cos(2θ) and sin(2θ) — θ and θ+π (and θ and θ+π/2 paired with Δn sign flip) produce physically identical matrices. PST extracts θ via atan2, producing output naturally bounded in [0, π). The ground-truth profile (random walk in layered_random mode) is folded to the same range via mod(θ, π) — the random walk occasionally crosses the π boundary, producing visible discontinuities at fold points.

Scope: This affects the visual appearance of the θ panel only. Δn inversion accuracy is unaffected (Δn RMSE ~10⁻⁵). Stokes propagation is unaffected (Jones matrix automatically handles π-periodicity).

Evaluation guidance: Use Δn RMSE as the primary quality metric for PST inversion. The θ panel validates the trend direction (whether the two traces roughly track each other), but deviations at folding points do not reflect algorithmic error.

Other Notes

• Set params.verbose = false to suppress console output
• All figures auto-saved to output/ by run_demo
• rng(42) seeds Part 3-4 for reproducible profiles
• Parts 2, 5, 6 are time-intensive and can be commented out for rapid iteration

---

License

MIT

Dependencies & References

• MATLAB R2020b or later, no toolboxes required
• Paths auto-added by init()

Key references:

• Stamnes K, et al. "Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media." Appl Opt, 1988
• Thomas GE, Stamnes K. Radiative Transfer in the Atmosphere and Ocean. Cambridge, 1999
• Egger H, Schlottbom M. "Analysis of the radiative transfer equation with almost isotropic scattering." M3AS, 2014
• Bosboom V, Schlottbom M, Schwenninger FL. "On the transition from radiative transfer to diffusion in the limit ω → 1." JMAA, 2024
• Lehner J, Wing GM. "On the spectrum of an unsymmetric operator arising in the transport theory of neutrons." CPAM, 1955
• Makita S, Yamanari M, Yasuno Y. "Polarization state tracing method for local birefringence imaging." BOE, 2020