Mahotas库:生物医学图像处理的Python利器

pip install mahotas

pip install mahotas numpy scipy matplotlib

import mahotas as mhfrom matplotlib import pyplot as plt# 读取图像image = mh.imread('cell_image.jpg')# 转换为灰度图gray = mh.colors.rgb2gray(image)# Otsu阈值分割T_otsu = mh.thresholding.otsu(gray)binary = gray > T_otsu# 显示结果plt.figure(figsize=(12, 4))plt.subplot(131)plt.imshow(gray, cmap='gray')plt.title('原始图像')plt.subplot(132)plt.imshow(binary, cmap='gray')plt.title(f'Otsu阈值={T_otsu:.2f}')plt.show()

# RC阈值分割T_rc = mh.thresholding.rc(gray)binary_rc = gray > T_rc

import numpy as npimport mahotas as mh# 读取图像image = mh.imread('touching_cells.jpg')gray = mh.colors.rgb2gray(image)# 使用Otsu阈值生成初始分割T_otsu = mh.thresholding.otsu(gray)seeds, _ = mh.label(gray > T_otsu)# 应用分水岭算法labeled = mh.cwatershed(gray.max() - gray, seeds)# 显示分割结果plt.figure(figsize=(10, 5))plt.subplot(121)plt.imshow(gray, cmap='gray')plt.title('原始图像')plt.subplot(122)plt.imshow(labeled, cmap='jet')plt.title('分水岭分割结果')plt.show()

# SLIC超像素分割labels = mh.segmentation.slic(gray, n_segments=100)plt.figure(figsize=(10, 5))plt.subplot(121)plt.imshow(gray, cmap='gray')plt.title('原始图像')plt.subplot(122)plt.imshow(mh.mark_boundaries(gray, labels))plt.title('SLIC超像素分割')plt.show()

# 计算Haralick纹理特征def compute_haralick_features(image):    """计算Haralick纹理特征"""    # 确保图像是uint8类型    if image.dtype != np.uint8:        image = (image * 255).astype(np.uint8)    # 计算Haralick特征（4个方向）    features = mh.features.haralick(image)    # 平均化（去除方向依赖性）    mean_features = features.mean(axis=0)    # 特征名称    feature_names = [        '对比度（Contrast）',        '相关性（Correlation）',        '差异性（Dissimilarity）',        '能量（Energy）',        '同质性（Homogeneity）',        'ASM（Angular Second Moment）',        '熵（Entropy）'    ]    return mean_features, feature_names# 使用示例gray = mh.colors.rgb2gray(mh.imread('tissue_image.jpg'))features, names = compute_haralick_features(gray)# 打印前7个主要特征for i, (name, value) in enumerate(zip(names, features[:7])):    print(f'{name}: {value:.4f}')

# 计算Zernike矩特征def compute_zernike_moments(image, degree=8):    """计算Zernike矩"""    # 二值化图像    T = mh.thresholding.otsu(image)    binary = image > T    # 计算Zernike矩    moments = mh.features.zernike_moments(binary, degree, radius=image.shape[0]//2)    return moments# 使用示例gray = mh.colors.rgb2gray(mh.imread('cell_shape.jpg'))zernike_moments = compute_zernike_moments(gray, degree=12)print(f"Zernike矩特征向量（{len(zernike_moments)}维）:")print(zernike_moments)

# 计算LBP特征def compute_lbp_features(image):    """计算LBP特征"""    if image.dtype != np.uint8:        image = (image * 255).astype(np.uint8)    # 计算LBP    lbp = mh.features.lbp(image, radius=3, points=24)    # 计算直方图    hist, _ = np.histogram(lbp, bins=np.arange(0, 27), range=(0, 26))    # 归一化    hist = hist.astype(float)    hist /= (hist.sum() + 1e-7)    return hist# 使用示例gray = mh.colors.rgb2gray(mh.imread('texture_image.jpg'))lbp_hist = compute_lbp_features(gray)print(f"LBP特征直方图（{len(lbp_hist)}维）")

