tckmap: eigen assertion fails on macOS

[bjeurissen]

CI for tckmap causes a failed assertion on macOS(other platforms appear to be fine):

➜  data git:(26b8e57) lldb ../../../bin/tckmap -- tracks.tck -vox 1 -
(lldb) target create "../../../bin/tckmap"
Current executable set to '/Users/ben/mrtrix3/bin/tckmap' (arm64).
(lldb) settings set -- target.run-args  "tracks.tck" "-vox" "1" "-"
(lldb) run
Process 97760 launched: '/Users/ben/mrtrix3/bin/tckmap' (arm64)
tckmap: [done] creating new template image
Assertion failed: ((!(RowsAtCompileTime!=Dynamic) || (rows==RowsAtCompileTime)) && (!(ColsAtCompileTime!=Dynamic) || (cols==ColsAtCompileTime)) && (!(RowsAtCompileTime==Dynamic && MaxRowsAtCompileTime!=Dynamic) || (rows<=MaxRowsAtCompileTime)) && (!(ColsAtCompileTime==Dynamic && MaxColsAtCompileTime!=Dynamic) || (cols<=MaxColsAtCompileTime)) && rows>=0 && cols>=0 && "Invalid sizes when resizing a matrix or array."), function resize, file PlainObjectBase.h, line 277.
Process 97760 stopped
* thread #1, queue = 'com.apple.main-thread', stop reason = hit program assert
    frame #4: 0x00000001000d1dfc tckmap`Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >::resize(this=0x0000000101505500, rows=-2, cols=4) at PlainObjectBase.h:273:7
   270      EIGEN_DEVICE_FUNC
   271      EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
   272      {
-> 273        eigen_assert(   EIGEN_IMPLIES(RowsAtCompileTime!=Dynamic,rows==RowsAtCompileTime)
   274                     && EIGEN_IMPLIES(ColsAtCompileTime!=Dynamic,cols==ColsAtCompileTime)
   275                     && EIGEN_IMPLIES(RowsAtCompileTime==Dynamic && MaxRowsAtCompileTime!=Dynamic,rows<=MaxRowsAtCompileTime)
   276                     && EIGEN_IMPLIES(ColsAtCompileTime==Dynamic && MaxColsAtCompileTime!=Dynamic,cols<=MaxColsAtCompileTime)
Target 0: (tckmap) stopped.
(lldb) thread backtrace
* thread #1, queue = 'com.apple.main-thread', stop reason = hit program assert
    frame #0: 0x00000001953be1b0 libsystem_kernel.dylib`__pthread_kill + 8
    frame #1: 0x00000001953f4cec libsystem_pthread.dylib`pthread_kill + 288
    frame #2: 0x000000019532e2c8 libsystem_c.dylib`abort + 180
    frame #3: 0x000000019532d620 libsystem_c.dylib`__assert_rtn + 272
  * frame #4: 0x00000001000d1dfc tckmap`Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >::resize(this=0x0000000101505500, rows=-2, cols=4) at PlainObjectBase.h:273:7
    frame #5: 0x0000000100289638 tckmap`MR::DWI::Tractography::Resampling::Upsampler::set_ratio(this=0x00000001015054f8, upsample_ratio=18446744073709551615) at upsampler.cpp:61:15
    frame #6: 0x000000010000b9fc tckmap`MR::DWI::Tractography::Mapping::TrackMapperBase::set_upsample_ratio(this=0x0000000101505390, i=18446744073709551615) at mapper.h:85:66
    frame #7: 0x00000001000067b8 tckmap`run() at tckmap.cpp:562:11
    frame #8: 0x0000000100002684 tckmap`main(cmdline_argc=5, cmdline_argv=0x000000016fdff410) at command.h:101:5
    frame #9: 0x00000001950cbe50 dyld`start + 2544

As raised in #2472 and #2566

More detail:

* thread #1, queue = 'com.apple.main-thread', stop reason = hit program assert
    frame #4: 0x000000010004a4f4 tckmap`Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >::resize(this=0x00000001015053b0, rows=-2, cols=4) at PlainObjectBase.h:281:7
   278      EIGEN_DEVICE_FUNC
   279      EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
   280      {
-> 281        eigen_assert(   EIGEN_IMPLIES(RowsAtCompileTime!=Dynamic,rows==RowsAtCompileTime)
   282                     && EIGEN_IMPLIES(ColsAtCompileTime!=Dynamic,cols==ColsAtCompileTime)
   283                     && EIGEN_IMPLIES(RowsAtCompileTime==Dynamic && MaxRowsAtCompileTime!=Dynamic,rows<=MaxRowsAtCompileTime)
   284                     && EIGEN_IMPLIES(ColsAtCompileTime==Dynamic && MaxColsAtCompileTime!=Dynamic,cols<=MaxColsAtCompileTime)
Target 0: (tckmap) stopped.
(lldb) frame variable
(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > *) this = 0x00000001015053b0
(Eigen::Index) rows = -2
(Eigen::Index) cols = 4

(lldb) up
frame #5: 0x000000010025e240 tckmap`MR::DWI::Tractography::Resampling::Upsampler::set_ratio(this=0x00000001015053a8, upsample_ratio=18446744073709551615) at upsampler.cpp:61:15
   58             if (upsample_ratio > 1) {
   59               const size_t dim = upsample_ratio - 1;
   60               Math::Hermite<value_type> interp (hermite_tension);
-> 61               M.resize (dim, 4);
   62               for (size_t i = 0; i != dim; ++i) {
   63                 interp.set ((i+1.0) / value_type(upsample_ratio));
   64                 for (size_t j = 0; j != 4; ++j)
(lldb) frame variable
(MR::DWI::Tractography::Resampling::Upsampler *) this = 0x00000001015053a8
(const size_t) upsample_ratio = 18446744073709551615
(const size_t) dim = 18446744073709551614
(MR::Math::Hermite<float>) interp = (w = (1.40129846E-45, 2.80259693E-45, 0, 3.82625779E-38), t = 0.0500000007)

(lldb) up
frame #6: 0x00000001000b8ca4 tckmap`MR::DWI::Tractography::Mapping::TrackMapperBase::set_upsample_ratio(this=0x0000000101505240, i=18446744073709551615) at mapper.h:85:66
   82               virtual ~TrackMapperBase() { }
   83
   84
-> 85               void set_upsample_ratio (const size_t i) { upsampler.set_ratio (i); }
   86               void set_map_zero (const bool i) { map_zero = i; }
   87               void set_use_precise_mapping (const bool i) {
   88                 if (i && ends_only)
(lldb) frame variable
(MR::DWI::Tractography::Mapping::TrackMapperTWI *) this = 0x0000000101505240
(const size_t) i = 18446744073709551615

(lldb) up
frame #7: 0x00000001000b5274 tckmap`run() at tckmap.cpp:562:11
   559    TrackLoader loader (file, num_tracks);
   560
   561    std::unique_ptr<TrackMapperTWI> mapper ((stat_tck == GAUSSIAN) ? (new Gaussian::TrackMapper (header, contrast)) : (new TrackMapperTWI (header, contrast, stat_tck)));
-> 562    mapper->set_upsample_ratio      (upsample_ratio);
   563    mapper->set_map_zero            (map_zero);
   564    mapper->set_use_precise_mapping (precise);
   565    mapper->set_map_ends_only       (ends_only);
(lldb) frame variable
(MR::DWI::Tractography::Properties) properties = {
  MR::KeyValues = size=25 {
    [0] = (first = "count", second = "500")
    [1] = (first = "downsample_factor", second = "3")
    [2] = (first = "fod_power", second = "0.25")
    [3] = (first = "init_threshold", second = "0.1")
    [4] = (first = "lmax", second = "8")
    [5] = (first = "max_angle", second = "45")
    [6] = (first = "max_dist", second = "250")
    [7] = (first = "max_num_attempts", second = "50000")
    [8] = (first = "max_num_tracks", second = "500")
    [9] = (first = "max_seed_attempts", second = "1")
    [10] = (first = "max_trials", second = "1000")
    [11] = (first = "method", second = "iFOD2")
    [12] = (first = "min_dist", second = "4")
    [13] = (first = "mrtrix_version", second = "0.3.12-325-gc203eda9")
    [14] = (first = "output_step_size", second = "1.25")
    [15] = (first = "rk4", second = "0")
    [16] = (first = "samples_per_step", second = "4")
    [17] = (first = "sh_precomputed", second = "1")
    [18] = (first = "source", second = "dwi2fod/out.mif")
    [19] = (first = "step_size", second = "1.25")
    [20] = (first = "stop_on_all_include", second = "0")
    [21] = (first = "threshold", second = "0.1")
    [22] = (first = "timestamp", second = "1443190374.822715044")
    [23] = (first = "total_count", second = "1748")
    [24] = (first = "unidirectional", second = "0")
  }
  include = {
    MR::DWI::Tractography::ROISetBase = {
      R = {
        std::__1::vector<MR::DWI::Tractography::ROI, std::__1::allocator<MR::DWI::Tractography::ROI> > = size=0 {}
      }
    }
  }
  exclude = {
    MR::DWI::Tractography::ROISetBase = {
      R = {
        std::__1::vector<MR::DWI::Tractography::ROI, std::__1::allocator<MR::DWI::Tractography::ROI> > = size=0 {}
      }
    }
  }
  mask = {
    MR::DWI::Tractography::ROISetBase = {
      R = {
        std::__1::vector<MR::DWI::Tractography::ROI, std::__1::allocator<MR::DWI::Tractography::ROI> > = size=0 {}
      }
    }
  }
  ordered_include = {
    MR::DWI::Tractography::ROISetBase = {
      R = {
        std::__1::vector<MR::DWI::Tractography::ROI, std::__1::allocator<MR::DWI::Tractography::ROI> > = size=0 {}
      }
    }
  }
  seeds = {
    seeders = {
      std::__1::vector<std::__1::unique_ptr<MR::DWI::Tractography::Seeding::Base, std::__1::default_delete<MR::DWI::Tractography::Seeding::Base> >, std::__1::allocator<std::__1::unique_ptr<MR::DWI::Tractography::Seeding::Base, std::__1::default_delete<MR::DWI::Tractography::Seeding::Base> > > > = size=0 {}
    }
    total_volume = 0
    total_count = 0
  }
  prior_rois = size=3 {
    [0] = (first = "mask", second = "mask.mif")
    [1] = (first = "mask", second = "mask.mif")
    [2] = (first = "seed", second = "mask.mif")
  }
  comments = {
    std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > = size=0 {}
  }
}
(MR::DWI::Tractography::Reader<float>) file = {
  MR::DWI::Tractography::__ReaderBase__ = {
    in = {
      std::__1::basic_istream<char, std::__1::char_traits<char> > = {
        std::__1::basic_ios<char, std::__1::char_traits<char> > = {
          std::__1::ios_base = {
            __fmtflags_ = 4098
            __precision_ = 6
            __width_ = 0
            __rdstate_ = 0
            __exceptions_ = 0
            __rdbuf_ = 0x000000016fdfee78{...}
            __loc_ = 0x00000001f0765800{...}
            __fn_ = 0x0000000000000000{...}
            __index_ = 0x0000000000000000{...}
            __event_size_ = 0
            __event_cap_ = 0
            __iarray_ = 0x0000000000000000{...}
            __iarray_size_ = 0
            __iarray_cap_ = 0
            __parray_ = 0x0000000000000000{...}
            __parray_size_ = 0
            __parray_cap_ = 0
          }
          __tie_ = nullptr
          __fill_ = -1
        }
        __gc_ = 0
      }
      __sb_ = {
        std::__1::basic_streambuf<char, std::__1::char_traits<char> > = {
          __loc_ = (__locale_ = 0x00000001f0765800)
          __binp_ = 0x0000000000000000
          __ninp_ = 0x0000000000000000
          __einp_ = 0x0000000000000000
          __bout_ = 0x0000000000000000
          __nout_ = 0x0000000000000000
          __eout_ = 0x0000000000000000
        }
        __extbuf_ = 0x0000000101808e00 ": 0\nc: 1\nMRViewDefaultTractGeomType: Lines\nMRViewDockFloating: 0\nMRViewWrapVolumes: 1\nMRViewShowComments: 0\n#BZeroThreshold: 110\nMRViewOrthoAsRow: 1\n#RealignTransform: 0\n"
        __extbufnext_ = 0x0000000000000000
        __extbufend_ = 0x0000000000000000
        __extbuf_min_ = "\xf0\xee\xdfo\U00000001\0\0"
        __ebs_ = 4096
        __intbuf_ = 0x0000000000000000
        __ibs_ = 0
        __file_ = 0x00000001f07642b8
        __cv_ = 0x00000001f0765618
        __st_ = (__mbstate8 = "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0", _mbstateL = 0)
        __st_last_ = (__mbstate8 = "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0", _mbstateL = 0)
        __om_ = 12
        __cm_ = 0
        __owns_eb_ = true
        __owns_ib_ = false
        __always_noconv_ = true
      }
    }
    dtype = (dt = 'F')
    current_index = 0
  }
  weights = {
    Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > = {
      m_storage = {
        m_data = 0x0000000000000000
        m_rows = 0
      }
    }
  }
}
(const size_t) num_tracks = 500
(MR::vector<double, 0>) voxel_size = {
  std::__1::vector<double, std::__1::allocator<double> > = size=3 {
    [0] = 1.0358282632233692E-309
    [1] = 1.0358282632233692E-309
    [2] = 1.0358282632233692E-309
  }
}
(MR::Header) header = {
  axes_ = {
    std::__1::vector<MR::Header::Axis, std::__1::allocator<MR::Header::Axis> > = size=3 {
      [0] = (size = 9223372036854775807, spacing = 1.0358282632233692E-309, stride = 1)
      [1] = (size = 9223372036854775807, spacing = 1.0358282632233692E-309, stride = 2)
      [2] = (size = 9223372036854775807, spacing = 1.0358282632233692E-309, stride = 3)
    }
  }
  transform_ = {
    m_matrix = {
      Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 4, 0, 3, 4> > = {
        m_storage = {
          m_data = {
            array ={...}
          }
        }
      }
    }
  }
  name_ = "tckmap image header"
  keyval_ = size=7 {
    [0] = (first = "comments", second = "track-weighted image")
    [1] = (first = "precise_mapping", second = "0")
    [2] = (first = "tck_source", second = "tracks.tck")
    [3] = (first = "twi_contrast", second = "tdi")
    [4] = (first = "twi_dimensionality", second = "greyscale")
    [5] = (first = "twi_tck_stat", second = "mean")
    [6] = (first = "twi_vox_stat", second = "sum")
  }
  format_ = 0x0000000000000000
  io = nullptr {
    __value_ = nullptr
  }
  datatype_ = (dt = 'D')
  offset_ = 0
  scale_ = 1
  realign_perm_ = {
    __elems_ = ([0] = 0, [1] = 1, [2] = 2)
  }
  realign_flip_ = {
    __elems_ = ([0] = false, [1] = false, [2] = false)
  }
}
(MR::vector<MR::App::ParsedOption, 0>) opt = {
  std::__1::vector<MR::App::ParsedOption, std::__1::allocator<MR::App::ParsedOption> > = size=0 {}
}
(const MR::DWI::Tractography::Mapping::contrast_t) contrast = TDI
(MR::DWI::Tractography::Mapping::vox_stat_t) stat_vox = V_SUM
(MR::DWI::Tractography::Mapping::tck_stat_t) stat_tck = T_MEAN
(float) gaussian_fwhm_tck = 0
(bool) backtrack = false
(MR::DWI::Tractography::Mapping::writer_dim) writer_type = GREYSCALE
(std::unique_ptr<MR::DWI::Directions::FastLookupSet, std::default_delete<MR::DWI::Directions::FastLookupSet> >) dirs = nullptr {
  __value_ = nullptr
}
(const bool) precise = false
(const bool) ends_only = false
(size_t) upsample_ratio = 18446744073709551615
(const bool) have_weights = false
(MR::DataType) default_datatype = (dt = '\x04')
(const bool) map_zero = false
(std::string) msg = "Generating greyscale image with density contrast and summed per-voxel statistic"
(MR::DWI::Tractography::Mapping::TrackLoader) loader = {
  reader = 0x000000016fdfee60
  tracks_to_load = 500
  progress = MR::ProgressBar @ 0x0000600002c00200 {
    __value_ = 0x0000600002c00200
  }
}
(std::unique_ptr<MR::DWI::Tractography::Mapping::TrackMapperTWI, std::default_delete<MR::DWI::Tractography::Mapping::TrackMapperTWI> >) mapper = MR::DWI::Tractography::Mapping::TrackMapperTWI @ 0x0000000101505240 {
  __value_ = 0x0000000101505240
}
(std::unique_ptr<MR::DWI::Tractography::Mapping::MapWriterBase, std::default_delete<MR::DWI::Tractography::Mapping::MapWriterBase> >) writer = MR::DWI::Tractography::Mapping::MapWriterBase @ 0x0000000000cbc000 {
  __value_ = 0x0000000000cbc000
}
*** Some of the displayed variables have a greater depth of members than the debugger will show by default. To increase the limit, use the --depth option to frame variable, or raise the limit by changing the target.max-children-depth setting.

Is this because macos CI uses a different eigen version than linux and windows? --> same problem with eigen 3.3 Or is this due to an update of macOS/xcode?

[bjeurissen]

From the above stack trace, the culprit appears to be (size_t) upsample_ratio = 18446744073709551615, which corresponds to the largest uint64. This leads to M.resize (18446744073709551614, 4). Clearly, something is going wrong.

[bjeurissen]

Found the culprit: https://github.com/MRtrix3/mrtrix3/blob/f633dfd7e9080f71877ea6a4619dabdde99a0fb6/cmd/tckmap.cpp#L275

For some reason this does not work as expected in combination with recent xcode... Instead of a vector of ones, this results in a vector of extremely small values.

[jdtournier]

OK, was looking into it, and really can't replicate any issues on Linux...

Can you try using:

voxel_size.resize (3, voxel_size.front()); 

instead of that line?

[bjeurissen]

indeed, I just tested resize and that works! Will push a fix

[jdtournier]

Hang on, I'm not sure it's a safe fix!

[bjeurissen]

Too late ;)

[jdtournier]

I think the issue is that std::vector::front() (and std::vector::operator[]) return a reference to the first element, not a copy. The first thing that the assign() (or resize()!) call might do is a memory re-allocation, which presumably involves deleting the original memory location, and so invalidate the data referenced in the earlier .front() call. The fact that it works is not an indication that it's safe, I think this is still technically undefined behaviour, whether we use .assign() or .resize().

The proper fix here would be something like:

274   if (voxel_size.size() == 1) {                                                                                           
275     auto v = voxel_size.front();                                                                                          
276     voxel_size.assign (3, v);                                                                                             
277   }   

[jdtournier]

... or potentially:

 voxel_size.assign (3, default_type (voxel_size.front())); 

since that will hopefully force a non-reference temporary value as the second argument...

[bjeurissen]

Thanks, I pushed the long option just to be sure it will work as expected