Skip to content
Permalink
Browse files

Capture and use slice encoding axis information

- When importing from DICOM, fill field 'SliceEncodingDirection', which is defined as part of BIDS.
- Whenever MRtrix3 performs a permutation of image axes (whether implicit or explicit), additionally adjust this field, just as is performed for DW and PE schemes.
- If attempting to run dwipreproc with slice-to-volume correction in eddy, ensure that if the slice encoding direction is specified, it corresponds to the third axis (which is what eddy will assume). Permutation of axes as part of dwipreproc may be addressed at a later date.
  • Loading branch information...
Lestropie committed Oct 23, 2017
1 parent b7a2807 commit eb447bc4b47458d20a4e4703447ed8886b517bf7
Showing with 65 additions and 58 deletions.
  1. +4 −1 bin/dwipreproc
  2. +6 −0 cmd/mrcalc.cpp
  3. +14 −0 cmd/mrconvert.cpp
  4. +2 −3 core/file/dicom/mapper.cpp
  5. +18 −7 core/file/json_utils.cpp
  6. +15 −3 core/header.cpp
  7. +2 −36 core/phase_encoding.cpp
  8. +4 −8 core/phase_encoding.h
@@ -147,6 +147,9 @@ if 'dw_scheme' not in dwi_header.keyval():
grad = dwi_header.keyval()['dw_scheme']
if not len(grad) == num_volumes:
app.error('Number of lines in gradient table (' + str(len(grad)) + ') does not match input image (' + str(num_volumes) + ' volumes); check your input data')
if if any(s.startswith('--mporder') for s in eddy_manual_options) and 'SliceEncodingDirection' in dwi_header.keyval():
if dwi_header.keyval()['SliceEncodingDirection'] != 'k':
app.error('DWI header indicates that 3rd spatial axis is not the slice axis; this is not yet compatible with --mporder option in eddy')


# Check the manual options being passed to eddy, ensure they make sense
@@ -714,7 +717,7 @@ else:
# output image, as they may have been useful for controlling pre-processing
# but are no longer required, and will just bloat the key-value listings of
# all subsequent derived images
keys_to_remove = [ 'EchoTime', 'FlipAngle', 'MultibandAccelerationFactor', 'PhaseEncodingDirection', 'RepetitionTime', 'SliceTiming', 'TotalReadoutTime', 'pe_scheme' ]
keys_to_remove = [ 'EchoTime', 'FlipAngle', 'MultibandAccelerationFactor', 'PhaseEncodingDirection', 'RepetitionTime', 'SliceEncodingDirection', 'SliceTiming', 'TotalReadoutTime', 'pe_scheme' ]
clear_property_options = ' ' + ' '.join(['-clear_property '+key for key in keys_to_remove if key in dwi_header.keyval() ])


@@ -611,6 +611,12 @@ void get_header (const StackEntry& entry, Header& header)
PhaseEncoding::clear_scheme (header);
}
}

auto slice_encoding_it = entry.image->keyval().find ("SliceEncodingDirection");
if (slice_encoding_it != entry.image->keyval().end()) {
if (header.keyval()["SliceEncodingDirection"] != slice_encoding_it->second)
header.keyval().erase (header.keyval().find ("SliceEncodingDirection"));
}
}


@@ -12,6 +12,7 @@
*/


#include "axes.h"
#include "command.h"
#include "header.h"
#include "image.h"
@@ -208,6 +209,18 @@ void permute_PE_scheme (Header& H, const vector<int>& axes)



void permute_slice_direction (Header& H, const vector<int>& axes)
{
auto it = H.keyval().find ("SliceEncodingDirection");
if (it == H.keyval().end())
return;
const Eigen::Vector3 orig_dir = Axes::id2dir (it->second);
const Eigen::Vector3 new_dir (orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]);
it->second = Axes::dir2id (new_dir);
}




template <class ImageType>
inline vector<int> set_header (Header& header, const ImageType& input)
@@ -233,6 +246,7 @@ inline vector<int> set_header (Header& header, const ImageType& input)
}
permute_DW_scheme (header, axes);
permute_PE_scheme (header, axes);
permute_slice_direction (header, axes);
} else {
header.ndim() = input.ndim();
axes.assign (input.ndim(), 0);
@@ -213,8 +213,6 @@ namespace MR {
// CSA mosaic defines these in ms; we want them in s
for (auto f : image.mosaic_slices_timing)
slices_timing.push_back (0.001 * f);
H.keyval()["SliceTiming"] = join (slices_timing, ",");
H.keyval()["MultibandAccelerationFactor"] = str (std::count (slices_timing.begin(), slices_timing.end(), 0.0f));
}
} else if (std::isfinite (frame.time_after_start)) {
DEBUG ("Taking slice timing information from CSA TimeAfterStart field");
@@ -237,8 +235,9 @@ namespace MR {
if (slices_acquired_at_zero < (image.images_in_mosaic ? image.images_in_mosaic : dim[1])) {
H.keyval()["SliceTiming"] = join (slices_timing, ",");
H.keyval()["MultibandAccelerationFactor"] = str (slices_acquired_at_zero);
H.keyval()["SliceEncodingDirection"] = "k";
} else {
DEBUG ("All slices acquired at same time; not writing slice timing information");
DEBUG ("All slices acquired at same time; not writing slice encoding information");
}
} else {
DEBUG ("No slice timing information obtained");
@@ -17,6 +17,7 @@
#include "file/json_utils.h"
#include "file/nifti_utils.h"

#include "axes.h"
#include "exception.h"
#include "header.h"
#include "mrtrix.h"
@@ -88,6 +89,13 @@ namespace MR
PhaseEncoding::set_scheme (H, pe_scheme);
INFO ("Phase encoding information read from JSON file modified according to expected header transform realignment");
}
auto slice_encoding_it = H.keyval().find ("SliceEncodingDirection");
if (slice_encoding_it != H.keyval().end() && (order[0] != 0 || order[1] != 1 || order[2] != 2)) {
const Eigen::Vector3 orig_dir (Axes::id2dir (slice_encoding_it->second));
const Eigen::Vector3 new_dir (orig_dir[order[0]], orig_dir[order[1]], orig_dir[order[2]]);
slice_encoding_it->second = Axes::dir2id (new_dir);
INFO ("Slice encoding direction read from JSON file modified according to expected header transform realignment");
}
}


@@ -98,6 +106,7 @@ namespace MR
auto pe_scheme = PhaseEncoding::get_scheme (H);
vector<size_t> order;
File::NIfTI::adjust_transform (H, order);
Header H_adj (H);
if (pe_scheme.rows() && (order[0] != 0 || order[1] != 1 || order[2] != 2 || H.stride(0) < 0 || H.stride(1) < 0 || H.stride(2) < 0)) {
// Assume that image being written to disk is going to have its transform adjusted,
// so modify the phase encoding scheme appropriately before writing to JSON
@@ -107,16 +116,18 @@ namespace MR
new_line[axis] = H.stride (order[axis]) > 0 ? pe_scheme(row, order[axis]) : -pe_scheme(row, order[axis]);
pe_scheme.row (row) = new_line;
}
Header H_adj (H);
PhaseEncoding::set_scheme (H_adj, pe_scheme);
for (const auto& kv : H_adj.keyval())
json[kv.first] = kv.second;
INFO ("Phase encoding information written to JSON file modified according to expected header transform realignment");
} else {
// Straight copy
for (const auto& kv : H.keyval())
json[kv.first] = kv.second;
}
auto slice_encoding_it = H_adj.keyval().find ("SliceEncodingDirection");
if (slice_encoding_it != H_adj.keyval().end() && (order[0] != 0 || order[1] != 1 || order[2] != 2)) {
const Eigen::Vector3 orig_dir (Axes::id2dir (slice_encoding_it->second));
const Eigen::Vector3 new_dir (orig_dir[order[0]], orig_dir[order[1]], orig_dir[order[2]]);
slice_encoding_it->second = Axes::dir2id (new_dir);
INFO ("Slice encoding direction written to JSON file modified according to expected header transform realignment");
}
for (const auto& kv : H_adj.keyval())
json[kv.first] = kv.second;
File::OFStream out (path);
out << json.dump(4);
}
@@ -12,8 +12,9 @@
*/


#include "mrtrix.h"
#include "axes.h"
#include "header.h"
#include "mrtrix.h"
#include "phase_encoding.h"
#include "stride.h"
#include "transform.h"
@@ -69,7 +70,7 @@ namespace MR

namespace {

std::string short_description (const Header& H)
std::string short_description (const Header& H)
{
vector<std::string> dims;
for (size_t n = 0; n < H.ndim(); ++n)
@@ -303,7 +304,7 @@ namespace MR
if (new_datatype != previous_datatype) {
new_datatype.unset_flag (DataType::BigEndian);
new_datatype.unset_flag (DataType::LittleEndian);
if (new_datatype != previous_datatype)
if (new_datatype != previous_datatype)
WARN (std::string ("requested datatype (") + previous_datatype.specifier() + ") not supported - substituting with " + H.datatype().specifier());
}

@@ -576,6 +577,17 @@ namespace MR
PhaseEncoding::set_scheme (*this, pe_scheme);
INFO ("Phase encoding scheme has been modified according to internal header transform realignment");
}

// If there's any slice encoding direction information present in the
// header, that's also necessary to update here
auto slice_encoding_it = keyval().find ("SliceEncodingDirection");
if (slice_encoding_it != keyval().end()) {
const Eigen::Vector3 orig_dir (Axes::id2dir (slice_encoding_it->second));
const Eigen::Vector3 new_dir (orig_dir[perm[0]], orig_dir[perm[1]], orig_dir[perm[2]]);
slice_encoding_it->second = Axes::dir2id (new_dir);
INFO ("Slice encoding direction has been modified according to internal header transform realignment");
}

}


@@ -45,41 +45,7 @@ namespace MR



std::string dir2id (const Eigen::Vector3& axis)
{
if (axis[0] == -1) {
assert (!axis[1]); assert (!axis[2]); return "i-";
} else if (axis[0] == 1) {
assert (!axis[1]); assert (!axis[2]); return "i";
} else if (axis[1] == -1) {
assert (!axis[0]); assert (!axis[2]); return "j-";
} else if (axis[1] == 1) {
assert (!axis[0]); assert (!axis[2]); return "j";
} else if (axis[2] == -1) {
assert (!axis[0]); assert (!axis[1]); return "k-";
} else if (axis[2] == 1) {
assert (!axis[0]); assert (!axis[1]); return "k";
} else {
throw Exception ("Malformed phase-encode direction: \"" + str(axis.transpose()) + "\"");
}
}
Eigen::Vector3 id2dir (const std::string& id)
{
if (id == "i-")
return { -1, 0, 0 };
else if (id == "i")
return { 1, 0, 0 };
else if (id == "j-")
return { 0, -1, 0 };
else if (id == "j")
return { 0, 1, 0 };
else if (id == "k-")
return { 0, 0, -1 };
else if (id == "k")
return { 0, 0, 1 };
else
throw Exception ("Malformed phase-encode identifier: \"" + id + "\"");
}




@@ -111,7 +77,7 @@ namespace MR
const auto it_time = header.keyval().find ("TotalReadoutTime");
if (it_dir != header.keyval().end() && it_time != header.keyval().end()) {
Eigen::Matrix<default_type, 4, 1> row;
row.head<3>() = id2dir (it_dir->second);
row.head<3>() = Axes::id2dir (it_dir->second);
row[3] = to<default_type>(it_time->second);
PE.resize ((header.ndim() > 3) ? header.size(3) : 1, 4);
PE.rowwise() = row.transpose();
@@ -19,10 +19,12 @@
#include <Eigen/Dense>

#include "app.h"
#include "axes.h"
#include "header.h"
#include "file/ofstream.h"



namespace MR
{
namespace PhaseEncoding
@@ -68,13 +70,7 @@ namespace MR



//! convert phase encoding direction between formats
/*! these helper functions convert the definition of
* phase-encoding direction between a 3-vector (e.g.
* [0 1 0] ) and a NIfTI axis identifier (e.g. 'i-')
*/
std::string dir2id (const Eigen::Vector3&);
Eigen::Vector3 id2dir (const std::string&);




@@ -118,7 +114,7 @@ namespace MR
} else {
erase ("pe_scheme");
const Eigen::Vector3 dir { PE(0, 0), PE(0, 1), PE(0, 2) };
header.keyval()["PhaseEncodingDirection"] = dir2id (dir);
header.keyval()["PhaseEncodingDirection"] = Axes::dir2id (dir);
if (PE.cols() >= 4)
header.keyval()["TotalReadoutTime"] = str(PE(0, 3), 3);
else

0 comments on commit eb447bc

Please sign in to comment.
You can’t perform that action at this time.