Source code for coseda.importers.import_mrc

from coseda.logging_utils import log_print
import struct
import mrcfile
import os
import datetime
import configparser
import platform
import psutil
import sys
import h5py
import numpy as np
from coseda.io import handle_input, parse_config, get_free_space_unix, get_free_space_windows, get_dir_size
from coseda.logging_utils import log_start, log_result, shoutout
from coseda.nexus.process import write_beam_incident_energy, write_nxprocess_import
from coseda.nexus.indices import ensure_image_key
from coseda.nexus.groups import ensure_nexus_parents
from coseda.nexus.detector import write_detector_geometry
from coseda.nexus.images import ensure_image_nxdata
from coseda.nexus.logs import ensure_dense_logs
from coseda.nexus.goniometer import ensure_goniometer_transforms, get_goniometer_transform_order


def _is_mrc_file(path):
    return os.path.isfile(path) and path.lower().endswith(".mrc")


def _mode_to_dtype(mode):
    mode = int(mode)
    if mode == 0:
        return "uint8"    # Unsigned 8-bit integer
    if mode == 1:
        return "int16"    # Signed 16-bit integer
    if mode == 2:
        return "float32"  # 32-bit float
    if mode == 6:
        return "uint16"   # Unsigned 16-bit integer
    return None


def _resolve_mrc_sources(input_paths):
    if isinstance(input_paths, str):
        input_paths = [input_paths]

    resolved = []
    for raw in input_paths:
        if raw is None:
            continue
        path = os.path.abspath(raw)
        if _is_mrc_file(path):
            resolved.append(path)
        elif os.path.isdir(path):
            resolved.append(path)
        else:
            log_print(f"Warning: {raw} is not a valid .mrc file or directory.")

    # Preserve order while deduplicating.
    return list(dict.fromkeys(resolved))


def _collect_mrc_files_from_source(source_path):
    source_path = os.path.abspath(source_path)
    if _is_mrc_file(source_path):
        return [source_path]
    if os.path.isdir(source_path):
        return sorted(
            os.path.join(source_path, f)
            for f in os.listdir(source_path)
            if f.lower().endswith(".mrc")
        )
    return []


def _extract_stage_tilt_from_header(mrc_obj):
    if mrc_obj.header.nsymbt <= 0:
        return np.nan, np.nan, np.nan, np.nan, np.nan
    extended_header_bytes = mrc_obj.extended_header.tobytes()
    try:
        alphatilt = struct.unpack('d', extended_header_bytes[100:108])[0]
        betatilt = struct.unpack('d', extended_header_bytes[108:116])[0]
        x_stage = struct.unpack('d', extended_header_bytes[116:124])[0]
        y_stage = struct.unpack('d', extended_header_bytes[124:132])[0]
        z_stage = struct.unpack('d', extended_header_bytes[132:140])[0]
        return alphatilt, betatilt, x_stage, y_stage, z_stage
    except Exception as exc:
        log_print(f"Error reading extended header stage/tilt values: {exc}")
        return np.nan, np.nan, np.nan, np.nan, np.nan


