Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merge SHARDRecon into main MRtrix3 codebase #2994

Open
wants to merge 434 commits into
base: dev
Choose a base branch
from
Open

Conversation

dchristiaens
Copy link
Member

@dchristiaens dchristiaens commented Sep 12, 2024

This PR aims to merge the SHARDRecon module into the main MRtrix3 repository. The recent changes to the build system and python interface make it challenging to maintain external projects that are not linked to one particular version of MRtrix3. Merging SHARDRecon into dev will ensure that the code remains compatible.

Key tasks:

  • merge code, preserving git history
  • fix compiler errors
  • update python scripts for new command interface
  • testing if output is correct
  • add documentation
  • add testing data

…niformly stored as yaw-pitch-roll Euler angles.
…tion.

The q-space basis is now evaluated per excitation, rather than per slice. This reduces the overall computational cost by effectively the multiband factor. Note that volume-to-volume reconstruction is treated as a special case where MB = Nslices, i.e., Nexcitations = 1.
This implementation uses the Lie-algebra se(3) of the rigid transformation group SE(3), coupled with a Levenberg-Marquardt solver for optimizing the asymmetric rigid registration problem.
github-actions[bot]

This comment was marked as outdated.

src/dwi/svr/mapping.h Outdated Show resolved Hide resolved
cmd/dwirecon.cpp Outdated
Comment on lines 179 to 181
// ----------------------- // TODO: Eddy uses a reverse LR axis for storing
PE.col(0) *= -1; // the PE table, akin to the gradient table.
// ----------------------- // Fix in the eddy import/export functions in core.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to determine whether this is an artifact of bad metadata handling in MRtrix3 (#2917), or whether it is in fact the case that the FSL text file formats for indicating the phase encoding table are expected to be based on a negative determinant image. My expectation was that it would be the former. If the comment is true and the issue is the latter, then there will need to be an explicit distinction between storage of the phase encoding table w.r.t. image axes and the import/export of that information from/to FSL.

Copy link
Member

@jdtournier jdtournier Sep 13, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've removed all of these on the version we're using at KCL now. You can see these changes on my fork: dchristiaens/shard-recon@master...jdtournier:shard-recon:master

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jdtournier Good point. I forgot to merge your PR before creating this one. Your last update is that you were going to confirm if the PE direction had been fixed in the main repo. Is this the case? If so, I don't understand why the test I did last week didn't flag a change in PE.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, all I can say is that the code we're currently running has these lines removed - but I did also have to modify the corresponding acqp.txt for it to work correctly... So that would suggest that it isn't fixed after all...

It terms of testing, I guess what's required is to first run eddy to verify the correctness of the acqp.txt file, and then run shardrecon with same acqp.txt, right? Must admit I haven't done that, and given the above, I suspect that will show that the problem persists...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  1. Is the derivation of your acqp.txt file going through MRtrix3 processing at all? I'm presuming it's not being loaded from DICOM headers, but are you loading it into an image header and then doing an export to text file? If so then the issues rectified in DWI metadata handling fixes #2917 could be playing a role.

  2. A useful test to do, with this line removed, would be to find the acqp.txt that results in topup doing the right thing, then convert the NIfTI between LAS and RAS and re-run with the same acqp.txt. If it then produces the wrong result, then the FSL interpretation of the phase encoding table is with respect to the image axes, as has always been my understanding. If it produces the correct result, then FSL expects the phase encoding table to be defined with respect to its internal left-handed coordinate system, and our command-line interface for phase encoding data handling has to be changed.

github-actions[bot]

This comment was marked as outdated.

github-actions[bot]

This comment was marked as outdated.

github-actions[bot]

This comment was marked as outdated.

pe_import_option = ['-import_pe_table', app.ARGS.pe_table]
elif app.ARGS.pe_eddy:
pe_import_option = ['-import_pe_eddy', app.ARGS.pe_eddy[0], app.ARGS.pe_eddy[1]]
# check output path - code no longer works; now done internally??
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct: arguments of type (Directory|File|Image)Out are now automatically checked, so devs are no longer responsible for adding these function calls.

