Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix MRC fileformat #445

Merged
merged 2 commits into from Apr 23, 2021
Merged

fix MRC fileformat #445

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8aaf9fe
Non regression test
kif Apr 23, 2021
fb1632a
Read missing "extended headers"
kif Apr 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
19 changes: 14 additions & 5 deletions fabio/mrcimage.py
View file
Open in desktop
Expand Up @@ -31,13 +31,16 @@

Specifications from:
http://ami.scripps.edu/software/mrctools/mrc_specification.php
New version on:
https://www.ccpem.ac.uk/mrc_format/mrc_format.php
https://www.fei-software-center.com/tem-apps/MRC-2014-Specifications/
"""

__authors__ = ["Jérôme Kieffer"]
__contact__ = "Jerome.Kieffer@terre-adelie.org"
__license__ = "MIT"
__copyright__ = "Jérôme Kieffer"
__date__ = "03/04/2020"
__date__ = "23/04/2021"

import logging
import numpy
Expand Down Expand Up @@ -85,12 +88,18 @@ def _readheader(self, infile):
int_block = numpy.frombuffer(infile.read(56 * 4), dtype=numpy.int32)
for key, value in zip(self.KEYS, int_block):
self.header[key] = value
if self.header["MAP"] != 542130509:
logger.info("Expected 'MAP ', got %s", self.header["MAP"].tobytes())
# convert some headers ...
self.header["MAP"] = self.header["MAP"].tobytes().decode()
if self.header["MAP"][:3] not in ('MAP ', 'FEI'):
logger.info("Expected 'MAP ', got %s", self.header["MAP"])

for i in range(10):
label = "LABEL_%02i" % i
self.header[label] = infile.read(80).strip()
self.header[label] = infile.read(80).decode().strip()

# Read extended header
self.header["extended"] = infile.read(self.header["NSYMBT"])

dim1 = int(self.header["NX"])
dim2 = int(self.header["NY"])
self._shape = dim2, dim1
Expand Down Expand Up @@ -126,7 +135,7 @@ def _calc_offset(self, frame):
:param frame: frame number
"""
assert frame < self.nframes
return 1024 + frame * self.imagesize
return 1024 + self.header["NSYMBT"] + frame * self.imagesize

def _makeframename(self):
self.filename = "%s$%04d" % (self.sequencefilename,
Expand Down
57 changes: 31 additions & 26 deletions fabio/test/codecs/test_mrcimage.py
View file
Open in desktop
Expand Up @@ -46,43 +46,48 @@

class TestMrc(unittest.TestCase):
"""basic test"""
mrcfilename = 'EMD-3001.map'
npyfilename = 'EMD-3001.npy'
results = """EMD-3001.map 73 43 -0.36814222 0.678705 0.006804995 0.1630334"""
mrcfilename = ('EMD-3001.map', "0p67_5s_0000.mrc")
npyfilename = ('EMD-3001.npy', "0p67_5s_0000.npy")
results = """EMD-3001.map 73 43 -0.36814222 0.678705 0.0062340125 0.16349247
0p67_5s_0000.mrc 2048 2048 -344.0 33553.0 82.653259 243.213061"""

def setUp(self):
"""Download files"""
self.fn = {}
for i in [self.mrcfilename, self.npyfilename]:
for i in self.mrcfilename + self.npyfilename:
self.fn[i] = UtilsTest.getimage(i + ".bz2")[:-4]
for i in self.fn:
assert os.path.exists(self.fn[i])

def test_read(self):
""" check we can read kcd images"""
vals = self.results.split()
dim1, dim2 = [int(x) for x in vals[1:3]]
shape = dim2, dim1
mini, maxi, mean, stddev = [float(x) for x in vals[3:]]
for ext in ["", ".gz", ".bz2"]:
try:
obj = openimage(self.fn[self.mrcfilename] + ext)
except Exception as err:
logger.error("unable to read: %s", self.fn[self.mrcfilename] + ext)
raise err
self.assertAlmostEqual(mini, obj.getmin(), 4, "getmin" + ext)
self.assertAlmostEqual(maxi, obj.getmax(), 4, "getmax" + ext)
self.assertAlmostEqual(mean, obj.getmean(), 4, "getmean" + ext)
self.assertAlmostEqual(stddev, obj.getstddev(), 4, "getstddev" + ext)
self.assertEqual(shape, obj.shape, "shape" + ext)
""" check we can read mrc images"""
for results in self.results.split("\n"):
vals = results.split()
dim1, dim2 = [int(x) for x in vals[1:3]]
shape = dim2, dim1
mini, maxi, mean, stddev = [float(x) for x in vals[3:]]
for ext in ["", ".gz", ".bz2"]:
filename = self.fn[vals[0]] + ext
try:
obj = openimage(filename)
except Exception as err:
logger.error("unable to read: %s", filename)
raise err
self.assertAlmostEqual(mini, obj.getmin(), 4, f"{filename} getmin")
self.assertAlmostEqual(maxi, obj.getmax(), 4, f"{filename} getmax")
self.assertAlmostEqual(mean, obj.getmean(), 4, f"{filename} getmean")
self.assertAlmostEqual(stddev, obj.getstddev(), 4, f"{filename} getstddev")
self.assertEqual(shape, obj.shape, f"{filename} shape")

def test_same(self):
""" see if we can read kcd images and if they are the same as the EDF """
mrc = MrcImage()
mrc.read(self.fn[self.mrcfilename])
npy = fabio.open(self.fn[self.npyfilename])
diff = (mrc.data - npy.data)
self.assertAlmostEqual(abs(diff).sum(), 0, 4)
""" see if we can read mrc images and if they are the same as the numpy dump"""
for mrcfilename, npyfilename in zip(self.mrcfilename, self.npyfilename):
logger.info("Comparing files %s and %s", mrcfilename, npyfilename)
mrc = MrcImage()
mrc.read(self.fn[mrcfilename])
npy = fabio.open(self.fn[npyfilename])
diff = (mrc.data - npy.data)
self.assertAlmostEqual(abs(diff).sum(), 0, 4, msg=mrcfilename)


def suite():
Expand Down