【Python】MiniAn 库

Octogen_hhhhmx

已于 2024-10-10 13:11:21 修改

阅读量506

点赞数 7

CC 4.0 BY-SA版权

分类专栏： python 文章标签： python 开发语言

于 2024-09-19 12:27:43 首次发布

本文链接：https://blog.youkuaiyun.com/hmxJim/article/details/142344643

python 专栏收录该内容

1 篇文章

订阅专栏

1 minian 库介绍

1.1 用途

显微数据处理，图像读取，可视化。

1.2 参考

1.2.1 Github 地址

Github: denisecailab / minian

1.2.2 安装教程

Installation Tutorial

2 虚拟环境配置

需在命令行或终端完成。

需要注意的是，官方 minian 库有千奇百怪的 bug。所以最好下载调好的包直接导入。

2.1 创建环境

conda create -y -n minian

这里也可以直接将配置好的包导入，跳过 2.2 与 2.3。
先将调好的环境下载。

minian.yaml

下载后，要将 conda 导入模式调为兼容。

conda config --set channel_priority flexible

然后直接导入包。值得注意的是，需要在 minian.yaml 所在的目录打开命令行。

conda env create -f minian.yaml -n minian

2.2 激活

conda activate minian

2.3 导入包

conda install -y -c conda-forge minian

2.4 Jupyter NoteBook

安装。

minian-install --notebooks

使用。

jupyter notebook

导入环境需要其他工具。

pip install --user ipykernel

将虚拟环境导入 Jupyter。

python -m ipykernel install --user --name=minian

查看 Jupyter 已有的环境。

jupyter kernelspec list

取消 Jupyter 中链接 conda 的环境。

jupyter kernelspec uninstall minian

创建批处理（bat）文件，将 jupyter notebook 添加到右键菜单。需要修改以下的 conda 安装路径，且要有管理员身份。

@echo off
setlocal

REM 创建注册表项，添加默认值
reg add "HKEY_CLASSES_ROOT\Directory\Background\shell\jupyter" /ve /t REG_SZ /d "Jupyter Notebook" /f

REM 添加图标
reg add "HKEY_CLASSES_ROOT\Directory\Background\shell\jupyter" /v "Icon" /t REG_SZ /d "D:\anaconda3\envs\minian\Menu\jupyter.ico" /f

REM 添加命令
reg add "HKEY_CLASSES_ROOT\Directory\Background\shell\jupyter\command" /ve /t REG_SZ /d "D:\anaconda3\python.exe D:\anaconda3\Scripts\jupyter-notebook-script.py %V" /f

endlocal

2.5 辅助工具

2.5.1 下载案例

minian-install --demo

2.5.2 快速跳转插件

可以借此创建节标题，快速跳转。

conda install -c conda-forge jupyter_contrib_nbextensions

2.5.3 加快导入

可以使用mamba加快包的导入。

conda install mamba

导入mamba后就可用mamba来导入其他包。

mamba install numpy

3 使用

注意，MiniAn 常在Jupyter Notebook 中运行。

3.1 参数详解

MiniAn documentation

3.2 准备工作

3.2.1 导入库

%%capture
%load_ext autoreload
%autoreload 2
# 导入所有库
import itertools as itt
import os
import sys

import holoviews as hv
import numpy as np
import xarray as xr
from dask.distributed import Client, LocalCluster
from holoviews.operation.datashader import datashade, regrid
from holoviews.util import Dynamic
from IPython.core.display import display

3.2.2 设置路径

# minian代码库位置
minian_path = "."
# 视频文件夹
dpath = "./demo_movies/"
minian_ds_path = os.path.join(dpath, "minian")
intpath = "./minian_intermediate"
subset = dict(frame=slice(0, None))
subset_mc = None
# 是否显示交互式绘图
interactive = True
# 控制所有图的相对大小（可>100）
output_size = 100
# 线程数
n_workers = int(os.getenv("MINIAN_NWORKERS", 4))
# 保存数据的目标文档夹和格式
param_save_minian = {
    # 保存数据的文档夹路径
    "dpath": minian_ds_path,
    # 为最终的标记数据结构构建元数据
    "meta_dict": dict(session=-1, animal=-2),
    # 是否覆盖数据
    "overwrite": True,
}

