walkthrough.qmd

---
title: "Classic vs TRX: A Complete Tractography Walkthrough"
subtitle: "HCP Data · Glasser 360 + Schaefer 4S456 · NODDI scalars · 500k streamlines"
execute:
  echo: true
  eval: true
  warning: false
  error: false
  freeze: auto
---

```{r}
#| label: r-setup
#| include: false

# Working directory for all bash chunks — created here, used throughout
wd <- file.path(dirname(normalizePath(knitr::current_input())), "walkthrough")
dir.create(wd, showWarnings = FALSE)
knitr::opts_knit$set(root.dir = wd)

# Prepend trx-mrtrix build binaries and miniforge3 so bash chunks find them
Sys.setenv(PATH = paste0(
  "/Users/mcieslak/projects/trx-mrtrix/build/bin",
  ":/Users/mcieslak/miniforge3/bin",
  ":", Sys.getenv("PATH")
))

# Number of streamlines — small for a fast demo render; scale up for production
Sys.setenv(N_TRACKS = "1000000")
```

::: {.callout-note}
## Rendering notes

This document executes during render (`eval: true`).  All MRtrix binaries come
from the `trx-mrtrix` build at `build/bin/`, not the system installation.
The demo uses **1M streamlines** for a long render; change `N_TRACKS`
in the R setup chunk to make it faster/slower.
:::

## Timing and resource measurement

Every major command is wrapped with `mtime` — a thin shell function around
`/usr/bin/time -l` (macOS) that emits a one-line summary after each step:

```
┌─ tckgen (TRX)              ─────────────────────────────────────────────┐
│  wall:    52.3s    user:   812.4s    sys:    6.1s    peak mem: 2.31 GB  │
└─────────────────────────────────────────────────────────────────────────┘
```

| Metric | What it measures |
|--------|-----------------|
| **wall** | Elapsed clock time — the number you actually wait |
| **user** | Total CPU time across all threads; `user >> wall` for multithreaded commands |
| **sys** | Kernel time (I/O, memory mapping) |
| **peak mem** | Physical memory footprint (`phys_footprint` via `TASK_VM_INFO`) — excludes clean read-only mmap pages, so TRX files are not double-counted |

::: {.callout-note collapse="true"}
## Why not RSS?

`maximum resident set size` (RSS) counts clean file-backed mmap pages — TRX
files are opened read-only via mmap, so RSS is inflated for TRX inputs even
though those pages impose no memory pressure.  **`peak memory footprint`**
counts only dirty/anonymous pages and matches what Activity Monitor shows.
:::

---

## Setup

```{bash}
#| label: setup

DERIV=/Users/mcieslak/data/qsirecon/derivatives/qsirecon-MRtrix3_act-HSVS/sub-100307/dwi
PARC=/Users/mcieslak/data/qsirecon/sub-100307/dwi

source ../timing_utils.sh

ln -sf "$DERIV/sub-100307_space-T1w_model-msmtcsd_param-fod_label-WM_dwimap.mif.gz" wm_fod.mif.gz
ln -sf "$DERIV/sub-100307_space-T1w_model-mtnorm_param-inliermask_dwimap.nii.gz"     mask.nii.gz
ln -sf "$PARC/sub-100307_space-T1w_model-noddi_param-icvf_dwimap.nii.gz"             icvf.nii.gz
ln -sf "$PARC/sub-100307_space-T1w_model-noddi_param-isovf_dwimap.nii.gz"            isovf.nii.gz
ln -sf "$PARC/sub-100307_space-T1w_seg-Glasser_dseg.mif.gz"                          glasser.mif.gz
ln -sf "$PARC/sub-100307_space-T1w_seg-Glasser_dseg.txt"                             glasser.txt
ln -sf "$PARC/sub-100307_space-T1w_seg-4S456Parcels_dseg.mif.gz"                     4S456.mif.gz
ln -sf "$PARC/sub-100307_space-T1w_seg-4S456Parcels_dseg.txt"                        4S456.txt

echo "N_TRACKS = $N_TRACKS"
echo "tckgen:   $(which tckgen)"
echo "trxlabel: $(which trxlabel)"
ls -lh *.mif.gz *.nii.gz *.txt
```

---

## Step 1 — Generate streamlines

Both pipelines start from a common `tckgen` run.  MRtrix selects the writer
from the output filename extension, so the only difference is `.tck` vs `.trx`.

:::: {.columns}

::: {.column width="49%"}
<span class="badge-classic">CLASSIC</span>

```{bash}
#| label: tckgen-classic
source ../timing_utils.sh
mtime "tckgen (classic)" \
tckgen wm_fod.mif.gz       \
    -algorithm iFOD2        \
    -seed_image mask.nii.gz \
    -mask       mask.nii.gz \
    -minlength  30          \
    -maxlength  250         \
    -select     $N_TRACKS   \
    -nthreads   8           \
    -force                  \
    tracks.tck
```

```{bash}
tckinfo tracks.tck
du -sh  tracks.tck
```
:::

::: {.column width="49%"}
<span class="badge-trx">TRX</span>

```{bash}
#| label: tckgen-trx
source ../timing_utils.sh
mtime "tckgen (TRX float16)" \
tckgen wm_fod.mif.gz       \
    -algorithm iFOD2        \
    -seed_image mask.nii.gz \
    -mask       mask.nii.gz \
    -minlength  30          \
    -maxlength  250         \
    -select     $N_TRACKS   \
    -nthreads   8           \
    -trx_float16            \
    -force                  \
    tracks_f16.trx
```

```{bash}
tckinfo tracks_f16.trx
du -sh  tracks_f16.trx
```
:::

::::

The float16 TRX is roughly **half the size** of the TCK at the same streamline
count.  For the remainder of this walkthrough we use `tckconvert` to produce a
float32 TRX from the same TCK, so both pipelines operate on **byte-for-byte
identical streamlines** and any output differences are purely due to format, not
stochastic variation.

```{bash}
#| label: tckconvert-trx
source ../timing_utils.sh
mtime "tckconvert → TRX" \
tckconvert tracks.tck tracks.trx -force
```

```{bash}
tckinfo tracks.trx
du -sh  tracks.trx
```

---

## Step 2 — SIFT2 weighting

SIFT2 calibrates per-streamline weights to match the WM FOD amplitudes.  In the
classic pipeline the weights land in a separate CSV.  In the TRX pipeline they
are embedded as a `data_per_streamline` field named `weights`.

:::: {.columns}

::: {.column width="49%"}
<span class="badge-classic">CLASSIC</span>

```{bash}
#| label: sift2-classic
source ../timing_utils.sh
mtime "tcksift2 (classic)" \
tcksift2                   \
    tracks.tck             \
    wm_fod.mif.gz          \
    weights_classic.csv    \
    -nthreads 8            \
    -force
```

```{bash}
wc -l weights_classic.csv
du -sh tracks.tck weights_classic.csv
```
:::

::: {.column width="49%"}
<span class="badge-trx">TRX</span>

```{bash}
#| label: sift2-trx
source ../timing_utils.sh
mtime "tcksift2 (TRX)" \
tcksift2               \
    tracks.trx         \
    wm_fod.mif.gz      \
    weights            \
    -nthreads 8        \
    -force
```

```{bash}
tckinfo tracks.trx
```
:::

::::

---

## Step 3 — Sample NODDI scalar maps

NODDI provides two microstructure metrics:

- **ICVF** (intracellular volume fraction) — neurite density proxy
- **ISOVF** (isotropic volume fraction) — free water content

`tcksample` interpolates a volumetric image at every streamline vertex, yielding:

- **dpv** (per-vertex) — full sampled profile; used for mrview coloring and along-tract statistics
- **dps mean** (per-streamline) — one number per streamline; useful as a connectivity weight or covariate

:::: {.columns}

::: {.column width="49%"}
<span class="badge-classic">CLASSIC</span>

