Skip to content

Extended header in SERI format #71

Description

@eisfabian

Hi Colin,

As discussed last week, here is my (Copilot supported) code for dealing with SERI extended headers. If you give me some pointers on where to start, I might give it a try to fit it into your code myself.

Thanks again and best,
Fabian

    def readExtendedHeader(self):
        """Read extended header."""

        with open(self.file_path, "rb") as mrc_file:
            mrc = mrc_file.read(4096)
            
        # MRC header format SERI (SERI is not covered by mrcfile)
        format = ""
        for char in struct.iter_unpack("s", mrc[104:108]):
            format += char[0].decode("utf-8")

        if format == "SERI":
            # Get number of sections
            section_number = struct.unpack("i", mrc[8: 12])[0]

            # Get bytes per section
            bytes_per_section = struct.unpack("h", mrc[128: 130])[0]

            # Bitflags
            self.extended_header_bitflags = struct.unpack("h", mrc[130: 132])[0]
            """
            https://bio3d.colorado.edu/imod/doc/mrc_format.txt
            1 = Tilt angle in degrees * 100  (2 bytes)
            2 = X, Y, Z piece coordinates for montage (6 bytes)
            4 = X, Y stage position in microns * 25   (4 bytes)
            8 = Magnification / 100 (2 bytes)
            16 = Intensity * 25000  (2 bytes)
            32 = Exposure dose in e-/A2, a float in 4 bytes
            128, 512: Reserved for 4-byte items
            64, 256, 1024: Reserved for 2-byte items
            If the number of bytes implied by these flags does
            not add up to the value in nint, then nint and nreal
            are interpreted as ints and reals per section
            """

            section_data = []
            for i in range(1024, 1024 + bytes_per_section * section_number, bytes_per_section): # extended header starts at byte 1024
                section_data.append(self.parseSectionData(mrc[i:i + bytes_per_section], self.extended_header_bitflags))

            self.extended_header = section_data

    @staticmethod
    def parseSectionData(section_data, bitflags):
        """Parse a single section of extended header data based on bitflags."""

        parsed_data = {}
        offset = 0

        # Check each bitflag and extract corresponding data
        if bitflags & 1:  # Tilt angle in degrees * 100 (2 bytes)
            parsed_data["TiltAngle"] = struct.unpack('h', section_data[offset:offset + 2])[0] / 100
            offset += 2

        if bitflags & 2:  # X, Y, Z piece coordinates for montage (6 bytes)
            parsed_data["PieceCoordinates"] = struct.unpack('hhh', section_data[offset:offset + 6])
            offset += 6

        if bitflags & 4:  # X, Y stage position in microns * 25 (4 bytes)
            parsed_data["StagePosition"] = tuple(coord / 25 for coord in struct.unpack('hh', section_data[offset:offset + 4]))
            offset += 4

        if bitflags & 8:  # Magnification / 100 (2 bytes)
            parsed_data["Magnification"] = struct.unpack('h', section_data[offset:offset + 2])[0] * 100
            offset += 2

        if bitflags & 16:  # Intensity * 25000 (2 bytes)
            parsed_data["Intensity"] = struct.unpack('h', section_data[offset:offset + 2])[0] / 25000
            offset += 2

        if bitflags & 32:  # Exposure dose in e-/A² (4 bytes, float)
            parsed_data["ExposureDose"] = struct.unpack('f', section_data[offset:offset + 4])[0]
            offset += 4

        # Reserved fields (optional, if needed)
        if bitflags & 128 or bitflags & 512:  # Reserved for 4-byte items
            offset += 4
        if bitflags & 64 or bitflags & 256 or bitflags & 1024:  # Reserved for 2-byte items
            offset += 2

        return parsed_data


    def writeExtendedHeader(self, file_path=None):
        """Write extended header."""

        file_path = Path(file_path) if file_path else self.file_path

        # If extended header comes from SERI format, it is a list
        if isinstance(self.extended_header, (list, tuple)):
            with open(file_path, "r+b") as mrc_file:
                # Write the extended header type (e.g., "SERI") to bytes 104–108
                mrc_file.seek(104)
                mrc_file.write(b"SERI")

                # Write the number of sections to bytes 8–12
                section_number = len(self.extended_header)
                mrc_file.seek(8)
                mrc_file.write(struct.pack("i", section_number))

                # Write the bytes per section to bytes 128–130
                bytes_per_section = len(self.serializeSectionData(self.extended_header[0], self.extended_header_bitflags))
                mrc_file.seek(128)
                mrc_file.write(struct.pack("h", bytes_per_section))

                # Write the bitflags to bytes 130–132
                mrc_file.seek(130)
                mrc_file.write(struct.pack("h", self.extended_header_bitflags))

                # Write the section data
                mrc_file.seek(1024)
                mrc_file.write(b"".join([self.serializeSectionData(section, self.extended_header_bitflags) for section in self.extended_header]))

        # If extended header comes from FEI format, use mrcfile built-in function
        else:
            with mrcfile.open(file_path, "r+") as mrc_file:
                mrc_file.set_extended_header(self.extended_header)

    @staticmethod
    def serializeSectionData(section_data, bitflags):
        """Convert parsed section data back to bytes based on bitflags."""
        
        section_bytes = b""

        # Serialize each field based on the bitflags
        if bitflags & 1:  # Tilt angle in degrees * 100 (2 bytes)
            tilt_angle = int(section_data.get("TiltAngle", 0) * 100)
            section_bytes += struct.pack('h', tilt_angle)

        if bitflags & 2:  # X, Y, Z piece coordinates for montage (6 bytes)
            piece_coordinates = section_data.get("PieceCoordinates", (0, 0, 0))
            section_bytes += struct.pack('hhh', *piece_coordinates)

        if bitflags & 4:  # X, Y stage position in microns * 25 (4 bytes)
            stage_position = tuple(int(coord * 25) for coord in section_data.get("StagePosition", (0, 0)))
            section_bytes += struct.pack('hh', *stage_position)

        if bitflags & 8:  # Magnification / 100 (2 bytes)
            magnification = int(section_data.get("Magnification", 0) / 100)
            section_bytes += struct.pack('h', magnification)

        if bitflags & 16:  # Intensity * 25000 (2 bytes)
            intensity = int(section_data.get("Intensity", 0) * 25000)
            section_bytes += struct.pack('h', intensity)

        if bitflags & 32:  # Exposure dose in e-/A² (4 bytes, float)
            exposure_dose = section_data.get("ExposureDose", 0.0)
            section_bytes += struct.pack('f', exposure_dose)

        # Reserved fields (optional, if needed)
        if bitflags & 128 or bitflags & 512:  # Reserved for 4-byte items
            section_bytes += b'\x00\x00\x00\x00'
        if bitflags & 64 or bitflags & 256 or bitflags & 1024:  # Reserved for 2-byte items
            section_bytes += b'\x00\x00'

        return section_bytes

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions