◆ FORGE Suite
GitHubMechanical Neuroimaging Lab · Univ. of Delaware
Sentinel.jl

Appearance

Your First Inversion ​

In this tutorial we will reconstruct a map of tissue stiffness from magnetic resonance elastography (MRE) displacement data, and print it out in pascals.

We will use a small dataset that ships inside Sentinel, so there is nothing to download. The whole reconstruction takes a few seconds.

Before you start ​

You need Julia 1.11 or newer. Check with:

sh
julia --version

If that command is not found, install Julia with juliaup first — see Getting Started for the installation steps.

Step 1: Get Sentinel ​

The dataset we will use lives inside the repository, so clone it:

sh
git clone https://github.com/mechneurolab/sentinel.git
cd sentinel

Step 2: Start Julia and load Sentinel ​

Start Julia on Sentinel's own environment — the --project=. is what makes Julia use the dependency versions Sentinel was tested with:

sh
julia --project=.

At the julia> prompt, load the package:

julia
using Sentinel

The first time you do this, Julia compiles Sentinel and its dependencies. This takes a minute or two and prints a running list as it goes:

Precompiling packages...
  27768.8 ms  ✓ Sentinel

How many packages appear depends on what your Julia installation has already compiled, so your list will be longer than this one. Wait for the julia> prompt to come back. Subsequent loads take a couple of seconds.

If it doesn't work

ERROR: ArgumentError: Package Sentinel not found means Julia is not using Sentinel's environment — quit and restart with julia --project=. from inside the sentinel directory.

A precompilation error saying some package is required but does not seem to be installed means the dependencies have not been downloaded yet. Run using Pkg; Pkg.instantiate() and then using Sentinel again.

Step 3: Find the bundled data ​

Sentinel ships with a small MRE dataset for testing. Locate it:

julia
data = joinpath(pkgdir(Sentinel), "test", "fixtures", "forge_mre_uiuc.mat")
isfile(data)

You should see:

true

If you see false, you are not running from a clone of the repository — go back to Step 1.

This .mat file holds a displacement field and the geometry it belongs to — the two things an MRE inversion needs. Sentinel will build the finite element mesh from it for us.

Step 4: Describe the inversion ​

Make a directory for the results, then describe the reconstruction we want:

julia
outdir = mkpath("first-inversion")

prob = InversionProblem(data = data,
                        frequency = 60.0,
                        model = :isotropic,
                        output = joinpath(outdir, "recon"),
                        max_global_iter = 3)

Those four settings are: the drive frequency this data was acquired at (60 Hz), the simplest of Sentinel's material models, where to write results, and a cap of three iterations to keep this tutorial quick.

Notice that nothing has been computed yet. InversionProblem only describes the reconstruction; it does not run it.

Step 5: Run the reconstruction ​

julia
result = solve(prob)

Sentinel now generates a mesh from the data and iterates. The output should look something like this:

Inverse solve: 3 global iterations max
  Zone edge: [0.0196, 0.0196, 0.0196], overlap: 15.0%
  3 iteration structures
  Grid: 147 nodes, 122 boundary nodes → numsig target
────────────────────────────────────────────────────────────────────────
  Global iter 1/3 — 1 zones (seed elem 2)...
    Zones: 1/1 (100%) — 3.2s elapsed, ETA 0s
  ── Iter 1/3 done in 3.6s ──
    epsilon=2.1449e-02 (tol=1.0e-03), obj=1.2323e-08, disp_err=1.2323e-08
    mean_shear=3340.9 Pa, zones=1/1 (seed=2) in 3.2s, CG_wt=0.040
  Global iter 2/3 — 1 zones (seed elem 7)...
    Zones: 1/1 (100%) — 0.1s elapsed, ETA 0s
  ── Iter 2/3 done in 0.1s ──
    epsilon=8.2444e-03 (tol=1.0e-03), obj=1.2296e-08, disp_err=1.2296e-08
    mean_shear=3366.4 Pa, zones=1/1 (seed=7) in 0.1s, CG_wt=0.025
  Global iter 3/3 — 1 zones (seed elem 5)...
    Zones: 1/1 (100%) — 0.1s elapsed, ETA 0s
  ── Iter 3/3 done in 0.1s ──
    epsilon=5.5457e-03 (tol=1.0e-03), obj=1.2421e-08, disp_err=1.2421e-08
    mean_shear=3378.0 Pa, zones=1/1 (seed=5) in 0.1s, CG_wt=0.010
────────────────────────────────────────────────────────────────────────
Inverse solve completed 3 iterations (not converged, eps=5.5457e-03)

Look at the mean_shear value on each iteration: 3340.9 Pa, then 3366.4 Pa, then 3378.0 Pa. That is the average stiffness of the reconstruction being refined as the solver works. epsilon is how much the estimate changed on that iteration, and it is falling.

Notice too that the first iteration took about three seconds and the rest took a tenth of a second each — the first one pays for compiling the solver.

Your numbers may differ in the last digit or two depending on your machine.

Stopping at the iteration cap is expected

The last line says not converged. The run stopped because it hit our three-iteration cap, not because anything went wrong — we set that cap in Step 4 to keep the tutorial fast. A production reconstruction runs for many more iterations.

Step 6: Check what came back ​

julia
result.state.converged
false

That is the same "not converged" from above, now as a value you can test in code. The reconstruction itself is perfectly usable.

Step 7: Look at the stiffness map ​

The reconstructed properties live on a mesh. Sample them onto a regular grid so we can look at them as an array:

julia
mesh = result.setup.meshes[result.setup.config.meshind[1, 1]]

grid = interpolate_to_grid(result.material, mesh,
                           collect(mesh.origin), collect(mesh.res), mesh.dims;
                           prop_indices = [1])

mu = grid.data["prop1_real"]
size(mu)
(5, 5, 2)

This dataset is deliberately tiny, so the map is a 5 × 5 × 2 block — small enough to print. Let's look at the first slice, rounded to whole pascals:

julia
round.(Int, mu[:, :, 1])
5×5 Matrix{Int64}:
 3217  3147  3070  3153  3263
 3249  3218  3180  3322  3368
 3432  3558  3618  3724  3587
 3490  3628  3636  3665  3487
 3440  3534  3454  3431  3339

This is the shear modulus map — the stiffness Sentinel recovered at each point, in pascals, roughly 3.1 to 3.7 kPa.

Check the range across the whole volume:

julia
extrema(mu)
(3070.213590426092, 3727.671873216156)

These values are already in pascals; the export helpers apply Sentinel's internal scaling for you. See How to Export Results to MATLAB and ParaView for the details and for the other export paths.

Don't read anything physical into these numbers

The bundled dataset is synthetic — its displacement field is a smooth numerical ramp, not a real scan — and the reconstruction starts from a uniform 3300 Pa guess. So the map you just produced stays within a few percent of that starting value, and its shape is not a measurement of anything. It tells you the pipeline ran correctly, which is what we came here for.

Step 8: Save it for ParaView ​

Write the same reconstruction to a VTK file you can open in ParaView:

julia
export_vtk(joinpath(outdir, "recon"), result.material, mesh; prop_indices = [1])

Check that it landed:

julia
readdir(outdir)
3-element Vector{String}:
 "generated_mesh"
 "recon.config.toml"
 "recon.vtu"

recon.vtu is the stiffness map, generated_mesh/ is the mesh Sentinel built from the .mat file in Step 5, and recon.config.toml is a full record of the settings that produced this run.

What you have done ​

You have run a complete MRE inversion: you took measured displacement data, reconstructed the shear modulus field that explains it, looked at the resulting stiffness map, and exported it for visualization.

That is the same pipeline used for real brain reconstructions. The difference is scale — this dataset had 147 nodes and one zone, where a brain dataset has hundreds of thousands of nodes spread over hundreds of zones.

Where to go next: