Transcranial ultrasound simulation in macaque
Subject zach; skull density from PETRA, simulation in BabelBrain
Introduction
Transcranial ultrasound is planned on a simulation, and the simulation needs a density at every voxel it passes through. The skull is where that matters: bone’s density sets the speed of sound and the attenuation, and so sets where the beam lands and how much of it arrives. CT measures density directly, and this animal has no CT. What it has is two MR images, a T1-weighted scan, and a PETRA, a sequence that reads the signal early enough to still see bone. This report is the route from those two images to a density map of the head, and from that map to a simulated pressure field at each target.
Processing
Input data
We start with two scans of the same animal, taken weeks apart. The PETRA came off the scanner as DICOM (one small file per slice) so the first job is to stack that folder into a single NIfTI volume and check that the geometry survived.
# Two images. The PETRA arrived as DICOM; the T1 as NIfTI.
petra = dicom_to_nifti(petra_dicom) # step 00
t1 = load_nifti(t1_file)Then a relabelling. The scanner assumes its default bore position (a person lying on their back in the tunnel) and writes the header to match. This animal is a macaque positioned in the sphinx position in a stereotaxic frame, with its head aligned along the scanner bore. Swapping the labels fixes it. No voxel value changes, and nothing is resampled.
# The scanner tagged the animal as a human lying supine. Relabel the axes so
# that "anterior" is the animal's nose. Nothing is resampled.
petra, t1 = reorient_to_animal(petra, t1, rotation) # step 01| PETRA | T1w | |
|---|---|---|
| Sequence | *Petra3d1 (noT1_petra_tra_FA1_S2_DIS3D) |
not recorded |
| Scanner | SIEMENS Prisma_fit, 3 T | not recorded |
| Echo time (TE) | 0.07 ms | not recorded |
| Repetition time (TR) | 3.61 ms | not recorded |
| Flip angle | 1° | not recorded |
| Matrix | 320 × 320 × 320 | 512 × 512 × 256 |
| Voxel size | 0.75 mm | 0.5 mm |
| Field of view | 240 × 240 × 240 mm | 256 × 256 × 128 mm |
| Coverage | whole head, clearing every face by at least 46 mm | a slab, 128 mm thick |
| Distortion correction | applied by the vendor (DIS3D/DIS2D) | unknown |
| Intensity range | 0–4095 | — |
The two volumes sit 59.8 mm apart in scanner coordinates — different table positions, weeks apart — so their stored headers cannot be used to bring them together, and step 03 has to start the registration from the centres of mass instead.
The T1’s acquisition parameters are blank because it was delivered as a reconstructed NIfTI with no DICOM tags. Its echo time is not needed anywhere in the pipeline — only the PETRA’s short one earns the bone signal — but whether gradient nonlinearity correction was applied to it is unknown, and no later step can find out.
Bias-field correction
The receive coil sits nearer some parts of the head than others, so the same tissue reads brighter where the coil is close. The map from brightness to density inside bone is a straight line, so bone under the coil would come out too light. We estimate the field with N4 over the whole head, skull included, and divide it out of both images.
# The receive coil is brighter near itself. Estimate that smooth field and
# divide it out, over the whole head including the skull.
petra, petra_field = n4_bias_correction(petra, head_mask(petra)) # step 02
t1, t1_field = n4_bias_correction(t1, head_mask(t1))| Shrink factor | 4 |
| Spline order | 3 |
| Initial spline distance | 120 mm |
| Levels | 5, finest spline distance 7.5 mm |
| Iterations per level | 50, 40, 30, 20, 10 |
| Histogram bins | 200 |
The spline’s knot spacing is the only free choice in the fit. The histogram sharpens the finer that spacing gets, so no measure of uniformity can pick one. We swept it and read the field inside the dark shell of the skull against the field in the tissue on either side of it: that ratio stays at one while the spline is too coarse to see the skull, and drops once it starts to follow it.
60 mm 1.026, 30 mm 1.055, 15 mm 1.047, 7.5 mm 1.023, 3.75 mm 0.975.
The finest spacing whose ratio is still at or above one is 7.5 mm, and that is the one the step uses.
The refit runs at a finest spacing of 15 mm, the coarsest that makes the two sides agree, on three classes divided by their own typical values before the fit (brain 1167, fluid 1173, scalp 1186):
- Scalp beside the vault bone, left over right: 1.108 before, 1.015 after.
- Bone over brain at the frontal pole: 0.618 before, 0.572 after.
Then a second pass. After the first, the brain beside the skull reads the same on both sides, but the scalp beside it is still brighter on the coil’s side. We refit N4 on the soft tissue alone, with bone, implant and air out of the fit and each tissue class divided by its own median first, and take the coarsest spline that brings the two sides together.
# The soft tissue beside the bone should read the same left and right. If it
# does not, the coil's field is still in the image: refit it on the soft
# tissue alone, with the bone excluded, and redo the segmentation once.
if not left_right_symmetric(petra, labels): # step 05
petra = n4_on_soft_tissue(petra, labels, exclude=bone | implant | air) # step 02b
skull, suspect, labels, bone, implant, air = (steps 04 and 05 again, once)The bone next to the formerly bright scalp comes down relative to the brain, which is correct if that brightness was the coil’s, and we read it so because the excess stood on both sides of the head, not one.
Below, the three images on one grid with each pass’s field beside them: what the first pass takes out of the head, and what the refit changes in the last millimetres under the scalp.
Registration
Then the T1 onto the PETRA’s grid. The two scans were taken weeks apart at different table positions, and their stored headers do not place one on the other. The skull does not deform between sessions, so the pose is six parameters. The contrasts differ, so we cannot match intensities: we maximise mutual information instead, rigid, and over the cranium alone, because the jaw and the neck did not hold still. The search starts from the two centres of mass, runs in two independent tools, and keeps the pose with the lower cost.
That is not yet enough. Mutual information’s optimum on this pair is flat, and the two tools stop a few millimetres apart. So we finish on the outer surface of the skull, the one boundary sharp in both contrasts and rigid between sessions: along each surface normal we read where the T1’s bone-to-scalp edge crosses, and move the pose until that crossing sits on the PETRA’s own. The T1 is then resampled once, onto the PETRA’s grid, and every later step works there. The PETRA is not touched.
# Put the T1 on the PETRA's grid. Mutual information gets within a few
# millimetres; the outer table of the skull, sharp in both, finishes it.
pose = register_rigid(t1, petra, cost="mutual information") # step 03
pose = refine_on_outer_table(t1, petra, pose)
t1_in_petra = resample(t1, onto=petra, pose)Do the skull’s surfaces coincide in the two contrasts?
Segmentation
BabelBrain reads one label volume, in SimNIBS’s numbering, and takes four things from it: the head, the skull’s outer and inner surfaces, and the bone between. No tool writes that volume for a macaque, so we write it.
We extract the brain with DeepBet, a network trained on macaque T1 scans, and bring its mask onto the PETRA’s grid with the pose the registration found. The classes inside the brain are the NMT template’s, warped in, because BabelBrain reads grey matter, white matter and fluid as one region and the source that is right by construction wins.
# Which voxels are what. Brain from a trained network on the T1; skull grown
# outward from the brain cavity through the PETRA's dark shell; tissue
# classes from a macaque template; everything dark that is not bone flagged.
brain = deepbet(t1) # step 04
brain = resample(brain, onto=petra, pose)
skull, suspect = grow_skull_from_cavity(petra, brain)
classes = warp_template_classes(nmt_template, t1_in_petra)
labels = assemble_labels(head_mask(petra), classes, skull)We find the skull by growing it outward from the brain cavity through the PETRA’s dark shell, so that position, and not brightness alone, is what separates bone from the material lying on it.
A voxel joins the skull when all four hold:
brightness 0.12 <= PETRA < 0.75 of the soft-tissue peak
side outside the intracranial space
depth no more than 12 mm from it
path a chain of voxels meeting the three above reaches
the surface of the intracranial spaceand then, once, the shell is closed and grown again:
closing closed by 3 mm; a voxel the closing encloses joins
if PETRA < 1.00, side and depth still holdingThe band is bounded at both ends. Below it is air, and an air pocket must stop the growth rather than join it. Above it is soft tissue: muscle, fat and skin are bright, so the flood cannot cross them.
Side and depth are measured from the intracranial space, which is the brain mask grown 4 mm through anything at or above 0.75 of the peak, so it stops on bone and its surface sits at the inner table. Where the skull is open, at the orbits and the foramen magnum, the reach is what bounds the leak.
The path condition is the one that makes this position rather than intensity. The cement over the vault is as dark as bone and touches it, so a chain does reach it: it is taken, and flagged for the next step. The fibrous layer over the temporalis is equally dark and is never reached, because no chain of dark voxels crosses bright muscle.
Marrow and open sutures read above the band and would halt the growth before the outer table. The closing lets the shell cross them, and accounts for 11.6% of the label; skull.t1_bright is what tests whether it overreached.
Bone no path from the cavity reaches, the face, the zygomatic arches and the jaw, is taken by a stricter threshold, barred from the vault. A threshold over the whole head was the first version of this rule: it took the cement and the fibrous layer for bone, bridged the temporalis between them, returned several times the skull the anatomy shows, and passed every numeric check.
The segmentation reads the PETRA the second bias pass refines, and that pass is fitted on these labels, so both run twice: the residual field had brightened the left vault past the value the growth stops at, leaving holes only the corrected image closes.
petra = n4_on_soft_tissue(petra, labels, exclude=bone | implant | air) # step 02b
skull, suspect, labels, bone, implant, air = (steps 04 and 05 again, once)The flagged material is then sorted by position: what sits immediately above the cavity is bone whatever it reads, beyond that material darker than bone is the cap, and only the largest piece counts. It is a sheet of cement over the back of the vault, exposed to air along its top and stopping short of the frontal bone, where the prefrontal trajectories enter; the post itself stands in air and reads as air, so only its bed is visible. Three classes and not two, because bone, cement and air read the same grey in a PETRA and go three different ways in the simulation.
bone, implant, air = delineate_implant(labels, petra, suspect) # step 05Does the bone label follow the thin dark shell between the brain and the scalp and nothing else, and where does the cap sit on it?
Density map
Every published mapping from PETRA to Hounsfield units is a straight line in normalised PETRA, the image divided by one number, the intensity of soft tissue. To calculate that number we divide by the mode of the PETRA’s histogram over the soft-tissue labels, brain and scalp with the implant excluded, which is the convention the mapping we apply was fitted with. A mode and not a mean, because the mask has a bright tail, voxels saturated beside the coil, and a dark tail, the head’s rim averaged with air.
We then applied the line inside the bone mask only.
# Brightness to density. Divide by the soft-tissue peak, map through a
# published line inside bone, then through a published table to density.
petra_n = petra / soft_tissue_peak(petra, labels) # step 06
hu = pseudo_ct(petra_n, bone, line="cph2025") # step 07We then turned those Hounsfield units into density with a published table of twelve points, and handed BabelBrain the density map rather than the pseudo-CT. Its CT input applies no calibration: it divides by the brightest bone voxel in the volume, so the same skull comes out differently from one of our reruns to the next. Its PETRA input redoes the normalisation itself. A density map is the only input that keeps our own calibration.
density = hu_to_density(hu, table="cph2025") # step 08
density[implant] = 1190 # acrylic, assumed
density[air] = 1.2BabelBrain converts that map back to Hounsfield units internally and takes speed of sound and attenuation from there, so the round trip is what reaches the physics, and we measure what it costs. Air carries the table’s own air value and not zero, because zero falls below the window BabelBrain looks for air in, and a cavity it does not detect is modelled as transparent tissue.
Normalisation
| The soft-tissue value | 1256 intensity units, width 161 (13% of it), from a histogram of 200 bins one of which is 1.0% of the value |
| One peak, or several? | one soft-tissue peak: the tallest other peak at least a tenth of the location away reaches 3% of it (limit 50%) |
| Where the other rules land | legacy_100_bins +0.0%, head_dominant_peak +0.3%, babelbrain_two_peaks +4.1%, soft_tissue_median +0.1% (limit 3%; one bin of the chosen histogram is 1.0%) |
| BabelBrain’s own rule | 134 on the whole volume (-89%) and 109 on the head’s bounding box (-91%), both of them background; 1308 (+4.1%) only with everything outside the head set to zero, which is the volume step 09 must hand it |
| Where the bone reads | over the vault, 0.52 of the value at the median and 0.27 to 0.85 between the 5th and 95th percentiles, inside the 0.1 to 0.9 the published mappings were fitted over |
Pseudo-CT
| The line applied | cph2025: HU = -2080 × normalised PETRA + 2133, never below soft tissue at 28.6 HU; background and enclosed air at -1000 HU |
| Bone in Hounsfield units | 1019 HU at the median and 296 to 1569 between the 5th and 95th percentiles, over 63.4 cm³ |
| A second published line exists | ucl: HU = -2930 × normalised PETRA + 3274, fitted on the same acquisition protocol and also on humans. At this animal’s bone 1019 HU (cph2025) and 1705 HU (ucl), 686 HU apart. No data of this animal’s can settle which is right, so the sensitivity sweep below prices the choice at the focus rather than this step resolving it |
| Written as air rather than mapped | 2.6 cm³ of the bone labels, 2.4 cm³ of it below the vault: 0.12 cm³ by the floor on the normalised value and 2.5 cm³ by the mask read on the raw image |
| How close that comes to an entry | frontal bone, over the orbits 0.00 cm³, temporal squama 0.01 cm³, temporal squama, mirrored to the left 0.00 cm³, within 10 mm (limit 0.05 cm³) |
| Where it does lie | 0.75 mm from the intracranial space at its nearest and 1.8 cm³ within 10 mm of it; the largest piece is 1.40 cm³ at the midline of the skull base, beneath the brain |
| Every bone voxel given a value? | yes: step 05’s mask (62.4 cm³) and step 04’s other bone, the face, the arches and the jaw (3.6 cm³), with 0 cm³ stranded |
Density
| The calibration | ct_to_density_calibration_cph2025_v1.csv, 12 points from -950 HU (1.2 kg/m³) to 1631 HU (2150 kg/m³), taken byte for byte from petra2density/maps rather than from BabelBrain, whose checkout ships only ct_to_density_calibration_cph2025_line_v1.csv (3 points), -57.5 kg/m³ at -1000 HU |
| Bone in density | 1669 kg/m³ at the median and 1220 to 2108 between the 5th and 95th percentiles |
| What the round trip costs | a median +3.8 HU and at most 28.5 (limit 100) over the 97.7% of bone the calibration interpolates; BabelBrain’s inverse carries 7 of this table’s own density anchors exactly and differs only in the Hounsfield value it puts on each |
| Air | 1.2 kg/m³, inside the 0.01 to 10 kg/m³ window BabelBrain looks for air in, along with 100% of the enclosed air cells; zero falls below that window and would be modelled as transparent tissue |
| The implant | 1190 kg/m³ from its own mask, 132 above the single value soft tissue carries here (1058) and below BabelBrain’s 1200 kg/m³ bone threshold |
| What the top-end clamp costs | 1.49 cm³ of bone (2.35%) reads above the calibration’s last point (1631 HU) and is held at 2150 kg/m³, which BabelBrain inverts to the single value 1660 HU: up to 223 HU, though 34% of it comes back higher than it went in |
| What BabelBrain’s CT input would have applied | density = 1000 + 1700 × HU / max(HU), normalised by the brightest bone voxel in the volume (1883 HU here): 1920 kg/m³ at the median against 1669, and a brightest voxel of 2100 HU instead would move it to 1825 |
The density map, and the two images it came from.
Simulation
Trajectories and BabelBrain’s inputs
We planned two sonications, one prefrontal target on each side. The targets were picked on the T1 after it had been resampled onto the PETRA’s grid, so nothing stands between the point we chose and the point the simulation focuses on. The two entry points on the skin were chosen as a mirror-symmetric pair, so that the two sides are insonated comparably.
| pfc-right | pfc-left | |
|---|---|---|
| Target (R, A, S) | -10.6, -0.5, 76.1 mm | -32.2, -0.5, 76.1 mm |
| Entry (R, A, S) | -1.9, 4.8, 92.7 mm | -40.5, 4.6, 92.1 mm |
| Bone crossed | 1.9 mm at 8° from normal | 1.8 mm at 6° from normal |
| Tissue to the target | 19.5 mm | 18.6 mm |
| Crosses the implant | no | no |
| Crosses enclosed air | no | no |
| Inside the transducer’s focal range | yes | yes |
BabelBrain reads four things, and each has to arrive in a particular shape.
| Input | What it carries | The form it has to be in |
|---|---|---|
| Label volume | which tissue each voxel is, in SimNIBS’s numbering (0, 1, 2, 3, 5, 7, 8, 9). BabelBrain meshes the head, the outer skull and the inner skull from it | four-dimensional, with a trailing axis of length one: 320 × 320 × 320 × 1 |
| T1 | the frame the simulation domain is built on. Its origin is unchanged, so the labels and the density map share it and no coregistration is asked for | one millimetre isotropic, which it asserts and will not work around: 1.0 × 1.0 × 1.0 mm |
| Density map | kilograms per cubic metre at every voxel. This, and not the label volume, is what makes a voxel bone | the window it looks for air in has to be given explicitly, 0.01 to 10 kg/m³; its own default is in Hounsfield units and matches nothing here |
| Trajectory | one per target, 2 in all | a four by four matrix, whose translation and third column both mean something particular; the last box in this section says what |
# What BabelBrain reads: the labels in its numbering, a 1 mm T1, the density
# map, and one trajectory per target.
inputs = write_babelbrain_inputs(labels, density, t1_in_petra, trajectories) # step 09
domain = babelbrain_step1(inputs) # its skull on its grid; air regions rewrittenBabelBrain does not simulate the map we hand it. It builds its own domain on a finer grid, takes the head’s geometry from the labels and the bone from the density map, and looks for enclosed air itself. Skull that its filter drops is not left as a hole: it is written as brain, so the beam passes through as if no bone were there. That error is invisible in the pressure field and it always makes the simulation optimistic, which is why we measure what survived rather than assume it.
The trajectory’s convention. A trajectory is a four by four matrix, and two things about it are easy to get backwards. Its translation is the target rather than the point on the skin, so putting the entry there would focus the beam at the scalp. Its third column points at the transducer rather than along the beam, and the wrong sign builds a domain that looks entirely healthy, with the target in the right place and plausible tissue around it, and the transducer on the far side of the head. Reading the matrix back cannot catch that, so we walk outward from the focal point instead and require the short way out of the head to be the entry we chose.
Its own grid. BabelBrain sets its voxel size from the slowest wave it will model, divided by the frequency and by the points it wants per wavelength. Here that is 0.49 mm, finer than the 0.75 mm grid everything upstream lives on. The consequence is not cosmetic: the filter it uses on the skull is seven voxels across, so it spans 3.4 mm on its grid where the same filter would span 5.25 mm on ours. Emulating that filter on a coarser grid, as the previous step did, overstates how much bone it removes.
The skull. It does not take our bone label. It thresholds the density map, applies that filter, closes the result by five millimetres and keeps the largest piece. What comes out holds 66.4 cm³ against the 63.4 cm³ we sent, which looks like a gain until the two are compared voxel by voxel: it keeps 87% of our bone, adds 11.1 cm³ and loses 8.1 cm³, for an overlap of 0.85. The bone it loses is thin bone at the skull base, which is what a majority vote in a cube does regardless of how dense that bone is.
The bone at the entries. This is the measurement the step exists for, because the entries were chosen through the thinnest bone that met the other constraints, and a filter of that size can erase a sheet thinner than itself. It did not happen: pfc-right survives at 2.3 mm where we sent 1.9; pfc-left survives at 2.4 mm where we sent 1.8.
The enclosed air. Its own detector found essentially nothing, 0.002 cm³ against the 2.60 cm³ we found. The reason is the grid again. Before anything looks for air, the density map is interpolated onto that finer grid, and a cavity a few voxels across is smoothed until it no longer reads as air at all: that accounts for 2.46 cm³ of the 2.61 cm³, with the detector’s own filtering taking most of the rest. Because the simulation treats those cavities as reflectors, we wrote that file ourselves from our own masks and kept BabelBrain’s answer beside it. This is the one difference here that we repaired rather than measured.
The implant. We wrote the cap as soft tissue, at a density just below the value BabelBrain calls bone, so the expectation was that none of it would be simulated as skull. In fact 2.34 cm³ of its 14.2 cm³ comes back as bone, because the five-millimetre closing reaches up into the cap from the skull beneath it. Small, and on the far side of the vault from both entries, but it means the density threshold is not the only thing standing between the cap and bone.
The focus. The point BabelBrain marks in its domain sits 0.27 mm and 0.33 mm from the targets we asked for, which is the size of one of its voxels and nothing else. Both land in brain.