基于Hyperspy的EELS谱图图像逐像素拟合并行处理需求

问题背景 #

使用Hyperspy处理EELS谱图图像（SI）时，双高斯模型直接拟合一维谱图效果不佳，因此采用双层循环逐像素拟合，但24万像素处理耗时约30分钟，速度过慢。尝试Python多进程并行处理未成功，需要保留像素对应结果及元数据的可行并行方案。

现有串行代码 #

SizeX = highlosscrop_2.axes_manager[0].size
SizeY = highlosscrop_2.axes_manager[1].size

highlosscrop_2_array = np.array(highlosscrop_2)

Ni2peak = np.zeros((SizeY, SizeX))
Ni4peak = np.zeros((SizeY, SizeX))

start = time.time()

# 补全原代码缺失的外层循环
for j in range(SizeY):
    for i in range(SizeX):
        NiEQspectrum = highlosscrop_2_array[j, i, :]

        NiEQspectrums = hs.signals.Signal1D(NiEQspectrum)

        NiEQspectrums.axes_manager[0].scale = highlosscrop_2.axes_manager[2].scale
        NiEQspectrums.axes_manager[0].units = highlosscrop_2.axes_manager[2].units
        NiEQspectrums.axes_manager[0].name = highlosscrop_2.axes_manager[2].name
        NiEQspectrums.axes_manager[0].offset = highlosscrop_2.axes_manager[2].offset


        m1=NiEQspectrums.create_model()

        g1 = hs.model.components1D.Gaussian() 
        g2 = hs.model.components1D.Gaussian() 

        # g1 fitting value（建议补充初始值加速拟合）
        # g2 fitting value（建议补充初始值加速拟合）

        m1.append(g1)
        m1.append(g2)


        m1.fit_component(g1, signal_range=(850., 856.), only_current=False)
        m1.fit_component(g2, signal_range=(854., 860.), only_current=False)


        m1.set_signal_range(845., 865.)
        m1.multifit()


        m1.reset_signal_range()
        # m1.plot(plot_components=True)

        #g1.print_current_values(fancy=True)
        #g2.print_current_values(fancy=True)
        
        Ni2peak[j,i] = g1.A.value
        Ni4peak[j,i] = g2.A.value
        
        
        del g1
        del g2

核心需求 #

• 保留Ni2peak、Ni4peak与原图像XY像素的一一对应关系
• 结果数据需继承原谱图的scale、offset、unit、name等元数据
• 实现并行处理，显著提升拟合速度

可行并行方案 #

方案1：Hyperspy内置并行拟合（推荐） #

Hyperspy原生支持谱图图像的模型并行拟合，无需手动编写循环，底层自动处理并行逻辑：

import hyperspy.api as hs

# 直接基于原信号创建全局模型
model = highlosscrop_2.create_model()

# 先通过平均谱获取拟合初始参数，大幅提升速度与稳定性
avg_spectrum = highlosscrop_2.mean((0,1))
avg_model = avg_spectrum.create_model()
avg_g1 = hs.model.components1D.Gaussian()
avg_g2 = hs.model.components1D.Gaussian()
avg_model.append(avg_g1)
avg_model.append(avg_g2)
avg_model.fit_component(avg_g1, signal_range=(850., 856.))
avg_model.fit_component(avg_g2, signal_range=(854., 860.))
avg_model.set_signal_range(845., 865.)
avg_model.fit()

# 初始化双高斯组件
g1 = hs.model.components1D.Gaussian()
g2 = hs.model.components1D.Gaussian()
# 用平均谱的拟合结果作为所有像素的初始值
g1.A.value = avg_g1.A.value
g1.centre.value = avg_g1.centre.value
g1.sigma.value = avg_g1.sigma.value
g2.A.value = avg_g2.A.value
g2.centre.value = avg_g2.centre.value
g2.sigma.value = avg_g2.sigma.value

# 添加组件到全局模型
model.append(g1)
model.append(g2)

# 启动并行拟合，n_jobs=-1表示使用全部CPU核心
model.multifit(n_jobs=-1, signal_range=(845., 865.))

# 提取每个像素的高斯振幅，直接生成带元数据的Signal2D
Ni2peak = model.components.Gaussian_0.A.as_signal2D((0,1))
Ni4peak = model.components.Gaussian_1.A.as_signal2D((0,1))

# 同步原谱图的元数据
Ni2peak.axes_manager.copy_axes_from(highlosscrop_2.axes_manager)
Ni4peak.axes_manager.copy_axes_from(highlosscrop_2.axes_manager)

方案2：手动用multiprocessing实现并行 #

若内置并行不满足需求，可手动基于进程池实现并行处理：

import multiprocessing as mp
import numpy as np
import hyperspy.api as hs

# 定义单个像素的拟合函数
def fit_single_pixel(args):
    spectrum_data, scale, units, name, offset = args
    # 创建信号并配置元数据
    s = hs.signals.Signal1D(spectrum_data)
    s.axes_manager[0].scale = scale
    s.axes_manager[0].units = units
    s.axes_manager[0].name = name
    s.axes_manager[0].offset = offset
    
    # 创建模型并添加组件（建议提前设置初始值）
    m = s.create_model()
    g1 = hs.model.components1D.Gaussian()
    g2 = hs.model.components1D.Gaussian()
    m.append(g1)
    m.append(g2)
    
    m.fit_component(g1, signal_range=(850., 856.), only_current=False)
    m.fit_component(g2, signal_range=(854., 860.), only_current=False)
    m.set_signal_range(845., 865.)
    m.multifit()
    
    return g1.A.value, g2.A.value

if __name__ == "__main__":
    # 提取原谱图的轴信息
    axis_info = highlosscrop_2.axes_manager[2]
    scale = axis_info.scale
    units = axis_info.units
    name = axis_info.name
    offset = axis_info.offset
    
    # 准备所有像素的拟合参数列表
    SizeX = highlosscrop_2.axes_manager[0].size
    SizeY = highlosscrop_2.axes_manager[1].size
    highloss_array = np.array(highlosscrop_2)
    args_list = []
    for j in range(SizeY):
        for i in range(SizeX):
            args_list.append((highloss_array[j,i,:], scale, units, name, offset))
    
    # 创建进程池并执行并行拟合
    pool = mp.Pool(processes=mp.cpu_count())
    results = pool.map(fit_single_pixel, args_list)
    pool.close()
    pool.join()
    
    # 将结果整理为二维数组
    Ni2peak = np.zeros((SizeY, SizeX))
    Ni4peak = np.zeros((SizeY, SizeX))
    idx = 0
    for j in range(SizeY):
        for i in range(SizeX):
            Ni2peak[j,i], Ni4peak[j,i] = results[idx]
            idx += 1
    
    # 转换为带元数据的Hyperspy信号
    Ni2peak_signal = hs.signals.Signal2D(Ni2peak)
    Ni4peak_signal = hs.signals.Signal2D(Ni4peak)
    Ni2peak_signal.axes_manager.copy_axes_from(highlosscrop_2.axes_manager)
    Ni4peak_signal.axes_manager.copy_axes_from(highlosscrop_2.axes_manager)

额外优化建议 #

• 必设初始参数：先对平均谱拟合得到初始参数，再应用到所有像素，可大幅降低拟合时间，提升稳定性
• 缩小拟合范围：始终限制在目标能量区间（845-865）内拟合，避免全谱计算浪费资源
• 减少重复初始化：串行代码中每次循环创建Signal1D和模型的开销极大，内置模型方案可复用组件，减少冗余操作

内容的提问来源于stack exchange，提问作者Ruiky94