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

Rename sh2power to sh2metric, new "entropy" operation #2985

Draft
wants to merge 1 commit into
base: dev
Choose a base branch
from

Conversation

Lestropie
Copy link
Member

@Lestropie Lestropie commented Sep 3, 2024

Had previously thought about this from the perspective of registration, but have also been looking at a collaborator dataset where doing satisfactory masking for FBA has been difficult, and using that to revised FBA documentation (#2685).

Here is the WM FOD l=0 image, which you can see doesn't do a very good job of separating what one would like to analyse vs. not:

Screenshot from 2024-09-03 16-13-43

Then this is the entropy of the WM FOD template (normalised & inverted):

Screenshot from 2024-09-03 16-14-08


  • Implement the Spherical Harmonic Entropy (SHE) metric as per:
    https://arxiv.org/pdf/1805.08084
    This may yield the same or similar result much faster

  • Add command-line option to modify the direction set that's used for amplitude sampling

  • Implement entropy operations for:

    • mrmath
    • fixel2voxel
  • Revisit templated entropy calculation code
    For some reason had trouble specifying the logarithm function as a template. The switch statement should get optimised out by the compiler, but maybe I'll have another go before merging.

Includes checking out a new test data commit that moves existing test data and adds new test data.
@Lestropie Lestropie self-assigned this Sep 3, 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 37. Check the log or trigger a new build to see more.


constexpr size_t default_direction_set = 1281;

std::vector<std::string> metrics = {"entropy", "power"};
Copy link

Choose a reason for hiding this comment

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

warning: variable 'metrics' is non-const and globally accessible, consider making it const [cppcoreguidelines-avoid-non-const-global-variables]

std::vector<std::string> metrics = {"entropy", "power"};
                         ^

}
// clang-format on

