neuroim2 0.19.1

Upgrading from CRAN 0.13.0

This release collects the changes recorded below for versions 0.14 through 0.19.1. It retains the R >= 4.3 requirement. Analyses and plotting code should account for these changes:

Documentation checks

The image(), as.raster() and ClusteredNeuroVol-to-DenseNeuroVol coercion help topics now include usage matching their documented arguments.

Decoded mapped reads and sequence orientation

MappedNeuroVec now applies source intensity slope/intercept exactly once, including per-volume scaling. Previously mapped access could return stored integers while dense and file-backed access returned decoded intensities. Re-read mapped inputs and rerun analyses that used scaled images.

series(NeuroVecSeq, ...) now preserves time-by-voxel orientation when the number of selected voxels equals the number of time points in a component. Previously those square blocks were silently transposed, also affecting voxel-group reductions. These two correctness fixes are shared with the separate 0.20.0 development work; lazy block iteration is not included here.

Sparse downsample(), scale_series(), arithmetic, concat() and coercion from ROIVec now preserve time-by-voxel orientation when the number of output voxels equals the number of time points. The former shape heuristic silently transposed those results. Internal sparse reads and dense-to-sparse conversions also declare their known matrix layout explicitly, avoiding spurious ambiguity warnings.

as.sparse() now preserves matrix dimensions when converting a DenseNeuroVec with one time point or a mask selecting one voxel. Numeric voxel indices and LogicalNeuroVol masks retain the same voxel-by-time matrix and time-by-voxel series contracts. Previously these singleton cases errored after dropping the selected data to a vector.

Numeric masks also retain each value at its original voxel when indices are unsorted or repeated. Previously unsorted indices could silently assign values to different voxels, and repeated indices could fail a cardinality check. Numeric and equivalent logical masks now agree, including negative exclusions, zeros and empty selections. Missing, non-finite and positive out-of-bounds indices are rejected. Recompute sparse conversions made with unsorted masks.

Reproducible soft-alpha overlays (#21)

soft_alpha_params() is now exported. It validates its inputs, accepts an explicit knee, and returns an alpha_floor, so the soft opacity curve can be previewed and reused. plot_overlay() gains alpha_knee, alpha_cap, alpha_floor, alpha_mid, gamma_min and gamma_max (appended to the signature), decoupling opacity from the colour limits. The resolved curve is recorded in attr(result, "soft_alpha"). Defaults are unchanged.

Bilateral filter bandwidths across maps (#20)

The bilateral_filter() and bilateral_filter_4d() help now explains that range_scale = NULL estimates the intensity bandwidth per input, so equal intensity_sigma values are not comparable across subjects or contrasts. For batch or group work, reuse one fixed range_scale.

Skull-stripped backgrounds (#19)

The background window of the redesigned plot functions is computed over head voxels, so bright tissue in skull-stripped images is no longer clipped (0% of brain voxels clipped on a skull-stripped MNI152 T1, versus 5.4% with the previous whole-image robust quantiles). bg_range = "data" remains available for a full-range window.

SparseNeuroVec matrix orientation (#31)

SparseNeuroVec() documents the matrix layout convention and accepts an explicit orientation argument ("auto", "voxels_x_time", "time_x_voxels"). When n_voxels == n_timepoints, auto mode keeps the historic voxels-by-time assumption but warns; pass orientation = "time_x_voxels" for a square series() result so the data are not silently transposed.

Redesigned slice figures (behavior change)

plot_overlay(), plot_montage(), plot_ortho(), plot_edge_overlay() and plot_checkerboard() now share one figure engine and look like one family: borderless black tiles cropped to the head, world-coordinate slice labels (z = -12 mm), haloed orientation letters, and a compact fixed-width colorbar placed beside the tiles.

Configurable colorbar titles

plot_overlay(), plot_montage(), and plot_ortho() gain a cbar_title argument naming the quantity drawn above the report-style colorbar. It defaults to "value", which is what the colorbar has always shown, so existing figures are unchanged. The new argument is appended after the existing parameters so positional callers of legend / crop / interpolate keep working. Set it to the statistic actually being displayed – cbar_title = "Semipartial r", "Delay coefficient (% signal change)" – so a figure does not assert a quantity it is not showing.

Sign-neutral overlay legend (behavior change)

The style = "report" legend strip previously labeled its swatches “Positive activation” / “Negative activation”, or “Activation” for unsigned maps. An overlay is not necessarily a BOLD activation map – it may be a correlation difference, a semipartial r, or a regression coefficient – so the labels are now “Positive” / “Negative” and “Suprathreshold”. The accompanying “higher than threshold” / “lower than threshold” subtitles are unchanged. Figures using style = "report" will render with the new wording.

neuroim2 0.19.0

Sparse NIfTI writing

write_vec() now writes SparseNeuroVec images directly through the public API, preserving voxel/time ordering, geometry, volume labels, and zeros outside the mask in .nii and .nii.gz files (#36). Writing still expands the sparse values into a full voxel buffer, but no explicit DenseNeuroVec conversion is required.

Index-only searchlight geometry

searchlight_indices(mask, radius, nonzero = TRUE) now compiles spherical searchlights directly to 1-based full-volume linear indices. It uses the same cached, anisotropic spherical-offset and compiled boundary-clipping machinery as the existing searchlight APIs, but constructs no ROI objects or coordinate matrices and never changes parallel execution state. The result records centre indices and spatial metadata so downstream sparse-measurement builders can use the geometry without taking on neuroim2’s analysis or execution policy.

Slice plots now preserve anatomical orientation

plot() for NeuroVol and NeuroSlice, plot_montage(), plot_ortho(), and the overlay/registration-QC helpers now share one affine-derived display transform. Native slice axes are permuted and flipped so right, anterior, and superior point toward increasing screen coordinates. This fixes axial plots whose anterior direction was rotated toward screen-right, corrects sagittal orientation labels from L/R to P/A, handles permuted voxel axes, and keeps oblique native slices on a regular raster without silently shifting pixels or resampling values.

Seven defects found by measuring the image-processing layer against nilearn 0.14.0 / scipy, which is the reference for the processing neuroim2 also does. The harness is in dev/bench/nilearn/ and the findings in dev/nilearn-gap-analysis.md.

Four of these change results; the sections that do say so explicitly.

gaussian_blur() now applies the smoothing it is asked for

window truncates the Gaussian at a fixed number of voxels while sigma is in millimetres, so the old default of window = 1 cut the kernel off at one voxel either side regardless of how wide the kernel was meant to be. Measured by smoothing an impulse and reading the width back out:

                                  delivered   requested   shortfall
gaussian_blur(sigma = 2, window = 1)  3.49 mm    4.71 mm     26%
gaussian_blur(sigma = 3, window = 1)  3.70 mm    7.06 mm     48%
gaussian_blur(sigma = 4, window = 1)  3.76 mm    9.42 mm     60%

Worse, the shortfall depended on the voxel size – sigma = 2, window = 1 delivered 3.49 mm FWHM on 2 mm voxels and 1.88 mm on 1 mm voxels, so the same call smoothed two acquisitions of the same brain differently – and nothing warned.

window now defaults to NULL, meaning “derive the half-width from sigma and the voxel size”, using scipy.ndimage.gaussian_filter’s rule (int(truncate * sigma / spacing + 0.5), truncate = 4). Requests are now honoured to better than 1% of FWHM on isotropic and anisotropic grids alike.

as_canonical() no longer discards data

Reorienting to RAS is a permutation and a flip of the array: the voxels are the same voxels, relabelled. It was implemented by building a target space and handing it to resample(), and reorient() rotated the affine while keeping the source dimensions. When the reorientation permuted axes – which is the whole point – the output grid did not contain the data:

                     in                 out            nonzero kept
nilearn reorder_img  AIR 60x72x56  ->   RAS 56x60x72      100%
as_canonical (was)   AIR 60x72x56  ->   RAS 60x72x56       76.4%

Nearly a quarter of the image was thrown away silently. Going through a registration library also meant values were interpolated where they should merely have moved, and that images smaller than four voxels on any axis were refused outright.

as_canonical() now permutes and flips the array and rebuilds the affine. It is exact (the output is a permutation of the input, checked with identical()), works at any size, matches nibabel::as_closest_canonical voxel for voxel, and is about 40x faster. reorient() on a NeuroSpace permutes dim and spacing along with the affine.

conn_comp()’s local maxima were never pruned

local_maxima_dist had no effect at any value. .pruneCoords() read the neighbour distances out of dbscan::kNN()’s result as ret$distances; that component is called dist, and $ does not partially match a name longer than the one it is looking for, so the lookup returned NULL, the comparison became logical(0), and the pruning loop exited on its first pass every time. The local_maxima table listed every in-mask voxel – 25,538 rows for 12 components on a 3 mm map, 25,527 of them from a single component, each one voxel from the next.

Underneath that, the rule did not deliver the minimum distance the argument documents even when it ran. It compared each voxel only against its single nearest neighbour, so three points at 0, 5 and 12 mm valued 5, 1 and 9 kept the 5 and the 9 at local_maxima_dist = 15: the 5 survived because its nearest neighbour was the 1. “The nearest neighbour” was not well defined either – on a voxel grid more than half of all points have several neighbours at exactly the same distance, and which one came back depended on the kd-tree’s traversal order.

The rule is now the one the argument names: a voxel is reported when no other voxel of its component within local_maxima_dist has a larger value. It is deterministic and the reported maxima are at least local_maxima_dist apart by construction.

                        local_maxima rows    closest pair of maxima
26-connect   was              25,538              3 mm (one voxel)
             now                  73              15.1 mm
6-connect    was              25,538              3 mm (one voxel)
             now               5,039              15.1 mm

Connected components in compiled code

conn_comp_3D() was a two-pass union-find written in interpreted R, building a 26x3 neighbour matrix per voxel and running find() as an R closure. src/conncomp.cpp existed but held only a commented-out sketch of it. Against scipy.ndimage.label() it was 690x slower on a small volume and 4,450x slower on a 1 mm brain – fifteen minutes for one labelling – and the ratio grew with size.

                                     scipy    before     after
40x48x38    (7,312 in-mask)         0.001 s    1.4 s     0.024 s
91x109x91  (90,357 in-mask)         0.016 s   38.4 s     0.31 s
182x218x182 (720,916 in-mask)       0.152 s  916 s       3.37 s
labelling alone, 91x109x91          0.016 s      --      0.027 s

The ratio to scipy no longer grows with size: it is flat at 18-22x across a 100x range in voxel count, and the labelling itself is at parity. What is left is the per-cluster voxel lists, the cluster table and the local maxima that conn_comp() returns and scipy.ndimage.label() does not; cluster_table = FALSE, local_maxima = FALSE halves the remainder.

The labelling is unchanged: the compiled version reproduces the R one’s labels, sizes and size-descending numbering including its tie-break, pinned in tests/testthat/test-conncomp-equivalence.R against the previous implementation, which is kept there as the reference.

conn_comp()’s wrapper was rewritten alongside it. It re-derived a coordinate matrix per component with as.matrix(), and then subset a data frame once per component to build the voxel lists, which cost 0.39 s at 4,989 components against 0.03 s for the labelling. ClusteredNeuroVol() no longer fills its cluster map with a name-by-name list lookup, which was quadratic in the cluster count and cost 2.2 s at 12,396 clusters.

automask() spent 92% of its time in the labeller and is 53x faster (14.0 s to 0.26 s on a 60-volume run) with an unchanged result.

The two searchlight iterators agree on their centres

searchlight_coords(mask, radius) returned one element per voxel in the grid where searchlight(mask, radius) returned one per nonzero mask voxel – 241,920 against 83,563 on an ordinary brain mask – although both documented the latter. searchlight_coords() now centres on nonzero mask voxels like every other iterator in the package, and nonzero decides only what each searchlight contains, as it does in searchlight(). It also validates its arguments.

A resample target that misses its source says so

NeuroSpace(dim, spacing, origin) builds a positive-diagonal affine, which mirrors the usual LAS-to-RAS source. Handing one to resample() as a target produced a near-empty image with no warning. resample() now compares the two world bounding boxes and warns when they barely overlap, naming the likely cause.

Reshapes that copied the payload

as.matrix() on a DenseNeuroVec duplicated the whole image to reshape it to voxels-by-time – 116 MB for a 60-volume run, and over half of automask()’s runtime. A reshape is a change of dim, not of memory, and R still copies on the first write, so the source is untouched either way. as.matrix() is now effectively free, mean() goes through .rowMeans() with no reshape at all (0.21 s to 0.12 s), and scale_series() uses the same route.

resample() no longer scans the image to infer an output datatype for headers whose geometry is all it uses.

Smoothing a 4-D image runs in compiled code

gaussian_blur() on a NeuroVec looped its volumes in R. Each iteration copied a volume out of the run and the result back in, re-allocated the kernel’s two scratch buffers, and – with a mask – recomputed the in-mask weight volume, which is the same for every volume. Smoothing 60 volumes cost 1.66 s against nilearn’s 0.51-0.71 s.

The separable kernel now has a 4-D driver that walks the run itself: nothing is copied per volume, the scratch buffers are allocated per worker, the in-mask weights are computed once, and volumes run in parallel (RcppParallel, so RcppParallel::setThreadOptions() controls it). 0.20 s, which is 2.6-3.6x faster than nilearn rather than 2.3x slower. The arithmetic per volume is unchanged and happens in the same order, so the result is bit-identical to smoothing the volumes one at a time – pinned in tests/testthat/test-processing-defects.R over both engines, both mask conventions and both normalize settings.

The dense kernel, which gaussian_blur() chooses when the mask is a small fraction of the volume, has no 4-D form and still loops; everything invariant across volumes is hoisted out of that loop.

neuroim2 0.18.0

Bug Fixes

NIfTI I/O rewritten on a compiled core

The read and write paths no longer go through readBin/writeBin. A single translation unit (src/nifti_data_io.cpp) converts between the file’s stored type and double in one pass, for plain and gzipped files alike, and both read_vol()/read_vec() and write_vol()/write_vec() now call it. On this machine, against nibabel 5.4.2 reading the same files:

                                   before     after    nibabel
read 3-D int16 .nii  (14 MB)        0.46 s    0.049 s    0.027 s
read 3-D float32 .nii.gz            0.45 s    0.220 s    0.247 s
read 4-D int16 .nii  (56 MB)        3.09 s    0.371 s    0.107 s
read 4-D int16 .nii.gz              2.81 s    0.750 s    0.585 s

Two things were behind the 4-D case in particular. read_mapped_vols() built a 236 MB vector of double indices via outer() simply to enumerate a contiguous range, gathered it element by element through the memory map, and handed back a time-by-voxel matrix that DenseNeuroVec() then transposed straight back. It now reads whole volumes sequentially and returns them voxels-by-volumes, the layout the object stores. And building the object no longer copies the payload: an S4 class extending array is that array with the slots as attributes, so the read path sets them in place, the same technique the ROI constructors already used. The equivalence to new() is pinned by identical() in the test suite.

read_mapped_vols() is internal, but its return layout changed from [volumes x voxels] to [voxels x volumes]. Two callers were also missing scl_slope/scl_inter scaling entirely; both are fixed.

Every NIfTI datatype the standard defines

int8, uint16, uint32, int64 and uint64 are read and written. They were rejected before with “Unsupported NIfTI data-type code” — uint16 in particular is not exotic.

int32 images no longer lose voxels holding -2147483648. That bit pattern is R’s NA_integer_, so routing the payload through an R integer vector turned it into NA and dropped it from every downstream statistic; the compiled reader lands in double and never sees an R integer.

Complex (DT_COMPLEX64/128/256), colour (DT_RGB24, DT_RGBA32), DT_FLOAT128 and 1-bit DT_BINARY are still not read, but the error now names the type and says why: a NeuroVol holds one real number per voxel, and inventing a projection — a modulus, a luminance — would be answering a question the caller did not ask.

Writing preserves the header

as_nifti_header() built every header from defaults plus the object’s geometry, so a read-write round trip silently discarded everything else. It now starts from the header the image was read from. Concretely, a round trip used to lose:

All of them survive now. The source header travels on a new header slot on NeuroObj, which every image class inherits; header(x) reads it, and that method now accepts images as its documentation always claimed.

write_vol()/write_vec() also default to the source file’s datatype rather than FLOAT, when the values are still exactly representable in it. An unmodified int16 image round-trips byte for byte instead of doubling in size; computed data that no longer fits is promoted to FLOAT rather than quietly quantised.

Integer output is scaled, not truncated

write_vol(x, f, data_type = "SHORT") went through writeBin(as.integer(x), ...), which truncates toward zero without a warning: data spanning ±3.7 came back with a maximum error of 0.998. scl_slope and scl_inter are now derived to fit the target type, giving 5.4e-05 on the same data — the same precision nibabel produces. Data that is already integral and in range is still written unscaled, so masks, label volumes and atlases stay exact.

NIfTI-2

Read and written, plain or gzipped. write_vol()/write_vec() take version = 2 or format = "NIFTI2", and select it automatically when a dimension exceeds 32767, which NIfTI-1’s 16-bit dim field cannot represent. Files written this way are read correctly by nibabel, including the affine, TR and units.

Conformance and diagnostics

Filtering

gaussian_blur() accepts a NeuroVec and smooths it volume by volume, with the argument checks, kernel and mask resolved once instead of per volume.

The separable kernel also got two exact optimisations. Its y and z passes accumulate over whole contiguous rows and slices rather than walking a column with stride d0 (or d0 * d1), which touched a fresh cache line per iteration; the additions for each output element still happen in the same order, so the result is bit-identical. And when the mask covers every voxel — what gaussian_blur(vol) with no mask means — the mask-insulation denominator is the zero-padded convolution of a constant, which factorises into three per-axis edge corrections, replacing three full-volume passes with d0 + d1 + d2 additions.

                                          before     after   nibabel/scipy
3-D, FWHM 6 mm, no mask                   0.60 s    0.21 s      0.13 s
4-D, FWHM 6 mm, 200 volumes               2.04 s    0.97 s      0.57 s

Equivalence to the dense kernel is checked over 432 configurations of dimension, window, sigma, spacing, mask and normalize in the test suite.

Documentation

Other changes

neuroim2 0.17.0

Testing

Documentation

Performance

Bug Fixes

Improvements

neuroim2 0.16.0

New Features

Improvements

Testing

neuroim2 0.15.0

New Features

Improvements

Testing

neuroim2 0.14.0

New Features

Testing

Documentation

Improvements

Testing

neuroim2 0.12.0

Bug Fixes

New Features

Testing

neuroim2 0.11.0

Bug Fixes

Dependency Changes

Improvements

Testing

Testing

Documentation

neuroim2 0.10.0

neuroim2 0.9.1

neuroim2 0.9.0

neuroim2 0.8.7

neuroim2 0.8.5

neuroim2 0.8.4

neuroim2 0.8.3

neuroim2 0.8.2

neuroim2 0.8.1