[docs] def extract_info_from_mrc_metadata(input_path): # Open the MRC file and extract relevant metadata with mrcfile.open(input_path, permissive=True) as mrc: # Extract resolution and pixel unit resolution_height, resolution_width = int(mrc.header.ny), int(mrc.header.nx) # Cast to int pixel_unit = "1/m" # Pixel unit as 1/m instrument_manufacturer = "Thermo Fisher Scientific" # Assuming it's always this for MRC # Determine bit depth based on mode (data type) mode = mrc.header.mode bit_depth = None if mode == 0: bit_depth = 8 # 8-bit unsigned integer elif mode == 1: bit_depth = 16 # 16-bit signed integer elif mode == 2: bit_depth = 32 # 32-bit float elif mode == 6: bit_depth = 16 # 16-bit unsigned integer else: bit_depth = "Unknown" # Initialize placeholders for extracted values binning_width = None binning_height = None exposure_time = None acceleration_voltage = None pixel_width = None pixel_height = None instrument_model = None instrument_id = None detector_used = None camera_length = None dose = None # Check if extended header exists and read it if mrc.header.nsymbt > 0: extended_header_bytes = mrc.extended_header.tobytes() try: # Extract parameters from extended header binning_width = struct.unpack('i', extended_header_bytes[427:431])[0] binning_height = struct.unpack('i', extended_header_bytes[431:435])[0] exposure_time = struct.unpack('d', extended_header_bytes[419:427])[0] acceleration_voltage = struct.unpack('d', extended_header_bytes[84:92])[0] pixel_width = struct.unpack('d', extended_header_bytes[156:164])[0] pixel_height = struct.unpack('d', extended_header_bytes[164:172])[0] instrument_model = extended_header_bytes[20:36].decode('utf-8').strip('\x00') instrument_id = extended_header_bytes[36:52].decode('utf-8').strip('\x00') detector_used = extended_header_bytes[435:451].decode('utf-8').strip('\x00') camera_length = struct.unpack('d', extended_header_bytes[301:309])[0] dose = struct.unpack('d', extended_header_bytes[92:100])[0] resolution_height = resolution_height * binning_height resolution_width = resolution_width * binning_width except Exception as e: log_print(f"Error extracting extended header metadata: {e}") # Return all the extracted parameters including bit depth return (exposure_time, resolution_height, resolution_width, binning_height, binning_width, acceleration_voltage, pixel_height, pixel_width, pixel_unit, camera_length, detector_used, instrument_id, instrument_manufacturer, instrument_model, bit_depth)
[docs] def create_insfile_mrc(input_path, target=None, source_path=None): metadata_file_path = os.path.abspath(input_path) source_path = os.path.abspath(source_path or metadata_file_path) if os.path.isdir(source_path): source_name = os.path.basename(source_path) source_kind = "folder" parent_directory = os.path.dirname(source_path) else: source_name = os.path.splitext(os.path.basename(source_path))[0] source_kind = "file" parent_directory = os.path.dirname(source_path) timestamp = datetime.datetime.now().strftime("%Y-%m-%d_%H-%M-%S") run_name = f"{source_name}_run_{timestamp}" if target is None: output_parent = parent_directory else: output_parent = os.path.abspath(target) outputfolder_path = os.path.join(output_parent, run_name) configfile = f"{run_name}.ini" configfile_path = os.path.join(output_parent, configfile) newfiles = configfile_path log_print(f'Configuration file created for source {source_path}') # Create output folder and plotting subfolder. os.makedirs(outputfolder_path, exist_ok=True) os.makedirs(os.path.join(outputfolder_path, "plots"), exist_ok=True) log_print(f"Results will be saved in {outputfolder_path}") logfile = f"{run_name}.log" logfile_path = os.path.join(outputfolder_path, logfile) config = configparser.ConfigParser() config.read(configfile_path) if not config.has_section('Paths'): config.add_section('Paths') config.set('Paths', 'configfile_path', configfile_path) config.set('Paths', 'outputfolder', os.path.basename(outputfolder_path)) config.set('Paths', 'originalfile', source_path) config.set('Paths', 'originalfile_path', source_path) config.set('Paths', 'originalfolder', os.path.basename(os.path.dirname(metadata_file_path))) config.set('Paths', 'logfile', os.path.basename(logfile_path)) config.set('Paths', 'framepath', 'entry/data/images') config.set('Paths', 'mrc_input_type', source_kind) if not config.has_section('AcquisitionDetails'): config.add_section('AcquisitionDetails') try: # Use the extract_info_from_mrc_metadata function to get the metadata (exposure_time, resolution_height, resolution_width, binning_height, binning_width, acceleration_voltage, pixel_height, pixel_width, pixel_unit, camera_length, detector_used, instrument_id, instrument_manufacturer, instrument_model, bit_depth) = extract_info_from_mrc_metadata(metadata_file_path) except Exception as e: log_print(f'Error: No valid metadata found in MRC file: {e}') return False, None # Write the metadata into the .ini file config.set('AcquisitionDetails', 'instrument_manufacturer', str(instrument_manufacturer)) config.set('AcquisitionDetails', 'instrument_model', str(instrument_model)) config.set('AcquisitionDetails', 'instrument_id', str(instrument_id)) config.set('AcquisitionDetails', 'detector_used', str(detector_used)) config.set('AcquisitionDetails', 'acceleration_voltage', str(acceleration_voltage)) config.set('AcquisitionDetails', 'camera_length', str(camera_length)) config.set('AcquisitionDetails', 'pixel_width', str(pixel_width)) config.set('AcquisitionDetails', 'pixel_height', str(pixel_height)) config.set('AcquisitionDetails', 'binning_width', str(binning_width)) config.set('AcquisitionDetails', 'binning_height', str(binning_height)) config.set('AcquisitionDetails', 'resolution_width', str(resolution_width)) config.set('AcquisitionDetails', 'resolution_height', str(resolution_height)) config.set('AcquisitionDetails', 'exposure_time', str(exposure_time)) config.set('AcquisitionDetails', 'pixel_unit', str(pixel_unit)) config.set('AcquisitionDetails', 'bit_depth', str(bit_depth)) if not config.has_section('Parameters'): config.add_section('Parameters') if not config.has_section('Output'): config.add_section('Output') # Write the config file to the parent folder with open(configfile_path, 'w') as configfile: config.write(configfile) # Write log file details log_start(logfile_path, '.ini file created') log_start(logfile_path, f'system: {platform.system()} {platform.architecture()}, version {platform.version()}') log_start(logfile_path, f'hardware: {psutil.cpu_count(logical=False)} physical cores, {psutil.cpu_count(logical=True)} logical cores, {round(psutil.virtual_memory().total / (1024 ** 3))} GB total memory') log_start(logfile_path, f'python version {sys.version}') log_start(logfile_path, f'MRC source type: {source_kind}, source path: {source_path}') log_start(logfile_path, f'data acquisition details: InstrumentManufacturer: {str(instrument_manufacturer)}, InstrumentModel: {str(instrument_model)}, InstrumentId: {str(instrument_id)}, DetectorUsed: {str(detector_used)}, AccelerationVoltage: {str(acceleration_voltage)}, CameraLength: {camera_length}, PixelWidth: {pixel_width}, PixelHeight: {pixel_height}, BinningWidth: {binning_width}, BinningHeight: {binning_height}, Resolution: ({resolution_width}, {resolution_height}), ExposureTime: {exposure_time}, PixelUnit: {pixel_unit}, BitDepth: {bit_depth}') return True, newfiles
[docs] def get_mrc_files(input_paths): mrc_files_in_sources = [] for source in _resolve_mrc_sources(input_paths): files = _collect_mrc_files_from_source(source) if files: mrc_files_in_sources.append(files) else: log_print(f"Warning: no .mrc files found in source '{source}'.") return mrc_files_in_sources
[docs] def create_insfiles_mrc_batch(input_path, target=None): ini_files = [] for source in _resolve_mrc_sources(input_path): source_files = _collect_mrc_files_from_source(source) if not source_files: log_print(f"Warning: no .mrc files found in source '{source}'.") continue metadata_file = source_files[0] success, newfiles = create_insfile_mrc(metadata_file, target=target, source_path=source) if success and newfiles: ini_files.append(newfiles) return ini_files
[docs] def mrc_folder_conversion( mrc_source, outputfolder, logfile_path, chunk_size=1000, acceleration_voltage=None, camera_length=None, camera_length_correction=1.0, goniometer_transform_order=None, ): try: h5file = outputfolder + ".h5" parentfolderpath = os.path.dirname(logfile_path) h5file_path = os.path.join(parentfolderpath, h5file) # Collect all source .mrc files: # - source is a folder of .mrc files # - or source is a single .mrc stack/movie all_files = _collect_mrc_files_from_source(mrc_source) if not all_files: return "No MRC files found in the provided source.", None # Scan inputs and validate they can be merged into a single frame stack. total_frames = 0 frame_height = None frame_width = None dtype = None mode = None per_file_frame_counts = [] for file_path in all_files: with mrcfile.open(file_path, permissive=True) as mrc: data_shape = tuple(int(v) for v in mrc.data.shape) file_mode = int(mrc.header.mode) file_dtype = _mode_to_dtype(file_mode) if file_dtype is None: return f"Unsupported MRC mode: {file_mode} in '{os.path.basename(file_path)}'", None if len(data_shape) == 2: file_frames = 1 this_h, this_w = data_shape elif len(data_shape) == 3: file_frames, this_h, this_w = data_shape else: return ( f"Unsupported MRC data shape {data_shape} in '{os.path.basename(file_path)}'. " "Expected 2D image or 3D stack." ), None if frame_height is None: frame_height, frame_width = int(this_h), int(this_w) elif (frame_height, frame_width) != (int(this_h), int(this_w)): return ( f"Inconsistent frame size in source files. Expected {(frame_height, frame_width)}, " f"got {(int(this_h), int(this_w))} in '{os.path.basename(file_path)}'." ), None if dtype is None: dtype = file_dtype mode = file_mode elif dtype != file_dtype: return ( f"Inconsistent MRC modes in source files ({mode} vs {file_mode}); " "mixed dtypes are not supported." ), None total_frames += int(file_frames) per_file_frame_counts.append((file_path, int(file_frames))) if total_frames <= 0: return "No frames found in the provided MRC source.", None # Build chunk shape as (frames, height, width). if isinstance(chunk_size, (tuple, list)): chunk_frames = int(chunk_size[0]) else: chunk_frames = int(chunk_size) chunk_frames = max(1, min(chunk_frames, total_frames)) chunk_shape = (chunk_frames, frame_height, frame_width) new_framepath = 'entry/data/images' new_indexpath = 'entry/data/index' stagepos_x_path = 'entry/data/stagepos_x' stagepos_x_refined_path = 'entry/data/stagepos_x_refined' stagepos_y_path = 'entry/data/stagepos_y' stagepos_y_refined_path = 'entry/data/stagepos_y_refined' stagepos_z_path = 'entry/data/stagepos_z' alphatilt_path = 'entry/data/alphatilt' betatilt_path = 'entry/data/betatilt' # Disk space checks if os.path.isdir(mrc_source): file_size = get_dir_size(mrc_source) free_space_base = os.path.dirname(mrc_source) else: file_size = os.path.getsize(mrc_source) free_space_base = os.path.dirname(mrc_source) if platform.system() == 'Windows': free_space = get_free_space_windows(free_space_base) else: free_space = get_free_space_unix(free_space_base) if free_space < 1.1 * file_size: return f"Insufficient disk space, ensure at least {1.1 * file_size} free.", None with h5py.File(h5file_path, 'w') as new_file: log_start(logfile_path, "HDF5 file created") new_dataset = new_file.create_dataset( new_framepath, shape=(total_frames, frame_height, frame_width), chunks=chunk_shape, dtype=dtype, ) index_dataset = new_file.create_dataset(new_indexpath, shape=(total_frames,), dtype='i4') # Buffer stage positions and only write datasets when finite values exist. stagepos_x_values = np.full(total_frames, np.nan, dtype=np.float32) stagepos_y_values = np.full(total_frames, np.nan, dtype=np.float32) stagepos_z_values = np.full(total_frames, np.nan, dtype=np.float32) # Keep tilt datasets for compatibility with existing workflows. alphatilt_dataset = new_file.create_dataset(alphatilt_path, shape=(total_frames,), dtype='float32') betatilt_dataset = new_file.create_dataset(betatilt_path, shape=(total_frames,), dtype='float32') frame_idx = 0 save_interval = 10000 for file_path, frames_in_file in per_file_frame_counts: with mrcfile.open(file_path, permissive=True) as mrc: data = mrc.data alphatilt, betatilt, x_stage, y_stage, z_stage = _extract_stage_tilt_from_header(mrc) if data.ndim == 2: data = data[np.newaxis, ...] if data.shape[1:] != (frame_height, frame_width): return ( f"Frame dimensions changed unexpectedly in '{os.path.basename(file_path)}': " f"{data.shape[1:]} != {(frame_height, frame_width)}" ), None for local_idx in range(frames_in_file): new_dataset[frame_idx] = np.asarray(data[local_idx], dtype=dtype) index_dataset[frame_idx] = frame_idx stagepos_x_values[frame_idx] = x_stage stagepos_y_values[frame_idx] = y_stage stagepos_z_values[frame_idx] = z_stage alphatilt_dataset[frame_idx] = alphatilt betatilt_dataset[frame_idx] = betatilt frame_idx += 1 if frame_idx % 1000 == 0: log_print(f'{frame_idx} frames processed') if frame_idx % save_interval == 0: new_file.flush() log_start(logfile_path, f"Saved after processing {frame_idx} frames") if np.isfinite(stagepos_x_values).any(): new_file.create_dataset(stagepos_x_path, data=stagepos_x_values, dtype='float32') new_file.create_dataset(stagepos_x_refined_path, data=stagepos_x_values, dtype='float32') else: log_start(logfile_path, "No finite stagepos_x values found; skipping stagepos_x(_refined) datasets.") if np.isfinite(stagepos_y_values).any(): new_file.create_dataset(stagepos_y_path, data=stagepos_y_values, dtype='float32') new_file.create_dataset(stagepos_y_refined_path, data=stagepos_y_values, dtype='float32') else: log_start(logfile_path, "No finite stagepos_y values found; skipping stagepos_y(_refined) datasets.") if np.isfinite(stagepos_z_values).any(): new_file.create_dataset(stagepos_z_path, data=stagepos_z_values, dtype='float32') else: log_start(logfile_path, "No finite stagepos_z values found; skipping stagepos_z dataset.") ensure_nexus_parents(new_file) ensure_image_nxdata(new_file) ensure_image_key(new_file) ensure_dense_logs(new_file) ensure_goniometer_transforms(new_file, goniometer_transform_order) write_nxprocess_import( new_file, program="coseda.importers.import_mrc", input_path=mrc_source, output_path=h5file_path, parameters={ "chunk_size": chunk_shape, "num_files": len(all_files), "num_frames": total_frames, }, ) write_beam_incident_energy(new_file, acceleration_voltage) write_detector_geometry(new_file, camera_length, camera_length_correction) log_start(logfile_path, "Dataset created and data copied") return None, h5file except Exception as e: return f"Error during file conversion: {str(e)}", None
[docs] def mrc_folder_conversion_batch(input_path, continue_on_error=False, return_errors=False): """ Convert one or many MRC import INIs to H5. Parameters ---------- input_path : str | list[str] Path(s) accepted by `handle_input`. continue_on_error : bool If True, process all config files and collect failures. If False, stop at first failure (legacy behavior). return_errors : bool If True, return `(success: bool, errors: list[str])`. If False, return only `success: bool` (legacy behavior). """ # Handle input path and configuration files configfiles, input_path = handle_input(input_path) all_successful = True error_messages = [] # Loop through all config files found in the input path for configfile in configfiles: logfile_path = None try: shoutout(configfile) # Parse the configuration file config_data = parse_config(configfile) if config_data is None: all_successful = False error_message = f"Failed to parse config file: {configfile}" error_messages.append(error_message) log_print(error_message) if not continue_on_error: break continue (config, outputfolder, originalfile, logfile, path, outputfolder_path, originalfile_path, logfile_path, framepath, _, _) = config_data chunk_size = 1000 # Define chunk size (modify if needed) log_start(logfile_path, f'Attempting to convert {os.path.basename(configfile)} with chunked frame stack') # Call the mrc_folder_conversion function try: acceleration_voltage = float(config.get('AcquisitionDetails', 'acceleration_voltage')) except Exception: acceleration_voltage = None try: camera_length = float(config.get('AcquisitionDetails', 'camera_length')) except Exception: camera_length = None try: camera_length_correction = float(config.get('AcquisitionDetails', 'camera_length_correction')) except Exception: camera_length_correction = 1.0 goniometer_transform_order = get_goniometer_transform_order(config) result, h5file = mrc_folder_conversion( originalfile_path, outputfolder, logfile_path, chunk_size, acceleration_voltage=acceleration_voltage, camera_length=camera_length, camera_length_correction=camera_length_correction, goniometer_transform_order=goniometer_transform_order, ) log_print(h5file) if result is None: # Conversion was successful config.set('Paths', 'h5file', h5file) # Save changes to the config file with open(configfile, 'w') as configfile_output: config.write(configfile_output) # Log the success log_start(logfile_path, 'Conversion completed without errors.') else: # Conversion failed, log the error log_result(logfile_path, 'Conversion failed', result) all_successful = False error_messages.append(f"{configfile}: {result}") if not continue_on_error: break except Exception as e: # Handle any unexpected exceptions error_message = f"Unexpected error: {str(e)}" if logfile_path: log_result(logfile_path, 'Conversion failed', error_message) else: log_print(error_message) all_successful = False error_messages.append(f"{configfile}: {error_message}") if not continue_on_error: break if return_errors: return all_successful, error_messages return all_successful