```{bash}
#| label: tcksample-classic
source ../timing_utils.sh
mtime "tcksample icvf dpv" \
tcksample tracks.tck icvf.nii.gz icvf.tsf -force

mtime "tcksample isovf dpv" \
tcksample tracks.tck isovf.nii.gz isovf.tsf -force

mtime "tcksample icvf mean" \
tcksample tracks.tck icvf.nii.gz icvf_mean.txt -stat_tck mean -force

mtime "tcksample isovf mean" \
tcksample tracks.tck isovf.nii.gz isovf_mean.txt -stat_tck mean -force

du -sh icvf.tsf isovf.tsf icvf_mean.txt isovf_mean.txt
```
:::

::: {.column width="49%"}
<span class="badge-trx">TRX</span>

```{bash}
#| label: tcksample-trx
source ../timing_utils.sh
mtime "tcksample icvf dpv" \
tcksample tracks.trx icvf.nii.gz icvf

mtime "tcksample isovf dpv" \
tcksample tracks.trx isovf.nii.gz isovf

mtime "tcksample icvf mean" \
tcksample tracks.trx icvf.nii.gz icvf_mean -stat_tck mean

mtime "tcksample isovf mean" \
tcksample tracks.trx isovf.nii.gz isovf_mean -stat_tck mean

tckinfo -prefix_depth 2 tracks.trx
```
:::

::::

---

## Step 4 — Visualize results in mrview

`mrview` loads TRX files directly.  The SIFT2 weights and NODDI scalars embedded
in the previous two steps are immediately available as colormap sources — no
sidecar files to locate or load separately.

Open `tracks.trx` in `mrview → Tools → Tractography → Colour → Scalar file → TRX field…`
and select `icvf`, `isovf`, or `weights` for per-vertex or per-streamline coloring.
If a field name exists in both dps and dpv, disambiguate using `dps:<name>` or `dpv:<name>`.

Command-line capture via `-capture.folder`, `-capture.prefix`, and `-capture.grab`
works with TRX inputs (requires a live GUI session):

```{bash}
mrview mask.nii.gz \
  -mode 2 \
  -imagevisible 0 \
  -tractography.load tracks.trx \
  -tractography.geometry lines \
  -tractography.opacity 1.0 \
  -tractography.thickness 0.25 \
  -tractography.slab -1 \
  -tractography.trx_scalar dps:weights \
  -tractography.tsf_range 0.15,2.5 \
  -capture.folder . \
  -capture.prefix trx_sift_overlay \
  -capture.grab \
  -exit
```

```{bash}
#| label: mrview-capture-and-embed
#| results: asis
source ../timing_utils.sh
rm -f trx_sift_overlay*.png mrview_trx_overlay.png mrview_trx_overlay_fallback.svg
rm -f ../mrview_trx_overlay.png ../mrview_trx_overlay_fallback.svg

set +e
mrview mask.nii.gz \
  -mode 2 \
  -imagevisible 0 \
  -tractography.load tracks.trx \
  -tractography.geometry lines \
  -tractography.opacity 1.0 \
  -tractography.thickness 0.25 \
  -tractography.slab -1 \
  -tractography.trx_scalar dps:weights \
  -tractography.tsf_range 0.15,2.5 \
  -capture.folder . \
  -capture.prefix trx_sift_overlay \
  -capture.grab \
  -exit
mrview_status=$?
set -e

capture_file=""
for f in trx_sift_overlay*.png; do
  if [ -f "$f" ]; then
    capture_file="$f"
    break
  fi
done

if [ -n "$capture_file" ]; then
  cp "$capture_file" mrview_trx_overlay.png
  cp "$capture_file" ../mrview_trx_overlay.png
  echo "Captured mrview image: \`$capture_file\`"
  echo ""
  echo "![TRX streamlines coloured by embedded SIFT2 weights](mrview_trx_overlay.png)"
fi
```

---

## Step 5 — Atlas labeling and connectome construction

`tck2connectome` assigns streamlines to node pairs on-the-fly — the assignment
is ephemeral and must be recomputed for every new atlas or metric.

`trxlabel` embeds the assignment as TRX groups permanently.  Multiple atlases
can be labeled in a **single pass** over the tractogram by repeating `-nodes`,
`-lut`, and `-prefix` — geometry is read exactly once regardless of how many
atlases are requested.

### 4a — Classic: two separate tck2connectome runs

The classic pipeline must reread all geometry and redo the radial search for each atlas.

```{bash}
#| label: tck2connectome-glasser-classic
source ../timing_utils.sh
mtime "tck2connectome Glasser" \
tck2connectome                  \
    tracks.tck                  \
    glasser.mif.gz              \
    connectome_glasser_classic.csv \
    -tck_weights_in weights_classic.csv \
    -assignment_radial_search 2 \
    -symmetric                  \
    -zero_diagonal              \
    -force
```

```{bash}
#| label: tck2connectome-4S456-classic
source ../timing_utils.sh
mtime "tck2connectome 4S456" \
tck2connectome                  \
    tracks.tck                  \
    4S456.mif.gz                \
    connectome_4S456_classic.csv \
    -tck_weights_in weights_classic.csv \
    -assignment_radial_search 2 \
    -symmetric                  \
    -zero_diagonal              \
    -force
```

### 4b — TRX: both atlases in a single trxlabel pass

`trxlabel` accepts multiple `-nodes`/`-lut`/`-prefix` triplets and labels all
atlases in one pass over the tractogram — no geometry is read twice.

```{bash}
#| label: trxlabel-both
source ../timing_utils.sh
mtime "trxlabel Glasser + 4S456" \
trxlabel                         \
    tracks.trx tracks.trx        \
    -nodes  glasser.mif.gz       \
    -nodes  4S456.mif.gz         \
    -lut    glasser.txt          \
    -lut    4S456.txt            \
    -prefix glasser              \
    -prefix 4S456                \
    -assignment_radial_search 2  \
    -force
```

```{bash}
#| label: tckinfo-after-label
tckinfo -prefix_depth 2 tracks.trx
```

```{bash}
#| label: trx2connectome-glasser
source ../timing_utils.sh
mtime "trx2connectome Glasser" \
trx2connectome                   \
    tracks.trx                   \
    connectome_glasser_trx.csv   \
    -tck_weights_in weights      \
    -out_node_names glasser_nodes.txt \
    -group_prefix glasser        \
    -lut glasser.txt             \
    -symmetric                   \
    -zero_diagonal               \
    -force
```

```{bash}
#| label: trx2connectome-4S456
source ../timing_utils.sh
mtime "trx2connectome 4S456" \
trx2connectome                 \
    tracks.trx                 \
    connectome_4S456_trx.csv   \
    -tck_weights_in weights    \
    -out_node_names 4S456_nodes.txt \
    -group_prefix 4S456        \
    -lut 4S456.txt             \
    -symmetric                 \
    -zero_diagonal             \
    -force
```

```{bash}
#| label: tckinfo-final
tckinfo -prefix_depth 2 tracks.trx
```

### 4c — TRX combined-atlas connectome in one matrix

Because both atlas group sets are embedded in one TRX file, `trx2connectome`
can emit a single matrix that includes:

- within-atlas edges (Glasser-Glasser, 4S456-4S456), and
- cross-atlas edges (Glasser-4S456),

without rerunning any node assignment.

```{bash}
#| label: trx2connectome-combined
source ../timing_utils.sh
mtime "trx2connectome combined-atlas matrix" \
trx2connectome \
    tracks.trx \
    connectome_combined_trx.csv \
    -tck_weights_in weights \
    -out_node_names combined_nodes.txt \
    -symmetric \
    -zero_diagonal \
    -force
```

```{bash}
#| label: combined-connectome-summary
python3 - <<'EOF'
import collections
import itertools
import numpy as np

mat = np.loadtxt("connectome_combined_trx.csv", delimiter=",")
names = [line.strip() for line in open("combined_nodes.txt") if line.strip()]

def atlas_prefix(name: str) -> str:
    return name.split("_", 1)[0] if "_" in name else "(no_prefix)"

if mat.shape[0] != len(names):
    raise RuntimeError(f"matrix shape {mat.shape} does not match node-name list length {len(names)}")

prefix_counts = collections.Counter(atlas_prefix(n) for n in names)
print(f"combined matrix shape: {mat.shape}")
print("node counts by atlas prefix:")
for k in sorted(prefix_counts):
    print(f"  {k}: {prefix_counts[k]}")

indices = {k: [i for i, n in enumerate(names) if atlas_prefix(n) == k] for k in sorted(prefix_counts)}
print("non-zero edge counts by atlas block:")
for a, b in itertools.combinations_with_replacement(sorted(indices), 2):
    block = mat[np.ix_(indices[a], indices[b])]
    nnz = int((block > 0).sum())
    total = float(block.sum())
    print(f"  {a:>8s} x {b:<8s}: nnz={nnz:6d}, sum={total:.3f}")
EOF
```

---

## Step 6 — Verify equivalence

With `-lut`, `trx2connectome` orders rows and columns by numeric node ID —
the same ordering `tck2connectome` uses.  The matrices are directly comparable
with no reordering.

```{bash}
#| label: diff-connectomes
python3 - <<'EOF'
import numpy as np, sys
ok = True
for atlas in [("Glasser", "connectome_glasser_classic.csv", "connectome_glasser_trx.csv"),
              ("4S456",   "connectome_4S456_classic.csv",   "connectome_4S456_trx.csv")]:
    name, cf, tf = atlas
    a = np.loadtxt(cf, delimiter=",")
    b = np.loadtxt(tf, delimiter=",")
    if a.shape != b.shape:
        print(f"{name}: SHAPE MISMATCH  classic={a.shape}  trx={b.shape}", file=sys.stderr)
        ok = False
        continue
    diff = np.abs(a - b)
    print(f"{name} {a.shape}  max|diff|={diff.max():.2e}  mean|diff|={diff.mean():.2e}")
    if diff.max() > 0.5:
        ok = False
        print(f"  FAIL: differences too large", file=sys.stderr)
if ok:
    print("All matrices match (within rounding).")
else:
    sys.exit(1)
EOF
```

```{bash}
#| label: heatmaps
python3 - <<'EOF'
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

pairs = [
    ("Glasser 360",    "connectome_glasser_classic.csv", "connectome_glasser_trx.csv"),
    ("Schaefer 4S456", "connectome_4S456_classic.csv",   "connectome_4S456_trx.csv"),
]

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
for row, (name, cf, tf) in enumerate(pairs):
    classic = np.loadtxt(cf, delimiter=",")
    trx     = np.loadtxt(tf, delimiter=",")
    for col, (mat, title) in enumerate([(classic, f"{name}\nClassic"), (trx, f"{name}\nTRX")]):
        ax = axes[row][col]
        im = ax.imshow(np.log1p(mat), cmap="inferno", aspect="auto")
        ax.set_title(title)
        ax.set_xlabel("Node")
        ax.set_ylabel("Node")
        plt.colorbar(im, ax=ax, label="log(1 + weight)")

plt.tight_layout()
plt.savefig("../connectome_comparison.png", dpi=150, bbox_inches="tight")
print("Saved connectome_comparison.png")
EOF
```

![Connectome comparison](connectome_comparison.png)

---

## Step 7 — File inventory

```{bash}
#| label: file-inventory
echo "=== Classic: files required to reproduce either connectome ==="
ls -lh tracks.tck weights_classic.csv icvf.tsf isovf.tsf icvf_mean.txt isovf_mean.txt \
        connectome_glasser_classic.csv connectome_4S456_classic.csv

echo ""
echo "=== TRX: everything in one file ==="
ls -lh tracks.trx connectome_glasser_trx.csv connectome_4S456_trx.csv
```

---

## What's next

::: {.callout-tip}
### Other commands that read/write TRX fields

| Command | TRX capability |
|---------|---------------|
| `tckedit` | Filters TRX by ROI, length, or count; remaps dps/dpv/groups to surviving streamlines |
| `tcksift` | Produces a TRX subset via `subset_streamlines` — all metadata preserved |
| `fixel2tsf` | Pass a bare field name as output to embed fixel scalars as dpv directly in TRX |
| `tsfinfo` / `tsfvalidate` | Inspect and validate TRX dpv fields |
| `mrview` | Loads TRX directly; colour by group, dps, or dpv field from the GUI |
:::