Conversation
New command that fills "invalid" voxels (non-finite by default, plus any flagged via -mask) with intensities predicted from surrounding valid data. Each 3D volume is processed independently; the set of voxels to impute is serialised via Voxel2Vector into the unknowns of a dense linear system solved per volume. Six methods are provided: laplacian / laplaciansq / biharmonic / spring (translated and generalised to 3D from the MatLab Inpaint_nans methods 0, 2, 3 and 4) and two new isotropic2 / isotropic4 methods that lift the 13-direction, spherical-harmonic-weighted finite-difference stencil of mrsense1fix into the same global solve. The numerical machinery lives in a new header-only module cpp/core/also/impute.h (MR::Impute): centred Eigen::Tensor stencils, a shared assembly/solve base class with per-method derived classes, and a duplicated SH-weighted stencil generator (mrsense1fix.cpp left untouched). Boundary handling reduces each operator to a form that remains exact for smooth fields (per-axis reduction for the Laplacian operators; skip-unless-fully-supported for directional stencils), so interior holes reproduce harmonic fields exactly while boundary voxels still extrapolate. Adds six binary regression tests (round-trip per method, -mask, 4D multi-volume, NaN/Inf default, degenerate guard, boundary extrapolation) and the generated command reference documentation. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add MR::Interp::Deform, a cubic interpolator for non-linear deformation fields that may contain non-finite voxels: at construction it uses the MR::Impute machinery to fill a scratch buffer (empirical data where finite, imputed values where not) and a bitwise validity mask governing voxel()/scanner() return values, so interpolation kernels never ingest NaNs. To make imputation viable as extrapolation across wide NaN borders rather than mere isolated-hole filling, extend the imputation library with a universal-kriging polynomial detrend (affine and quadratic), a biharmonic operator with natural boundary conditions, and a Hessian (Frobenius) energy method whose free boundaries extrapolate boundary trends without bias. Deform defaults to the Hessian method with a quadratic detrend; standalone mrimpute usage defaults to no detrending. Relocate the imputation header from core/also to core/algo. Regression tests and regenerated command documentation are included. Session prompts: 1. > Build an MR::Interp adapter dedicated to the application of non-linear deformation fields with cubic interpolation. The adapter takes as input a non-linear deformation field image, which may contain within it voxels with non-finite values. During construction the adapter will make use of the imputation capabilities introduced with new command mrimpute to fill a scratch buffer image that contains a combination of the input empirical warp field data where finite and imputed values where not (this facilitates the use of cubic interpolation of the deformation field without incurring data loss where non-finite values enter the interpolation kernel). The class should also generate a bitwise mask scratch image where the value reflects the presence of non-finite values in the input deformation field. This should be used to determine the return value of the .voxel() and .scanner() functions: they should return true only if the voxel in which the vertex resides is true in the mask. 2. > Class MR::Interp::Deform has been added to improve the efficacy of non-linear deformation field application, by permitting cubic interpolation even where the boundary of the image is filled with NaNs due to the absence of valid image data from the other image participanting in the registration. However many of the imputation algorithms devised during prior addition of the mrimpute command are principally designed for imputing isolated pixels / voxels, not extrapolating beyond the borders of valid image data. Investigate online any algorithms for extrapolation of image data beyond the borders of the image, which would therefore also be applicable in the situation where there are many NaN values around the borders of the input image. Place particular emphasis on algorithms that can be implemented using linear algebra, or would otherwise require minimal code modification without taking on additional third-party libraries. Ignore any methods based on deep learning. > Implement recommendation 2 (universal-kriging: detrend -> impute -> re-add trend), then recommendation 1 (modify biharmonic with natural boundary conditions). > Are you still there? > Implement regression test. Make the default behaviour for MR::Interp::Deform to include an affine detrend, but the default for any other usage of imputation should be no detrending. Run syntax and code format checks. Regenerate user documentation. 3. > In MR::Interp::Deform, a global affine de-trend has been added to improve the quality of imputation particularly beyond the bounds of the valid deformation region where the image becomes NaN-filled. This global affine may however be only a poor estimate of the local pattern of deformation at that particular edge of the valid data region. Comment on the feasibility and likely computational expense of performing a local affine detrending, where for each voxel to be imputed, the detrend is based on a local affine estimate where the weighted contribution of each prioximal voxel with valid image data is suitably proportional to the distance of that proximal voxel to the one being imputed. Also search the online literatore for any alternative proposals for solving this issue. Also present any alternative approaches that Claude may devise in this regard. 4. > 1. Implement the Hessian-energy natural-BC formulation from Stein et al, and make this the default imputation method for MR::Interp::Deform. Ensure that additional tests are added, and that documentation pages are re-generated. > 2. Extend the enumeration that currently triggers presence / absence of global affine detrending to detrend a quadratic rather than linear affine as a third option; default to this option for MR::Interp::Deform. 5. > Move cpp/core/also/impute.h to cpp/core/algo/impute.h and re-establish compilation. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
mrimpute now processes multi-volume images with one thread per 3D volume, looping over all non-spatial axes via a ThreadedLoop so that independent volumes are imputed concurrently. It additionally tallies, per volume, the number of voxels for which imputation was performed, reporting the full statistics through INFO and warning when no imputation occurred in any volume. New regression tests cover the per-volume statistics reporting and the no-imputation warning. Session prompts: 1. > In cpp/cmd/mrimpute.cpp, utilise a ThreadedLoop that causes the > command to utilise one thread per volume (ie. threads loop over > all non-spatial axes). 2. > mrimpute should additionally generate statistics on the number of > voxels for which imputation took place within each 3D volume. The > full statistics should be reported using the INFO() macro. If no > imputation occurred in any volume, then WARN() should be invoked. 3. > Add tests for both per-volume test statistics, and to ensure that > the warning message in the absence of imputation is emitted. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
Adds an optional 26-voxel-neighbourhood connectivity mode to the mask dilation and erosion filters, selectable through a new maskfilter -connectivity option with accompanying tests. Reworks MR::Interp::Deform so that construction imputes only the two-voxel halo required to support cubic interpolation around valid voxels rather than every non-finite voxel: a validity mask is computed and its grid padded where necessary so valid voxels retain a two-voxel margin, the mask is dilated twice at 26-connectivity, and the resulting halo is filled by the imputation machinery using the valid data alone, with the solve's field of view restricted to that halo. Session prompts: 1. > 1. Modify the mask dilation / erosion filters to accept a Boolean toggle between only selecting the 6 shared-face adjacent voxels (current behaviour) and the 26 shared-corner adjacent voxels (new behaviour). > 2. Enable access to this functionality in the maskfilter command through option -connectivity. Ensure that new tests are added; do not attempt to run the tests as reference data will not yet exist; the user will generate this. > 3. In MR::Interp::Deform, requisite modifications at construction will be performed to guarantee that cubic interpolation can be performed in all voxels with valid input values, but imputation computation is not performed for voxels where it is not necessary. > 3.1. A mask of all voxels with valid data will be generated. > 3.2. If necessary, a new voxel grid will be generated, which guarantees that all voxels in this mask are at least two voxels away from the image field of view edge. The initial mask must then be resampled onto this new grid. > 3.3. The mask dilation filter, modified in step 1 to have a mode invoking 26-nearest-neighbour connectivity, is then executed with two passes. > 3.4. Generate a mask of voxels that are present in the output of step 3.3 but absent in its input. > 3.5. Use the imputation oeprations added on this code branch to generate a new scratch image. Its voxel grid is that determined in step 3.2. For voxels with valid input image data, the values in this scratch image are taken from the input image (shifted in index if necessary to conform to step 3.2). For voxels in the mask determined in step 3.4, the values are determined through imputation. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
Rework the image comparison machinery in cpp/core/image_diff.h so that differences are detected not only by supra-threshold value deviations but also by inconsistency in the presence or absence of NaN values, correctly handling complex-valued images. Comparison arithmetic is now selected via SFINAE at the highest width and precision of the two input types: integers compared as integers, values promoted to complex floats only when an input is truly complex, and NaN-location checks restricted to floating-point data. A new bitwise comparison performs exact equality without any tolerance. The testing_diff_image command now reads both headers, determines the appropriate comparison data type, instantiates the matching MR::Image<> classes, and rejects tolerance options when the inputs are bitwise. Session prompts: 1. > 1. Complete porting of cpp/core/image_diff.h to ensure that it is > not merely supra-threshold differences between two images that > are detected, but also an inconsistency in presence / absence of > NaN values. This must currently be updated to properly handle > complex-valued images. > 2. Perform SFINAE template specialisation of image comparison > operations in cpp/core/image_diff.h. Numerical comparisons should > be performed using only the highest width / precision of the two > input image data types. Eg. If both are integers, then perform > integer comparison. Only cast image values to complex-valued > floats if at least one of the images is truly of complex > floating-point type. Only perform check that either both or > neither image is NaN if data type is floating-point. > 3. Add bitwise image comparison to cpp/core/image_diff.h, which by > definition cannot have a comparison tolerance. > 4. Update command testing/tools/testing_diff_image.cpp to read the > headers of both images, determine the appropriate data type to > use for the comparison, and instantiate MR::Image<> classes and > the comparison function accordingly. Ensure that tolerance > command-line options are not used if input images are bitwise. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Reimplement non-linear deformation-field inversion (warpinvert / Registration::Warp::invert_deformation) on top of the imputation-aware Interp::Deform interpolator. Per-iteration optimisation now uses a Newton search over the cubic value and Jacobian with step-halving to stay within the valid input region, and a consecutive-blocked heuristic that marks voxels NaN where the warp cannot be inverted or its origin lies in extrapolated data. The inverse is seeded by a scatter of non-extrapolated input voxels followed by a greedy 26-neighbour gap fill, with debug image dumps of its intermediate state. Interp::Deform gains a processing-type template parameter for value-and-gradient access and an image-based constructor; a small identity round-trip regression test is added. Session prompts: 1. > Modify the code responsible for inversion of non-linear deformation fields. > 1. Make use of new class MR::Interp::Deform. The presence of NaNs in the output image should be based primarily on where the warp can be inverted, and secondarily on where the origin point is reported by the interpolator as residing in a voxel containing valid input image data without extrapolation; ie. if a voxel of the inverted warp field maps to an input image reconstructed via extrapolation, its value in the inverted warp should be NaN. During inversion of the warp field, per-iteration optimisation should be based on a Newton search rather than the existing projection. Use the cubic interpolator that yields both values and gradients to access this information. Where the Newton search projects to a location outside of the valid input image region, the step size should be iteratively halved until the projection no longer exits the valid input data region. A heuristic is required to detect when successive iterations continue to attempt to drive the projected point outside of the valid region of the input image, at which point it should be assumed not possible to map the inverted warp field voxel to the origin point in the input warp image; propose suitable heuristics and prompt user for confirmation. > 2. Implement a better-informed initialisation of the inverted field. On the inverted field voxel grid, generate two scratch images: one is the minimal distance between the centre of that voxel and the location from which a non-extrapolated input warp field voxel samples, and the other is the scanner-space position of that nearest input warp field voxel. After looping over all non-extrapolated input warp field voxels, it is necessary to fill gaps in this image so that all output inverted warp field voxels have an appropriate initialisation location. The user proposal is: rank-order uninitialised voxels by the number of voxels in their 26-voxel shared-corner neighbourhood that already possess an initial value; for all voxels with the equal largest number of such voxels, compute the average value as the initial value for that output voxel---unless this maps to a location in the input image that cannot be interpolated, in which case choose the input voxel to which the distance between the mapped location in the output image and the voxel under consideration is minimal; update the rank-ordering of all adjacent uninitialised voxels; iterate until all voxels possess an initial value. Include within this code invocations of MR::Image::dump_to_mrtrix_file() so that the internal operation of the algorithm can be debugged. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Interp::Deform previously filled the halo around the valid region of a deformation field with a dense global imputation solve, which becomes computationally infeasible for large non-linear warps as the halo's surface-area voxel count grows. Replace it with rank-ordered local polynomial extrapolation: halo voxels are filled outward from the valid region in onion-peel layers, each voxel fitted by a small least-squares affine or quadratic model over nearby finite samples, exploiting the field being locally low-order. The padded grid and scratch buffer are retained, and the imputation PDE machinery remains available to mrimpute untouched. This makes large-field extrapolation tractable while keeping the cubic interpolator well-posed at the boundary. Session prompts: 1. > A parallel Claude session has identified that the approach currently in > use in MR::Interp::Deform, which involves padding the region of valid > image data and utilising a dense imputation solver, becomes > computationally infeasible for large non-linear deformation fields. > Propose an alternative solution that will improve the computational > feasibility of the approach. The generation of a padded voxel grid and > corresponding scratch image storing the extrapolated values will remain, > since it minimises re-computation and will be better-posed even for warp > fields without missing data than will be the default clamp operation in > the cubic interpolator. Comment on whether use of a sparse matrix solver > would improve computational feasibility. Otherwise, generate proposals > for separable extrapolation of the deformation field. Extrapolation of > data for voxels requiring it could be rank-ordered on the availability of > valid data in their vicinity. The deformation field can be assumed to be > locally affine or quadratic, and so could be solved primarily based on > this premise rather than this being merely a de-trending operation for > PDE-based imputation. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Replace the dense coefficient matrix and dense factorisations in the mrimpute imputation machinery with Eigen sparse equivalents. The global system, assembled per 3D volume from finite-difference stencils, is inherently sparse since each equation couples only a handful of neighbouring unknowns. The matrix is now built from triplets into an Eigen::SparseMatrix, the least-squares methods are solved with SparseQR and the square method with SparseLU, and the underdetermined-column guard becomes a column-major inner-iterator scan. Because each system is full column rank, the solution is unique and the results match the former dense solvers to numerical precision; all command tests and the extrapolate unit tests pass unchanged. Session prompts: 1. > Refactor the machinery behind command mrimpute to make use of > Eigen::SparseMatrix<>. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Refine the non-linear deformation-field inversion so that the valid-voxel set produced by inverting a structured warp twice is governed predictably rather than scattering across validity boundaries, adding diagnostic controls over halo extrapolation and the validity policy. Extend the testing image-comparison machinery so that differences can be assessed within a binary mask: the per-voxel comparison is factored into shared helpers and gains mask-aware overloads of every comparison mode, with a 3D mask broadcast through Adapter::Replicate so it can modulate a 4D / 5D comparison, exposed via a new testing_diff_image -mask option. A regression test double-inverts the NaN-slab warp.mif, confirms the round-trip agrees in value within the region valid in both fields, and bounds the drift in the valid-voxel set. Session prompts: 1. > When a non-linear deformation field is inverted, and then inverted again, it is expected that the set of voxels containing valid values may not be precisely equivalent. There is however greater discrepancy than expected given the current implementation. For instance, for a warp field where one slice contains all valid values and the next slice contains all invalid values, this double-inversion yields scattered voxels with valid and invalid values for both slices on either side of the transition in the input image. But with a combination of 2 x dilation to extrapolate values suitable for cubic interpolation, and retaining data in all output voxels for which the voxel inside which the projected location resides in the input image contains valid data, unambiguous acceptance vs. rejection was expected for these two slices. Evaluate the implementation and propose hypothetical mechanisms by which the presence or absence of valid data in a twice-inverted non-linear deformation field may be either biased or imprecise relative to user expectation. > > Image "warp.mif" in BinariesTestData can be used to manifest the fault without requiring construction of a new test input image. > > Test data has been updated to contain a version of warp.mif that contains the requisite NaN slab structure. 2. > Round-trip test on warpinvert should not be passing given that the set of voxels containing finite image data does change in response to the round trip, and the library functions for comparing two images to detect regressions were recently updated to flag discrepancies in the presence of NaN values between the two images. Diagnose. 3. > Add a test that involves a double-inversion of image warp.mif. > 1. In cpp/core/image_diff.h, add alternative versions of relevant functions that provide the capability to specify a bitwise mask image, with comparisons performed only within that mask. The mask image should itself go through Adapter::Replicate, so that a 3D mask volume can be used to modulate which voxels contribute to the operation of a 4D / 5D image comparison. > 2. Generate a mask of voxels that are finite in both warp.mif and the twice-inverted version of such. > 3. Compute the maximal discrepancy between these two images within this mask; report that value here, round up to the nearest one-signifant-figure value and apply that as the tolerance within the test. > 4. Within the test, compute the number of voxels that are finite in one image and non-finite in the other; round up to the nearest one-significant-figure value and apply that as a threshold on the number of such voxels found within the test. 4. > For image_based differences, use testing_diff_image. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Rename the imputation-aware cubic deformation-field interpolator from Interp::Deform to Interp::Warp and generalize it to serve both deformation fields (absolute positions) and displacement fields (offsets), since its validity-masking and halo-extrapolation logic operates only on the raw field components. Give it a single templated field constructor and add a WarpReslice adapter so the generic reslice machinery can drive it, then route the deformation-field interpolation paths — the warp regrid in Filter::warp and transformcompose — through it so they gain robustness to non-finite voxels and cubic accuracy. The halo-extrapolation helper is renamed to match. Conversion of ComposeHalfwayKernel and tcktransform is deferred with in-code reminders, and the regression tests whose reference data is expected to drift at warp boundaries are annotated. Session prompts: 1. > Examine the session log in "~/OneDrive/Documents/MRtrix3/warp_interp_audit.md". For all code paths that involve performing interpolation of a non-linear deformation warp image, the goal of this session is for MR::Interp::Deform to be invoked when reading from the input warp. Propose a refactoring plan to achieve this. 2. > 1. Defer modification of ComposeHalfwayKernel, but remind user to investigate the interaction of this code with the intent behind MR::Interp::Deform. > 2. Add comments to those test files where regression of tests due to change of behaviour is expected given the planned code modifications. > 3. Augment plan with the refactoring and renaming of code relating to interpolation of deformation fields to also work for displacement fields. Most code should work for both forms of data. > 4. Proceed with rest of plan. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
tcktransform now guarantees valid warp information at every streamline vertex, however far it lies beyond the region over which the input deformation field was estimated. The field is cubically interpolated against a buffer whose invalid voxels are filled by rank-ordered local polynomial extrapolation, and that buffer is shared across worker threads and grown on demand (under a shared mutex) whenever a streamline reaches beyond it. As a result streamlines are transformed in full rather than being truncated at the field-of-view edge, and the command warns when any such growth beyond the input data occurred. This replaces the previous linear-interpolation path that silently dropped out-of-bounds vertices. Session prompts: 1. > Command tcktransform is to be augmented with a modified version of > MR::Interp::Warp. Currently that class extrapolates non-linear warp > data only for a 2-voxel halo around the region in which the input > warp image is valid. For tcktransform, it will be asserted that valid > warp information be available for all input streamline vertices, no > matter how far beyond the valid input data they are. > 1. On construction, instead of only initialising warp field data for > a 2-voxel halo, after generating a padded voxel grid (if > necessary), all voxels that do not contain valid non-linear warp > field data will have an extrapolated value derived, following the > existing iterative rank-ordering strategy. > 2. Any time when the cubic interpolation of the warp field would > require relying on the cubic interpolator's clamp() member function > to prevent attempting to access image data beyond the field of view > (ie. not at least 2 voxels between that in which the vertex resides > and the last voxel of the image), the following will occur: > 2.1. The thread will look at all vertices of the streamline > currently being processed, to discover its axis-aligned > bounding box. > 2.2. The thread will lock a std::shared_mutex. > 2.3. The padding for the scratch image in the MR::Interp::Warp class > will be expanded beyond its current state so as to support > cubic interpolation of the axis-aligned bounding box of the > streamline. > 2.4. Values for all new voxels on this voxel grid are to be computed > using the existing iterative rank-ordered process. > 3. If at least one such additional padding occurred during > transformation of streamline data, issue a warning to the user that > some streamline transformation was based on an extrapolation of the > non-linear warp field beyond the input data. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
Continued the warp-field interpolation refactor by routing the two remaining registration sites the audit flagged onto the imputation-aware cubic interpolator. The multi-resolution displacement-field upsampling in non-linear registration now uses a new Extrapolated validity policy on Interp::Warp so cubic interpolation is well-posed and finite at every output voxel, while the demons hot-path inversion is deliberately kept linear and documented as such. The warpinvert command's displacement path now inverts through the same cubic Newton machinery as its deformation path, with the -extrapolate and -validity options applying to both, backed by a new regression test. The mrregister non-linear reference data is stale as a result, and the affected tests carry comments describing how to regenerate it. Session prompts: 1. > In one or more previous sessions, the interpolation of non-linear warp > field image data was modified in order to make use of cubic interpolation > without incurring loss of valid data near image field of view edges / > non-finite input values. There were some components of the registration > code that were identified as being relevant to the context of > interpolation of non-linear warp field images, but for which the > appropriate refactoring plan was uncertain. The goal of this session is to > revisit this topic given recent modifications to the code, identify those > sections that have not yet been updated, and perform a deeper dive. > The results of the initial audit can be found in > "~/OneDrive/Documents/MRtrix3/warp_interp_audit.md". Cross-reference that > report with the current state of the code. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Continue migrating non-linear registration from linear to cubic interpolation, targeting the remaining "full" warp field format and the internal Demons field-inversion machinery. Introduce a PaddedField abstraction that encapsulates the padded voxel grid and its extrapolation so the composition and inversion code can sample the warp at sub-voxel locations cubically. Update Interp::Warp, the warp compose/invert/ extrapolate helpers and the nonlinear registration driver to route through this path, and add targeted unit tests for the new component. Reference expectations for the rigidnonlinear FOD registration tests are updated to reflect the cubic sampling. Session prompts: 1. > This session will continue progress on the refactoring of code relating to non-linear registration utilising cubic rather than linear interpolation. Prior session logs can be found in ~/OneDrive/Documents/MRtrix3/warp_interp_audit.md and ~/OneDrive/Documents/MRtrix3/warp_interp_progress.md. The remaining code that is still utilising linear interpolation relates to the so-called "full" warp field format and the internal Demons calculations including field inversion. > 1. Confirm and tabulate those sections of code that are still utilising linear interpolation to sample a non-linear warp field at sub-voxel locations. > 2. Determine and report the deepest code loop that invokes field inversion and/or interpolation, and report the relevant code location and the approximate number of times that code would be executed for a typical mrregister invocation with default parameters. > 3. Generate proposals for the refactoring of this code to utilise cubic rather than linear interpolation. One proposal from the user is: instead of generating a padded voxel grid and then extrapolating, one could instead generate the padded voxel grid at the very first iteration only, with subsequent iterations initially performing resampling of the whole image using cubic interpolation with the cubic interpolator's built-in image field of view clamping, then selectively NaN-filling and extrapolating the slabs corresponding to the two excess voxels at the border of the image. Generate other hypothetical candidate strategies if possible. Generated-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The migration of the non-linear registration warp composition and inversion path to tri-cubic Catmull-Rom interpolation can violate the diffeomorphic-demons requirement that every composed displacement field have an everywhere-positive Jacobian determinant, because cubic interpolation overshoots near steep transitions in a way linear interpolation cannot. This work characterised the failure mode, derived the exact 3D Catmull-Rom injectivity bound (B < 32/225, a per-sub-step divisor of 225/32), and tightened the scaling-and-squaring step bound accordingly so the algorithm never feeds a folding-amplitude update to the interpolator. A debug-only check reinstates the per-iteration Jacobian-determinant assertion and additionally logs the realised minimum determinant under -debug, so the conservatism of the bound can be assessed against real image data. The deliberately-failing unit-test demonstrator that first exposed the violation is retained as commented reference documentation rather than an executed test. Session prompts: 1. > In the previous commit, deep internal operation of the image registration framework was modified to utilise cubic rather than linear interpolation. A potential consequence of this modification flagged by pre-existing code comments and by Claude's analysis is that the use of higher-order interpolation could result in violation of internal requirements of the Demons algorithm. > - Scan the logs of prior sessions executed from this working directory. > - Identify and concisely state the nature of the potential violation. > - Hypothesize what features of the image data content could potentially lead to displacement field content that could violate that requirement. > - Describe how one might synthesize sample data that could expose the violation, and whether the synthesis of such data could be achieved by Claude. > - Generate a plan for construction of a test utilising such data that could expose failure with the current code whereas execution would succeed for prior code utilising linear interpolation. 2. > Proceed with implementation. 3. > Evaluate prospective mechanisms by which the diffeomorphism violation within the diffeomorphic demons registration implementation may be mitigated. It is necessary to reduce the magnitude of the field update in such a way that the diffeomorphism is not violated, for any possible sub-voxel position sampled by the tri-cubic interpolator. Determine, either from literature or from first principles, whether it is more suitable to seek a local modulation of field updates to prevent the violation, or whether the field update should be modulated using a single scalar across the whole image. If the Catmull-Rom spline has an algebraic form, such that the location and magnitude of extrema are known, then perhaps that logic can be expanded to three dimensions, such that for every 4x4x4 cubic interpolation window, one can compute algebraically the extreme value within the centre of that window. It should then be possible to derive the extent to which one would need to modulate the magnitude of the displacement field updates in order to keep the minimum Jacobian determinant within the window greater than some minimum threshold. From there would would either need to e.g. solve a sparse system to find the minimal sum of deviation of voxel intensities that yields a diffeomorphic update, or find the maximal requisite scalar modulation for any window and apply that to the whole field update. 4. > Dispatch two sub-agents: > 1. Implement Tier 1, with code comments describing what each candidate value of K corresponds to, and reinstate the check for diffeomorphism, running only when compiled in debug mode; finally re-run tests. > 2. Determine the 3D Catmull-Rom injectivity bound. 5. > 1. Set the default K to 225/32. > 2. Over and above checking for non-diffeomorphism, also log the minimum Jacobian determinant anywhere in the image for each iteration. When compiled in debug mode and `-debug` is specified, write to the terminal this information at an appropriate level of detail, so that if the rigorous injectivity bounds are unnecessarily pessimistic with respect to real image data this may be detected and modulated. 6. > Comment out the content of the failing Demons test, adding a comment on what purpose it served and why it should not be executed, and remove addition of the test within the relevant CMakeFiles.txt. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
Audited every location where a Jacobian or its determinant is derived from a non-linear warp field, establishing that at voxel centres the tri-cubic Catmull-Rom derivative is algebraically identical to the existing finite-difference scheme, so the voxel-centre consumers need no change. Corrected the cubic and linear interpolators' scanner-space gradient transform, whose operand order was wrong for anisotropic voxels while coincidentally correct for isotropic ones. Where a deformation field is resliced onto a new grid before image warping, its Jacobian determinant is now evaluated analytically from the source spline during that interpolation rather than finite-differencing the resampled field afterward. The warp adapter consumes this precomputed determinant when available and retains the finite-difference path as a fallback. Session prompts: 1. > Perform an audit of all locations where the Jacobian (including the > determinant of such) is computed from a non-linear warp field. Given > that previously interpolation of such images was always performed > using linear interpolation, many locations of such will be based on > the finite difference approach. Now that cubic interpolation is used > in all instances, the derivatives of the warp field should be > accessible from the interpolator; and now that extrapolation is > performed, those gradients should be well-formed at all > (non-extrapolated) image voxel centres. For each source code > location, report on the current algebraic approach for gradient > quantification in use, what the gradient is utilised for, and whether > the computation is always exclusively at voxel centres or may involve > sub-voxel positions. While sub-voxel cubic interpolation involves a > 4x4x4 neighbourhood, user expectation is that Catmull-Rom derivatives > at voxel centres can likely be computed from the 3x3x3 neighbourhood; > comment. Finally, generate a plan for porting all discovered code > locations such that any interpretation of derivatives of a non-linear > warp field image is based on Catmull-Rom. 2. > Does the Jacobian determinant at a sub-voxel interpolated position > not require more than just the three axis-aligned derivatives > provided by the cubic interpolator with the ValueAndDerivative flag? > It may be possible to get the Jacobian and therefore the determinant > algebraically given knowledge of Catmull-Rom coefficients, but I > would expect there to be interaction terms between axes. If so, the > enum that specifies ValueAndDerivative should be expanded to include > ValueAndDeterminant and ValueAndJacobian. 3. > Implement parts C and E. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
This was referenced Jul 26, 2026
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
1 participant
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.
A second Claude tinker-box that I need to push off my system to focus on more pertinent pursuits. So will likely be sitting in draft form for a while.
The original idea behind this branch is as follows.
I've been concerned for some time about the quality of non-linear registration outcomes with MRtrix3 compared to ANTs.
There's a lot in there that I definitely have not yet fully explored or understood.
One thing that jumped out at me straight away was the interpolation of non-linear warp fields.
This is currently done exclusively using linear interpolation.
Given all we know about interpolation I figured this was likely to be detrimental to some greater or lesser extent.
I believe one reason why linear interpolation is currently used is because interpolation may have to occur: 1) near the FoV edges (e.g. for multi-resolution pyramid); 2) near voxels where no warp field data are present. Cubic interpolation may therefore erode the volume in which interpolation can be performed. For the FoV edge at least there's tricks to get around it, but for NaN-padded warp fields the damage is worse. Another reason is that one must provide guarantees of preserving diffeomorphism during the Demons update, and this is trickier if it must be guaranteed for all sub-voxel trilinearly-interpolated positions; but it was investigated thoroughly here. Definitely let me know if anyone knows of / recalls other reasons.
I decided to try to enable cubic interpolation of warp fields. If a non-linear warp field could be reasonably padded with extrapolated data by two voxels in all directions (whether at the outer FoV edge, or in the NaN-filled voxels arising from registration of images on different grids), then any sub-voxel location in any voxel with finite data in the input warp (displacement or deformation) could undergo cubic interpolation without NaN corruption. And if this could be implemented as a wrapper it could be utilised wherever warp fields are interpreted, with no modification to images on the filesystem.
I initially tried to approach this using imputation, which I'd experimented with years ago. So there's a new command
mrimputehere for that purpose; it doesn't make any attempt to identify voxels for which data should be imputed, it just generates predicted values for voxels in a user-specified mask. This is however the wrong tool to use here: all of the voxels requiring imputed values are at the extremities of the valid data, whereas all of the imputing operators are based on being enveloped by adequate available data. Despite not being used in the end, I will nevertheless leave this command here, as it has other potential applications.What I landed on instead is computing a local linear or quadratic affine based on the surrounding voxels that do have valid data. This is iteratively extended out to produce the requisite two-voxel buffer.
This comes with a particular benefit for
tcktransform. We've historically had issues with streamlines where vertices lie outside the valid warp field region, leading to NaN-filled vertices that were originally interpreted as streamline delimiters by the.tckformat. We've now got some handling in place, but it's a band-aid; streamlines end up fragmented. Here, when such a vertex is found, the extent of the warp field padding is just increased. So every vertex of every streamline gets transformed to a suitable extrapolated position. Indeed even a.tsfwould retain correspondence with the output.tckfile.Testing via double-inversion of a warp field, NaN values no longer infiltrate deeper into the image, and imprecision of the round-trip is halved.
Where progress on this branch ceased was the discovery of #3392. That is more deserving of imminent attention. Yet it may not have been discovered for some time were it not for this train of experimentation.
Unrelated items implemented here: