Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions cpp/mrd-tools/converters.cc
Original file line number Diff line number Diff line change
Expand Up @@ -964,6 +964,7 @@ ISMRMRD::Waveform convert(mrd::Waveform<uint32_t>& wfm) {
// Convert mrd::Image<T> to ISMRMRD::Image<T>
template <class T>
ISMRMRD::Image<T> convert(mrd::Image<T>& image) {
// ISMRMRD Image takes MRD Image without the frequency dimension (x, y, z, channels)
ISMRMRD::Image<T> im(image.Cols(), image.Rows(), image.Slices(), image.Channels());

mrd::ImageHeader& head = image.head;
Expand Down Expand Up @@ -1053,7 +1054,9 @@ ISMRMRD::Image<T> convert(mrd::Image<T>& image) {
for (int z = 0; z < im.getMatrixSizeZ(); z++) {
for (int y = 0; y < im.getMatrixSizeY(); y++) {
for (int x = 0; x < im.getMatrixSizeX(); x++) {
im(x, y, z, c) = image.data(c, z, y, x);
// ISMRMRD does not support frequency dimension from MRD image
// Take 0th image on frequency dimension from MRD to ISMRMRD
im(x, y, z, c) = image.data(c, z, y, x, 0);
}
}
}
Expand Down Expand Up @@ -2060,18 +2063,19 @@ mrd::Image<T> convert(ISMRMRD::Image<T>& im) {
image.head.user_float.push_back(im.getUserFloat(i));
}

mrd::ImageData<T> data({im.getNumberOfChannels(), im.getMatrixSizeZ(), im.getMatrixSizeY(), im.getMatrixSizeX()});
// Add an extra singleton dimension from ISMRMRD to represent (channel, z, y, x, frequency) dimensions for MRD image
mrd::ImageData<T> data({im.getNumberOfChannels(), im.getMatrixSizeZ(), im.getMatrixSizeY(), im.getMatrixSizeX(), 0});
for (int c = 0; c < im.getNumberOfChannels(); c++) {
for (int z = 0; z < im.getMatrixSizeZ(); z++) {
for (int y = 0; y < im.getMatrixSizeY(); y++) {
for (int x = 0; x < im.getMatrixSizeX(); x++) {
data(c, z, y, x) = im(x, y, z, c);
data(c, z, y, x, 0) = im(x, y, z, c);
}
}
}
}

image.data = xt::view(data, xt::all(), xt::all(), xt::all(), xt::all());
image.data = xt::view(data, xt::all(), xt::all(), xt::all(), xt::all(), xt::all());

std::string attrib;
im.getAttributeString(attrib);
Expand Down
16 changes: 8 additions & 8 deletions cpp/mrd-tools/mrd_phantom.cc
Original file line number Diff line number Diff line change
Expand Up @@ -22,24 +22,24 @@ xt::xtensor<std::complex<float>, N> generate_noise(std::array<size_t, N> shape,
}

mrd::ImageData<std::complex<float>> generate_coil_kspace(size_t matrix, size_t ncoils, uint32_t oversampling) {
xt::xtensor<std::complex<float>, 4> phan = shepp_logan_phantom(matrix);
xt::xtensor<std::complex<float>, 4> coils = generate_birdcage_sensitivities(matrix, ncoils, 1.5);
xt::xtensor<std::complex<float>, 5> phan = shepp_logan_phantom(matrix);
xt::xtensor<std::complex<float>, 5> coils = generate_birdcage_sensitivities(matrix, ncoils, 1.5);
coils = phan * coils;

if (oversampling > 1) {
std::array<size_t, 4> padded_shape = coils.shape();
std::array<size_t, 5> padded_shape = coils.shape();
padded_shape[3] *= oversampling;
xt::xtensor<std::complex<float>, 4> padded = xt::zeros<std::complex<float>>(padded_shape);
xt::xtensor<std::complex<float>, 5> padded = xt::zeros<std::complex<float>>(padded_shape);
auto pad = (oversampling - 1) * matrix / 2;
xt::view(padded, xt::all(), xt::all(), xt::all(), xt::range(pad, pad + matrix)) = coils;
xt::view(padded, xt::all(), xt::all(), xt::all(), xt::range(pad, pad + matrix), 0) = coils;
coils = padded;
}

coils = fftshift(coils);
for (unsigned int c = 0; c < ncoils; c++) {
auto tmp1 = xt::fftw::fft2(xt::xarray<std::complex<float>>(xt::view(coils, c, 0, xt::all(), xt::all())));
auto tmp1 = xt::fftw::fft2(xt::xarray<std::complex<float>>(xt::view(coils, c, 0, xt::all(), xt::all(), 0)));
tmp1 /= std::sqrt(1.0f * tmp1.size());
xt::view(coils, c, 0, xt::all(), xt::all()) = tmp1;
xt::view(coils, c, 0, xt::all(), xt::all(), 0) = tmp1;
}
return fftshift(coils);
}
Expand Down Expand Up @@ -241,7 +241,7 @@ int main(int argc, char** argv) {
acq.head.idx.kspace_encode_step_2 = 0;
acq.head.idx.slice = 0;
acq.head.idx.repetition = r;
acq.data = xt::view(kspace, xt::all(), 0, line, xt::all());
acq.data = xt::view(kspace, xt::all(), 0, line, xt::all(), 0);
w->WriteData(acq);
}
}
Expand Down
8 changes: 4 additions & 4 deletions cpp/mrd-tools/shepp_logan_phantom.cc
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,15 @@
#include "mrd/types.h"

mrd::ImageData<std::complex<float>> phantom(std::vector<PhantomEllipse>& ellipses, unsigned int matrix_size) {
std::array<size_t, 4> shape = {1, 1, matrix_size, matrix_size};
std::array<size_t, 5> shape = {1, 1, matrix_size, matrix_size, 1};
mrd::ImageData<std::complex<float>> out = xt::zeros<std::complex<float>>(shape);
for (std::vector<PhantomEllipse>::iterator it = ellipses.begin(); it != ellipses.end(); it++) {
for (unsigned int y = 0; y < matrix_size; y++) {
float y_co = (1.0f * y - (matrix_size >> 1)) / (matrix_size >> 1);
for (unsigned int x = 0; x < matrix_size; x++) {
float x_co = (1.0f * x - (matrix_size >> 1)) / (matrix_size >> 1);
if (it->isInside(x_co, y_co)) {
out(0, 0, y, x) += std::complex<float>(it->getAmplitude(), 0.0);
out(0, 0, y, x, 0) += std::complex<float>(it->getAmplitude(), 0.0);
}
}
}
Expand Down Expand Up @@ -56,7 +56,7 @@ std::vector<PhantomEllipse> modified_shepp_logan_ellipses() {
mrd::ImageData<std::complex<float>> generate_birdcage_sensitivities(unsigned int matrix_size, unsigned int ncoils, float relative_radius) {
// This function is heavily inspired by the mri_birdcage.m Matlab script in Jeff Fessler's IRT packake
// http://web.eecs.umich.edu/~fessler/code/
std::array<size_t, 4> shape = {ncoils, 1, matrix_size, matrix_size};
std::array<size_t, 5> shape = {ncoils, 1, matrix_size, matrix_size, 1};
mrd::ImageData<std::complex<float>> out = xt::zeros<std::complex<float>>(shape);

for (unsigned int c = 0; c < ncoils; c++) {
Expand All @@ -69,7 +69,7 @@ mrd::ImageData<std::complex<float>> generate_birdcage_sensitivities(unsigned int
float x_co = (1.0f * x - (matrix_size >> 1)) / (matrix_size >> 1) - coilx;
float rr = std::sqrt(x_co * x_co + y_co * y_co);
float phi = atan2(x_co, -y_co) + coil_phase;
out(c, 0, y, x) = std::polar(1 / rr, phi);
out(c, 0, y, x, 0) = std::polar(1 / rr, phi);
}
}
}
Expand Down
14 changes: 10 additions & 4 deletions cpp/mrd/binary/protocols.cc
Original file line number Diff line number Diff line change
Expand Up @@ -485,6 +485,8 @@ struct IsTriviallySerializable<mrd::ImageHeader> {
std::is_standard_layout_v<__T__> &&
IsTriviallySerializable<decltype(__T__::flags)>::value &&
IsTriviallySerializable<decltype(__T__::measurement_uid)>::value &&
IsTriviallySerializable<decltype(__T__::measurement_frequency)>::value &&
IsTriviallySerializable<decltype(__T__::measurement_frequency_label)>::value &&
IsTriviallySerializable<decltype(__T__::field_of_view)>::value &&
IsTriviallySerializable<decltype(__T__::position)>::value &&
IsTriviallySerializable<decltype(__T__::col_dir)>::value &&
Expand All @@ -504,8 +506,8 @@ struct IsTriviallySerializable<mrd::ImageHeader> {
IsTriviallySerializable<decltype(__T__::image_series_index)>::value &&
IsTriviallySerializable<decltype(__T__::user_int)>::value &&
IsTriviallySerializable<decltype(__T__::user_float)>::value &&
(sizeof(__T__) == (sizeof(__T__::flags) + sizeof(__T__::measurement_uid) + sizeof(__T__::field_of_view) + sizeof(__T__::position) + sizeof(__T__::col_dir) + sizeof(__T__::line_dir) + sizeof(__T__::slice_dir) + sizeof(__T__::patient_table_position) + sizeof(__T__::average) + sizeof(__T__::slice) + sizeof(__T__::contrast) + sizeof(__T__::phase) + sizeof(__T__::repetition) + sizeof(__T__::set) + sizeof(__T__::acquisition_time_stamp_ns) + sizeof(__T__::physiology_time_stamp_ns) + sizeof(__T__::image_type) + sizeof(__T__::image_index) + sizeof(__T__::image_series_index) + sizeof(__T__::user_int) + sizeof(__T__::user_float))) &&
offsetof(__T__, flags) < offsetof(__T__, measurement_uid) && offsetof(__T__, measurement_uid) < offsetof(__T__, field_of_view) && offsetof(__T__, field_of_view) < offsetof(__T__, position) && offsetof(__T__, position) < offsetof(__T__, col_dir) && offsetof(__T__, col_dir) < offsetof(__T__, line_dir) && offsetof(__T__, line_dir) < offsetof(__T__, slice_dir) && offsetof(__T__, slice_dir) < offsetof(__T__, patient_table_position) && offsetof(__T__, patient_table_position) < offsetof(__T__, average) && offsetof(__T__, average) < offsetof(__T__, slice) && offsetof(__T__, slice) < offsetof(__T__, contrast) && offsetof(__T__, contrast) < offsetof(__T__, phase) && offsetof(__T__, phase) < offsetof(__T__, repetition) && offsetof(__T__, repetition) < offsetof(__T__, set) && offsetof(__T__, set) < offsetof(__T__, acquisition_time_stamp_ns) && offsetof(__T__, acquisition_time_stamp_ns) < offsetof(__T__, physiology_time_stamp_ns) && offsetof(__T__, physiology_time_stamp_ns) < offsetof(__T__, image_type) && offsetof(__T__, image_type) < offsetof(__T__, image_index) && offsetof(__T__, image_index) < offsetof(__T__, image_series_index) && offsetof(__T__, image_series_index) < offsetof(__T__, user_int) && offsetof(__T__, user_int) < offsetof(__T__, user_float);
(sizeof(__T__) == (sizeof(__T__::flags) + sizeof(__T__::measurement_uid) + sizeof(__T__::measurement_frequency) + sizeof(__T__::measurement_frequency_label) + sizeof(__T__::field_of_view) + sizeof(__T__::position) + sizeof(__T__::col_dir) + sizeof(__T__::line_dir) + sizeof(__T__::slice_dir) + sizeof(__T__::patient_table_position) + sizeof(__T__::average) + sizeof(__T__::slice) + sizeof(__T__::contrast) + sizeof(__T__::phase) + sizeof(__T__::repetition) + sizeof(__T__::set) + sizeof(__T__::acquisition_time_stamp_ns) + sizeof(__T__::physiology_time_stamp_ns) + sizeof(__T__::image_type) + sizeof(__T__::image_index) + sizeof(__T__::image_series_index) + sizeof(__T__::user_int) + sizeof(__T__::user_float))) &&
offsetof(__T__, flags) < offsetof(__T__, measurement_uid) && offsetof(__T__, measurement_uid) < offsetof(__T__, measurement_frequency) && offsetof(__T__, measurement_frequency) < offsetof(__T__, measurement_frequency_label) && offsetof(__T__, measurement_frequency_label) < offsetof(__T__, field_of_view) && offsetof(__T__, field_of_view) < offsetof(__T__, position) && offsetof(__T__, position) < offsetof(__T__, col_dir) && offsetof(__T__, col_dir) < offsetof(__T__, line_dir) && offsetof(__T__, line_dir) < offsetof(__T__, slice_dir) && offsetof(__T__, slice_dir) < offsetof(__T__, patient_table_position) && offsetof(__T__, patient_table_position) < offsetof(__T__, average) && offsetof(__T__, average) < offsetof(__T__, slice) && offsetof(__T__, slice) < offsetof(__T__, contrast) && offsetof(__T__, contrast) < offsetof(__T__, phase) && offsetof(__T__, phase) < offsetof(__T__, repetition) && offsetof(__T__, repetition) < offsetof(__T__, set) && offsetof(__T__, set) < offsetof(__T__, acquisition_time_stamp_ns) && offsetof(__T__, acquisition_time_stamp_ns) < offsetof(__T__, physiology_time_stamp_ns) && offsetof(__T__, physiology_time_stamp_ns) < offsetof(__T__, image_type) && offsetof(__T__, image_type) < offsetof(__T__, image_index) && offsetof(__T__, image_index) < offsetof(__T__, image_series_index) && offsetof(__T__, image_series_index) < offsetof(__T__, user_int) && offsetof(__T__, user_int) < offsetof(__T__, user_float);
};

template <typename T>
Expand Down Expand Up @@ -1883,7 +1885,7 @@ template<typename Y, yardl::binary::Writer<Y> WriteY>
return;
}

yardl::binary::WriteNDArray<Y, WriteY, 4>(stream, value);
yardl::binary::WriteNDArray<Y, WriteY, 5>(stream, value);
}

template<typename Y, yardl::binary::Reader<Y> ReadY>
Expand All @@ -1893,7 +1895,7 @@ template<typename Y, yardl::binary::Reader<Y> ReadY>
return;
}

yardl::binary::ReadNDArray<Y, ReadY, 4>(stream, value);
yardl::binary::ReadNDArray<Y, ReadY, 5>(stream, value);
}

[[maybe_unused]] void WriteImageHeader(yardl::binary::CodedOutputStream& stream, mrd::ImageHeader const& value) {
Expand All @@ -1904,6 +1906,8 @@ template<typename Y, yardl::binary::Reader<Y> ReadY>

yardl::binary::WriteFlags<mrd::ImageFlags>(stream, value.flags);
yardl::binary::WriteInteger(stream, value.measurement_uid);
yardl::binary::WriteOptional<yardl::DynamicNDArray<uint32_t>, yardl::binary::WriteDynamicNDArray<uint32_t, yardl::binary::WriteInteger>>(stream, value.measurement_frequency);
yardl::binary::WriteOptional<yardl::DynamicNDArray<std::string>, yardl::binary::WriteDynamicNDArray<std::string, yardl::binary::WriteString>>(stream, value.measurement_frequency_label);
yardl::binary::WriteFixedNDArray<float, yardl::binary::WriteFloatingPoint, 3>(stream, value.field_of_view);
yardl::binary::WriteFixedNDArray<float, yardl::binary::WriteFloatingPoint, 3>(stream, value.position);
yardl::binary::WriteFixedNDArray<float, yardl::binary::WriteFloatingPoint, 3>(stream, value.col_dir);
Expand Down Expand Up @@ -1933,6 +1937,8 @@ template<typename Y, yardl::binary::Reader<Y> ReadY>

yardl::binary::ReadFlags<mrd::ImageFlags>(stream, value.flags);
yardl::binary::ReadInteger(stream, value.measurement_uid);
yardl::binary::ReadOptional<yardl::DynamicNDArray<uint32_t>, yardl::binary::ReadDynamicNDArray<uint32_t, yardl::binary::ReadInteger>>(stream, value.measurement_frequency);
yardl::binary::ReadOptional<yardl::DynamicNDArray<std::string>, yardl::binary::ReadDynamicNDArray<std::string, yardl::binary::ReadString>>(stream, value.measurement_frequency_label);
yardl::binary::ReadFixedNDArray<float, yardl::binary::ReadFloatingPoint, 3>(stream, value.field_of_view);
yardl::binary::ReadFixedNDArray<float, yardl::binary::ReadFloatingPoint, 3>(stream, value.position);
yardl::binary::ReadFixedNDArray<float, yardl::binary::ReadFloatingPoint, 3>(stream, value.col_dir);
Expand Down
12 changes: 10 additions & 2 deletions cpp/mrd/hdf5/protocols.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1441,6 +1441,8 @@ struct _Inner_ImageHeader {
_Inner_ImageHeader(mrd::ImageHeader const& o)
: flags(o.flags),
measurement_uid(o.measurement_uid),
measurement_frequency(o.measurement_frequency),
measurement_frequency_label(o.measurement_frequency_label),
field_of_view(o.field_of_view),
position(o.position),
col_dir(o.col_dir),
Expand All @@ -1465,6 +1467,8 @@ struct _Inner_ImageHeader {
void ToOuter (mrd::ImageHeader& o) const {
yardl::hdf5::ToOuter(flags, o.flags);
yardl::hdf5::ToOuter(measurement_uid, o.measurement_uid);
yardl::hdf5::ToOuter(measurement_frequency, o.measurement_frequency);
yardl::hdf5::ToOuter(measurement_frequency_label, o.measurement_frequency_label);
yardl::hdf5::ToOuter(field_of_view, o.field_of_view);
yardl::hdf5::ToOuter(position, o.position);
yardl::hdf5::ToOuter(col_dir, o.col_dir);
Expand All @@ -1488,6 +1492,8 @@ struct _Inner_ImageHeader {

mrd::ImageFlags flags;
uint32_t measurement_uid;
yardl::hdf5::InnerOptional<yardl::hdf5::InnerDynamicNdArray<uint32_t, uint32_t>, yardl::DynamicNDArray<uint32_t>> measurement_frequency;
yardl::hdf5::InnerOptional<yardl::hdf5::InnerDynamicNdArray<yardl::hdf5::InnerVlenString, std::string>, yardl::DynamicNDArray<std::string>> measurement_frequency_label;
yardl::FixedNDArray<float, 3> field_of_view;
yardl::FixedNDArray<float, 3> position;
yardl::FixedNDArray<float, 3> col_dir;
Expand Down Expand Up @@ -1525,7 +1531,7 @@ struct _Inner_Image {
}

mrd::hdf5::_Inner_ImageHeader head;
yardl::hdf5::InnerNdArray<_T_Inner, T, 4> data;
yardl::hdf5::InnerNdArray<_T_Inner, T, 5> data;
yardl::hdf5::InnerMap<yardl::hdf5::InnerVlenString, std::string, yardl::hdf5::InnerVlen<::InnerUnion3<yardl::hdf5::InnerVlenString, std::string, int64_t, int64_t, double, double>, mrd::ImageMetaValue>, std::vector<mrd::ImageMetaValue>> meta;
};

Expand Down Expand Up @@ -2060,6 +2066,8 @@ struct _Inner_ImageArray {
H5::CompType t(sizeof(RecordType));
t.insertMember("flags", HOFFSET(RecordType, flags), H5::PredType::NATIVE_UINT64);
t.insertMember("measurementUid", HOFFSET(RecordType, measurement_uid), H5::PredType::NATIVE_UINT32);
t.insertMember("measurementFrequency", HOFFSET(RecordType, measurement_frequency), yardl::hdf5::OptionalTypeDdl<yardl::hdf5::InnerDynamicNdArray<uint32_t, uint32_t>, yardl::DynamicNDArray<uint32_t>>(yardl::hdf5::DynamicNDArrayDdl<uint32_t, uint32_t>(H5::PredType::NATIVE_UINT32)));
t.insertMember("measurementFrequencyLabel", HOFFSET(RecordType, measurement_frequency_label), yardl::hdf5::OptionalTypeDdl<yardl::hdf5::InnerDynamicNdArray<yardl::hdf5::InnerVlenString, std::string>, yardl::DynamicNDArray<std::string>>(yardl::hdf5::DynamicNDArrayDdl<yardl::hdf5::InnerVlenString, std::string>(yardl::hdf5::InnerVlenStringDdl())));
t.insertMember("fieldOfView", HOFFSET(RecordType, field_of_view), yardl::hdf5::FixedNDArrayDdl(H5::PredType::NATIVE_FLOAT, {3}));
t.insertMember("position", HOFFSET(RecordType, position), yardl::hdf5::FixedNDArrayDdl(H5::PredType::NATIVE_FLOAT, {3}));
t.insertMember("colDir", HOFFSET(RecordType, col_dir), yardl::hdf5::FixedNDArrayDdl(H5::PredType::NATIVE_FLOAT, {3}));
Expand Down Expand Up @@ -2087,7 +2095,7 @@ template <typename _T_Inner, typename T>
using RecordType = mrd::hdf5::_Inner_Image<_T_Inner, T>;
H5::CompType t(sizeof(RecordType));
t.insertMember("head", HOFFSET(RecordType, head), mrd::hdf5::GetImageHeaderHdf5Ddl());
t.insertMember("data", HOFFSET(RecordType, data), yardl::hdf5::NDArrayDdl<_T_Inner, T, 4>(T_type));
t.insertMember("data", HOFFSET(RecordType, data), yardl::hdf5::NDArrayDdl<_T_Inner, T, 5>(T_type));
t.insertMember("meta", HOFFSET(RecordType, meta), yardl::hdf5::InnerMapDdl<yardl::hdf5::InnerVlenString, yardl::hdf5::InnerVlen<::InnerUnion3<yardl::hdf5::InnerVlenString, std::string, int64_t, int64_t, double, double>, mrd::ImageMetaValue>>(yardl::hdf5::InnerVlenStringDdl(), yardl::hdf5::InnerVlenDdl(::InnerUnion3Ddl<yardl::hdf5::InnerVlenString, std::string, int64_t, int64_t, double, double>(false, yardl::hdf5::InnerVlenStringDdl(), "string", H5::PredType::NATIVE_INT64, "int64", H5::PredType::NATIVE_DOUBLE, "float64"))));
return t;
}
Expand Down
12 changes: 12 additions & 0 deletions cpp/mrd/ndjson/protocols.cc
Original file line number Diff line number Diff line change
Expand Up @@ -2590,6 +2590,12 @@ void to_json(ordered_json& j, mrd::ImageHeader const& value) {
if (yardl::ndjson::ShouldSerializeFieldValue(value.measurement_uid)) {
j.push_back({"measurementUid", value.measurement_uid});
}
if (yardl::ndjson::ShouldSerializeFieldValue(value.measurement_frequency)) {
j.push_back({"measurementFrequency", value.measurement_frequency});
}
if (yardl::ndjson::ShouldSerializeFieldValue(value.measurement_frequency_label)) {
j.push_back({"measurementFrequencyLabel", value.measurement_frequency_label});
}
if (yardl::ndjson::ShouldSerializeFieldValue(value.field_of_view)) {
j.push_back({"fieldOfView", value.field_of_view});
}
Expand Down Expand Up @@ -2656,6 +2662,12 @@ void from_json(ordered_json const& j, mrd::ImageHeader& value) {
if (auto it = j.find("measurementUid"); it != j.end()) {
it->get_to(value.measurement_uid);
}
if (auto it = j.find("measurementFrequency"); it != j.end()) {
it->get_to(value.measurement_frequency);
}
if (auto it = j.find("measurementFrequencyLabel"); it != j.end()) {
it->get_to(value.measurement_frequency_label);
}
if (auto it = j.find("fieldOfView"); it != j.end()) {
it->get_to(value.field_of_view);
}
Expand Down
2 changes: 1 addition & 1 deletion cpp/mrd/protocols.cc

Large diffs are not rendered by default.

Loading
Loading