# Now that FORCE_OVERWRITE has been set,
# check any user-specified output paths
try:
for key in vars(ARGS):
value = getattr(ARGS, key)
if isinstance(value, Parser._UserOutPathExtras): # pylint: disable=protected-access
value.check_output()
except FileExistsError as exception:
sys.stderr.write('\n')
sys.stderr.write(f'{EXEC_NAME}: {ANSI.error}[ERROR] {exception}{ANSI.clear}\n')
sys.stderr.flush()
sys.exit(1)

if not image.match(header, 'mask.mif', up_to_dim=3):
raise MRtrixError('Provided mask image does not match input DWI')
else:
run.command('dwi2mask in.mif mask.mif')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Needs to be updated to new dwi2mask interface


# __________ Copy outputs __________

run.command(['mrconvert', f'recon-{nepochs}.mif', app.ARGS.output], mrconvert_keyval=f'recon-{nepochs}.mif', force=app.FORCE_OVERWRITE)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
run.command(['mrconvert', f'recon-{nepochs}.mif', app.ARGS.output], mrconvert_keyval=f'recon-{nepochs}.mif', force=app.FORCE_OVERWRITE)
run.command(['mrconvert', f'recon-{nepochs}.mif', app.ARGS.output],
mrconvert_keyval=app.ARGS.input,
force=app.FORCE_OVERWRITE)

Purpose here is to have command_history in the output image reflect the invocation of dwimotioncorrect, rather than listing all of its internal steps.

Comment on lines +327 to +330
if app.ARGS.export_motion:
run.command(f'cp motion.txt {app.ARGS.export_motion}')
if app.ARGS.export_weights:
run.command(f'cp sliceweights.txt {app.ARGS.export_weights}')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In addition to using Python functions rather than Unix cp, this could potentially be improved by adding command_history to the header comments of these files.

@Lestropie
Copy link
Member

