fod2fixel: Improved quadrature - #3419
Merged
Merged
Conversation
Provides reduced variance in dixel sample weights, superior conditioning which provides robustness against floating-point truncation, and greater rotational invariance for AFD.
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| } | ||
|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]
const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4);
^| } | ||
|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); |
There was a problem hiding this comment.
warning: no header providing "size_t" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:20:
+ #include <cstddef>|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); | ||
| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); |
There was a problem hiding this comment.
warning: no header providing "Eigen::Dynamic" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:20:
+ #include <Eigen/src/Core/util/Constants.h>|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); | ||
| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); |
There was a problem hiding this comment.
warning: no header providing "Eigen::Matrix" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:20:
+ #include <Eigen/src/Core/Matrix.h>| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); | ||
| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); | ||
| for (size_t row = 0; row != dirs.size(); ++row) | ||
| cartesian.row(row) = dirs[row]; |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]
cartesian.row(row) = dirs[row];
^| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); | ||
| for (size_t row = 0; row != dirs.size(); ++row) | ||
| cartesian.row(row) = dirs[row]; | ||
| const Eigen::MatrixXd SH2At = |
There was a problem hiding this comment.
warning: no header providing "Eigen::MatrixXd" is directly included [misc-include-cleaner]
const Eigen::MatrixXd SH2At =
^| for (size_t row = 0; row != dirs.size(); ++row) | ||
| cartesian.row(row) = dirs[row]; | ||
| const Eigen::MatrixXd SH2At = | ||
| Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose(); |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]
Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose();
^| cartesian.row(row) = dirs[row]; | ||
| const Eigen::MatrixXd SH2At = | ||
| Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose(); | ||
| Eigen::VectorXd integrals = Eigen::VectorXd::Zero(SH2At.rows()); |
There was a problem hiding this comment.
warning: no header providing "Eigen::VectorXd" is directly included [misc-include-cleaner]
Eigen::VectorXd integrals = Eigen::VectorXd::Zero(SH2At.rows());
^| const Eigen::MatrixXd SH2At = | ||
| Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose(); | ||
| Eigen::VectorXd integrals = Eigen::VectorXd::Zero(SH2At.rows()); | ||
| integrals[0] = 2.0 * std::sqrt(Math::pi); |
There was a problem hiding this comment.
warning: no header providing "MR::Math::pi" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:19:
- #include "math/sphere.h"
+ #include "math/math.h"
+ #include "math/sphere.h"| const Eigen::MatrixXd SH2At = | ||
| Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose(); | ||
| Eigen::VectorXd integrals = Eigen::VectorXd::Zero(SH2At.rows()); | ||
| integrals[0] = 2.0 * std::sqrt(Math::pi); |
There was a problem hiding this comment.
warning: no header providing "std::sqrt" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:20:
+ #include <cmath>
| @@ -1,4 +1,4 @@ | |||
| /* Copyright (c) 2008-2026 the MRtrix3 contributors. | |||
| q/* Copyright (c) 2008-2026 the MRtrix3 contributors. | |||
There was a problem hiding this comment.
warning: unknown type name 'q' [clang-diagnostic-error]
q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
^| } | ||
|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); |
There was a problem hiding this comment.
warning: no header providing "size_t" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:0:
- q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
+ #include <cstddef>
+ q/* Copyright (c) 2008-2026 the MRtrix3 contributors.|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); | ||
| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); |
There was a problem hiding this comment.
warning: no header providing "Eigen::Dynamic" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:0:
- q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
+ #include <Eigen/src/Core/util/Constants.h>
+ q/* Copyright (c) 2008-2026 the MRtrix3 contributors.|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); | ||
| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); |
There was a problem hiding this comment.
warning: no header providing "Eigen::Matrix" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:0:
- q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
+ #include <Eigen/src/Core/Matrix.h>
+ q/* Copyright (c) 2008-2026 the MRtrix3 contributors.| const Eigen::MatrixXd SH2At = | ||
| Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose(); | ||
| Eigen::VectorXd integrals = Eigen::VectorXd::Zero(SH2At.rows()); | ||
| integrals[0] = 2.0 * std::sqrt(Math::pi); |
There was a problem hiding this comment.
warning: no header providing "MR::Math::pi" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:0:
- q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
+ #include "math/math.h"
+ q/* Copyright (c) 2008-2026 the MRtrix3 contributors.| const Eigen::MatrixXd SH2At = | ||
| Math::SH::init_transform(Math::Sphere::cartesian2spherical(cartesian), calibration_lmax).transpose(); | ||
| Eigen::VectorXd integrals = Eigen::VectorXd::Zero(SH2At.rows()); | ||
| integrals[0] = 2.0 * std::sqrt(Math::pi); |
There was a problem hiding this comment.
warning: no header providing "std::sqrt" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:0:
- q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
+ #include <cmath>
+ q/* Copyright (c) 2008-2026 the MRtrix3 contributors.
Lestropie
added a commit
that referenced
this pull request
Jul 11, 2026
|
|
||
| data = A.householderQr().solve(integral_results); | ||
| const size_t calibration_lmax = Math::SH::LforN((3 * dirs.size()) / 4); | ||
| Eigen::Matrix<default_type, Eigen::Dynamic, 3> cartesian(dirs.size(), 3); |
There was a problem hiding this comment.
warning: no header providing "Eigen::Dynamic" is directly included [misc-include-cleaner]
cpp/core/dwi/fmls.cpp:19:
- #include <cmath>
+ #include <Eigen/src/Core/util/Constants.h>
+ #include <cmath>GCC -O3 value-range analysis incorrectly flags Eigen::Array<bool> allocations in the null-lobe construction block as potential overflows, triggering -Walloc-size-larger-than= under -Werror. The warning is a false positive caused by inlining of the FOD_lobe(mask_type) constructor. A GCC-only diagnostic suppression pragma wraps the affected block; the __clang__ guard prevents Clang from rejecting the unknown option. Prompt: > Compilation in Linux GCC fails with compilation report posted below. Diagnose and report / fix. Generated-by: Claude Opus 4.8 <noreply@anthropic.com>
Lestropie
force-pushed
the
fod2fixel_quadrature_fix
branch
from
July 11, 2026 11:12
bfaf490 to
e78377c
Compare
Lestropie
added a commit
that referenced
this pull request
Aug 3, 2026
The FOD FMLS segmenter previously constructed a hardcoded 1281-direction set in every command that uses it, dominating the runtime of the commands that segment FODs. A new FMLSDirectionSet configuration file entry, together with an -fmls_directions command-line option on fod2fixel, now selects either a built-in set by integer count or an external file of directions; resolution follows the precedence of explicit option, then configuration entry, then the unchanged hard-coded default. All commands that segment FODs — afdconnectivity, fod2fixel, tckgen, tcksift, tcksift2, and the dynamic seeding path — obtain the set through a shared accessor rather than constructing it directly, and the reference documentation is regenerated accordingly. Control of the direction set outside fod2fixel is deferred exclusively to the configuration entry, so no other command gains a new option. The tests of those commands now request 321 directions to cut execution time, which required the regression reference data of only fod2fixel tests to be regenerated. Other slow tests were modified to use reduced data. Prompts: 1. > Build the software in debug mode with tests enabled. Run all tests with the "binary" label and no parallel execution. Sort tests by execution time, longest to shortest. Find a suitable bifurcation point in this distribution to select a group of tests for which modifications to reduce execution time is warranted. For each test, analyse the input data consumed by the test. Within the clone of the test data repository that resides within the build directory, propose a modification to the test data that would reduce the computational expense of the test. Some input test data may be consumed by multiple tests; where this occurs, the test itself can be modified to select a subset of the input data before executing the command under evaluation. Regenerate any reference output data for regression testing within the clone of the test data repository within the build directory. 2. > Recently (GitHub #3419) the quadrature for FOD integration was changed from householderQr() to completeOrthogonalDecomposition(), with a reduced lmax. This may be contributing to slow execution of [afdconnectivity, tcksift, tcksift2, tckgen_seed_dynamic]. Testing during the Claude session that produced that change proved that QR decomposition was the wrong choice in that context, as was the use of a large lmax. This worktree has now been rebased to include those changes. Investigate (from prior Claude session content & online) whether an alternative decomposition would yield comparable results to completeOrthogonalDecomposition() but be substantially faster. Speed test the rebased code against prior results in this session. If this does not have more than a 50% effect on test execution time, implement a new MRtrix config file entry that controls which of the built-in direction sets is utilised by default by the FOD FMLS segmenter, and update command documentation of those commands in this list offering an explicit command-line option to choose this set / provide an external text file containing a set to state that in the absence of that command-line option it falls back to that config file entry, or the hard-coded default if that is also absent. 3. > Only add the -fmls_directions option to the fod2fixel command; for other relevant commands, control over this direction set will be deferred exclusively to the config file entry. For relevant tests, use either this option where available or the -config option where not available to utilise the 513 direction set. Indicate which existing tests now fail and therefore require regeneration of reference data for regression testing. Generated-by: Claude Opus 5 <noreply@anthropic.com>
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.
While attempting to solve another problem it was discovered that a couple of my prior decisions around improving AFD integration (#699) were very sub-optimal: if directions were read as truncated strings from a header key-value rather than from the in-built double-precision tables, the resulting weights would be drastically different.
Mostly the instability washed out because the FOD features are so coarse relative to the sampling, but they can nevertheless be improved.
LforN(0.75 x N)where benefits in rotation invariance roll off and the condition number ramps up.Differences aren't really visible in
mrviewso not worth showing.devbranch ontest_datarepo upon merge.