diff --git a/gridData/mrc.py b/gridData/mrc.py index e4551bf..803a969 100644 --- a/gridData/mrc.py +++ b/gridData/mrc.py @@ -107,12 +107,15 @@ def read(self, filename): "alpha={0}, beta={1}, gamma={2}".format( h.cellb.alpha, h.cellb.beta, h.cellb.gamma)) # mrc.data[z, y, x] indexed: convert to x,y,z as used in GridDataFormats - # mapc, mapr, maps = 1, 2, 3 for Fortran-ordering and 3, 2, 1 for C-ordering. - # Other combinations are possible. We reorder the data for the general case - # by sorting mapc, mapr, maps in descending order. + # together with the axes orientation information in mapc/mapr/maps. + # mapc, mapr, maps = 1, 2, 3 for Fortran-ordering and 3, 2, 1 for C-ordering. + # Other combinations are possible. We reorder the data for the general case + # by sorting mapc, mapr, maps in ascending order, i.e., to obtain x,y,z. + # mrcfile provides the data in zyx shape (without regard to map*) so we first + # transpose it to xyz and then reorient with axes_c_order. axes_order = np.hstack([mrc.header.mapc, mrc.header.mapr, mrc.header.maps]) - axes_c_order = np.argsort(axes_order)[::-1] - self.array = np.transpose(mrc.data, axes = axes_c_order) + axes_c_order = np.argsort(axes_order) + self.array = np.transpose(np.transpose(mrc.data), axes=axes_c_order) self.delta = np.diag([mrc.voxel_size.x, mrc.voxel_size.y, mrc.voxel_size.z]) self.origin = np.array([h.origin.x, h.origin.y, h.origin.z]) self.rank = 3 diff --git a/gridData/tests/test_mrc.py b/gridData/tests/test_mrc.py index 161a8a3..a5fc49f 100644 --- a/gridData/tests/test_mrc.py +++ b/gridData/tests/test_mrc.py @@ -90,8 +90,15 @@ def test_ccp4_read_header(ccp4data, name, value): else: assert_equal(ccp4data.header[name], value) +def test_axes_orientation(ccp4data): + # correctly interpret mapc, mapr, maps = 2, 1, 3 + # for nx, ny, nz = 96, 76, 70. + # see also #76 + assert_equal(ccp4data.array.shape, (76, 96, 70)) + def test_triclinic_ValueError(): with pytest.raises(ValueError, match="Only orthorhombic unitcells are currently " "supported, not"): Grid(datafiles.MRC_EMD3001, file_format="MRC") +