17#ifndef __image_helpers_h__
18#define __image_helpers_h__
31 template <
class AxesType>
32 FORCE_INLINE auto __ndim (
const AxesType& axes) ->
decltype (axes.size(), size_t()) {
return axes.size(); }
34 template <
class AxesType>
35 FORCE_INLINE auto __ndim (
const AxesType& axes) ->
decltype (axes.ndim(), size_t()) {
return axes.ndim(); }
37 template <
class AxesType>
38 FORCE_INLINE auto __get_index (
const AxesType& axes,
size_t axis) ->
decltype (axes.size(), ssize_t())
39 {
return axes[
axis]; }
41 template <
class AxesType>
42 FORCE_INLINE auto __get_index (
const AxesType& axes,
size_t axis) ->
decltype (axes.ndim(), ssize_t())
43 {
return axes.index(
axis); }
45 template <
class AxesType>
46 FORCE_INLINE auto __set_index (AxesType& axes,
size_t axis, ssize_t
index) ->
decltype (axes.size(), void())
49 template <
class AxesType>
50 FORCE_INLINE auto __set_index (AxesType& axes,
size_t axis, ssize_t
index) ->
decltype (axes.ndim(), void())
54 template <
class... DestImageType>
59 template <
class ImageType>
63 template <
class... DestImageType>
68 template <
class ImageType>
72 template <
class... DestImageType>
76 template <
class ImageType>
80 template <
class... DestImageType>
84 template <
class ImageType>
88 template <
class ImageType>
90 template <
class... DestImageType>
92 size_t last_axis = to_axis;
93 apply (__max_axis<DestImageType...> (last_axis), std::tie (ref, dest...));
94 for (
size_t n = from_axis; n < last_axis; ++n)
95 apply (__assign<DestImageType...> (n, __get_index (ref, n)), std::tie (dest...));
98 const size_t from_axis, to_axis;
102 template <
class ImageType,
typename IntType>
104 template <
class... DestImageType>
107 apply (__assign<DestImageType...> (a, __get_index (ref, a)), std::tie (dest...));
109 const ImageType& ref;
110 const vector<IntType> axes;
114 template <
typename ValueType>
class Image;
119 template <
class HeaderType,
typename ReturnType>
121 typedef decltype ((void) (
122 std::declval<HeaderType>().ndim() +
123 std::declval<HeaderType>().size(0) +
124 std::declval<HeaderType>().name().size()
125 ), std::declval<ReturnType>()) type;
129 template<
typename HeaderType>
131 typedef char yes[1], no[2];
132 template<
typename C>
static yes& test(
typename enable_if_header_type<HeaderType,int>::type);
133 template<
typename C>
static no& test(...);
135 static bool const value =
sizeof(test<HeaderType>(0)) ==
sizeof(yes);
142 template <
class ImageType,
typename ReturnType>
144 typedef decltype ((void) (
145 std::declval<ImageType>().ndim() +
146 std::declval<ImageType>().size(0) +
147 std::declval<ImageType>().name().size() +
148 std::declval<ImageType>().value() +
149 std::declval<ImageType>().index(0)
150 ), std::declval<ReturnType>()) type;
155 template<
typename ImageType>
157 typedef char yes[1], no[2];
158 template<
typename C>
static yes& test(
typename enable_if_image_type<ImageType,int>::type);
159 template<
typename C>
static no& test(...);
161 static bool const value =
sizeof(test<ImageType>(0)) ==
sizeof(yes);
170 template<
class ImageType>
172 static bool const value = std::is_same<ImageType, ::MR::Image<typename ImageType::value_type>>
::value;
176 template<
class ImageType>
193 template <
class ImageType>
195 assign_pos_of (
const ImageType& reference,
size_t from_axis = 0,
size_t to_axis = std::numeric_limits<size_t>::max())
197 return { reference, from_axis, to_axis };
210 template <
class ImageType,
typename IntType>
212 assign_pos_of (
const ImageType& reference,
const vector<IntType>& axes)
214 return { reference, axes };
217 template <
class ImageType,
typename IntType>
219 assign_pos_of (
const ImageType& reference,
const vector<IntType>&& axes)
221 return assign_pos_of (reference, axes);
226 template <
class ImageType>
227 FORCE_INLINE bool is_out_of_bounds (
const ImageType& image,
228 size_t from_axis = 0,
size_t to_axis = std::numeric_limits<size_t>::max())
230 for (
size_t n = from_axis; n < std::min<size_t> (to_axis, image.ndim()); ++n)
231 if (image.index(n) < 0 || image.index(n) >= image.size(n))
236 template <
class HeaderType,
class VectorType>
238 size_t from_axis = 0,
size_t to_axis = std::numeric_limits<size_t>::max())
240 for (
size_t n = from_axis; n < std::min<size_t> (to_axis, header.ndim()); ++n)
241 if (pos[n] < 0 || pos[n] >= header.size(n))
248 template <
class HeaderType>
249 FORCE_INLINE void check_3D_nonunity (
const HeaderType& in)
252 throw Exception (
"Image \"" + in.name() +
"\" does not represent spatial data (less than 3 dimensions)");
253 if (std::min ({ in.size(0), in.size(1), in.size(2) }) == 1)
254 throw Exception (
"Image \"" + in.name() +
"\" does not represent spatial data (has axis with size 1)");
260 template <
class HeaderType>
261 FORCE_INLINE void check_effective_dimensionality (
const HeaderType& in,
size_t N)
264 throw Exception (
"Image \"" + in.name() +
"\" does not represent " +
str(N) +
"D data (too few dimensions)");
265 for (
size_t i = N; i < in.ndim(); ++i)
267 throw Exception (
"Image \"" + in.name() +
"\" does not represent " +
str(N) +
"D data (axis " +
str(i) +
" has size " +
str(in.size(i)) +
")");
271 template <
class HeaderType>
272 inline size_t voxel_count (
const HeaderType& in,
size_t from_axis = 0,
size_t to_axis = std::numeric_limits<size_t>::max())
274 if (to_axis > in.ndim()) to_axis = in.ndim();
275 assert (from_axis < to_axis);
277 for (
size_t n = from_axis; n < to_axis; ++n)
283 template <
class HeaderType>
284 inline size_t voxel_count (
const HeaderType& in,
const char* specifier)
287 for (
size_t n = 0; n < in.ndim(); ++n)
288 if (specifier[n] !=
' ') fp *= in.size (n);
293 template <
class HeaderType>
294 inline size_t voxel_count (
const HeaderType& in,
const std::initializer_list<size_t> axes)
303 template <
class HeaderType>
304 inline int64_t voxel_count (
const HeaderType& in,
const vector<size_t>& axes)
307 for (
size_t n = 0; n < axes.size(); ++n) {
308 assert (axes[n] < in.ndim());
309 fp *= in.size (axes[n]);
314 template <
typename ValueType>
315 inline int64_t footprint (int64_t count) {
316 return count *
sizeof(ValueType);
320 inline int64_t footprint<bool> (int64_t count) {
324 inline int64_t footprint (int64_t count, DataType dtype) {
325 return dtype ==
DataType::Bit ? (count+7)/8 : count*dtype.bytes();
329 template <
class HeaderType>
331 return footprint (voxel_count (in, from_dim, up_to_dim), in.datatype());
335 template <
class HeaderType>
337 return footprint (voxel_count (in, specifier), in.datatype());
342 template <
class HeaderType1,
class HeaderType2>
343 inline bool spacings_match (
const HeaderType1& in1,
const HeaderType2& in2,
const double tol=0.0)
345 if (in1.ndim() != in2.ndim())
return false;
346 for (
size_t n = 0; n < in1.ndim(); ++n)
347 if (
abs(in1.spacing (n) - in2.spacing (n)) > tol * 0.5 * (in1.spacing (n) + in2.spacing (n)))
return false;
351 template <
class HeaderType1,
class HeaderType2>
352 inline bool spacings_match (
const HeaderType1& in1,
const HeaderType2& in2,
size_t from_axis,
size_t to_axis,
const double tol=0.0)
354 assert (from_axis < to_axis);
355 if (to_axis > in1.ndim() || to_axis > in2.ndim())
return false;
356 for (
size_t n = from_axis; n < to_axis; ++n)
357 if (
abs(in1.spacing (n) - in2.spacing (n)) > tol * 0.5 * (in1.spacing (n) + in2.spacing (n)))
return false;
361 template <
class HeaderType1,
class HeaderType2>
362 inline bool spacings_match (
const HeaderType1& in1,
const HeaderType2& in2,
const vector<size_t>& axes,
const double tol=0.0)
364 for (
size_t n = 0; n < axes.size(); ++n) {
365 if (in1.ndim() <= axes[n] || in2.ndim() <= axes[n])
return false;
366 if (
abs(in1.spacing (axes[n]) - in2.spacing(axes[n])) > tol * 0.5 * (in1.spacing (axes[n]) + in2.spacing(axes[n])))
return false;
373 template <
class HeaderType1,
class HeaderType2>
374 inline bool dimensions_match (
const HeaderType1& in1,
const HeaderType2& in2)
376 if (in1.ndim() != in2.ndim())
return false;
377 for (
size_t n = 0; n < in1.ndim(); ++n)
378 if (in1.size (n) != in2.size (n))
return false;
382 template <
class HeaderType1,
class HeaderType2>
383 inline bool dimensions_match (
const HeaderType1& in1,
const HeaderType2& in2,
size_t from_axis,
size_t to_axis)
385 assert (from_axis < to_axis);
386 if (to_axis > in1.ndim() || to_axis > in2.ndim())
return false;
387 for (
size_t n = from_axis; n < to_axis; ++n)
388 if (in1.size (n) != in2.size (n))
return false;
392 template <
class HeaderType1,
class HeaderType2>
393 inline bool dimensions_match (
const HeaderType1& in1,
const HeaderType2& in2,
const vector<size_t>& axes)
395 for (
size_t n = 0; n < axes.size(); ++n) {
396 if (in1.ndim() <= axes[n] || in2.ndim() <= axes[n])
return false;
397 if (in1.size (axes[n]) != in2.size (axes[n]))
return false;
404 template <
class HeaderType>
405 std::string dim2str (
const HeaderType& in)
407 std::string msg =
str(in.size(0));
409 msg +=
"," +
str(in.size(
axis));
414 template <
class HeaderType1,
class HeaderType2>
415 inline void check_dimensions (
const HeaderType1& in1,
const HeaderType2& in2)
417 if (!dimensions_match (in1, in2))
418 throw Exception (
"dimension mismatch between \"" + in1.name() +
"\" and \"" + in2.name() +
"\"" +
419 " (" + dim2str(in1) +
" vs. " + dim2str(in2) +
")");
422 template <
class HeaderType1,
class HeaderType2>
423 inline void check_dimensions (
const HeaderType1& in1,
const HeaderType2& in2,
size_t from_axis,
size_t to_axis)
425 if (!dimensions_match (in1, in2, from_axis, to_axis))
426 throw Exception (
"dimension mismatch between \"" + in1.name() +
"\" and \"" + in2.name() +
"\" between axes " +
str(from_axis) +
" and " +
str(to_axis-1) +
427 " (" + dim2str(in1) +
" vs. " + dim2str(in2) +
")");
430 template <
class HeaderType1,
class HeaderType2>
431 inline void check_dimensions (
const HeaderType1& in1,
const HeaderType2& in2,
const vector<size_t>& axes)
433 if (!dimensions_match (in1, in2, axes))
434 throw Exception (
"dimension mismatch between \"" + in1.name() +
"\" and \"" + in2.name() +
"\" for axes [" +
join(axes,
",") +
"]" +
435 " (" + dim2str(in1) +
" vs. " + dim2str(in2) +
")");
438 template <
class HeaderType1,
class HeaderType2>
439 inline void check_voxel_grids_match_in_scanner_space (
const HeaderType1& in1,
const HeaderType2& in2,
const double tol = 1.0e-3) {
440 Eigen::IOFormat FullPrecFmt(Eigen::FullPrecision, 0,
", ",
"\n",
"[",
"]");
441 if (!voxel_grids_match_in_scanner_space (in1, in2, tol))
442 throw Exception (
"images \"" + in1.name() +
"\" and \"" + in2.name() +
"\" do not have matching header transforms "
443 +
"\n" +
str(in1.transform().matrix().format(FullPrecFmt))
444 +
"\nvs\n" +
str(in2.transform().matrix().format(FullPrecFmt)) +
")");
449 template <
class HeaderType1,
class HeaderType2>
450 inline bool voxel_grids_match_in_scanner_space (
const HeaderType1 in1,
const HeaderType2 in2,
451 const double tol = 1.0e-3) {
452 if (!dimensions_match(in1, in2, 0, 3))
455 const Eigen::Vector3d vs1 (in1.spacing(0), in1.spacing(1), in1.spacing(2));
456 const Eigen::Vector3d vs2 (in2.spacing(0), in2.spacing(1), in2.spacing(2));
458 Eigen::MatrixXd voxel_coord = Eigen::MatrixXd::Zero(4,4);
459 voxel_coord.row(3).fill(1.0);
460 voxel_coord(0,1) = voxel_coord(0,2) = 0.5 * (in1.size(0) + in2.size(0));
461 voxel_coord(1,1) = voxel_coord(1,3) = 0.5 * (in1.size(1) + in2.size(1));
462 voxel_coord(2,2) = voxel_coord(2,3) = 0.5 * (in1.size(2) + in2.size(2));
464 double diff_in_scannercoord = std::sqrt((vs1.asDiagonal() * in1.transform().matrix() * voxel_coord -
465 vs2.asDiagonal() * in2.transform().matrix() * voxel_coord).colwise().squaredNorm().maxCoeff());
466 DEBUG (
"transforms_match: FOV difference in scanner coordinates: "+
str(diff_in_scannercoord));
467 return diff_in_scannercoord < (0.5*(vs1+vs2)).minCoeff() * tol;
470 template <
class HeaderType>
471 inline void squeeze_dim (HeaderType& in,
size_t from_axis = 3)
473 size_t n = in.ndim();
474 while (in.size(n-1) <= 1 && n > from_axis) --n;
483 template <
class ImageType>
488 Index (
const Index&) =
delete;
494 FORCE_INLINE ssize_t operator++ (
int) {
auto p =
get(); move( 1);
return p; }
495 FORCE_INLINE ssize_t operator-- (
int) {
auto p =
get(); move(-1);
return p; }
500 friend std::ostream&
operator<< (std::ostream& stream,
const Index& p) { stream << p.get();
return stream; }
505 FORCE_INLINE void move (ssize_t amount) { image.move_index (
axis, amount); }
509 template <
class ImageType>
520 template <
typename OtherType>
526 friend std::ostream&
operator<< (std::ostream& stream,
const Value& V) { stream << V.get();
return stream; }
534 template <
class ImageType>
537 ConstRow (ImageType& image,
size_t axis) :
axis (
axis), image (image) { assert (
axis >= 0 &&
axis < image.ndim()); }
538 ssize_t size ()
const {
return image.size (
axis); }
543 template <
typename,
int,
int,
int,
int,
int>
friend class Eigen::Matrix;
544 template <
class Derived>
friend class Eigen::MatrixBase;
545 template <
class OtherImageType>
friend class Row;
549 template <
class ImageType>
551 public ConstRow<ImageType>
557 Row (ImageType& image,
size_t axis) : ConstRow<ImageType> (image,
axis) { }
559 template <
class OtherImageType>
560 Row (ConstRow<OtherImageType>&& other) {
561 assert (image.size(
axis) == other.image.size(other.axis));
562 for (image.index(
axis) = 0, other.image.index(other.axis);
563 image.index(
axis) < image.size(
axis);
564 ++image.index(
axis), ++other.image.index(other.axis))
568 using ConstRow<ImageType>::image;
571#define MRTRIX_OP(ARG) \
572 template <class Derived> \
573 FORCE_INLINE void operator ARG (const Eigen::MatrixBase<Derived>& vec) { \
574 assert (vec.rows() == image.size(axis)); \
575 assert (vec.cols() == 1); \
576 for (image.index(axis) = 0; image.index(axis) < image.size(axis); ++image.index(axis)) \
577 image.value() ARG vec[image.index(axis)]; \
584#define MRTRIX_OP(ARG) \
585 FORCE_INLINE void operator ARG (value_type val) { \
586 for (image.index(axis) = 0; image.index(axis) < image.size(axis); ++image.index(axis)) \
587 image.value() ARG val; \
597 assert (image.size(
axis) == other.image.size(other.axis));
598 for (image.index(
axis) = 0, other.image.index(other.axis) = 0;
599 image.index(
axis) < image.size(
axis);
600 ++image.index(
axis), ++other.image.index(other.axis))
601 image.value() = other.image.value();
604#define MRTRIX_OP(ARG) \
605 template <class OtherImageType> \
606 FORCE_INLINE void operator ARG (ConstRow<OtherImageType>&& other) { \
607 assert (image.size(axis) == other.image.size(other.axis)); \
608 for (image.index(axis) = 0, other.image.index(other.axis) = 0; \
609 image.index(axis) < image.size(axis); \
610 ++image.index(axis), ++other.image.index(other.axis)) \
611 image.value() ARG typename OtherImageType::value_type (other.image.value()); \
625 template <
class Derived,
typename ValueType>
627 {
MEMALIGN (ImageBase<Derived,ValueType>)
634 FORCE_INLINE Helper::Value<Derived>
value () {
return {
static_cast<Derived&
> (*this) }; }
635 FORCE_INLINE ValueType
value ()
const {
return static_cast<const Derived*
>(
this)->get_value(); }
651 FORCE_INLINE Helper::ConstRow<Derived> row (
size_t axis)
const {
return {
static_cast<Derived&
> (*this),
axis }; }
652 FORCE_INLINE Helper::Row<Derived> row (
size_t axis) {
return {
static_cast<Derived&
> (*this),
axis }; }
Array & operator=(const MR::Helper::ConstRow< ImageType > &row)
static constexpr uint8_t Bit
Matrix(const MR::Helper::ConstRow< ImageType > &row)
Derived & operator-=(const MR::Helper::ConstRow< ImageType > &row)
Derived & operator+=(const MR::Helper::ConstRow< ImageType > &row)
VectorType::Scalar value(const VectorType &coefs, typename VectorType::Scalar cos_elevation, typename VectorType::Scalar cos_azimuth, typename VectorType::Scalar sin_azimuth, int lmax)
std::ostream & operator<<(std::ostream &stream, const App::ParsedArgument &arg)
MR::default_type value_type
void set(HeaderType &header, const List &stride)
set the strides of header from a vector<ssize_t>
List get(const HeaderType &header)
return the strides of header as a vector<ssize_t>
void apply(F &&f, T &&t)
invoke f(x) for each entry in t
std::string join(const vector< std::string > &V, const std::string &delimiter)
constexpr std::enable_if< std::is_arithmetic< X >::value &&std::is_unsigned< X >::value, X >::type abs(X x)
std::string str(const T &value, int precision=0)
T to(const std::string &string)