I'm going to generate a list here of things that come to mind; nothing major, just thoughts that I think are worth contemplating. Don't want to block progress by demanding massive pieces of work, just starting dialogues. Some I can help out with.

  • dwirecon command (New command: dwirecon #2915 relevant):
    At the moment this is hard-coded to yield multi-shell spherical harmonics data, however from the command name one would intuitively expect DWI data. There's a few different options here, and it might need a combination of things to get fully resolved with New command: dwirecon #2915:

    • Change name of this command to msshrecon. Completely unblocks New command: dwirecon #2915.
    • Make dwirecon an operation-based command. Essentially take what's in New command: dwirecon #2915 and make what is contained within cmd/dwirecon.cpp in this PR another operation within that command, with a set of command-line options that are applicable only to that operation.
    • Have a toggle that changes the output of this code between multi-shell spherical harmonics and reconstructed DWI data.
      (I would possibly lean towards a toggle here, rather than having a default output and then a command-line option that outputs the data in the other form, since it's large data and users are likely to want one or the other)
  • Data type for command naming:
    Currently here "mssh" is used for multi-shell spherical harmonics. I am however reminded that dwi2fod msmt_csd started its life as msdwi2fod. Here there is a stronger argument for a different datatype indicator: mssh and sh are less interchangeable from a data description perspective that single-shell vs. multi-shell DWI data. It may however be the case that many operations intended for SH data can be extended to multi-shell SH data.

    • Possibility of just using sh from a command name perspective, and anything that expects multi-shell SH but receives a single SH, or vice versa, throws an Exception? Whether this might be reasonable would require looking across the package for SH interfaces.
    • Is mssh the best indicator? Would "multi-spherical-harmonics (msh)" work? Are there any contexts in which one may have multiple spherical harmonics functions, where the groups of coefficients do not correspond to shells? Multiple ODFs within a single 5D image comes to mind. I suppose shard isn't correct from a data description perspective?
  • dwirecon needs a more comprehensive Description.

  • I'd like for there to be more explanation of how the slice profile is derived / interfaced so that I can try to make proper use of it for my own data.

  • I know that for my own data, having the ability to correct for eddy current distortions is going to be important. I tried to read up on the Lie algebra used for the motion modeling, but it's well outside my current expertise to know how it might be extended from rigid-body to an affine model, or what the consequences of doing so on other code might be. @dchristiaens Do you have an intuition of how involved such a change might be? I'm happy to help on this one where I can from a code perspective , I'm moreso doubtful on my capacity to derive / follow the math.

  • Some commands involve specification of multi-band factors to influence slice-to-volume registration. This could be streamlined by instead making use of SliceTiming.

  • With mrfieldunwarp, I'm wondering if this could be handled differently. If one has a susceptibility field in Hz, then that can be converted to a deformation field based on phase encoding direction & readout time. That's just a non-linear warp field that can be composed with any other transformation, whether it be rigid-body motion (as `mrfieldunwarp currently handles on a per-slice-stack basis) or anything else (eg. I'm imagining composing this with the gradient non-linearity geometric distortion field). Therefore, is this something that could be adequately handled by:

    • Giving warpconvert the ability to convert inhomogeneity field to deformation field based on phase encoding information
    • Having some other way to compose this deformation with per-slice-stack motion parameters
      ? It's not immediately clicking what the best form of the second point would be. Explicitly composing transforms and then subsequently applying those transforms would be the most modular, but highly inefficient. Maybe if mrtransform could accept both the deformation field and the rigid-body motion parameters and would then be responsible for the composition?
  • @dchristiaens Are you able to share some exemplar data for which you know the reconstruction pipeline behaves properly and yields an acceptable result, whether privately or publicly? Just if I'm going to muck around with the code I'd like to know that I'm not breaking anything.
    (If I do come up with any changes, I won't commit directly to this branch, I'll make a PR)

  • Some Python commands need to have their typed arguments updated.

  • Python commands using 4 spaces for indentation, whereas the internal convention is 2.

  • Some Python code could make better use of F-strings, which are now used throughout the Python code base.

  • File copyright headers currently differ from the main code base.

@MRtrix3 MRtrix3 deleted a comment from github-actions bot Sep 18, 2024
Copy link

@github-actions github-actions bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 342. Check the log or trigger a new build to see more.

cmd/dwirecon.cpp Show resolved Hide resolved
cmd/dwirecon.cpp Outdated Show resolved Hide resolved
cmd/dwirecon.cpp Outdated Show resolved Hide resolved
cmd/dwirecon.cpp Outdated Show resolved Hide resolved
cmd/dwisliceoutliergmm.cpp Outdated Show resolved Hide resolved
cmd/dwisliceoutliergmm.cpp Outdated Show resolved Hide resolved
cmd/dwisliceoutliergmm.cpp Outdated Show resolved Hide resolved
cmd/dwisliceoutliergmm.cpp Show resolved Hide resolved
cmd/dwisliceoutliergmm.cpp Show resolved Hide resolved
cmd/dwisliceoutliergmm.cpp Outdated Show resolved Hide resolved
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Copy link

@github-actions github-actions bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 321. Check the log or trigger a new build to see more.

y.setZero();
ssize_t j = 0;
for (auto lv = Loop("loading image data", {0, 1, 2, 3})(dwisub); lv; lv++, j++) {
float const w = Wsub(dwisub.index(2), dwisub.index(3)) * Wvox[j];

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: no viable conversion from 'const CwiseBinaryOp<internal::scalar_product_op<typename internal::traits<IndexedView<Matrix<float, -1, -1, 0, -1, -1>, Index<Extract1D<Image>>, Index<Extract1D<Image>>>>::Scalar, typename internal::promote_scalar_arg<Scalar, float, (Eigen::internal::has_ReturnType<Eigen::ScalarBinaryOpTraits<Scalar, float, Eigen::internal::scalar_product_op<Scalar, float>>>::value)>::type>, const IndexedView<Matrix<float, -1, -1, 0, -1, -1>, Index<Extract1D<Image>>, Index<Extract1D<Image>>>, const typename internal::plain_constant_type<IndexedView<Matrix<float, -1, -1, 0, -1, -1>, Index<Extract1D<Image>>, Index<Extract1D<Image>>>, typename internal::promote_scalar_arg<Scalar, float, (Eigen::internal::has_ReturnType<Eigen::ScalarBinaryOpTraits<Scalar, float, Eigen::internal::scalar_product_op<Scalar, float>>>::value)>::type>::type>' (aka 'const CwiseBinaryOp<scalar_product_op<float, float>, const Eigen::IndexedView<Eigen::Matrix<float, -1, -1, 0, -1, -1>, MR::Helper::Index<MR::Adapter::Extract1D<MR::Image>>, MR::Helper::Index<MR::Adapter::Extract1D<MR::Image>>>, const CwiseNullaryOp<scalar_constant_op, const Eigen::Matrix<float, -1, -1, 0, -1, -1>>>') to 'const float' [clang-diagnostic-error]

    float const w = Wsub(dwisub.index(2), dwisub.index(3)) * Wvox[j];
                ^

}

private:
const ssize_t nv, nz, ne;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'nv' of type 'const ssize_t' (aka 'const long') is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

  const ssize_t nv, nz, ne;
                ^

}

private:
const ssize_t nv, nz, ne;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'nz' of type 'const ssize_t' (aka 'const long') is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

  const ssize_t nv, nz, ne;
                    ^

}

private:
const ssize_t nv, nz, ne;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'ne' of type 'const ssize_t' (aka 'const long') is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

  const ssize_t nv, nz, ne;
                        ^

VecType posterior() const { return Rin.array().exp(); }

private:
const int niter;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: member 'niter' of type 'const int' is const qualified [cppcoreguidelines-avoid-const-or-ref-data-members]

  const int niter;
            ^

// Log-residuals
int k = 0;
Eigen::VectorXf res(E.rows() * shells[s].count());
for (size_t v : shells[s].get_volumes())

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: variable 'v' of type 'size_t' (aka 'unsigned long') can be declared 'const' [misc-const-correctness]

Suggested change
for (size_t v : shells[s].get_volumes())
for (size_t const v : shells[s].get_volumes())

int k = 0;
Eigen::VectorXf res(E.rows() * shells[s].count());
for (size_t v : shells[s].get_volumes())
res.segment(E.rows() * (k++), E.rows()) = E.col(v);

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]

      res.segment(E.rows() * (k++), E.rows()) = E.col(v);
                                                      ^

for (size_t v : shells[s].get_volumes())
res.segment(E.rows() * (k++), E.rows()) = E.col(v);
// clip at non-zero minimum
value_type nzmin = res.redux([](value_type a, value_type b) {

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: variable 'nzmin' of type 'value_type' (aka 'float') can be declared 'const' [misc-const-correctness]

Suggested change
value_type nzmin = res.redux([](value_type a, value_type b) {
value_type const nzmin = res.redux([](value_type a, value_type b) {

Comment on lines +278 to +279
else
return (b > 0) ? b : std::numeric_limits<value_type>::infinity();

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: do not use 'else' after 'return' [readability-else-after-return]

Suggested change
else
return (b > 0) ? b : std::numeric_limits<value_type>::infinity();
return (b > 0) ? b : std::numeric_limits<value_type>::infinity();

else
return (b > 0) ? b : std::numeric_limits<value_type>::infinity();
});
Eigen::VectorXf logres = res.array().max(nzmin).log();

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

warning: variable 'logres' of type 'Eigen::VectorXf' (aka 'Matrix<float, Dynamic, 1>') can be declared 'const' [misc-const-correctness]

Suggested change
Eigen::VectorXf logres = res.array().max(nzmin).log();
Eigen::VectorXf const logres = res.array().max(nzmin).log();

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants