Skip to content

Commit

Permalink
correct orientation of axes from MRC file
Browse files Browse the repository at this point in the history
- mrcfile only reshapes the data to nz,ny,nx but does not take
  the mapc, mapr, maps into account. For gridData we need to
  do the coordinate system reorientation after transposing the
  mrcfile.data to nx, ny, nz.
- add test (using CCP4_1JZV from the test data which contains
  mapc, mapr, maps = 2, 1, 3
- manually checked 1cbs.ccp4 from issue #76 (also a 2, 1, 3)
  and orientation is now correct but there is still an offset error
  that needs to be fixed
  • Loading branch information
orbeckst committed Feb 19, 2022
1 parent 6b326b9 commit 933cf9d
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 5 deletions.
13 changes: 8 additions & 5 deletions gridData/mrc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 7 additions & 0 deletions gridData/tests/test_mrc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")

0 comments on commit 933cf9d

Please sign in to comment.