# 计算TAS特征def compute_tas_features(image):    """计算TAS特征"""    if image.dtype != np.uint8:        image = (image * 255).astype(np.uint8)    # 计算TAS特征    tas = mh.features.tas(image)    return tas# 使用示例gray = mh.colors.rgb2gray(mh.imread('image.jpg'))tas_features = compute_tas_features(gray)print(f"TAS特征: {tas_features}")

# Sobel边缘检测def sobel_edge_detection(image):    """Sobel边缘检测"""    if image.dtype != np.uint8:        image = (image * 255).astype(np.uint8)    # 计算Sobel边缘    edges = mh.sobel(image)    return edges# 使用示例gray = mh.colors.rgb2gray(mh.imread('cell_image.jpg'))edges = sobel_edge_detection(gray)plt.figure(figsize=(12, 4))plt.subplot(131)plt.imshow(gray, cmap='gray')plt.title('原始图像')plt.subplot(132)plt.imshow(edges, cmap='gray')plt.title('Sobel边缘检测')plt.show()

# 距离变换def distance_transform(binary_image):    """计算距离变换"""    # 确保输入是二值图像    binary = binary_image.astype(bool)    # 计算距离变换    dmap = mh.distance(binary)    return dmap# 使用示例binary = mh.imread('binary_shape.jpg', as_grey=True) > 128dmap = distance_transform(binary)plt.figure(figsize=(12, 4))plt.subplot(131)plt.imshow(binary, cmap='gray')plt.title('二值图像')plt.subplot(132)plt.imshow(dmap, cmap='hot')plt.title('距离变换')plt.show()

# 形态学操作import mahotas as mh# 腐蚀操作eroded = mh.erode(binary, np.ones((3, 3)))# 膨胀操作dilated = mh.dilate(binary, np.ones((3, 3)))# 开运算opened = mh.open(binary, np.ones((3, 3)))# 闭运算closed = mh.close(binary, np.ones((3, 3)))plt.figure(figsize=(15, 4))plt.subplot(151)plt.imshow(binary, cmap='gray')plt.title('原始二值图像')plt.subplot(152)plt.imshow(eroded, cmap='gray')plt.title('腐蚀')plt.subplot(153)plt.imshow(dilated, cmap='gray')plt.title('膨胀')plt.subplot(154)plt.imshow(opened, cmap='gray')plt.title('开运算')plt.subplot(155)plt.imshow(closed, cmap='gray')plt.title('闭运算')plt.show()

import numpy as npimport mahotas as mhfrom matplotlib import pyplot as pltdef count_cells(image_path, min_size=100):    """    细胞计数与分割    参数:        image_path: 图像路径        min_size: 最小细胞面积（像素）    返回:        cell_count: 细胞数量        labeled: 标记图像        props: 细胞属性列表    """    # 1. 读取并预处理图像    image = mh.imread(image_path)    gray = mh.colors.rgb2gray(image)    # 2. 高斯模糊去噪    blurred = mh.gaussian_filter(gray, 2)    # 3. Otsu阈值分割    T = mh.thresholding.otsu(blurred)    binary = blurred > T    # 4. 形态学操作    # 先闭运算填充细胞内部孔洞，再开运算去除小噪点    cleaned = mh.close(binary, np.ones((5, 5)))    cleaned = mh.open(cleaned, np.ones((3, 3)))    # 5. 标记连通区域    labeled, n_objects = mh.label(cleaned)    # 6. 过滤小区域    sizes = mh.labeled.labeled_size(labeled)    mask = sizes >= min_size    # 重新标记    labeled = mh.labeled.remove_regions_where(labeled, ~mask)    labeled, cell_count = mh.labeled.relabel(labeled)    # 7. 提取细胞属性    props = []    for i in range(1, cell_count + 1):        region = labeled == i        size = region.sum()        bbox = mh.labeled.bbox(region)        centroid = mh.center_of_mass(region)        props.append({            'id': i,            'size': size,            'bbox': bbox,            'centroid': centroid        })    return cell_count, labeled, props# 使用示例cell_count, labeled, props = count_cells('cells.jpg')print(f"检测到 {cell_count} 个细胞")print("\n细胞属性:")for p in props[:5]:  # 显示前5个细胞    print(f"细胞 {p['id']}: 面积={p['size']}, 中心={p['centroid']}")# 可视化plt.figure(figsize=(15, 5))image = mh.imread('cells.jpg')gray = mh.colors.rgb2gray(image)plt.subplot(131)plt.imshow(gray, cmap='gray')plt.title('原始图像')plt.subplot(132)plt.imshow(labeled, cmap='jet')plt.title(f'细胞分割结果 (共{cell_count}个)')plt.subplot(133)plt.imshow(mh.overlay(gray, labeled))plt.title('叠加显示')plt.show()

import numpy as npimport mahotas as mhfrom sklearn.ensemble import RandomForestClassifierfrom sklearn.model_selection import train_test_splitdef extract_texture_features(image, patch_size=64, stride=32):    """    从图像中提取纹理特征    参数:        image: 输入图像        patch_size: 图像块大小        stride: 步长    返回:        features: 特征矩阵 (n_samples, n_features)        positions: 位置信息    """    height, width = image.shape    features = []    positions = []    # 滑动窗口提取特征    for y in range(0, height - patch_size + 1, stride):        for x in range(0, width - patch_size + 1, stride):            patch = image[y:y+patch_size, x:x+patch_size]            # 计算Haralick特征            haralick = mh.features.haralick(patch.astype(np.uint8))            feature_vector = haralick.mean(axis=0)  # 平均化            features.append(feature_vector)            positions.append((x, y))    return np.array(features), positionsdef classify_tissue_regions(image_path, model=None):    """    组织区域分类    参数:        image_path: 图像路径        model: 预训练模型（可选）    返回:        classified_image: 分类结果    """    # 读取图像    image = mh.imread(image_path)    gray = mh.colors.rgb2gray(image)    # 提取纹理特征    features, positions = extract_texture_features(gray)    # 如果没有提供模型，使用随机森林（实际应用中应使用预训练模型）    if model is None:        model = RandomForestClassifier(n_estimators=100, random_state=42)        # 这里应该使用已标注数据进行训练        # 示例中随机分配标签（实际应用请删除）        labels = np.random.randint(0, 3, len(features))        model.fit(features, labels)    # 预测    predictions = model.predict(features)    # 创建分类图像    classified = np.zeros(gray.shape, dtype=np.int32)    for (x, y), pred in zip(positions, predictions):        classified[y:y+64, x:x+64] = pred    return classified, model# 使用示例classified, model = classify_tissue_regions('pathology_image.jpg')plt.figure(figsize=(12, 5))image = mh.imread('pathology_image.jpg')gray = mh.colors.rgb2gray(image)plt.subplot(121)plt.imshow(gray, cmap='gray')plt.title('原始病理图像')plt.subplot(122)plt.imshow(classified, cmap='jet')plt.title('区域分类结果')plt.show()

import numpy as npimport mahotas as mhfrom sklearn.svm import SVCfrom sklearn.preprocessing import StandardScalerdef extract_zernike_features(image, degree=12):    """    提取Zernike矩特征    参数:        image: 输入图像        degree: Zernike矩的阶数    返回:        features: Zernike特征向量    """    # 二值化    T = mh.thresholding.otsu(image)    binary = image > T    # 计算Zernike矩    radius = min(image.shape) // 2    moments = mh.features.zernike_moments(binary, degree, radius=radius)    return momentsdef classify_cells_by_shape(images, labels):    """    基于形状特征对细胞进行分类    参数:        images: 图像列表        labels: 对应的标签    返回:        model: 训练好的模型        scaler: 标准化器        accuracy: 分类准确率    """    # 提取特征    features = []    for img in images:        if img.ndim == 3:            img = mh.colors.rgb2gray(img)        feat = extract_zernike_features(img)        features.append(feat)    features = np.array(features)    # 标准化    scaler = StandardScaler()    features_scaled = scaler.fit_transform(features)    # 训练SVM分类器    model = SVC(kernel='rbf', random_state=42)    model.fit(features_scaled, labels)    # 计算准确率    predictions = model.predict(features_scaled)    accuracy = np.mean(predictions == labels)    return model, scaler, accuracy# 示例使用# 假设有两组细胞图像normal_cells = [mh.imread(f'normal_{i}.jpg') for i in range(10)]abnormal_cells = [mh.imread(f'abnormal_{i}.jpg') for i in range(10)]images = normal_cells + abnormal_cellslabels = [0] * 10 + [1] * 10  # 0=正常, 1=异常# 训练模型model, scaler, accuracy = classify_cells_by_shape(images, labels)print(f"分类准确率: {accuracy*100:.2f}%")# 测试新图像test_image = mh.imread('test_cell.jpg')test_gray = mh.colors.rgb2gray(test_image)test_features = extract_zernike_features(test_gray)test_scaled = scaler.transform([test_features])prediction = model.predict(test_scaled)print(f"预测结果: {'正常'if prediction[0] == 0else'异常'}")

import numpy as npimport mahotas as mhfrom matplotlib import pyplot as pltdef preprocess_microscopy_image(image_path):    """    显微镜图像预处理    参数:        image_path: 图像路径    返回:        preprocessed: 预处理后的图像        steps: 处理步骤结果    """    # 1. 读取图像    image = mh.imread(image_path)    gray = mh.colors.rgb2gray(image)    steps = {'original': gray.copy()}    # 2. 去噪 - 高斯模糊    denoised = mh.gaussian_filter(gray, 1.5)    steps['denoised'] = denoised.copy()    # 3. 对比度增强 - 使用自适应直方图均衡化    # 简单版本：使用全局直方图均衡化    from skimage import exposure    equalized = exposure.equalize_hist(denoised)    steps['equalized'] = equalized.copy()    # 4. 背景减除（如果需要）    # 估计背景    background = mh.gaussian_filter(equalized, 50)    background_subtracted = equalized - background + 0.5    background_subtracted = np.clip(background_subtracted, 0, 1)    steps['background_subtracted'] = background_subtracted.copy()    # 5. 锐化    from scipy.ndimage import gaussian_filter    blurred = gaussian_filter(background_subtracted, sigma=1)    sharpened = background_subtracted + (background_subtracted - blurred) * 0.5    sharpened = np.clip(sharpened, 0, 1)    steps['sharpened'] = sharpened.copy()    return sharpened, steps# 使用示例preprocessed, steps = preprocess_microscopy_image('microscopy.jpg')# 可视化处理步骤fig, axes = plt.subplots(2, 3, figsize=(15, 10))axes = axes.ravel()for i, (name, img) in enumerate(steps.items()):    axes[i].imshow(img, cmap='gray')    axes[i].set_title(name)    axes[i].axis('off')axes[5].imshow(preprocessed, cmap='gray')axes[5].set_title('最终结果')axes[5].axis('off')plt.tight_layout()plt.show()

import osimport numpy as npimport mahotas as mhimport pandas as pddef batch_process_images(input_dir, output_dir, processing_func):    """    批量处理图像    参数:        input_dir: 输入目录        output_dir: 输出目录        processing_func: 处理函数    """    # 创建输出目录    os.makedirs(output_dir, exist_ok=True)    results = []    # 遍历输入目录    for filename in os.listdir(input_dir):        if filename.lower().endswith(('.jpg', '.png', '.tif', '.tiff')):            input_path = os.path.join(input_dir, filename)            output_path = os.path.join(output_dir, filename)            try:                # 处理图像                result = processing_func(input_path)                # 保存结果                if isinstance(result, np.ndarray):                    mh.imsave(output_path, result)                results.append({                    'filename': filename,                    'status': 'success',                    'output_path': output_path                })                print(f"处理完成: {filename}")            except Exception as e:                results.append({                    'filename': filename,                    'status': f'error: {str(e)}',                    'output_path': None                })                print(f"处理失败 {filename}: {str(e)}")    # 保存结果报告    df = pd.DataFrame(results)    df.to_csv(os.path.join(output_dir, 'processing_report.csv'), index=False)    return df# 示例处理函数def count_cells_in_image(image_path):    """统计细胞数量"""    image = mh.imread(image_path)    gray = mh.colors.rgb2gray(image)    T = mh.thresholding.otsu(gray)    binary = gray > T    labeled, n_cells = mh.label(binary)    # 在图像上标注细胞数量    # ...（添加标注代码）    return labeled# 使用示例report = batch_process_images(    input_dir='raw_images/',    output_dir='processed_images/',    processing_func=count_cells_in_image)print("\n处理报告:")print(report)

# 不好的做法：使用float64image = mh.imread('large_image.tif').astype(np.float64)# 好的做法：使用uint8image = mh.imread('large_image.tif').astype(np.uint8)# Haralick特征需要uint8features = mh.features.haralick(image)

from concurrent.futures import ProcessPoolExecutorimport numpy as npdef process_single_image(image_path):    """处理单张图像"""    image = mh.imread(image_path)    # ... 处理逻辑    return resultdef parallel_process_images(image_paths, n_workers=4):    """并行处理多张图像"""    with ProcessPoolExecutor(max_workers=n_workers) as executor:        results = list(executor.map(process_single_image, image_paths))    return results# 使用image_paths = ['img1.jpg', 'img2.jpg', 'img3.jpg', 'img4.jpg']results = parallel_process_images(image_paths)

# 不好的做法：多次转换gray = mh.colors.rgb2gray(image)uint8_img = (gray * 255).astype(np.uint8)features = mh.features.haralick(uint8_img)# 好的做法：一次性转换gray = mh.colors.rgb2gray(image).astype(np.uint8)features = mh.features.haralick(gray)

import mahotas as mhimport numpy as npdef prepare_for_deep_learning(image_path, target_size=(256, 256)):    """    为深度学习模型准备数据    参数:        image_path: 图像路径        target_size: 目标尺寸    返回:        processed_image: 处理后的图像        features: 提取的传统特征（可选）    """    # 读取图像    image = mh.imread(image_path)    # 转换为灰度    if image.ndim == 3:        gray = mh.colors.rgb2gray(image)    else:        gray = image    # 调整大小    from skimage.transform import resize    resized = resize(gray, target_size, anti_aliasing=True)    # 归一化到[0, 1]    normalized = (resized - resized.min()) / (resized.max() - resized.min())    # 可选：提取传统特征作为额外输入    features = mh.features.haralick((normalized * 255).astype(np.uint8))    features_mean = features.mean(axis=0)    return normalized, features_mean# 使用示例image, features = prepare_for_deep_learning('cell_image.jpg')# image可以直接输入深度学习模型# features可以作为额外特征拼接combined_input = np.concatenate([image.flatten(), features])

# 解决方案：使用国内镜像pip install mahotas -i https://pypi.tuna.tsinghua.edu.cn/simple# 或使用condaconda install -c conda-forge mahotas

# 检查依赖import numpyimport scipyprint(f"NumPy version: {numpy.__version__}")print(f"SciPy version: {scipy.__version__}")# 重新安装pip install --upgrade mahotas numpy scipy

