Fix image gradients w.r.t. scanner for non-1mm & anisotropic voxels - #3392
Draft
Lestropie wants to merge 4 commits into
Draft
Fix image gradients w.r.t. scanner for non-1mm & anisotropic voxels#3392Lestropie wants to merge 4 commits into
Lestropie wants to merge 4 commits into
Conversation
Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
Erroneously initialised Gradient1D::axis using Gradient3D::wrt_scanner, and left Gradient1D::wrt_scanner to take default value. Also simplified code for interpolators that provide gradients to use scanner2voxel transforms rather than manually composing them, which was erroneous prior to parent commit.
Member
Author
|
Thought just hit me: Did the reported empirical improvement of non-linear registration when first upsampling from 2.5mm to 1.25mm have anything to do with reducing the mismatch in magnitude between voxel-space and scanner-space gradients, thereby mitigating the effect of this fault? |
This was referenced Jun 20, 2026
Lestropie
changed the title
A new `-diagnostics_stats` option in `mrregister` writes one tab-delimited row of per-iteration statistics per non-linear (Demons) stage, covering cost, image gradient magnitudes, the proportion of voxels in the gradient-dominated regime, raw and applied update magnitudes, and accumulated displacement field norms. The instrumentation lives in a new `demons_diagnostics.h` header and is inert unless the option is supplied. This infrastructure was used to run a controlled experiment comparing the pre- and post-fix gradient adapter on 2.5 mm test data, where the voxel-size gradient inflation factor is 2.65×. The pre-fix build produced 90.6 % of voxels in the gradient-dominated regime versus 58.6 % with the fix, stalled after 16 instead of 166 productive iterations, and reached a final cost 9.7 % higher while recovering only two-thirds of the deformation magnitude. Results and methodology are documented in two companion Markdown files, including attribution confirming the entire effect traces to the `Gradient3D` construction fix rather than the interpolator reordering. Prompt: > Note that this session is executing on a branch derived from the "master" branch of MRtrix3, which possess a different filesystem structure and build system to that for which the project-level Claude configuration is written. > This PR involves a bug fix to the computation of image gradients in 3D. The goal of the session is to perform a data-driven investigation into the effects of this bug and provide evidence for the benefits of the fix, particularly in the context of non-linear image registration where these gradients are utilised. > In ~/src/test_data/, images moving.mif.gz and template.mif.gz are intended for the testing of image registration software features. They possess a voxel spacing of 2.5mm, which should expose the erroneous scaling of gradients prior to the fix. Using file moving2template.txt as an initial affine transformation to align the two images, perform non-linear registration between the two images, with and without the fix in core/adapter/gradient3D.h. For both software versions, export statistics that provide relevant information about the registration process; consider e.g. maximum displacement per iteration, mean displacement per iteration, cost function value, any other parameters of the Demons formulation that may be of utility in assessment. Write these data to Markdown files in the current working directory. Derive a conclusion about the magnitude of the effect of this bug on non-linear registration performance. Generated-by: Claude Opus <noreply@anthropic.com>
The nonlinear registration multiplied the Demons displacement increment by the mean voxel spacing of the current pyramid level before applying it. This was the correct companion to the pre-fix per-index gradient but becomes spurious once the gradient is expressed in intensity per millimetre: the update is already a displacement in millimetres, and the multiplication produces units of mm² and scales the applied increment by the downsampling factor at each level. Removing it resolves the coarse-level overshoot, the dimensional inconsistency, and the premature iteration-cap binding identified as follow-up items in the companion investigation. The default gradient step is re-derived from 91 registrations across six test cases and updated from 0.5 to 1.5; the permitted range is widened to (0.0001, 10.0]. Under the previous default the coarsest pyramid level delivered no net cost reduction; it now achieves 43 % in that stage, and results for the same problem posed at 1.0, 2.5, and 5.0 mm are bit-identical across gradient steps, which they were not before. Prompt: > Read prior investigation in file "gradient3D_fix_investigation.md". Session goal is to derive suitable modifications to the non-linear registration implementation to accompany the bug fix on this worktree. It is likely that heuristic parameters were previously tuned based on the erroneous implementation and must therefore now be re-derived. > For the multi-resolution pyramid, it is possible that dividing the step size by the downsampling factor may be necessary to prevent overshoot at coarse resolutions. Comment and evaluate on test data. Contrast this description against previous session; unclear whether some tests hinted at this. > If performance at step size of 0.75 is superior across test cases after the fix above, change the software default to this value. > Revise maximum steps per iteration based on performance with the changes above. > Produce preserved derivative images of outcomes of non-linear registration of test data for user communication and reporting. Generated-by: Claude Opus <noreply@anthropic.com>
Lestropie
changed the title
Member
Author
|
Changed issue name: might be easy to disregard if the error is communicated as only manifesting for anisotropic voxels, which would be very rare in this context. The bigger problem is the distinction between magnitude of gradients being w.r.t. voxels or realspace mm. From synthetic testing, the fix outperforms across the board. Behaviour goes from regularly getting trapped on the second iteration of the first resolution to being able to triple the step size and not get trapped. Maximal number of iterations shouldn't need to change. It's possible that current default
|
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.
Discovered by Claude while working on another task.
Gradient3D, erroneous construction ofGradient1Dwas passingbool Gradient3D::wrt_scannertoGradient1D::axis, and leavingGradient1D::wrt_scanneras the default value offalse. In conjunction with use ofTransform::image2scannerto transform the gradient to be with respect to scanner space, this resulted in erroneous scaling for non-1mm voxels and erroneous orientations for anisotropic voxels.check_syntaxboolean usage as function arguments, instead forcing use of class enumerations; that would have caught this mistake as it precludes implicit cast to integer.wrt_spacing" forGradient1Dand "wrt_scanner" for interpolators &Gradient3Dmake sense. Throughout transform handling code we have the distinctions between "voxel", "image" and "scanner" space. To me it would make more sense if these classes conformed to those definitions (suitably banning "scanner" forGradient1D), and the same choices should be available for the interpolators.Gradient3Dis utilised by non-linear registration. Going to take some time to diagnose what the consequences are here; eg. was the bug only detrimental for anisotropic images or was the erroneous scaling by voxel size detrimental also. Certainly somemrregistertests are failing...mrregistertest failures