# 视频参数
param_load_videos = {
    "pattern": "msCam[0-9]+\.avi$",
    "dtype": np.uint8,
    "downsample": dict(frame=1, height=1, width=1),
    "downsample_strategy": "subset",
}
# 降噪参数，最好看过效果后手动调试
param_denoise = {"method": "median", "ksize": 7}
# 去除背景参数，最好看过效果后手动调试
param_background_removal = {"method": "tophat", "wnd": 15}

# Motion Correction Parameters#
subset_mc = None
param_estimate_motion = {"dim": "frame"}

# Initialization Parameters#
param_seeds_init = {
    "wnd_size": 1000,
    "method": "rolling",
    "stp_size": 500,
    "max_wnd": 15,
    "diff_thres": 3,
}
param_pnr_refine = {"noise_freq": 0.06, "thres": 1}
param_ks_refine = {"sig": 0.05}
param_seeds_merge = {"thres_dist": 10, "thres_corr": 0.8, "noise_freq": 0.06}
param_initialize = {"thres_corr": 0.8, "wnd": 10, "noise_freq": 0.06}
param_init_merge = {"thres_corr": 0.8}

# CNMF Parameters#
param_get_noise = {"noise_range": (0.06, 0.5)}
param_first_spatial = {
    "dl_wnd": 10,
    "sparse_penal": 0.01,
    "size_thres": (25, None),
}
param_first_temporal = {
    "noise_freq": 0.06,
    "sparse_penal": 1,
    "p": 1,
    "add_lag": 20,
    "jac_thres": 0.2,
}
param_first_merge = {"thres_corr": 0.8}
param_second_spatial = {
    "dl_wnd": 10,
    "sparse_penal": 0.01,
    "size_thres": (25, None),
}
param_second_temporal = {
    "noise_freq": 0.06,
    "sparse_penal": 1,
    "p": 1,
    "add_lag": 20,
    "jac_thres": 0.4,
}

os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MINIAN_INTERMEDIATE"] = intpath

3.2.3 导入MiniAn

%%capture
sys.path.append(minian_path)
from minian.cnmf import (
    compute_AtC,
    compute_trace,
    get_noise_fft,
    smooth_sig,
    unit_merge,
    update_spatial,
    update_temporal,
    update_background,
)
from minian.initialization import (
    gmm_refine,
    initA,
    initC,
    intensity_refine,
    ks_refine,
    pnr_refine,
    seeds_init,
    seeds_merge,
)
from minian.motion_correction import apply_transform, estimate_motion
from minian.preprocessing import denoise, remove_background
from minian.utilities import (
    TaskAnnotation,
    get_optimal_chk,
    load_videos,
    open_minian,
    save_minian,
)
from minian.visualization import (
    CNMFViewer,
    VArrayViewer,
    generate_videos,
    visualize_gmm_fit,
    visualize_motion,
    visualize_preprocess,
    visualize_seeds,
    visualize_spatial_update,
    visualize_temporal_update,
    write_video,
)
dpath = os.path.abspath(dpath)
hv.notebook_extension("bokeh", width=100)
cluster = LocalCluster(
    n_workers=n_workers,
    # 分配内存
    memory_limit="2GB",
    resources={"MEM": 1},
    # 线程
    threads_per_worker=2,
    # 进程PID
    dashboard_address=":8787",
)
annt_plugin = TaskAnnotation()
cluster.scheduler.add_plugin(annt_plugin)
client = Client(cluster)

3.3 预处理

3.3.1 导入视频

# 查看视频导入配置
param_load_videos
# 导入视频
varr = load_videos(dpath, **param_load_videos)
# 确定合适分块数
chk, _ = get_optimal_chk(varr, dtype=float)
# 储存导入好的视频数据
varr = save_minian(
    varr.chunk({"frame": chk["frame"], "height": -1, "width": -1}).rename("varr"),
    intpath,
    overwrite=True,
)
# 查看储存结果
varr

值得注意的是，如果在 3.2.3 中 cluster 的内存限制 memory_limit 设置得太小，那么这里会报错 KilledWorker。

3.3.2 可视化结果

hv.output(size=output_size)
if interactive:
    # 可以使用4种方法："mean", "min", "max", "diff"
    vaviewer = VArrayViewer(varr, framerate=5, summary=["mean", "max"])
    display(vaviewer.show())

# 设置蒙版
if interactive:
    try:
        subset_mc = list(vaviewer.mask.values())[0]
    except IndexError:
        pass

# 检查视频问题
varr.sel(height=slice(0, 99), width=slice(0, 99))

3.3.3 去除背景光

# 选取样例
varr_ref = varr.sel(subset)
# 去除背景光
varr_min = varr_ref.min("frame").compute()
varr_ref = varr_ref - varr_min
# 查看效果
hv.output(size=int(output_size * 0.7))
if interactive:
    vaviewer = VArrayViewer(
        [varr.rename("original"), varr_ref.rename("glow_removed")],
        framerate=5,
        summary=None,
        layout=True,
    )
    display(vaviewer.show())

3.3.4 降噪

# 查看降噪参数
param_denoise
# 查看不同级别降噪效果
hv.output(size=int(output_size * 0.6))
if interactive:
    display(
        visualize_preprocess(
            varr_ref.isel(frame=0).compute(),
            denoise,
            method=["median"],
            ksize=[5, 7, 9],
        )
    )

确定最好效果后，开始降噪。

param_denoise['ksize'] = 5  # 需要自行修改
varr_ref = denoise(varr_ref, **param_denoise)

3.3.5 去除背景

# 查看去除背景参数
param_background_removal
# 查看不同级别去背景效果
hv.output(size=int(output_size * 0.6))
if interactive:
    display(
        visualize_preprocess(
            varr_ref.isel(frame=0).compute(),
            remove_background,
            method=["tophat"],
            wnd=[10, 15, 20],
        )
    )

开始去背景。

param_background_removal['wnd'] = 15  # 需要自行修改
varr_ref = remove_background(varr_ref, **param_background_removal)

3.3.6 保存结果

值得注意的是，3.3.1 中我们已经保存了原视频导入结果，并将结果重命名为 varr。如果这里仍然用 varr，则会覆盖之前的结果。应命名为 varr_ref 较为稳妥。

%%time
varr_ref = save_minian(varr_ref.rename("varr_ref"), dpath=intpath, overwrite=True)

3.4 动作校正

3.4.1 动作估计

查看参数。

param_estimate_motion
param_save_minian

开始估计动作并保存。

%%time
motion = estimate_motion(varr_ref.sel(subset_mc), **param_estimate_motion)

motion = save_minian(
    motion.rename("motion").chunk({"frame": chk["frame"]}), **param_save_minian
)

3.4.2 动作可视化

hv.output(size=output_size)
visualize_motion(motion)

3.4.3 应用更改并保存

%%time
Y = apply_transform(varr_ref, motion, fill=0)

Y_fm_chk = save_minian(Y.astype(float).rename("Y_fm_chk"), intpath, overwrite=True)
Y_hw_chk = save_minian(
    Y_fm_chk.rename("Y_hw_chk"),
    intpath,
    overwrite=True,
    chunks={"frame": -1, "height": chk["height"], "width": chk["width"]},
)

3.4.4 校正结果可视化

hv.output(size=int(output_size * 0.7))
if interactive:
    vaviewer = VArrayViewer(
        [varr_ref.rename("before_mc"), Y_fm_chk.rename("after_mc")],
        framerate=5,
        summary=None,
        layout=True,
    )
    display(vaviewer.show())

观察最大投影。

im_opts = dict(
    frame_width=500,
    aspect=varr_ref.sizes["width"] / varr_ref.sizes["height"],
    cmap="Viridis",
    colorbar=True,
)
(
    regrid(
        hv.Image(
            varr_ref.max("frame").compute().astype(np.float32),
            ["width", "height"],
            label="before_mc",
        ).opts(**im_opts)
    )
    + regrid(
        hv.Image(
            Y_hw_chk.max("frame").compute().astype(np.float32),
            ["width", "height"],
            label="after_mc",
        ).opts(**im_opts)
    )
)

3.4.5 生成校正后视频

%%time
vid_arr = xr.concat([varr_ref, Y_fm_chk], "width").chunk({"width": -1})
write_video(vid_arr, "minian_mc.mp4", dpath)