Bumpy flatmap redesign - #720
Open
alexhuth wants to merge 18 commits into
Open
Conversation
The bumpy flatmap gives a flatmap relief that survives full data coverage. It was computed in javascript, at viewer load time, and it was noisy enough that the result was treated with iterated Laplacian smoothing. Reading the shipped code, two things stood out. The height was `computeFlatVolumeHeight(wmareas, vertexvolumes)` -- the denominator was the *folded* white matter vertex area, not the flattened one -- so with `A_pial = r**2 * A_wm` it collapsed to `thickness * (1 + r + r**2) / 3` and never looked at the flatmap at all. The version that does use flat areas was commented out, because it spikes: `V / A_flat` has a denominator that goes to zero wherever flattening crushed a triangle, and smoothing a field of ratios afterwards is dominated by exactly those outliers. Replace it with a physical model, computed in python and cached. cortex.polyutils.FlatSlab treats the tissue between the white matter and pial surfaces as a soft elastic solid, pins the white matter side to the flatmap and lets the pial side settle, minimising a stable Neo-Hookean energy (Smith, De Goes & Kim 2018) over wedge elements cut into tetrahedra. What keeps real tissue from spiking is not gravity -- over a 3mm slab `rho*g*t/mu` is about 0.026, a couple of percent -- but shear: a tall narrow column is expensive because it shears against its neighbours and relieves that by spreading sideways. So the pial vertices get three degrees of freedom rather than one, gyral crowns end up wider at the top than at the bottom and sulcal fundi narrower, and the relief comes out smooth without anything being smoothed. With no body force and a uniform material the shear modulus factors out of the minimiser entirely, so Poisson's ratio is the only material parameter that affects the result and it is the only one the API takes. On S1's left hemisphere this drops the height's p99.9-to-median ratio from 2.4 (naive) and 3.8 (the old javascript) to 1.9, with the tallest point at 5.9mm rather than 13.9mm, and the height correlates with mean curvature at +0.52. Also move the equivolume depth sampling's vertex areas out of javascript, which lets them use the cotangent-weighted operator in Surface.smooth rather than a uniform umbrella, and removes the last of the per-load numeric kernels. Plumbing: - surfinfo.bumpy_flatmap and surfinfo.equivolume_areas cache the results. Both store more than one value per vertex, so they avoid the `left`/`right` keys that get_surfinfo turns into a Vertex by concatenation. - get_surfinfo takes `generate=False` for info too expensive to compute behind someone's back. The bumpy flatmap uses it: it is generated when a flatmap is imported, and a subject imported before this existed gets a flat flatmap and a message saying how to fill it in, rather than a viewer that disappears for a quarter of an hour. - import_flat generates it, and now also clears the surface info derived from the flatmap it just replaced -- clear_cache only empties cache/, so distortion and flat_border could go stale there too. - The offsets ride to the viewer in a `flatoffset` ctm attribute. The surface shaders use all 16 of the vertex attribute slots webgl guarantees (see #715), so javascript unpacks them into `flatbump.w` and the unused `w` of `wm` and of the flat morph target rather than adding a seventeenth. - bumpy_flatmap_scale exaggerates the relief for display without touching the cached geometry. - The ctm cache token goes to v4, since older packs have neither attribute. Tests check the relaxation against cases with known answers rather than against stored output: an identity flat map has to leave the slab exactly alone, a uniformly stretched thin slab has to reproduce the homogeneous solution in its interior, and the requested Poisson ratio has to come back out of the small strain elasticity tensor. A headless test renders the flatmap with and without the relief and fails if the two images match. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
Three measured changes, on S1's left hemisphere (874k tetrahedra): Work in the transposed convention throughout. A column slice of a C-ordered (n, 3, 3) array is strided while a row slice is contiguous, and the objective spent much of its time in cross products and matrix multiplies over strided views. Transposing the deformation gradient, the cofactor and the per-node gradient -- and folding the element volume into a precomputed `vol0 * Dm^-1` rather than scaling afterwards -- takes an energy-and-gradient evaluation from 1.31s to 0.72s in isolation. The gradients agree with the old code to 1.2e-11 relative. Scatter the gradient in one pass. Three bincounts over strided component slices each had to copy their weights out first; one bincount over a precomputed flattened index takes the per-element array's own buffer. Raise the L-BFGS memory from the default 10 correction pairs to 60. At a fixed 150 iterations this reaches a 29% lower energy (1.16e5 against 1.64e5) for about 10% more time per iteration. Relax the two hemispheres in parallel, which is worth much less than the factor of two it looks like -- the objective is memory-bandwidth-bound, so two copies compete -- but is still about 15% off the wall clock. Together: both hemispheres of S1 at 400 iterations now take 15.6 minutes rather than roughly 25. Two things that did not work and are recorded here so they are not retried: single precision is unusable, because a cortical mesh's element volumes span enough orders of magnitude that `vol0 * Dm^-1` loses its dynamic range in float32 and the energy comes out 12x wrong; and reusing preallocated buffers is slower than allocating, because writing into `(n, 3, 3)` slices is strided. Cutting the iteration count is not a safe optimisation: at 150 iterations the heights still differ from the 400-iteration answer by 8% at the median. Also store the cached offsets compressed and in single precision -- they end up in a float32 vertex attribute regardless -- which takes the file for a subject from about 7MB to 3.3MB, and rewrite the module docstring as technical prose with references. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
Relaxing the cortical slab takes about a quarter of an hour for both hemispheres, which is a long time to wait to see what the feature does on the subject that ships with pycortex. Committing the result means S1's flatmap has relief straight out of the box, and the gallery example and the headless viewer test do not each have to compute it. Generated with the shipped defaults (poisson_ratio=0.45, 400 iterations), compressed and stored as float32, which is 3.3MB. The rest of surface-info/ stays ignored, as before; this one file is checked in deliberately. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
Now that S1's relaxation is in the filestore the example can look it up instead of recomputing a shortened version of it, which takes the page from 259 seconds to build to 5, and shows the converged result for both hemispheres rather than one truncated hemisphere. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
A quasi-Newton method takes far longer to resolve the long wavelength part of this problem -- how the cortical sheet as a whole slides as it settles -- than the local detail, and on a mesh of 150,000 vertices that is what dominates the cost. Solving a smaller mesh first finds exactly that part cheaply. Getting a coarser mesh is easy here only because a flatmap is planar. Take a maximal independent set of the vertices, so no two survivors were neighbours, and retriangulate them with a two dimensional Delaunay triangulation; no surface-aware decimation is needed, and none is available anyway since polyutils.decimate needs tvtk. Delaunay triangulates the convex hull, though, so it bridges the medial wall and every relaxation cut. Those bridges join vertices that are close together in the plane but far apart across the surface, so a coarse triangle is kept only if its vertices are within a few steps of each other in the fine mesh's own graph. How many steps that should be depends on how regular the mesh is, so it is not fixed but measured: the number chosen is the one whose coarse mesh covers the same total area as the fine one. On S1 that lands on five steps and matches the fine area to 0.01%, against 17% too much with no rejection at all. Each level is about three and a half times smaller than the one above, so on S1 the hierarchy is 146k, 42k and 11k vertices. A coarse solution is prolonged onto the next level by barycentric interpolation, with the few vertices that fall outside the coarse triangulation taking their nearest coarse vertex's value. Every level gets the same iteration count, which means the hierarchy adds about 40% to the cost of one pass over the full mesh. Giving the coarse levels more than that was measurably a waste: they were converged well past the point where the full mesh could still tell the difference, and spent 65% of the total time doing it. The tests check the properties this rests on rather than a stored result: that the surviving vertices really are an independent set, that a hole punched in a mesh is still a hole after coarsening (which fails if the bridge rejection is removed), that prolongation reproduces affine fields exactly and is a partition of unity, and that a coarse-to-fine solve reaches the same energy and the same geometry as a single-level solve on a patch small enough for both to converge. Also from a docs audit: add the new public API to the sphinx reference, make polyutils re-export everything bumpy.py declares public, say precisely why equivolume_areas cannot be a Vertex (it holds four scalar maps, not a multi-component one), and delete a commented-out javascript fallback that referred to attributes which no longer exist. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
With the coarse-to-fine solve in place the full mesh no longer needs anything
like as many iterations, because it no longer starts from scratch. Measured on
S1's left hemisphere, energy against wall clock:
single level, 150 iterations 250s 1.16e5
single level, 250 iterations 344s 5.49e4
single level, 400 iterations 555s 2.10e4 <- the old default
coarse to fine, 60 iterations 128s 7.20e4
coarse to fine, 100 iterations 209s 3.86e4
coarse to fine, 150 iterations 301s 1.61e4
So 150 iterations coarse to fine reaches a *lower* energy than 400 did on its
own, in a little over half the time. Comparing at matched energy rather than
matched iterations, the hierarchy is worth between two and two and a half times.
Default max_iter accordingly drops from 400 to 150.
Relaxing the two hemispheres in separate processes now buys nothing and stops
being the default. The solve is memory-bandwidth-bound rather than
compute-bound, so two at once mostly compete for the same bandwidth: on S1 that
is 611 seconds in parallel against 602 sequential, for twice the peak memory.
The parameter stays for machines with bandwidth to spare.
Regenerate the shipped S1 offsets with the new defaults. Both hemispheres now
take 611 seconds rather than 938, and land at a lower energy: 1.61e4 against
2.10e4 on the left. Volume loss rises from about 9% to about 12%, which is the
solution getting closer to the true minimum rather than a regression -- at a
Poisson ratio of 0.45 the material is compressible and the minimum genuinely
sheds that much; raising the ratio towards 0.5 recovers it, which is what
test_incompressible_slab_conserves_volume checks.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
bumpy_flatmap_scale was only reachable by editing the config file, which is an awkward way to adjust something whose whole purpose is to be looked at. Put it in the surface controls next to the bumpy_flatmap toggle, bound straight to the shader uniform the way depth and opacity are, so dragging it redraws immediately. The slider starts wherever the configuration file put it. Its range runs from 0 to five times true scale, or to twice the configured value when that is already higher, so someone who deliberately set 8 does not find it clamped to 5 the first time the menu is drawn. A missing or unparseable setting falls back to 1; that cannot be written as `parseFloat(...) || 1`, because 0 is a legitimate value and has to survive. Adding it to the surface menu is enough to make it settable from python and saved with a view: JSMixer.view_props is built from the menu's own controls. The headless test now also drives the slider and checks the render changes, which is what catches the uniform not reaching the shader. It deliberately does not read the value back afterwards: reading any surface menu property through the javascript proxy returns an empty dict, for unfold, depth and opacity just as much as for this one, so a read-back assertion would be testing that pre-existing limitation rather than this feature. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
Two independent defects, both mine, and neither to do with Poisson's ratio -- which is why sweeping it made no difference. A slab one element thick. The whole design argument was that shear across the thickness of the slab supplies the regularisation, and then the slab was built with exactly one element layer through it. A single linear element cannot represent a shear profile across a layer, so the coupling being relied on did not exist: neighbouring columns of tissue exchanged material far too freely and the relief flattened towards a slab of uniform thickness. On a patch of S1 the correlation between relief and mean curvature is 0.77 for the unrelaxed prism estimate, 0.52 after the one-layer relaxation, and 0.79 with four layers. No regularisation at all. The javascript this replaced smoothed vertex areas over five iterations and thickness over twenty. This smoothed nothing, so millimetre-scale noise in the segmentation went straight into the relief. Smoothing the white-to-pial displacement -- the folding is left alone, only the local thickness is regularised -- takes mesh-scale roughness from 0.21 to 0.10 relative to the signal. Neither substitutes for the other. Smoothing alone leaves the correlation at 0.57; layers alone leave the roughness at 0.21. Four layers is four times the elements, so `resolution` now stops the hierarchy before the finest mesh and interpolates the rest of the way. The relief varies over millimetres and solving it at the scale of a single triangle mostly buys the opportunity to wrinkle; solving S1's 1.6mm level rather than its 0.9mm one costs a fifth as much and band-limits the result by construction. Solving every level at four layers takes upwards of forty minutes a hemisphere and is not worth it. Both hemispheres of S1 now take 729 seconds against 611, and volume loss drops from about 12% to 4-6% because the relaxation can no longer squeeze material as freely. The test that should have caught this did not exist: the relaxation was only ever checked for energy, volume and convergence, all of which improved while the picture got worse. test_relief_follows_the_folding_and_is_not_mesh_scale_noise now asserts both properties on a real patch and pins the cause by requiring a one-layer unsmoothed slab to do worse on both. test_coarse_to_fine_matches_a_single_level_solve needed resolution=0 to keep meaning what it says, since at the default the hierarchy deliberately stops before the finest mesh. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JG9XE4ZAitPgsQYRDwtwVT
…scale Three problems showed up once the relief was actually usable. It was far too slow. Both hemispheres of S1 took 729 seconds at import. Solving at 3.1 mm rather than 1.6 mm, three element layers rather than four, and 60 iterations rather than 150 bring that to 57. This is not free -- the correlation between the relief and high-passed mean curvature falls from 0.70 to 0.57, so real gyral-scale detail is being given up -- but a subject import that pauses for a minute is a different thing from one that pauses for twelve. Two of those three are nearly free on their own: a fourth element layer costs 60% more time and moves the correlation by 0.002 in the wrong direction, and 100 iterations against 60 costs a third more time for 0.001. Almost all of the loss is the resolution, and almost all of the speed is too. It was not smooth. The heights looked fine; the lighting did not. A shading normal is the derivative of the height field, and every level finer than `resolution` is reached by barycentric interpolation, which is only C0 -- so the field had creases along the coarse triangle edges that were invisible in the heights and glaring in the shading. Smoothing the finished offsets (`polish`) takes the RMS angle between the normals of neighbouring triangles from 8.8 degrees to 2.4, and *improves* the correlation with curvature, so at this strength it removes noise rather than signal. The shading did not track the scale slider. The geometry grew and the lighting stayed where it was. Exaggeration is now vertical only, like on a topographic map, so turning the slider up does not smear the relief sideways relative to the data underneath it -- and because it is vertical only, the corrected normal is exact rather than approximate: scaling a height field by s scales the normal's in-plane components by s and leaves the out-of-plane one alone. It needs no recomputation when the slider moves and degenerates correctly to the flat normal at zero. Two further defects found while tracing the normals. The displacement ramped over anatomical-to-inflated while the shading normal ramped over inflated-to-flat, so the two disagreed across the whole first half of the unfold -- harmless when the offset was a height along the current normal, wrong once I made it a fixed offset in flatmap coordinates, which mean nothing on a folded surface. Both now ramp over the last segment. And `computeNormal` gave vertices in no triangle a zero normal, which the fragment shader then divided by; they fall back to the sheet's own normal. Two speedups that cost nothing: `Surface.smooth` factorises its operator on every call, so smoothing a three-component field paid for three factorisations of the same matrix; and the coarsening's hop search multiplied full n-by-n reachability matrices when only the surviving third of the rows was ever read. Not fixed here, but found: the pick and depth shaders have no HASFLAT branch, so they render undisplaced geometry while the colour pass renders displaced geometry. Picking drifts from what is drawn, and the drift grows with the scale slider.
The bumpy flatmap read as fine texture rather than hills and valleys, for the third round running. The reason it kept coming back is that every metric it had ever been tuned against was blind to the thing being complained about: `silk` compares neighbouring face normals, `roughness` compares a vertex to its one-ring, and the correlation guard high-passed the field before measuring it. All three are local. Nothing in the repo could see a relief that was beautifully smooth and anatomically empty, which is exactly what those three parameters had been driven to produce. Measuring it by spatial band shows the mistake. `Surface.smooth`'s `factor` is a diffusion time, and one backward-Euler step has transfer function 1/(1 + k^2 t), so its half-power point sits at a wavelength of 2*pi*sqrt(t) -- a factor of 4.4 longer than the Gaussian-equivalent sigma that made these values look mild when I picked them. smooth = 1.0 cuts at 6.3 mm polish = 4.0 cuts at 12.6 mm correlation_length = 2.5 mm cuts at 15.7 mm Gyri are 8-16 mm across. All three were filtering the signal, not the noise. `smooth` was the worst of them, and it is now off by default. It was added to keep segmentation noise out of the relief and it does, but at 6.3 mm it was removing the thickness variation that makes gyri thicker than sulci: on S1 it roughly halved the correlation with mean curvature, 0.47 to 0.25 in the 8-16 mm band and 0.62 to 0.42 in 16-32 mm. Even 0.05 measured slightly worse than zero. `polish` removes the same noise afterwards for a fraction of the signal, so there is no regime where it helps. `polish` drops 4.0 -> 2.0 and `correlation_length` 2.5 -> 0.5 mm. The second of those should not matter at all -- it sets only the starting point -- and it matters because the solve stops at `max_iter` rather than at a tolerance, which nothing printed and I had never checked. It does not converge: `info['converged']` is False with a final gradient sup-norm around 3e3 against scipy's 1e-5. An unconverged solve stays near its starting point, so over-smoothing that point costs real amplitude. `_relax_hemisphere` now says so out loud. On S1, against what shipped before, and slightly faster because `smooth` was itself a sparse solve: runtime 56.7 s -> 49.1 s hp-corr 0.567 -> 0.739 height std 0.536 -> 0.703 mm silk 2.41 -> 2.80 deg, but 4.52 -> 4.00 per unit of relief That beats the old resolution=1.5 configuration (0.700) while running four times faster than it, so the resolution trade made last round was never the real cost. Two things this rules out. Poisson's ratio is irrelevant -- 0.499 against 0.45 changes amplitude by 1.02x and makes the correlation worse -- which settles the earlier experiments with 0.45/0.49/0.499. And the flattening area change correlates 0.00-0.05 with both mean and Gaussian curvature in every band, so it carries no folding signal at all; the relief is essentially cortical thickness, and a theory I had that flattening compression is Gaussian-curvature-driven and therefore intrinsically speckled is simply wrong. Adds a band-resolved test, since the absence of one is what allowed this. Still open, and now visible rather than hidden: the solve does not converge, and at the default `resolution` the floor clamp leaves only the coarsest level actually solved, so the coarse-to-fine cascade is not running. Adding levels coarser than the floor would cost very little.
The relief still reads as texture rather than terrain, and no amount of solver work will change that, because shaded relief is itself a high-pass filter: brightness follows the surface gradient, so a band at wavelength L and amplitude A contributes contrast in proportion to A/L. A 30 mm swell needs roughly ten times the amplitude of a 3 mm wrinkle to read as strongly. Terrain looks like hills and valleys because its amplitude spectrum falls steeply with frequency; cortical relief's does not. `polish` cannot fix it -- a low-pass shrinks the whole relief and has to be chased with the viewer's scale slider, which rescales every band equally and so changes nothing about the balance. `tilt` boosts the long wavelengths and leaves the short ones alone. On S1, the share of shading above 32 mm against the share in 4-16 mm: tilt=1 25% vs 40% silk 2.80 relief std 0.70 mm tilt=2 33% vs 33% silk 3.00 relief std 1.22 mm tilt=3 38% vs 28% silk 3.24 relief std 1.75 mm tilt=5 43% vs 24% silk 3.75 relief std 2.81 mm Note how little the normal continuity moves for how much relief is added -- tilt=3 nearly triples it for half a degree -- because a long wavelength carries very little gradient per unit amplitude. That is the same fact that makes the relief look high-pass, used in reverse. `tilt_scale` defaults to 24 mm, chosen to sit over the bands that actually carry folding: the relief's correlation with mean curvature is 0.62 at 16-32 mm and 0.66 at 32-64 mm, but only 0.30 above 64 mm, so exaggerating much past 40 mm looks smoother while spending the exaggeration on the one band that is barely anatomical. This is a display choice and it is documented as one. `tilt=1` is off and is the default, so nothing changes unless it is asked for. Also plumbs FlatSlab's keyword arguments through `surfinfo.bumpy_flatmap`, which previously accepted only `poisson_ratio` -- the parameter that turns out to matter least. `db.get_surfinfo` folds them into the cache filename, but `brainctm` looks for the unsuffixed file, so the docstring gives the recipe for trying a setting in the viewer.
`tilt` boosted the relief's long wavelengths, on the reasoning that shaded relief is a high-pass filter so the long wavelengths need help. Measured on the real field, it makes things worse, and the report that it "just increases the overall height" is exactly right. The relief has a whole-map swell in it. Per-band standard deviation on S1, with each band's own mean removed: <2 2-4 4-8 8-16 16-32 32-64 >64 total 0.008 0.019 0.054 0.112 0.157 0.147 0.421 0.700 That last column is 36% of the variance sitting in one mode broad enough to span the flatmap. A long-wavelength shelf boosts it hardest -- of all the variance `tilt=3` added, 53% went there -- so the map gets taller and domed while the gyri stay where they were. Normalising the gyral gain by the height increase: shelf tilt=3 1.94 / 2.49 = 0.78 <- worse than doing nothing bandpass 8-48 g=3 2.32 / 1.89 = 1.23 dome removed 0.88 / 0.66 = 1.33 dome + bandpass 2.20 / 1.64 = 1.34 The shelf also lowered the correlation with mean curvature in every band (0.62 -> 0.58 at 16-32 mm). So the right move is the opposite one: stop spending amplitude on structure larger than folding. `detrend` removes it, and the band-pass then adds nothing that the viewer's scale slider does not already provide, so it is not worth a parameter. The swell is real -- cortex is regionally thicker in some lobes than others -- but it is not folding: it correlates with mean curvature at 0.30 where the gyral bands manage 0.62-0.66, and `bumpy_flatmap_scale` multiplies it along with everything else, so it eats the visual dynamic range that should be showing the folding. Removing it takes the relief from 0.70 mm to 0.46 mm and leaves the gyral bands within 10-15% of where they were, so the same slider setting times 1.5 gets back to the same height with the folding half again as strong. Correlation improves in every band, most above 64 mm (0.30 -> 0.56). On S1: 55.3 s, silk 2.76 (was 2.80), hp-corr 0.735 (was 0.739), height std 0.460 (was 0.703). The two quality numbers are unchanged because both are blind to this band by construction -- hp-corr high-passes at 31 mm -- which is the same blind spot that let the earlier defaults drift. Found while checking the previous commit's own claim, and it took two attempts: the first measurement subtracted the mean over interior vertices only, leaving a step at the flatmap boundary that put 70% of the apparent shading below 2 mm. Bands are now taken from the untouched field and reported as per-band standard deviations.
…ot fix
The relief reads as round knobs rather than gyral ridges. Measured cause:
with resolution=3.2 only the coarsest ~3.1 mm level is solved, so the
finished relief lives entirely in the span of that mesh's hat functions,
and on a near-equilateral Delaunay mesh each hat is a circular disc about
6.2 mm across -- the width of a gyral ridge. A 6 x 20 mm ridge holds
about 14 coarse vertices: 6-7 along its length and 2 across it. Across
the ridge that is one linear tent and no shape; along it, a chain of
independent round hats.
Two changes here.
`detail` puts ridge-scale structure back. Shear only couples the slab
over roughly its own thickness, and a ridge is several times wider than
that -- the regime where the cheap volume-preserving prism height is
already close to right. So the relaxation supplies the long wavelengths,
the full-resolution prism height supplies the crossover down to
`detail_floor`, and nothing comes from below that.
`anisotropy` makes the polish follow the gyri. Isotropic diffusion
applied to a chain of knobs rounds each one off; it cannot join them up.
The direction comes from the structure tensor of mean curvature, damped
in proportion to coherence so an isotropic patch is not combed into
streaks from noise. The stiffness is the standard FEM generalisation of
the cotangent formula, reusing the rotated edge vectors already cached on
Surface; with an identity tensor it reproduces the existing operator to
1.8e-15, which is asserted in the tests because a subtly wrong stiffness
would still produce a plausible-looking smooth field.
Honest accounting of what this buys. Ridgeness -- the inequality of the
principal curvatures at crests -- goes 0.677 to 0.699, and correlation
with thickness 0.668 to 0.714. That is real and it is not enough, because
the ceiling is low for the quantity being computed:
mean curvature 0.833
pia/wm area ratio r 0.776
legacy javascript height 0.736
thickness 0.705
this relief 0.695
naive prism 0.688
The legacy javascript version really is more ridge-like than the
relaxation, and the decomposition says why. It computes
thickness * (1 + r + r**2) / 3; thickness alone scores only 0.705, and
what carries the ridges is r, the pial-to-white area ratio, which tracks
folding. The relaxation instead computes a thickness-and-flattening
quantity, and flattening compression correlates with curvature at 0.00 to
0.05 -- essentially noise -- which leaves it as smoothed thickness. Every
field available to rebalance sits at 0.69-0.71, so no reweighting of them
reaches the ridges; pushing the crossover to 48 mm hands almost
everything to the prism and asymptotes at 0.710.
Getting ridges means putting the folding term into the elastic reference,
which is a change to what is being relaxed rather than to how it is
filtered, and is not attempted here.
Also adds the first test in this file with anisotropic geometry. Every
existing synthetic case is a square grid stretched equally in both
directions, so not one of them could fail on a relief made of round
knobs.
The relief gave almost no sense of folding, and the reason was the quantity rather than the machinery. Both candidate heights are the same folded volume over an area; they differ only in which area. The old one divided by the *flattened* triangle area, the honest model of a slab laid flat. But the flatmap's own area distortion measures essentially uncorrelated with both mean and Gaussian curvature -- as a denominator it contributes no anatomy and injects the flattening algorithm's artifacts in its place. What is left is close to a map of cortical thickness, and thickness is a blobby field that reads as round knobs. Dividing by the *folded* white matter area instead gives thickness * (1 + r + r^2) / 3 with r = sqrt(A_pia/A_wm), and r tracks folding directly: the pia carries more area than the white matter beneath a gyral crown and less in a fundus. Measured on S1 with every field smoothed alike, by anisotropy of the Hessian at crests: r scores 0.776 where thickness alone scores 0.705 and mean curvature -- the folding itself -- scores 0.833. This is algebraically what the old javascript implementation computed, which is why that one looked more folded despite being cruder in every other respect. It scores 0.734; the relief now scores 0.759, and its correlation with smoothed curvature goes from 0.739 to 0.894/0.854. It has to be computed on the full-resolution mesh. Every coarse level subsamples the white and pial surfaces and retriangulates them with the *flat* Delaunay connectivity, so each coarse edge is a chord through the fold rather than an arc along it. The pial flare is second order -- curvature times thickness -- while the surfaces are zeroth order, so it is the first thing that deficit removes: the coarse area ratio collapses towards one and the coarse solve minimises an energy whose reference has already had the gyri taken out of it. Both the detail splice and the coarse warm start now take the folding term from the finest mesh. Fixing the coarse reference properly is recorded in the module docstring. Two consequences worth being plain about. The relief no longer shows flattening compression at all -- a patch that is flat but whose flatmap is distorted now produces no relief, by design. And the relaxation contributes much less to what is on screen: with max_iter=60 the solve stays near its starting point, and the starting point is now the folding height. The one-layer control in the folding test has been removed rather than adjusted, because it no longer isolates anything. Costs 68.5s for both hemispheres against a 60s target, and the relief is about twice as tall, so the raw inter-face angle rises from 2.8 to 3.8 degrees -- but per unit relief the surface is smoother than before. New tests: folding at constant thickness must still produce relief; folding height must be invariant to flatmap distortion where the prism height is not; and the frustum algebra is pinned.
The relief is right, so this removes the machinery that was used to find it rather than to express it. bumpy.py goes from 1452 lines to 415 and test_bumpy.py from 866 to 354; the feature as a whole from +3104 lines to +1413, and from 68.5s to 11.3s for both hemispheres. What goes is the elastic solver and its multigrid: the Neo-Hookean energy, the tetrahedral element assembly, L-BFGS, the maximal-independent -set coarsening, the barycentric prolongation, the hierarchy, and the eight parameters that only fed them. Also the anisotropic ridge filter -- 130 lines for 0.017 of ridgeness. Deleting the solver was decided by measurement against a threshold fixed beforehand, and the honest account is that the threshold said keep. The solver moves the answer by 31% RMS, far over the 5% agreed. But correlation between the two is 0.991 and ridgeness 0.741 against 0.738: the difference is almost entirely a gain factor of about 0.85, and the in-plane slip only the solver can produce is 0.047mm rms against a relief of 1.0mm. So it was acting as a constant the viewer's scale setting absorbs, at 34 of the 38 seconds. The gain is given back by setting bumpy_flatmap_scale to 1.55 rather than the 1.8 that was asked for, so what reaches the screen is what was approved. On S1 the relief comes out at hp-corr 0.886/0.880 against 0.894/0.854 before. Silk rises from 3.84 to 4.66 degrees, part from the taller relief and part from losing the anisotropic filter; scaled by the new default that is about 4% rougher on screen than before. With no relaxation the relief is purely vertical, so the shader carries it in one component instead of three: flatbump.w alone, giving back the spare w lanes of `wm` and the flat morph target. And at 5s a hemisphere the plumbing built for a slow computation goes too -- get_surfinfo's `generate` flag, brainctm's not-generated-on-demand contract, the import_flat hook, and surfinfo's parallel path. Kept: legacy_js_height and naive_prism_height, the two reference heights the folding quantity was judged against, and which the gallery example plots. The coarse-reference chord deficit that motivated computing the folding term at full resolution is recorded in the module docstring; it is no longer reachable, but it is why the answer is not simply relaxed.
Docstrings that recounted how the answer was found rather than what the
code does. The scorecards ("0.776 where thickness scores 0.705"), the
account of the elastic solver's removal, and the notes on which metric
misled which round of tuning are all history, and history belongs in the
log. bumpy.py 415 -> 388, test_bumpy.py 354 -> 345, and the feature from
+1413 to +1345.
Kept the facts that stop someone making a mistake: that a diffusion time
`t` is half power at a wavelength of 2*pi*sqrt(t) and not sqrt(t), which
is a factor of 4.4 and easy to get backwards; why the denominator is the
folded area; and why the log-space regularisation is the right averaging
for a ratio.
Two real errors fixed along the way. Every top-level def in bumpy.py had
lost its blank-line separation -- the previous commit normalised
whitespace with a regex that collapsed them. And the shader still carried
a paragraph explaining that the pial surface slides sideways so the
offset is a full vector rather than a height, immediately above the line
that now writes two zeros; the relief has been purely vertical since the
relaxation went.
…tils You were right about the duplication -- there was more of it than the helpers I had noticed. bumpy.py goes 388 -> 251, and the feature from +1345 to +1201. Four duplicates, three of them of things already in polyutils: `_face_areas(pts, polys)` was `misc.face_area(pts[polys])` written out again. `_vertex_areas` was `_lumped` of that. `face_prism_volumes` was `misc.face_volume` -- the same three-tetrahedron decomposition, except that misc's version loops in python and prints its progress every thousand faces, which is presumably why it got rewritten rather than used. Vectorised it in place instead, checked against `brick_vol`, and deleted the copy. `_smooth_vectors` was the fourth and the least obvious: it was `Surface.smooth`'s operator assembly copied out so that one factorisation could serve three columns. The relief has been a single column since the relaxation went, so the reason for it had already gone -- `Surface.smooth` is now called directly. `_regularise_log_height` was assembling that same operator a third time; since M - t(W-V) is M + tL, the screened-Poisson solve it wants is exactly what Surface.smooth does, so it is now a log, a smooth and an exp. Also removes `legacy_js_height` and `_umbrella_smooth`. It served its purpose -- it is what showed the relief should be driven by the pia/white area ratio -- but the answer is in `folding_height` now and the example does not need three curves to make the point. `misc.face_volume` is the one behaviour change outside the feature: same signature and same numbers, no longer O(n) in python and no longer printing. It was exported but unused.
…e slip Checking for stale references turned up a real bug rather than only prose. The example's plotting loops used `height` as the loop variable while `height` was also the folding relief being plotted. So after the first loop the name was bound to the naive field, and the histogram's "folding height" curve and the final curvature hexbin both drew that instead. The shadowing came in when `relaxed` was renamed to `height` two commits ago; the example still ran and still produced plausible figures, which is why nothing caught it. The relief is now `relief` and the loop variable `field`. Removed with it: a whole section plotting "how far the pial surface slides sideways", which is identically zero now -- it was described as "the part a vertical-prism model cannot represent at all". A sentence about the javascript height that trailed off mid-clause when that curve was deleted. And the headless test's claim that javascript packs the offset into `flatbump`, `wm` and the flat morph target, when it has been `flatbump.w` alone since the relief went vertical. Also "relaxation cuts" -> "flattening cuts" in the example, which was the freesurfer term for the cuts made to flatten a surface but now reads as a reference to a relaxation that no longer exists.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This change revamps how the bumpy flatmap and equivolume transform are computed. Previously, they were both computed on the client side in javascript. This was the source of a lot of complexity on the javascript side and was inconsistent with how most other things in pycortex work. Now, both the bumpy flatmap offsets and equivolume are computed in python and then passed to the web viewer, similar to curvature, etc. Moving these computations to python let us do things that are a little fancier and end up looking better. In particular, the bumpy flatmap computation now uses better smoothing than before, making it look prettier. (see screenshot)
The python version of these algorithms is a bit more succinct than the javascript version, but this change ends up adding a decent number of lines of code due to additional code to transport the data from python to javascript, test code, examples, and documentation.