void run_entropy() {
Copy link

Choose a reason for hiding this comment

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

warning: function 'run_entropy' has cognitive complexity of 30 (threshold 25) [readability-function-cognitive-complexity]

void run_entropy() {
     ^
Additional context

cmd/sh2metric.cpp:83: +1, including nesting penalty of 0, nesting level increased to 1

  for (size_t i = 0; i != argument.size() - 2; ++i) {
  ^

cmd/sh2metric.cpp:87: +2, including nesting penalty of 1, nesting level increased to 2

    if (i == 0) {
    ^

cmd/sh2metric.cpp:89: +1, nesting level increased to 2

    } else {
      ^

cmd/sh2metric.cpp:90: +3, including nesting penalty of 2, nesting level increased to 3

      if (!voxel_grids_match_in_scanner_space(header, H_out))
      ^

cmd/sh2metric.cpp:105: nesting level increased to 1

    Processor(const std::vector<Image<float>> &SH_images,
    ^

cmd/sh2metric.cpp:112: +2, including nesting penalty of 1, nesting level increased to 2

      for (const auto &i : SH_images)
      ^

cmd/sh2metric.cpp:116: nesting level increased to 1

    Processor(const Processor &that)
    ^

cmd/sh2metric.cpp:118: +2, including nesting penalty of 1, nesting level increased to 2

      for (const auto &i : that.images)
      ^

cmd/sh2metric.cpp:122: nesting level increased to 1

    bool operator()(Iterator &pos) {
    ^

cmd/sh2metric.cpp:124: +2, including nesting penalty of 1, nesting level increased to 2

      for (auto &image : images) {
      ^

cmd/sh2metric.cpp:128: +3, including nesting penalty of 2, nesting level increased to 3

        for (auto l = Loop(3)(image); l; ++l)
        ^

cmd/sh2metric.cpp:136: +2, including nesting penalty of 1, nesting level increased to 2

      } catch (Exception &) {
        ^

cmd/sh2metric.cpp:154: nesting level increased to 1

      Shared(const std::vector<Image<float>> &SH_images, const DWI::Directions::Set &dirs, const bool normalise)
      ^

cmd/sh2metric.cpp:157: +2, including nesting penalty of 1, nesting level increased to 2

        for (ssize_t row = 0; row != dirs.size(); ++row)
        ^

cmd/sh2metric.cpp:160: +2, including nesting penalty of 1, nesting level increased to 2

        for (const auto &i : SH_images) {
        ^

cmd/sh2metric.cpp:162: +3, including nesting penalty of 2, nesting level increased to 3

          if (transforms.find(lmax) == transforms.end())
          ^

cmd/sh2metric.cpp:166: +2, including nesting penalty of 1, nesting level increased to 2

        if (normalise)
        ^

cmd/sh2metric.cpp:169: nesting level increased to 1

      vector_type operator()(vector_type &SH_coefs) const {
      ^

cmd/sh2metric.cpp:175: nesting level increased to 1

      size_t get_num_dirs() const { return num_dirs; }
      ^

cmd/sh2metric.cpp:176: nesting level increased to 1

      default_type normalise(const default_type in) const { return normalisation(in); }
      ^

cmd/sh2metric.cpp:184: nesting level increased to 1

        Normalisation()
        ^

cmd/sh2metric.cpp:187: nesting level increased to 1

        void initialise(const size_t num_images, const DWI::Directions::Set &dirs, const transform_type transform) {
        ^

cmd/sh2metric.cpp:200: nesting level increased to 1

        default_type operator()(const default_type in) const {
        ^

cmd/sh2metric.cpp:201: +2, including nesting penalty of 1, nesting level increased to 2

          if (!std::isfinite(lower) || !std::isfinite(upper))
          ^

cmd/sh2metric.cpp:201: +1

          if (!std::isfinite(lower) || !std::isfinite(upper))
                                    ^

cmd/sh2metric.cpp:182: nesting level increased to 1

      class Normalisation {
            ^

cmd/sh2metric.cpp:182: nesting level increased to 1

      class Normalisation {
            ^

cmd/sh2metric.cpp:182: nesting level increased to 1

      class Normalisation {
            ^

cmd/sh2metric.cpp:151: nesting level increased to 1

    class Shared {
          ^

cmd/sh2metric.cpp:151: nesting level increased to 1

    class Shared {
          ^

cmd/sh2metric.cpp:151: nesting level increased to 1

    class Shared {
          ^

cmd/sh2metric.cpp:151: nesting level increased to 1

    class Shared {
          ^

cmd/sh2metric.cpp:151: nesting level increased to 1

    class Shared {
          ^

cmd/sh2metric.cpp:102: nesting level increased to 1

  class Processor {
        ^

cmd/sh2metric.cpp:102: nesting level increased to 1

  class Processor {
        ^

for (size_t i = 0; i != argument.size() - 2; ++i) {
Header header(Header::open(argument[i]));
Math::SH::check(header);
max_lmax = std::max(max_lmax, Math::SH::LforN(header.size(3)));
Copy link

Choose a reason for hiding this comment

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

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

    max_lmax = std::max(max_lmax, Math::SH::LforN(header.size(3)));
                                                  ^

}
H_out.ndim() = 3;

DWI::Directions::Set dirs(default_direction_set);
Copy link

Choose a reason for hiding this comment

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

warning: variable 'dirs' of type 'DWI::Directions::Set' can be declared 'const' [misc-const-correctness]

Suggested change
DWI::Directions::Set dirs(default_direction_set);
DWI::Directions::Set const dirs(default_direction_set);

H_out.ndim() = 3;

DWI::Directions::Set dirs(default_direction_set);
Image<float> image_out(Image<float>::create(argument[argument.size() - 1], H_out));
Copy link

Choose a reason for hiding this comment

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

warning: variable 'image_out' of type 'Image' can be declared 'const' [misc-const-correctness]

Suggested change
Image<float> image_out(Image<float>::create(argument[argument.size() - 1], H_out));
Image<float> const image_out(Image<float>::create(argument[argument.size() - 1], H_out));

Normalisation()
: lower(std::numeric_limits<default_type>::quiet_NaN()),
upper(std::numeric_limits<default_type>::quiet_NaN()) {}
void initialise(const size_t num_images, const DWI::Directions::Set &dirs, const transform_type transform) {
Copy link

Choose a reason for hiding this comment

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

warning: the const qualified parameter 'transform' is copied for each invocation; consider making it a reference [performance-unnecessary-value-param]

Suggested change
void initialise(const size_t num_images, const DWI::Directions::Set &dirs, const transform_type transform) {
void initialise(const size_t num_images, const DWI::Directions::Set &dirs, const transform_type& transform) {

void initialise(const size_t num_images, const DWI::Directions::Set &dirs, const transform_type transform) {
Eigen::Matrix<default_type, Eigen::Dynamic, 1> delta_coefs;
Math::SH::delta(
delta_coefs, Eigen::Matrix<default_type, 3, 1>({0.0, 0.0, 1.0}), Math::SH::LforN(transform.cols()));
Copy link

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 'int' is implementation-defined [bugprone-narrowing-conversions]

              delta_coefs, Eigen::Matrix<default_type, 3, 1>({0.0, 0.0, 1.0}), Math::SH::LforN(transform.cols()));
                                                                               ^

void initialise(const size_t num_images, const DWI::Directions::Set &dirs, const transform_type transform) {
Eigen::Matrix<default_type, Eigen::Dynamic, 1> delta_coefs;
Math::SH::delta(
delta_coefs, Eigen::Matrix<default_type, 3, 1>({0.0, 0.0, 1.0}), Math::SH::LforN(transform.cols()));
Copy link

Choose a reason for hiding this comment

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

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

              delta_coefs, Eigen::Matrix<default_type, 3, 1>({0.0, 0.0, 1.0}), Math::SH::LforN(transform.cols()));
                                                                                               ^

// and pad out zeroes to simulate all other input SH functions being zero
Eigen::Matrix<default_type, Eigen::Dynamic, 1> amps = transform * delta_coefs;
amps.conservativeResizeLike(
Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero(num_images * transform.rows()));
Copy link

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]

              Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero(num_images * transform.rows()));
                                                                   ^

Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Zero(num_images * transform.rows()));
lower = Math::Entropy::nits(amps);
upper = Math::Entropy::nits(
Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Constant(num_images * transform.rows(), 1.0));
Copy link

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]

              Eigen::Matrix<default_type, Eigen::Dynamic, 1>::Constant(num_images * transform.rows(), 1.0));
                                                                       ^

@Lestropie
Copy link
Member Author

Was potentially a little too idealistic about this metric being useful as a scalar surrogate for multi-contrast registration (as multiple SH images can be provided as input to the entropy calculation). Kind of interesting, and maybe useful in specific use cases such as the unusual template above, but not a game-changer.

Single-subject HCP data:
MSMT_tissues
MSMT_entropy

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

Successfully merging this pull request may close these issues.

1 participant