Showing revision 1.1
Redbird

3-D Model-based Tomographic Reconstruction Platform

for Diffuse Optical Tomography (DOT), functional Near-Infrared Spectroscopy (fNIRS), Microwave Imaging, and more

MATLAB / Octave
Python
20+
Years Research
45+
Rich Forward+Inverse Functions
FEM
Fast&Deterministic
Download MATLABDownload PythonDownload Manual

Core Features

🔬

FEM Forward Solver

Diffusion equation solver validated against Monte Carlo

🎯

Gauss-Newton Recon

Iterative nonlinear reconstruction with Tikhonov

🌈

Multi-Spectral

Chromophore estimation: HbO, HbR, water, lipids

🧠

Structure Priors

Laplacian, Helmholtz, compositional priors

📡

Wide-Field Sources

Planar, pattern, Fourier-basis illumination

⚡

Dual-Mesh

Separate forward/inverse meshes for speed

Quick Start: cfg = rbmeshprep(cfg); [det, phi] = rbrun(cfg);

Forward Modeling

Photon Propagation

Solve diffusion equation for fluence using FEM. CW and frequency-domain supported.

rbrun(cfg)One-liner forward solver
rbfemlhs()Build FEM stiffness matrix
rbfemrhs()Build RHS source vectors
rbfemgetdet()Extract detector values
cfg.omega = 2*pi*100e6 enables 100 MHz FD mode
% Create mesh and set optical properties
[cfg.node, cfg.face, cfg.elem] = meshabox([0 0 0], [60 60 30], 5);
cfg.prop = [0 0 1 1; 0.01 1 0 1.37];
cfg.seg = ones(size(cfg.elem,1),1);

% Source & detector configuration
cfg.srcpos = [30 30 0]; cfg.srcdir = [0 0 1];
cfg.detpos = [30 40 0; 40 30 0];
cfg.omega = 0;  % CW mode

% Run forward simulation
cfg = rbmeshprep(cfg);
[detval, phi] = rbrun(cfg);
import redbirdpy as rb
from iso2mesh import meshabox
import numpy as np

# Create mesh and set optical properties
node, face, elem = meshabox([0,0,0], [60,60,30], 5)
cfg = {'node': node, 'elem': elem,
  'prop': np.array([[0,0,1,1], [0.01,1,0,1.37]]),
  'srcpos': [[30,30,0]], 'omega': 0}

cfg, sd = rb.meshprep(cfg)
detval, phi = rb.run(cfg)

Image Reconstruction

Iterative Gauss-Newton with regularization, dual-mesh, and structure priors.

Bulk Fitting

Single property set for entire domain

Segmented (Hard-Prior)

One property per labeled tissue segment

Soft-Prior

Spatial priors as soft constraints from MRI/CT

Unconstrained

Independent per-node properties with Tikhonov

Mesh Generation

iso2mesh Integration

Seamless tetrahedral mesh generation from surfaces, volumes, or analytical shapes.

meshabox()Create box mesh
meshasphere()Create sphere mesh
mergemesh()Combine meshes
s2m()Surface to tetrahedral
rbmeshprep(cfg) computes face, area, evol, nvol, reff, deldotdel
% Create box with spherical inclusion
[nobbx, fcbbx] = meshabox([0 0 0], [60 60 30], 5);
[nosp, fcsp] = meshasphere([30 30 15], 8, 1);
[no, fc] = mergemesh(nobbx, fcbbx, nosp, fcsp);

% Generate tetrahedral mesh
[cfg.node, cfg.elem] = s2m(no, fc(:,1:3), 1, 5);
cfg = rbmeshprep(cfg);
import iso2mesh as i2m

nobbx, fcbbx, _ = i2m.meshabox([0,0,0], [60,60,30], 5)
nosp, fcsp, _ = i2m.meshasphere([30,30,15], 8, 1)
no, fc = i2m.mergemesh(nobbx, fcbbx, nosp, fcsp)

node, elem, _ = i2m.s2m(no, fc[:,:3], 1, 5)
cfg, sd = rb.meshprep({'node': node, 'elem': elem})

Linear Solvers

Multiple solver backends for sparse FEM systems. Block QMR for multiple RHS.

Direct Solvers

PARDISO, UMFPACK, CHOLMOD, SuperLU

Iterative Solvers

CG, GMRES, BiCGSTAB, MINRES, QMR

Block QMR

Multiple RHS — up to 10× faster

Examples

🔬 Basic Forward Simulation

Simple CW forward simulation comparing with analytical solution

% Create mesh and setup simulation
[cfg.node, cfg.face, cfg.elem] = meshabox([0 0 0], [60 60 30], 1);
cfg.seg = ones(size(cfg.elem, 1), 1);
cfg.srcpos = [30 30 0];
cfg.srcdir = [0 0 1];
cfg.detpos = [30 30 30];
cfg.detdir = [0 0 -1];

% Optical properties [mua, mus', g, n]
cfg.prop = [0 0 1 1; 0.005 1 0 1.37];
cfg.omega = 0;  % CW mode

% Prepare mesh and run
cfg = rbmeshprep(cfg);
[detphi, phi] = rbrun(cfg);

% Analytical solution for comparison
phicw = cwdiffusion(cfg.prop(2,1), cfg.prop(2,2)*(1-cfg.prop(2,3)), ...
                    cfg.reff, [30 30 0], cfg.node);

% Visualization
[cutpos, cutval] = qmeshcut(cfg.elem, cfg.node, phi(:,1), 'x=29.5');
contour(xi, yi, log10(vphi), 0:-0.5:-5, 'r-');
import redbirdpy as rb
import numpy as np
from iso2mesh import meshabox

# Create mesh
node, face, elem = meshabox([0,0,0], [60,60,30], 1)

cfg = {
    'node': node, 'elem': elem,
    'seg': np.ones(elem.shape[0], dtype=int),
    'srcpos': np.array([[30,30,0]]),
    'srcdir': np.array([[0,0,1]]),
    'detpos': np.array([[30,30,30]]),
    'detdir': np.array([[0,0,-1]]),
    'prop': np.array([[0,0,1,1], [0.005,1,0,1.37]]),
    'omega': 0
}

# Run simulation
cfg, sd = rb.meshprep(cfg)
detphi, phi = rb.run(cfg)

# Analytical comparison
from redbirdpy.analytical import semi_infinite_cw
phicw = semi_infinite_cw(cfg['prop'][1,0], 
    cfg['prop'][1,1]*(1-cfg['prop'][1,2]),
    cfg['prop'][1,-1], 1.0, cfg['srcpos'][0], cfg['node'])

🧊 Heterogeneous Domain

Forward simulation with box inclusion (different optical properties)

% Create mesh with box inclusion
boxsize = [60, 50, 40];
[nbox1, fbox1] = meshabox([0,0,0], boxsize, 4);
[nbox2, fbox2] = meshabox([10,10,10], [30,30,30], 4);
[nbox1, fbox1] = removeisolatednode(nbox1, fbox1(:,1:3));
[nbox2, fbox2] = removeisolatednode(nbox2, fbox2(:,1:3));
[no1, fc1] = mergemesh(nbox1, fbox1, nbox2, fbox2);

% Seed points for region labels
regionseed = [1,1,1; 11,11,11];
[cfg.node, cfg.elem] = s2m(no1, fc1, 1, 4, 'tetgen', regionseed);
cfg.seg = cfg.elem(:,5);
cfg.elem(:,5) = [];

% Properties: background (label 1) and inclusion (label 2)
cfg.prop = [0 0 1 1;              % label 0
            0.006 0.8 0 1.37;     % label 1 (background)
            0.02  1   0 1.37];    % label 2 (inclusion)

cfg.srcpos = [25, 25, 0];
cfg.srcdir = [0, 0, 1];
cfg.detpos = [35, 25, max(cfg.node(:,3))];
cfg.detdir = [0, 0, -1];

cfg = rbmeshprep(cfg);
[detphi, phi] = rbrunforward(cfg);

% phi(:,1) is from source, phi(:,2) is from detector
plotmesh([cfg.node log10(phi(:,1))], cfg.elem, 'y > 25');
import redbirdpy as rb
from redbirdpy import forward
import iso2mesh as i2m
import numpy as np

# Create mesh with box inclusion
boxsize = [60, 50, 40]
nbox1, fbox1, _ = i2m.meshabox([0,0,0], boxsize, 4)
nbox2, fbox2, _ = i2m.meshabox([10,10,10], [30,30,30], 4)
nbox1, fbox1 = i2m.removeisolatednode(nbox1, fbox1[:,:3])[:2]
nbox2, fbox2 = i2m.removeisolatednode(nbox2, fbox2[:,:3])[:2]
no1, fc1 = i2m.mergemesh(nbox1, fbox1[:,:3], nbox2, fbox2[:,:3])[:2]

regionseed = [[1,1,1], [11,11,11]]
node, elem, _ = i2m.s2m(no1, fc1[:,:3], 1, 4, 'tetgen', regionseed)

cfg = {
    'node': node, 'elem': elem[:,:4], 'seg': elem[:,4].astype(int),
    'prop': np.array([[0,0,1,1], [0.006,0.8,0,1.37], [0.02,1,0,1.37]]),
    'srcpos': np.array([[25,25,0]]),
    'srcdir': [0,0,1],
    'detpos': np.array([[35,25,np.max(node[:,2])]]),
    'detdir': [0,0,-1]
}

cfg, sd = rb.meshprep(cfg)
detphi, phi = forward.runforward(cfg, sd=sd)
# phi[:,0] from source, phi[:,1] from detector

🎯 Image Reconstruction

CW reconstruction of absorption target with spherical inclusion

% Create phantom with spherical inclusion
s0 = [70, 50, 20];
[nobbx, fcbbx] = meshabox([40 0 0], [160 120 60], 10);
[nosp, fcsp] = meshasphere(s0, 5, 1);
[no, fc] = mergemesh(nobbx, fcbbx, nosp, fcsp);
[cfg0.node, cfg0.elem] = s2m(no, fc(:,1:3), 1, 40, 'tetgen', [41 1 1; s0]);
cfg0.seg = cfg0.elem(:,5);

% Source/detector grid
[xi, yi] = meshgrid(60:20:140, 20:20:100);
cfg0.srcpos = [xi(:), yi(:), zeros(numel(yi),1)];
cfg0.detpos = [xi(:), yi(:), 60*ones(numel(yi),1)];
cfg0.prop = [0 0 1 1; 0.008 1 0 1.37; 0.016 1 0 1.37];
cfg0 = rbmeshprep(cfg0);

% Generate data from heterogeneous phantom
detphi0 = rbrun(cfg0);

% Setup homogeneous forward mesh
[node, face, elem] = meshabox([40 0 0], [160 120 60], 10);
cfg = rbsetmesh(cfg, node, elem, cfg.prop, ones(size(node,1),1));
sd = rbsdmap(cfg);

% Coarse reconstruction mesh
[recon.node, ~, recon.elem] = meshabox([40 0 0], [160 120 60], 20);
[recon.mapid, recon.mapweight] = tsearchn(recon.node, recon.elem, cfg.node);

% Run reconstruction
[newrecon, resid, newcfg] = rbrun(cfg, recon, detphi0, sd, 'lambda', 1e-4);

% Plot results
plotmesh([newcfg.node, newcfg.prop(:,1)], newcfg.elem, 'z=20');
import redbirdpy as rb
import iso2mesh as i2m
import numpy as np

# Create phantom with spherical inclusion
s0 = [70, 50, 20]
nobbx, fcbbx, _ = i2m.meshabox([40,0,0], [160,120,60], 10)
nosp, fcsp, _ = i2m.meshasphere(s0, 5, 1)
no, fc = i2m.mergemesh(nobbx, fcbbx, nosp, fcsp[:,:3])
node, elem, _ = i2m.s2m(no, fc[:,:3], 1, 40, 'tetgen', [[41,1,1], s0])

xi, yi = np.meshgrid(np.arange(60,141,20), np.arange(20,101,20))
cfg0 = {
    'node': node, 'elem': elem[:,:4], 'seg': elem[:,4].astype(int),
    'srcpos': np.c_[xi.flat, yi.flat, np.zeros(xi.size)],
    'detpos': np.c_[xi.flat, yi.flat, 60*np.ones(xi.size)],
    'prop': np.array([[0,0,1,1], [0.008,1,0,1.37], [0.016,1,0,1.37]]),
    'omega': 0
}
cfg0 = rb.meshprep(cfg0)[0]
detphi0 = rb.run(cfg0)[0]

# Setup recon mesh
recon = {}
recon['node'], _, recon['elem'] = i2m.meshabox([40,0,0], [160,120,60], 20)
recon['mapid'], recon['mapweight'] = i2m.tsearchn(
    recon['node'], recon['elem'], cfg['node'])

# Run reconstruction
newrecon, resid, newcfg = rb.run(cfg, recon, detphi0, sd, lambda_=1e-4)[:3]

🌈 Multi-Spectral Reconstruction

Reconstruct chromophore concentrations (HbO, HbR) from multi-wavelength data

% Create phantom with inclusion
s0 = [70, 50, 20]; rs = 5;
[nobbx, fcbbx] = meshabox([40 0 0], [160 120 60], 10);
[nosp, fcsp] = meshasphere(s0, rs, 1);
[no, fc] = mergemesh(nobbx, fcbbx, nosp, fcsp);
[cfg0.node, cfg0.elem] = s2m(no, fc(:,1:3), 1, 40, 'tetgen', [41 1 1; s0]);
cfg0.seg = cfg0.elem(:,5);

% Define chromophore concentrations (wavelength-independent)
cfg0.param.hbo = [15 30];  % HbO in uM [background, inclusion]
cfg0.param.hbr = [4  8];   % HbR in uM [background, inclusion]

% Multi-wavelength optical properties (to be updated from param)
cfg0.prop = containers.Map();
cfg0.prop('690') = [0 0 1 1; 0 1 0 1.37; 0 1 0 1.37];
cfg0.prop('830') = [0 0 1 1; 0 0.8 0 1.37; 0 0.8 0 1.37];

cfg0.omega = 0;
cfg0 = rbmeshprep(cfg0);
detphi0 = rbrun(cfg0);  % Forward for all wavelengths

% Setup reconstruction mesh
[node, face, elem] = meshabox([40 0 0], [160 120 60], 10);
cfg = rbsetmesh(cfg, node, elem, cfg.prop, ones(size(node,1),1));
[recon.node, ~, recon.elem] = meshabox([40 0 0], [160 120 60], 25);
[recon.mapid, recon.mapweight] = tsearchn(recon.node, recon.elem, cfg.node);
sd = rbsdmap(cfg);

% Initial guesses for chromophores
recon.bulk = struct('hbo', 8, 'hbr', 2);
recon.param = struct('hbo', 8, 'hbr', 2);
recon.prop = containers.Map({'690','830'}, {[], []});

% Bulk fitting first
[newrecon, resid] = rbrun(cfg, recon, detphi0, sd, 'mode', 'bulk', 'lambda', 1e-3);

% Full image reconstruction
recon.bulk.hbo = newrecon.param.hbo;
recon.bulk.hbr = newrecon.param.hbr;
[newrecon, resid, newcfg] = rbrun(cfg, recon, detphi0, sd, 'mode', 'image', 'lambda', 1e-3);

% Plot reconstructed HbO
plotmesh([newrecon.node, newrecon.param.hbo(:)], newrecon.elem, 'z=20');
import redbirdpy as rb
import iso2mesh as i2m
import numpy as np

# Create phantom with inclusion
s0 = [70, 50, 20]; rs = 5
nobbx, fcbbx, _ = i2m.meshabox([40,0,0], [160,120,60], 10)
nosp, fcsp, _ = i2m.meshasphere(s0, rs, 1)
no, fc = i2m.mergemesh(nobbx, fcbbx, nosp, fcsp[:,:3])
node, elem, _ = i2m.s2m(no, fc[:,:3], 1, 40, 'tetgen', [[41,1,1], s0])

# Chromophore concentrations [background, inclusion]
cfg0 = {
    'node': node, 'elem': elem[:,:4], 'seg': elem[:,4].astype(int),
    'param': {'hbo': [15, 30], 'hbr': [4, 8]},  # in uM
    'prop': {
        '690': np.array([[0,0,1,1], [0,1,0,1.37], [0,1,0,1.37]]),
        '830': np.array([[0,0,1,1], [0,0.8,0,1.37], [0,0.8,0,1.37]])
    },
    'omega': 0
}

cfg0 = rb.meshprep(cfg0)[0]
detphi0 = rb.run(cfg0)[0]

# Setup reconstruction
recon = {
    'bulk': {'hbo': 8, 'hbr': 2},
    'param': {'hbo': 8, 'hbr': 2},
    'prop': {'690': [], '830': []}
}

# Bulk fitting then image reconstruction
newrecon, _ = rb.run(cfg, recon, detphi0, sd, mode='bulk', lambda_=1e-3)
recon['bulk'] = newrecon['param']
newrecon, resid, newcfg = rb.run(cfg, recon, detphi0, sd, mode='image', lambda_=1e-3)[:3]

📡 Wide-Field Sources

Planar source/detector for spatial frequency domain imaging

% Create regular grid mesh
[cfg.node, cfg.elem] = meshgrid6(0:2:60, 0:2:60, 0:2:30);
cfg.elem(:,1:4) = meshreorient(cfg.node(:,1:3), cfg.elem(:,1:4));
cfg.face = volface(cfg.elem);
cfg.seg = ones(size(cfg.elem,1), 1);

% Planar source configuration
cfg.srctype = 'planar';
cfg.srcpos = [9.5 9.5 0];
cfg.srcparam1 = [40 0 0 0];
cfg.srcparam2 = [0 40 0 0];
cfg.srcdir = [0 0 1];

% Planar detector configuration
cfg.dettype = 'planar';
cfg.detpos = [10 10 30];
cfg.detparam1 = [40 0 0 0];
cfg.detparam2 = [0 40 0 0];
cfg.detdir = [0 0 -1];

cfg.prop = [0 0 1 1; 0.005 1 0 1.37];
cfg.omega = 0;
cfg.srcweight = 2;

cfg = rbmeshprep(cfg);

% Solve
[Amat, deldotdel] = rbfemlhs(cfg, cfg.deldotdel);
[rhs, loc, bary] = rbfemrhs(cfg);
phi = rbfemsolve(Amat, rhs);
phi(phi < 0) = 0;
detval = rbfemgetdet(phi, cfg, rhs);

% Visualization
plotmesh([cfg.node log10(phi(:,1))], cfg.elem, 'x>30');
import redbirdpy as rb
from redbirdpy import forward
from redbirdpy.solver import femsolve
import iso2mesh as i2m
import numpy as np

# Create mesh
node, face, elem = i2m.meshabox([0,0,0], [60,60,30], 2)

cfg = {
    'node': node, 'elem': elem[:,:4], 'face': face[:,:3],
    'seg': np.ones(elem.shape[0], dtype=int),
    # Planar source
    'srctype': 'planar',
    'srcpos': [9.5, 9.5, 0],
    'srcparam1': [40, 0, 0, 0],
    'srcparam2': [0, 40, 0, 0],
    'srcdir': [0, 0, 1],
    # Planar detector
    'dettype': 'planar',
    'detpos': [10, 10, 30],
    'detparam1': [40, 0, 0, 0],
    'detparam2': [0, 40, 0, 0],
    'detdir': [0, 0, -1],
    'prop': np.array([[0,0,1,1], [0.005,1,0,1.37]]),
    'omega': 0, 'srcweight': 2
}

cfg = rb.meshprep(cfg)[0]

# Build and solve
deldotdel, _ = forward.deldotdel(cfg)
Amat = forward.femlhs(cfg, deldotdel)
rhs, loc, bary, _ = forward.femrhs(cfg)
phi, flag = femsolve(Amat, rhs)
phi = np.real(phi); phi[phi < 0] = 0
detval = forward.femgetdet(phi, cfg, rhs, loc, bary)

Get Started

1

Register & Download

Redbird for MATLAB/Octave (stable) Redbird for Python (stable)

2

MATLAB / Octave

addpath('/path/to/redbird/matlab')

3

Python

import redbirdpy as rb

4

Install iso2mesh

Required for mesh generation: iso2mesh.sf.net

5

For developers

Redbird Github (unstable) RedbirdPy Github (unstable)

📥 Redbird Github🐍 RedbirdPy Github
Created by Qianqian Fang · Supported by NIH R01-CA204443
Powered by Habitat