Core Data Types

This section covers the fundamental data structures used for 3D tomography reconstruction, outlining the distinct representations needed for phantoms, grids, and the mathematical coefficients used in both mean and covariance estimations.

3D Phantom

For a 3D phantom consisting of several Gaussian blobs, we use a specialized wrapper:

HeteroTomo3D.KernelPhantom3D — Type

KernelPhantom3D{T<:Real}

A 3D phantom represented functionally as a linear combination of Gaussian kernels. This avoids heavy 3D voxel grid allocations and provides exact analytic X-ray transforms.

\[f_i (\mathbf{z}) = \sum_{l=1}^L a_{l,i} \exp(-\gamma_l \| \mathbf{z} - \mathbf{c}_l \|^2)\]

Fields

weights::Matrix{T}: The weight coefficients of size L × n for n phantoms.
centers::Vector{NTuple{3,T}}: The center coordinates of each Gaussian kernel in 3D space.
gammas::Vector{T}: The bandwidth parameters of each 3D Gaussian kernel.

source

HeteroTomo3D.rand_center_grid — Function

rand_center_grid(L::Int; seed::Union{Nothing,Int}=nothing) -> Vector{NTuple{3,T}}

Create L random 3D centers within a unit sphere.

Arguments

L::Int: Number of random 3D centers to generate.

Keyword Argument

seed::Union{Nothing, Int}: Random seed for reproducibility. Defaults to nothing.

Examples

julia> using HeteroTomo3D;

julia> using Random; Random.seed!(42);

julia> L = 2;

julia> centers = rand_center_grid(L; seed=42)
2-element Vector{Tuple{Float64, Float64, Float64}}:
 (0.25869024628521786, -0.09932211880761277, -0.04518571313436448)
 (0.4062596980064028, 0.3466922912789925, -0.6682111304137319)

source

using HeteroTomo3D, BlockArrays, LinearAlgebra, Random

centers = [
    (0.3, 0.3, 0.3),
    (-0.3, -0.3, 0.3),
    (-0.4, 0.4, -0.4),
    (0.3, -0.3, -0.3)
]
weights = reshape([1.0, 0.8, 0.6, 0.4], 4, 1)
gammas = [5.0, 4.0, 6.0, 4.0]

phantom = KernelPhantom3D(weights, centers, gammas)

KernelPhantom3D{Float64}([1.0; 0.8; 0.6; 0.4;;], [(0.3, 0.3, 0.3), (-0.3, -0.3, 0.3), (-0.4, 0.4, -0.4), (0.3, -0.3, -0.3)], [5.0, 4.0, 6.0, 4.0])

Observational Grid

The collection of 2D evaluation points and viewing angles are represented by arrays of type NTuples{2, I} and UnitQuaternion{T}, respectively:

HeteroTomo3D.EvaluationGrid — Type

EvaluationGrid{I<:Integer}

3D array of size (2, s, r, n) storing the discrete 2D detector voxel coordinates for X-ray evaluation, i.e., $X[i, j, k] \in \mathbb{Z}^{2}$.

Fields

blocks::Array{NTuple{2,I}, 3}: 3D array of size (s, r, n) containing the coordinates.
s::Int64: Number of evaluation points per view.
r::Int64: Number of projection angles.
n::Int64: Number of 3D functions.
m::Int64: Voxel grid resolution/dimension parameter.

julia> using HeteroTomo3D;

julia> s, r, n, m = 2, 3, 2, 64;

julia> blocks = fill((1, 1), s, r, n);

julia> eval_grid = EvaluationGrid(blocks, s, r, n, m);

julia> size(eval_grid) == (s, r, n)
true

Mean Estimation

In mean estimation, we solve for a single global function. The Tikhonov solution of the representer theorem is defined by a standard vector of coefficients.

The continuous function is expanded as:

\[\hat{f} = \sum_{i=1}^{n} \sum_{j=1}^{r} \sum_{k=1}^{s} \mathbf{a}_{ijk} \varphi_{\gamma} (\mathbf{R}_{\mathbf{q}_{ij}}, \mathbf{x}_{ijk})\]

Here, the coefficient is conceptualized as a vector of length s * r * n and is stored using column-major indexing at [k + (j-1) * s + (i-1) * s * r].

Covariance Estimation

Unlike mean estimation, which outputs a single volume, covariance estimation models the spatial relationship between points. This requires a much larger tensor-based expansion.

The covariance solution is spanned as follows:

\[\hat{f} = \sum_{i=1}^{n} \sum_{j, j'=1}^{r} \sum_{k, k'=1}^{s} \mathbf{A}_{i,jk, j'k'} \varphi_{ijk} \otimes \varphi_{ij'k'},\]

where $\varphi_{ijk} = \varphi_{\gamma} (\mathbf{R}_{\mathbf{q}_{ij}}, \mathbf{x}_{ijk})$.

Coefficient Data Types

Because the covariance coefficients form a block-diagonal tensor rather than a simple vector, we provide custom types to manage them:

HeteroTomo3D.BlockDiagonal — Type

BlockDiagonal{T, R<:(AbstractArray{<:AbstractArray{T, 2}, 1})}

Blocked array of square matrices of type T arranged in a block diagonal structure.

Fields

blocks::R: stores the vector of matrices of type T being wrapped.

Supertype Hierarchy

BlockDiagonal{T, R<:(AbstractArray{<:AbstractArray{T, 2}, 1})} <: AbstractArray{T, 1} <: Any

Examples

julia> using HeteroTomo3D;

julia> A = [1 2; 3 4];

julia> B = [5 6; 8 9];

julia> bd = BlockDiagonal([A, B]);

julia> bd.blocks[1]
2×2 Matrix{Int64}:
 1  2
 3  4

source

HeteroTomo3D.block_outer — Function

block_outer(y)
y ⊙ y

Computes the block-wise outer product of a BlockVector, returning a BlockDiagonal matrix.

For a BlockVector y with blocks y₁, y₂, ..., this function computes a BlockDiagonal matrix where the i-th block is the outer product yᵢ * yᵢ'.

The infix operator ⊙ is an alias for this function. It must be used on the same object (e.g., y ⊙ y).

Arguments

y::AbstractBlockVector: The input block vector.

Examples

julia> using HeteroTomo3D, BlockArrays;

julia> v = [1.0, 2.0, 3.0];

julia> y = BlockVector(v, [2, 1]);

julia> (y ⊙ y).blocks[1]
2×2 Matrix{Float64}:
 1.0  2.0
 2.0  4.0

source

HeteroTomo3D.:⊙ — Function

y ⊙ y

Computes the block-wise outer product of a BlockVector, returning a BlockDiagonal matrix.

The infix operator ⊙ is an alias for block_outer(y), and must be used on the same object (e.g., y ⊙ y).

Lazy Tensor Operations

To avoid heap memory allocations when working with the Khatri-Rao product of the mean Gram matrix, we wrap these operations in a lazy fashion. This ensures space complexity remains optimized during continuous solver iterations:

HeteroTomo3D.AbstractBlockTensor — Type

AbstractBlockTensor{T}

Element-free tensor type for block-wise operations.

source

HeteroTomo3D.BlockOuter — Type

BlockOuter(blocks, [workspace])
K ⊙ K

Represents a linear operator that performs a block-wise outer product. Specifically, for a BlockDiagonal matrix A, the output B = (K ⊙ K)(A) is computed as B.blocks[j] = ∑ᵢ K.blocks[j, i] * A.blocks[i] * K.blocks[i, j].

The infix operator ⊙ is an alias for the constructor BlockOuter(K). It requires that both operands are the same object (e.g., K ⊙ K).

Fields

K::AbstractBlockMatrix: The block matrix that defines the operator.
workspace::Matrix: Pre-allocated matrix buffer for intermediate iterations.

Examples

julia> using HeteroTomo3D, BlockArrays;

julia> K = BlockMatrix(rand(3, 3), [2, 1], [2, 1]);

julia> L = K ⊙ K;

julia> L isa BlockOuter
true

julia> size(L)
(5, 5)

Functional PCA

Once the covariance coefficients are found, you can extract the primary modes of variation (eigenfunctions) using the Conjugate Lanczos Algorithm:

HeteroTomo3D.conj_lanczos — Function

conj_lanczos(b0, D, K; ...)

Performs the C-Lanczos algorithm for the generalized eigenvalue problem D*C*q = λ*q.

This iterative method generates a K-orthonormal basis Q for the Krylov subspace and a symmetric tridiagonal matrix T that represents the projection of the operator D * K onto that subspace.

Arguments

b0::AbstractBlockVector{T}: The starting vector for the iteration.
D::BlockDiagonal{T}: A block-diagonal matrix to factorize.
K::AbstractBlockMatrix{T}: A PSD matrix defining the inner product.

Keyword Arguments

itmax::Int=100: Maximum number of iterations.
tol::Float64=1e-8: Tolerance for the K-norm of the residual to determine convergence.
history::Bool=false: If true, the history of residual K-norms is stored and returned.
reortho_level::Symbol=:full: Level of reorthogonalization. Use :full for numerical stability, or :none to observe loss of orthogonality.

Returns

NamedTuple: A named tuple (T, Q, history) containing:
- T: A SymTridiagonal matrix.
- Q: A matrix whose columns are the K-orthonormal Lanczos basis vectors.
- history: A vector of residual K-norms, or nothing if history=false.

source

HeteroTomo3D.fpca — Function

fpca(k, A, K; itmax=50) -> Tuple{Vector, Matrix}

Solves the k leading eigenvalue problem A*K*v = λ*v under K-orthogonality. It uses the Conjugate Lanczos algorithm (conj_lanczos) to iteratively find the eigenvalues.

Arguments

k::Int: The number of principal components to compute.
A::BlockDiagonal: A block-diagonal matrix representing the prior covariance of the coefficients (the A matrix in the eigenproblem).
K::BlockMatrix: A Gram matrix representing the inner product for the eigenproblem.

Keyword Arguments

itmax::Int=50: The maximum number of Krylov subspace iterations.

source