# 解决方案：分块处理def process_large_image(image_path, block_size=1024):    """分块处理大图像"""    # 获取图像尺寸    import PIL.Image as PILImage    with PILImage.open(image_path) as img:        width, height = img.size    # 分块处理    results = []    for y in range(0, height, block_size):        for x in range(0, width, block_size):            block = mh.imread(image_path, as_grey=True)            block = block[y:y+block_size, x:x+block_size]            # 处理块            processed = process_block(block)            results.append(processed)    # 合并结果    return merge_results(results, width, height)

# 解决方案1：使用适当的数据类型image = image.astype(np.uint8)  # 而非float64# 解决方案2：减少不必要的计算# 预先缓存结果T = mh.thresholding.otsu(gray)  # 计算一次binary = gray > T  # 多次使用# 解决方案3：使用并行处理from joblib import Parallel, delayedresults = Parallel(n_jobs=4)(delayed(process)(img) for img in images)

# 解决方案：尝试多种阈值方法def auto_threshold(gray):    """自动选择最佳阈值方法"""    # Otsu    T_otsu = mh.thresholding.otsu(gray)    # Riddler-Calvard    T_rc = mh.thresholding.rc(gray)    # 选择效果更好的方法    # 可以基于图像质量指标评估    return T_otsu  # 或T_rc# 或者使用自适应阈值from skimage.filters import threshold_localadaptive_thresh = threshold_local(gray, block_size=35)

# 数据处理import numpy as npimport pandas as pd# 图像处理import mahotas as mhimport cv2from skimage import io, filters, morphology# 机器学习from sklearn import svm, ensembleimport tensorflow as tf# 可视化import matplotlib.pyplot as pltimport seaborn as sns# 科学计算from scipy import ndimage

project/├── data/│   ├── raw/│   ├── processed/│   └── results/├── src/│   ├── preprocessing.py│   ├── segmentation.py│   ├── feature_extraction.py│   └── utils.py├── notebooks/│   └── analysis.ipynb├── tests/│   └── test_mahotas.py└── requirements.txt

# 好的函数设计def analyze_cells(image_path, params=None):    """    分析细胞图像    参数:        image_path: 图像路径        params: 可选参数字典            - min_size: 最小细胞面积            - threshold_method: 阈值方法            - output_dir: 输出目录    返回:        dict: 包含分析结果的字典            - cell_count: 细胞数量            - features: 特征数组            - labeled_image: 标记图像    """    # 设置默认参数    if params is None:        params = {}    min_size = params.get('min_size', 100)    threshold_method = params.get('threshold_method', 'otsu')    # 实现    # ...    return results

def validate_results(segmented_image, ground_truth=None):    """    验证分割结果    参数:        segmented_image: 分割结果        ground_truth: 真值（可选）    返回:        dict: 验证指标    """    if ground_truth is not None:        # 计算准确率、召回率等指标        from sklearn.metrics import accuracy_score        accuracy = accuracy_score(ground_truth.flatten(),                                   segmented_image.flatten())    else:        # 无真值时的启发式评估        # 例如：检查连通区域的合理性        n_regions = mh.labeled.labeled_size(segmented_image > 0)        print(f"检测到 {n_regions} 个区域")    return validation_metrics

# 常用函数速查# 图像读取image = mh.imread('image.jpg')gray = mh.colors.rgb2gray(image)# 阈值分割T = mh.thresholding.otsu(gray)T_rc = mh.thresholding.rc(gray)# 特征提取haralick = mh.features.haralick(gray.astype(np.uint8))zernike = mh.features.zernike_moments(binary, degree=12)lbp = mh.features.lbp(gray.astype(np.uint8), radius=3, points=24)# 分割labeled, n_objects = mh.label(binary)watershed = mh.cwatershed(gray.max() - gray, seeds)# 形态学eroded = mh.erode(binary, np.ones((3,3)))dilated = mh.dilate(binary, np.ones((3,3)))# 边缘检测edges = mh.sobel(gray)