commit b8a8ff2bc662c89e102f22a3fefcd154b6bc32cd Author: 刘航宇 <3364451258@qq.com> Date: Wed May 6 19:41:26 2026 +0800 feat: cDNA微阵列图像处理作业 - Python实现 实现内容： - 网格划分：投影分析 + 自相关估周期 + 白顶帽去背景 + 质心提取 - 三种阈值分割：人工阈值、Otsu自动阈值、迭代阈值 - TV去噪（Chambolle投影算法） - 后处理：去小连通域 + 保留最大连通域 - 完整可视化：网格叠加、阈值对比、收敛曲线、分割结果 参考MATLAB代码：NewGridAndCV/demo_GriddingAndCV.m diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..aeb3680 --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +# 大型数据文件（>10MB） +*.tif +*.TIF + +# MATLAB临时文件 +*.asv + +# Python +__pycache__/ +*.pyc + +# OS +.DS_Store +Thumbs.db diff --git a/README.md b/README.md new file mode 100644 index 0000000..574af15 --- /dev/null +++ b/README.md @@ -0,0 +1,118 @@ +# cDNA 微阵列图像处理作业资料 + +## 作业概述 + +本次作业涉及 **cDNA 微阵列（基因芯片）图像处理**，主要研究两个核心问题： + +1. **网格划分（Gridding）** - 定位微阵列图像中每个点的精确位置 +2. **图像分割（Segmentation）** - 区分前景（基因点）与背景 + +--- + +## 资料清单 + +### 1. 参考论文 + +| 文件名 | 作者 | 内容简介 | +|--------|------|----------| +| `2019-3-高污染基因芯片图像的网格划分_芦碧波.pdf` | 芦碧波 | 提出针对高污染基因芯片的网格划分方法，利用图像增强、分块处理和自动阈值检测来提高鲁棒性 | +| `封面+低对比度cDNA图像分割的局部水平集方法_芦碧波.pdf` | 芦碧波、刘利群、张霄宏、林忠华 | 提出基于局部信息的水平集方法，解决低对比度cDNA图像分割问题，引入局部图像拟合能量 | +| `显微图像分割.pdf` | - | 显微图像分割相关资料（密码保护，需另存为可读版本） | + +### 2. MATLAB 代码 + +#### `NewGridAndCV/` - 网格划分与C-V分割实现 + +| 文件名 | 功能说明 | +|--------|----------| +| `demo_GriddingAndCV.m` | **主程序** - 演示网格划分和Chan-Vese分割的完整流程 | +| `GriddingAndCV.m` | 网格划分核心算法，使用自相关估计点间距 | +| `cvseg.m` | Chan-Vese 水平集分割算法 | +| `chenvese.m` | C-V 模型的另一种实现 | +| `tvdenoise.m` | TV（Total Variation）去噪算法 | +| `choice.m` | 剔除面积过小的连通区域 | +| `choosemaxobj.m` | 保留最大连通区域 | +| `contour_bw.m` | 轮廓提取与二值化 | +| `fillingholes.m` | 填充孔洞 | +| `kappa.m` | 曲率计算 | +| `maskcircle2.m` | 圆形掩膜生成 | +| `redcolorcontour.m` | 红色轮廓显示 | +| `showphi.m` | 显示水平集函数 | + +#### `cDNA图像处理实例/` - 基础示例 + +| 文件名 | 功能说明 | +|--------|----------| +| `图像处理实例.pptx` | 实例讲解PPT | +| `数据/cDNA/Demo_cdna.m` | 基础演示代码 | +| `数据/cDNA/*.tif` | 测试图像数据 | + +### 3. 图像数据 + +| 文件名 | 说明 | +|--------|------| +| `GSM16390_CH1.tif` | 通道1原始图像（~26MB） | +| `GSM16390_CH2.tif` | 通道2原始图像（~26MB） | +| `GSM16390_CH2color.tif` | 彩色合成图像（~79MB） | +| `*_small.tif` | 上述图像的缩小版本 | +| `cDNA.png` | 测试用cDNA图像 | +| `I_bw.jpg` | 二值化结果示例 | +| `I_griddingout.tif` | 网格划分输出示例 | + +--- + +## 核心算法简介 + +### 网格划分算法流程 + +``` +1. 读取图像 → 转灰度 +2. 计算行/列投影（均值） +3. 自相关分析 → 估计点间距 +4. 形态学滤波 → 去除背景 +5. 阈值分割 → 标记峰值区域 +6. 提取质心 → 确定网格点位置 +``` + +### Chan-Vese 水平集分割 + +- **核心思想**：不依赖图像梯度，基于区域统计信息分割 +- **能量函数**：$E = \mu \cdot Length(C) + \nu \cdot Area(C) + \lambda_1 \int_{inside(C)} |I - c_1|^2 + \lambda_2 \int_{outside(C)} |I - c_2|^2$ +- **优势**：对模糊边缘和低对比度图像效果好 + +### 局部水平集方法（论文贡献） + +- 引入局部图像拟合能量 $E^{LIF}$ +- 使用高斯核 $K_\sigma$ 提取局部信息 +- 能有效处理灰度不均匀的cDNA图像 + +--- + +## 运行环境 + +- **MATLAB** R2016b 或更高版本 +- **需要工具箱**： + - Image Processing Toolbox + - Signal Processing Toolbox（用于自相关计算） + +--- + +## 快速开始 + +```matlab +% 1. 进入代码目录 +cd NewGridAndCV + +% 2. 运行演示 +demo_GriddingAndCV + +% 3. 或单独运行网格划分 +GriddingAndCV +``` + +--- + +## 参考文献 + +1. 芦碧波. 高污染基因芯片图像的网格划分[J]. 河南理工大学学报(自然科学版), 2019. +2. 芦碧波, 刘利群, 张霄宏, 林忠华. 低对比度cDNA图像分割的局部水平集方法[J]. 中国图象图形学报, 2014. diff --git a/cDNA图像处理实例/图像处理实例.pptx b/cDNA图像处理实例/图像处理实例.pptx new file mode 100644 index 0000000..a5f0c46 Binary files /dev/null and b/cDNA图像处理实例/图像处理实例.pptx differ diff --git a/cDNA图像处理实例/数据/cDNA/Demo_cdna.m b/cDNA图像处理实例/数据/cDNA/Demo_cdna.m new file mode 100644 index 0000000..ee0c161 --- /dev/null +++ b/cDNA图像处理实例/数据/cDNA/Demo_cdna.m @@ -0,0 +1,14 @@ +clc +clear +close all + +im = imread('cDNA.png'); +imgray = rgb2gray(im); + +%% sum of row , sum of col +sum_col = sum(imgray); +sum_row = sum(imgray'); + +figure,imshow(imgray) +figure,plot(sum_col),title('sum of col'); +figure,plot(sum_row),title('sum of row'); diff --git a/cDNA图像处理实例/数据/cDNA/I_bw.jpg b/cDNA图像处理实例/数据/cDNA/I_bw.jpg new file mode 100644 index 0000000..302e499 Binary files /dev/null and b/cDNA图像处理实例/数据/cDNA/I_bw.jpg differ diff --git a/cDNA图像处理实例/数据/cDNA/cDNA.png b/cDNA图像处理实例/数据/cDNA/cDNA.png new file mode 100644 index 0000000..494542a Binary files /dev/null and b/cDNA图像处理实例/数据/cDNA/cDNA.png differ diff --git a/results/result_I_bw.png b/results/result_I_bw.png new file mode 100644 index 0000000..a150615 Binary files /dev/null and b/results/result_I_bw.png differ diff --git a/results/result_full_segmentation.png b/results/result_full_segmentation.png new file mode 100644 index 0000000..40657af Binary files /dev/null and b/results/result_full_segmentation.png differ diff --git a/results/result_gridding.png b/results/result_gridding.png new file mode 100644 index 0000000..bac3ff5 Binary files /dev/null and b/results/result_gridding.png differ diff --git a/results/result_gridding_overlay.png b/results/result_gridding_overlay.png new file mode 100644 index 0000000..1adb487 Binary files /dev/null and b/results/result_gridding_overlay.png differ diff --git a/results/result_iterative_convergence.png b/results/result_iterative_convergence.png new file mode 100644 index 0000000..91b9992 Binary files /dev/null and b/results/result_iterative_convergence.png differ diff --git a/results/result_threshold_compare.png b/results/result_threshold_compare.png new file mode 100644 index 0000000..25eacfa Binary files /dev/null and b/results/result_threshold_compare.png differ diff --git a/src/cDNA_segmentation.py b/src/cDNA_segmentation.py new file mode 100644 index 0000000..8d9dc37 --- /dev/null +++ b/src/cDNA_segmentation.py @@ -0,0 +1,432 @@ +""" +cDNA微阵列图像处理 - 网格划分与阈值分割 +========================================== +作业要求： +1. 分析涉及的图像处理技术 +2. 编写阈值分割代码（人工阈值 / 迭代阈值 / Otsu自动阈值） +3. 选取cDNA图像进行分割 + +参考MATLAB代码：NewGridAndCV/demo_GriddingAndCV.m, GriddingAndCV.m +""" + +import os +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import rcParams +from PIL import Image +from scipy import ndimage +from skimage import filters, morphology, color + +# 中文字体设置 +rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans'] +rcParams['axes.unicode_minus'] = False + +# 路径配置（使用脚本位置向上两级作为基准） +_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) +_BASE_DIR = os.path.dirname(_SCRIPT_DIR) # 项目根目录 +DATA_DIR = os.path.join(_BASE_DIR, 'cDNA图像处理实例', '数据', 'cDNA') +OUTPUT_DIR = os.path.join(_BASE_DIR, 'results') + + +# ============================================================ +# 第一部分：阈值分割算法 +# ============================================================ + +def manual_threshold(img_gray: np.ndarray, T: int) -> np.ndarray: + """人工固定阈值分割：灰度 > T 为前景""" + return (img_gray > T).astype(np.uint8) + + +def otsu_threshold(img_gray: np.ndarray) -> tuple[np.ndarray, float]: + """Otsu自动阈值分割：自动寻找最佳阈值""" + T = filters.threshold_otsu(img_gray) + binary = (img_gray > T).astype(np.uint8) + return binary, T + + +def iterative_threshold(img_gray: np.ndarray, tol: float = 1.0, max_iter: int = 100) -> tuple[np.ndarray, float, list]: + """ + 迭代阈值分割 + 算法流程： + 1. 初始阈值 T = 图像均值 + 2. 用T将像素分为前景和背景 + 3. 计算前景均值 μ_fg 和背景均值 μ_bg + 4. 更新 T_new = (μ_fg + μ_bg) / 2 + 5. 若 |T_new - T| < tol，停止；否则重复 + """ + T = float(np.mean(img_gray)) + history = [T] + for _ in range(max_iter): + fg = img_gray[img_gray > T] + bg = img_gray[img_gray <= T] + if len(fg) == 0 or len(bg) == 0: + break + T_new = (float(np.mean(fg)) + float(np.mean(bg))) / 2.0 + history.append(T_new) + if abs(T_new - T) < tol: + T = T_new + break + T = T_new + binary = (img_gray > T).astype(np.uint8) + return binary, T, history + + +# ============================================================ +# 第二部分：TV去噪（参考tvdenoise.m，Chambolle投影算法） +# ============================================================ + +def tv_denoise(f: np.ndarray, lam: float = 12.0, tol: float = 1e-2, max_iter: int = 200) -> np.ndarray: + """ + 全变分去噪 - Rudin-Osher-Fatemi模型 + min TV(u) + λ/2 ||f - u||² + 使用Chambolle投影算法求解 + """ + f = f.astype(np.float64) + dt = 0.25 + p1 = np.zeros_like(f) + p2 = np.zeros_like(f) + divp = np.zeros_like(f) + + for _ in range(max_iter): + lastdivp = divp.copy() + z = divp - f * lam + z1 = np.roll(z, -1, axis=1) - z + z2 = np.roll(z, -1, axis=0) - z + denom = 1.0 + dt * np.sqrt(z1**2 + z2**2) + p1 = (p1 + dt * z1) / denom + p2 = (p2 + dt * z2) / denom + divp = p1 - np.roll(p1, 1, axis=1) + p2 - np.roll(p2, 1, axis=0) + if np.max(np.abs(divp - lastdivp)) < tol: + break + + return f - divp / lam + + +# ============================================================ +# 第三部分：网格划分（参考GriddingAndCV.m） +# ============================================================ + +def estimate_spacing(profile: np.ndarray) -> int: + """ + 通过自相关分析估计点间距 + 参考MATLAB: ac = xcov(xProfile); maxima = find(s1>0 & s2<0); estPeriod = median(diff(maxima)) + """ + p = profile - np.mean(profile) + ac = np.correlate(p, p, mode='full') + ac = ac[len(p) - 1:] # 只取正半部分（lag >= 0） + + # 检测峰值：左斜率>0 且 右斜率<0（与MATLAB一致） + s1 = np.diff(ac, prepend=ac[0]) # 左斜率 + s2 = np.diff(ac, append=ac[-1]) # 右斜率 + maxima = np.where((s1 > 0) & (s2 < 0))[0] + + # 过滤掉lag=0附近的峰和过小的峰 + min_gap = len(profile) // 30 # 最小间距保护 + maxima = maxima[maxima > min_gap] + + if len(maxima) < 2: + return max(len(profile) // 16, 10) + + spacing = int(np.round(np.median(np.diff(maxima)))) + return max(spacing, 10) + + +def find_grid_centers(profile: np.ndarray, est_period: int) -> np.ndarray: + """ + 从投影曲线中提取网格中心位置 + 参考MATLAB流程: + imtophat(profile, strel('line', estPeriod, 0)) → 去背景 + graythresh → bwlabel → regionprops.Centroid → 提取质心 + """ + # 白顶帽变换：用长度为est_period的线性结构元素去除背景 + # MATLAB: seLine = strel('line', estPeriod, 0); xProfile2 = imtophat(xProfile, seLine) + # MATLAB中xProfile是1×N向量，这里用scipy的1D白顶帽实现 + profile_f = profile.astype(np.float64) + selem = np.ones(est_period) # 1D线性结构元素 + enhanced = ndimage.white_tophat(profile_f, structure=selem) + + # 自动阈值二值化（与MATLAB graythresh一致） + if enhanced.max() == 0: + return np.array([]) + T = filters.threshold_otsu(enhanced) + bw = (enhanced > T).astype(int) + + # 标记连通域并提取质心（与MATLAB bwlabel + regionprops一致） + labeled, num = ndimage.label(bw) + if num == 0: + return np.array([]) + + centers = [] + for i in range(1, num + 1): + indices = np.where(labeled == i)[0] + centers.append(float(np.mean(indices))) + return np.array(sorted(centers)) + + +def compute_grid_lines(centers: np.ndarray) -> np.ndarray: + """从中心点计算网格分割线（相邻中心的中点）""" + if len(centers) < 2: + return np.array([]) + gaps = np.diff(centers) / 2.0 + first = centers[0] - gaps[0] + last = centers[-1] + gaps[-1] + lines = [first] + for i in range(len(centers)): + if i < len(gaps): + lines.append(centers[i] + gaps[i]) + else: + lines.append(last) + return np.round(lines).astype(int) + + +def gridding(gray: np.ndarray) -> tuple: + """ + 完整网格划分流程 + 返回: x_grid, y_grid, col_profile, row_profile, x_period, y_period + """ + col_profile = np.mean(gray, axis=0).astype(np.float64) + row_profile = np.mean(gray, axis=1).astype(np.float64) + + x_period = estimate_spacing(col_profile) + y_period = estimate_spacing(row_profile) + + x_centers = find_grid_centers(col_profile, x_period) + y_centers = find_grid_centers(row_profile, y_period) + + x_grid = compute_grid_lines(x_centers) + y_grid = compute_grid_lines(y_centers) + + return x_grid, y_grid, col_profile, row_profile, x_period, y_period + + +# ============================================================ +# 第四部分：后处理（参考choice.m, choosemaxobj.m） +# ============================================================ + +def remove_small_objects(binary: np.ndarray, min_size: int = 20) -> np.ndarray: + """去除面积小于min_size的连通域""" + labeled, num = ndimage.label(binary) + result = binary.copy() + for i in range(1, num + 1): + if np.sum(labeled == i) < min_size: + result[labeled == i] = 0 + return result + + +def keep_largest_object(binary: np.ndarray, min_size: int = 20) -> np.ndarray: + """只保留最大连通域""" + labeled, num = ndimage.label(binary) + if num == 0: + return binary + areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)] + max_area = max(areas) + if max_area < min_size: + return np.zeros_like(binary) + max_idx = int(np.argmax(areas)) + 1 + return (labeled == max_idx).astype(np.uint8) + + +# ============================================================ +# 第五部分：可视化 +# ============================================================ + +def plot_threshold_comparison(gray_block: np.ndarray, block_name: str = "子块") -> plt.Figure: + """对比三种阈值方法的分割结果""" + manual_T = int(np.mean(gray_block) * 0.6) + bw_manual = manual_threshold(gray_block, manual_T) + bw_otsu, T_otsu = otsu_threshold(gray_block) + bw_iter, T_iter, _ = iterative_threshold(gray_block) + + fig, axes = plt.subplots(2, 2, figsize=(10, 10)) + fig.suptitle(f'{block_name} - 三种阈值分割方法对比', fontsize=14) + + axes[0, 0].imshow(gray_block, cmap='gray') + axes[0, 0].set_title('原始灰度图') + axes[0, 0].axis('off') + + axes[0, 1].imshow(bw_manual, cmap='gray') + axes[0, 1].set_title(f'人工阈值 (T={manual_T})') + axes[0, 1].axis('off') + + axes[1, 0].imshow(bw_otsu, cmap='gray') + axes[1, 0].set_title(f'Otsu阈值 (T={T_otsu:.1f})') + axes[1, 0].axis('off') + + axes[1, 1].imshow(bw_iter, cmap='gray') + axes[1, 1].set_title(f'迭代阈值 (T={T_iter:.1f})') + axes[1, 1].axis('off') + + plt.tight_layout() + return fig + + +def plot_gridding_result(gray: np.ndarray, x_grid: np.ndarray, y_grid: np.ndarray, + col_profile: np.ndarray, row_profile: np.ndarray) -> plt.Figure: + """绘制网格划分结果""" + fig = plt.figure(figsize=(14, 10)) + fig.suptitle('cDNA微阵列网格划分结果', fontsize=14) + + ax1 = fig.add_subplot(2, 2, 1) + ax1.imshow(gray, cmap='gray') + for x in x_grid: + ax1.axvline(x=x, color='cyan', linewidth=0.5, alpha=0.7) + for y in y_grid: + ax1.axhline(y=y, color='cyan', linewidth=0.5, alpha=0.7) + ax1.set_title(f'网格叠加 ({len(x_grid)-1}列 x {len(y_grid)-1}行)') + ax1.axis('off') + + ax2 = fig.add_subplot(2, 2, 2) + ax2.plot(col_profile, 'b-', linewidth=0.5) + for x in x_grid: + ax2.axvline(x=x, color='r', linewidth=0.5, alpha=0.5) + ax2.set_title('列投影 (列均值)') + ax2.set_xlabel('列') + + ax3 = fig.add_subplot(2, 2, 3) + ax3.plot(row_profile, 'b-', linewidth=0.5) + for y in y_grid: + ax3.axvline(x=y, color='r', linewidth=0.5, alpha=0.5) + ax3.set_title('行投影 (行均值)') + ax3.set_xlabel('行') + + ax4 = fig.add_subplot(2, 2, 4) + ax4.hist(gray.ravel(), bins=50, color='gray', edgecolor='black', linewidth=0.3) + T_otsu = filters.threshold_otsu(gray) + ax4.axvline(x=T_otsu, color='r', linestyle='--', label=f'Otsu T={T_otsu:.0f}') + ax4.set_title('灰度直方图') + ax4.set_xlabel('灰度值') + ax4.legend() + + plt.tight_layout() + return fig + + +def plot_full_segmentation(gray: np.ndarray, bw_result: np.ndarray, title: str = "完整分割结果") -> plt.Figure: + """绘制完整分割结果""" + fig, axes = plt.subplots(1, 2, figsize=(12, 6)) + fig.suptitle(title, fontsize=14) + axes[0].imshow(gray, cmap='gray') + axes[0].set_title('原始灰度图') + axes[0].axis('off') + axes[1].imshow(bw_result, cmap='gray') + axes[1].set_title('二值分割结果') + axes[1].axis('off') + plt.tight_layout() + return fig + + +# ============================================================ +# 第六部分：主流程 +# ============================================================ + +def main() -> None: + os.makedirs(OUTPUT_DIR, exist_ok=True) + + # ---- 读取图像 ---- + img_path = os.path.join(DATA_DIR, 'cDNA.png') + print(f"读取图像: {img_path}") + img = np.array(Image.open(img_path)) + if img.ndim == 3: + gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8) + else: + gray = img + print(f"图像尺寸: {gray.shape}, 灰度范围: [{gray.min()}, {gray.max()}]") + + # ---- 步骤1: 网格划分 ---- + print("\n[步骤1] 网格划分...") + x_grid, y_grid, col_prof, row_prof, x_per, y_per = gridding(gray) + print(f" 列间距估计: {x_per}, 行间距估计: {y_per}") + print(f" 网格: {len(x_grid)-1} 列 x {len(y_grid)-1} 行") + + fig_grid = plot_gridding_result(gray, x_grid, y_grid, col_prof, row_prof) + fig_grid.savefig(os.path.join(OUTPUT_DIR, 'result_gridding.png'), dpi=150, bbox_inches='tight') + print(" 保存: result_gridding.png") + + # ---- 步骤2: 提取单个子块，三种阈值方法对比 ---- + print("\n[步骤2] 单块阈值分割对比...") + bi, bj = len(y_grid) // 2, len(x_grid) // 2 + if len(x_grid) >= 3 and len(y_grid) >= 3: + r1, r2 = y_grid[bi], y_grid[bi + 1] + c1, c2 = x_grid[bj], x_grid[bj + 1] + block = gray[r1:r2, c1:c2] + else: + h, w = gray.shape + bi, bj = 0, 0 + block = gray[h // 4:3 * h // 4, w // 4:3 * w // 4] + print(f" 选取子块 [{bi},{bj}]: 尺寸{block.shape}") + + fig_compare = plot_threshold_comparison(block, f"子块[{bi},{bj}]") + fig_compare.savefig(os.path.join(OUTPUT_DIR, 'result_threshold_compare.png'), dpi=150, bbox_inches='tight') + print(" 保存: result_threshold_compare.png") + + # 迭代阈值收敛过程 + _, T_iter, history = iterative_threshold(block) + fig_conv, ax_conv = plt.subplots(figsize=(6, 4)) + ax_conv.plot(history, 'bo-', markersize=3) + ax_conv.set_title('迭代阈值收敛过程') + ax_conv.set_xlabel('迭代次数') + ax_conv.set_ylabel('阈值T') + ax_conv.axhline(y=T_iter, color='r', linestyle='--', label=f'最终T={T_iter:.1f}') + ax_conv.legend() + fig_conv.savefig(os.path.join(OUTPUT_DIR, 'result_iterative_convergence.png'), dpi=100, bbox_inches='tight') + print(" 保存: result_iterative_convergence.png") + + # ---- 步骤3: 全图逐块分割（Otsu + TV去噪） ---- + print("\n[步骤3] 全图逐块分割...") + bw_full = np.zeros_like(gray) + if len(x_grid) >= 2 and len(y_grid) >= 2: + for i in range(len(y_grid) - 1): + for j in range(len(x_grid) - 1): + r1, r2 = y_grid[i], y_grid[i + 1] + c1, c2 = x_grid[j], x_grid[j + 1] + blk = gray[r1:r2, c1:c2].copy() + if blk.size == 0: + continue + # 暗图像增强（参考MATLAB: 均值<5则×5，<30则×1.5） + blk_mean = float(np.mean(blk)) + if blk_mean < 5: + blk = np.clip(blk.astype(float) * 5, 0, 255).astype(np.uint8) + elif blk_mean < 30: + blk = np.clip(blk.astype(float) * 1.5, 0, 255).astype(np.uint8) + # TV去噪 + blk_denoised = tv_denoise(blk.astype(float), lam=12.0) + # Otsu分割 + try: + T = filters.threshold_otsu(blk_denoised.astype(np.uint8)) + bw_blk = (blk_denoised > T).astype(np.uint8) + except ValueError: + bw_blk = np.zeros(blk.shape, dtype=np.uint8) + # 后处理：保留最大连通域 + bw_blk = keep_largest_object(bw_blk, min_size=8) + bw_full[r1:r2, c1:c2] = bw_blk + + bw_full = remove_small_objects(bw_full, min_size=20) + + fig_full = plot_full_segmentation(gray, bw_full, "全图逐块Otsu分割结果") + fig_full.savefig(os.path.join(OUTPUT_DIR, 'result_full_segmentation.png'), dpi=150, bbox_inches='tight') + print(" 保存: result_full_segmentation.png") + + bw_img = Image.fromarray((bw_full * 255).astype(np.uint8)) + bw_img.save(os.path.join(OUTPUT_DIR, 'result_I_bw.png')) + print(" 保存: result_I_bw.png") + + # ---- 步骤4: 网格叠加图 ---- + if img.ndim == 3: + overlay = img[:, :, :3].copy() + else: + overlay = np.stack([gray] * 3, axis=-1) + overlay = overlay.astype(np.uint8) + for x in x_grid: + if 0 <= x < overlay.shape[1]: + overlay[:, max(x - 1, 0):min(x + 2, overlay.shape[1]), :] = [0, 255, 255] + for y in y_grid: + if 0 <= y < overlay.shape[0]: + overlay[max(y - 1, 0):min(y + 2, overlay.shape[0]), :, :] = [0, 255, 255] + Image.fromarray(overlay).save(os.path.join(OUTPUT_DIR, 'result_gridding_overlay.png')) + print(" 保存: result_gridding_overlay.png") + + print(f"\n全部完成！输出文件在: {OUTPUT_DIR}") + + +if __name__ == '__main__': + main() diff --git a/参考资料/2019-3-高污染基因芯片图像的网格划分_芦碧波.pdf b/参考资料/2019-3-高污染基因芯片图像的网格划分_芦碧波.pdf new file mode 100644 index 0000000..861e718 Binary files /dev/null and b/参考资料/2019-3-高污染基因芯片图像的网格划分_芦碧波.pdf differ diff --git a/参考资料/NewGridAndCV/GriddingAndCV.m b/参考资料/NewGridAndCV/GriddingAndCV.m new file mode 100644 index 0000000..cf0e6fa --- /dev/null +++ b/参考资料/NewGridAndCV/GriddingAndCV.m @@ -0,0 +1,304 @@ +%% Microarray Spot Finding Example +% This example shows a simple method for locating spots on a microarray and +% extracting the intensties of the spots. It can be downloaded from *MATLAB +% Central*. +% http://www.mathworks.com/matlabcentral +% +% Copyright 2004-2010 RBemis The MathWorks, Inc. + +%% Start with clean slate +clear %empty workspace (no variables) +close all %no figures +clc %empty command window + +%% Read image file +% MATLAB can read many standard image formats including TIFF, GIF and BMP +% using the |imread| command. In addition, the *Image Procesing Toolbox* +% provides support for working with specialized image file formats such as +% DICOM. This microarray image was stored as a J-PEG file. The image is +% much larger than the screen size, so |imshow| scales it down to fit and +% let's you know with a warning message. +% x = imread('MicroArraySlide.JPG'); +% imageSize = size(x) +% screenSize = get(0,'ScreenSize') +% +% iptsetpref('ImshowBorder','tight') +% imshow(x) +% title('original image') + +%% Crop specified region +% Next we use |imcrop| to extract a region of interest. You can repeat this +% for all print-tip blocks for a full microarray study. +% y = imcrop(x,[622 2467 220 227]); + + +y=imread('test.tif'); +% y=imread('cDNA.png'); +% y=imread('test001.tif'); +% y = imread('C:\Users\Administrator\Desktop\йҽѧѧʵԱȣ\ʵ\test001.tif'); +% y = imread('C:\Users\Administrator\Desktop\йҽѧѧʵԱȣ\ʵ\test003.tif'); +f1 = figure('position',[40 46 285 280]); +imshow(y) + +%% Display red & green layers +% This image was stored in RGB format. We are only interested in the red +% and green planes. To extract the red plane, simply index layer 1. For the +% green plane, layer 2. Custom colormaps make visualization more intuitive. +% Notice that spot shapes are not necessarily the same in both colors. +f2 = figure('position',[265 163 647 327]); +subplot(121) +redMap = gray(256); +redMap(:,[2 3]) = 0; +subimage(y(:,:,1),redMap) +axis off +title('red (layer 1)') +subplot(122) +greenMap = gray(256); +greenMap(:,[1 3]) = 0; +subimage(y(:,:,2),greenMap) +axis off +title('green (layer 2)') + +%% Convert RGB image to grayscale for spot finding +% Initially we care more about where the spots are located than their red +% and green intensities. Converting from RGB color to grayscale allows us +% to focus first on spot locations. +z = rgb2gray(y); +figure(f1) +imshow(z) + +%% Create horizontal profile +% We are looking for a regular grid of spots so we start by looking at the +% mean intensity for each column of the image. This will help us identify +% where the centres of the spots are and where the gaps between the spots +% can be found. +xProfile = mean(z); +f2 = figure('position',[39 346 284 73]); +plot(xProfile) +title('horizontal profile') +axis tight + +%% Estimate spot spacing by autocorrelation +% Ideally the spots would be periodicaly spaced consistently printed, but +% in practice they tend to have different sizes and intensities, so the +% horizontal profile is irregular. We can use autocorrelation to enhance +% the self similarity of the profile. The smooth result promotes peak +% finding and estimation of spot spacing. The *Signal Processing Toolbox* +% allows easy computation of the autocorrelation function using the |xcov| +% command. + +ac = xcov(xProfile); %unbiased autocorrelation +f3 = figure('position',[-3 427 569 94]); +plot(ac) +s1 = diff(ac([1 1:end])); %left slopes +s2 = diff(ac([1:end end])); %right slopes +maxima = find(s1>0 & s2<0); %peaks +estPeriod = round(median(diff(maxima))) %nominal spacing +hold on +plot(maxima,ac(maxima),'r^') +hold off +title('autocorrelation of profile') +axis tight + +%% Remove background morphologically +% We can use the spacing estimate to help design a filter to remove the +% background noise from the intensity profile. We do this with the +% |imtophat| function from the *Image Processing Toolbox*. The |strel| +% command creates a simple rectangular 1D window or line shaped structuring +% element. +seLine = strel('line',estPeriod,0); +xProfile2 = imtophat(xProfile,seLine); +f4 = figure('position',[40 443 285 76]); +plot(xProfile2) +title('enhanced horizontal profile') +axis tight + +%% Segment peaks +% Now that we have clean and anchored gaps between the peaks, we can number +% each peak region with the |bwlabel| command. These regions were segmented +% by thresholding with |im2bw|. The threshold value was automatically +% determined by statistical properties of the data using |graythresh|. This +% is a good example of image processing techniques are often useful for 1D +% data analysis. +level = graythresh(xProfile2/255)*255 +bw = im2bw(xProfile2/255,level/255); +L = bwlabel(bw); +f5 = figure('position',[40 540 285 70]); +plot(L) +axis tight +title('labelled regions') + +%% Locate centers +% We can extract the centroids of the peaks. These correspond to the +% horizontal centres of the spots. This is a common blob analysis or +% feature extraction task that can be done with |regionprops|. +stats = regionprops(L); +centroids = [stats.Centroid]; +xCenters = centroids(1:2:end) +figure(f5) +hold on +plot(xCenters,1:max(L),'ro') +hold off +title('region centers') + +%% Determine divisions between spots +% The midpoints between adjacent peaks provides grid point locations. +gap = diff(xCenters)/2; +first = xCenters(1)-gap(1); +xGrid = round([first xCenters(1:end)+gap([1:end end])]) +figure(f2) +for i=1:length(xGrid) + line(xGrid(i)*[1 1],ylim,'color','m') +end +title('vertical separators') + +%% Transpose and repeat +% We just did the analysis on the vertical grid. Now we want to do the same +% for the horizontal spacing. To do this, we simply transpose the image and +% repeat all the steps used above. This time without intermediate graphics +% display commands in order to summarize the mathematical steps of this +% algorithm. + +yProfile = mean(z'); %peak profile +ac = xcov(yProfile); %cross correlation +p1 = diff(ac([1 1:end])); +p2 = diff(ac([1:end end])); +maxima = find(p1>0 & p2<0); %peak locations +estPeriod = round(median(diff(maxima))) %spacing estimate +seLine = strel('line',estPeriod,0); +yProfile2 = imtophat(yProfile,seLine); %background removed +level = graythresh(yProfile2/255); %automatic threshold level +bw = im2bw(yProfile2/255,level); %binarized peak regions +L = bwlabel(bw); %labeled regions +stats = regionprops(L); +centroids = [stats.Centroid]; %centroids +yCenters = centroids(1:2:end) %Y parts only +gap = diff(yCenters)/2; %inner region half widths +first = yCenters(1)-gap(1); + +% list defining vertical boundaries between spot regions +yGrid = round([first yCenters(1:end)+gap([1:end end])]) + +%% Put bounding boxes around each spot +% We have now found the rectangular grid. Using pairs of neighboring grid +% points we can form bounding box regions to address each spot +% individually. The position and size coordinates of each bounding box were +% tabulated for convenience into a 4-column matrix called |ROI|, which +% stands for regions of interest. +% figure(f1) +figure() +imshow(z) +line(xGrid'*[1 1],yGrid([1 end]),'color','b') +line(xGrid([1 end]),yGrid'*[1 1],'color','b') +[X,Y] = meshgrid(xGrid(1:end-1),yGrid(1:end-1)); %xGridyGridд洢˷ָӵ +[dX,dY] = meshgrid(diff(xGrid),diff(yGrid)); +ROI = [X(:) Y(:) dX(:) dY(:)]; +% first few rows of ROI table +ROI(1:5,:) +%% +% +I=y; +% location_row=xGrid; +% location_col=yGrid; + +location_row=yGrid; +location_col=xGrid; + +for ii = 1:length(location_row)-1 + for jj = 1: length(location_col)-1 + subimage = I([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:);%Iȡÿgrid + subimage_double = double(subimage); + + %%%% תݵ[0,1]֮䣬Ȼת[0,255]֮ + sumimage_01 = (subimage_double-min(subimage_double(:)))./max(subimage_double(:))-min(subimage_double(:)); %ÿgridת[0,1]֮ +% subimage_norm = uint8(255*sumimage_01);%ԭ + subimage_norm = uint8(255*sumimage_01); + + %%%% TVȥ + subimage_norm(:,:,1)= tvdenoise(double(subimage_norm(:,:,1)),0.01); + subimage_norm(:,:,2)= tvdenoise(double(subimage_norm(:,:,2)),0.01); + subimage_norm(:,:,3)= tvdenoise(double(subimage_norm(:,:,3)),0.01); + + I_norm([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:)=subimage_norm;%ÿgrid + + + +%% + %CVзָ + + if mean(mean(subimage_norm(:,:,1)))<5 %%ֵС5ǿ5ͼ + subimage_norm=5*subimage_norm; + else if mean(mean(subimage_norm(:,:,1)))<30 %%ֵС30ǿ1.5ͼ + subimage_norm=1.5*subimage_norm; + end + end + + iter=500; + dt=0.1; + u0=cvseg(subimage_norm,iter,dt); %C-Vˮƽָ +% u1=im2bw(u0); + + %ֵĵΪĬֵ0 + [m_u0,n_u0]=size(u0); + for i=1:m_u0 + for j=1:n_u0 + if(isnan(u0(i,j))) + u0(i,j)=0; + end + end + end +% subimage_bw=u0; + %һеĵһһΪɫȡ + u1=u0; + if(u1(end,1)==1 || u1(end,end)) + u2=1-u1; + else + u2=u1; + end + subimage_bw=u2; + +% subimage_bwΪõĶֵͼ񣬼Ϊݵһͼ + %% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ָ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + I_bw([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:)=subimage_bw;%ÿgrid; + + + end +end +figure,imshow(I_bw) +%ʾԭͼѧȵĶֵͼ +all_bw=zeros(size(y(:,:,1))); +all_bw(1:location_row(end),1:location_col(end))=I_bw; +figure(),imshow(all_bw) +%% +I_bw1=all_bw; +I_bw_01=choice(I_bw1,100); %%޳100Ŀ꣨100Ϊ趨ֵ,ֱͼķȷֵ +figure(),imshow(I_bw_01); +%% +num_sub_I_bw_02=0; +for ii = 1:length(location_row)-1 + for jj = 1: length(location_col)-1 + sub_I_bw_01 = I_bw_01([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:); + sub_I_bw_02=choosemaxobj(sub_I_bw_01,8);%޳С20ĵ + + [m,n]=size(sub_I_bw_02); + for i=1:m + for j=1:n + if(sub_I_bw_02(1,j)==1 || sub_I_bw_02(i,1)==1|| sub_I_bw_02(m,j)==1 || sub_I_bw_02(i,n)==1) + num_sub_I_bw_02=num_sub_I_bw_02+1; + end + end + end + + if(num_sub_I_bw_02>10) + sub_I_bw_02=zeros(size(sub_I_bw_02)); + end + + num_sub_I_bw_02=0;%ÿμ + I_bw_Last_01([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:)=sub_I_bw_02;%ÿgrid; + end +end +all_bw_01=zeros(size(y(:,:,1))); +all_bw_01(1:location_row(end),1:location_col(end))=I_bw_Last_01; +figure(),imshow(all_bw_01); diff --git a/参考资料/NewGridAndCV/Heaviside.m b/参考资料/NewGridAndCV/Heaviside.m new file mode 100644 index 0000000..90f7d46 --- /dev/null +++ b/参考资料/NewGridAndCV/Heaviside.m @@ -0,0 +1,14 @@ +function H=Heaviside(z) +% Heaviside step function (smoothed version) +% Copyright (c) 2009, +% Yue Wu @ ECE Department, Tufts University +% All Rights Reserved + +Epsilon=10^(-5); +H=zeros(size(z,1),size(z,2)); +idx1=find(z>Epsilon); +idx2=find(z-Epsilon); +H(idx1)=1; +for i=1:length(idx2) + H(idx2(i))=1/2*(1+z(idx2(i))/Epsilon+1/pi*sin(pi*z(idx2(i))/Epsilon)); +end; \ No newline at end of file diff --git a/参考资料/NewGridAndCV/MicroArraySlide.JPG b/参考资料/NewGridAndCV/MicroArraySlide.JPG new file mode 100644 index 0000000..93819d5 Binary files /dev/null and b/参考资料/NewGridAndCV/MicroArraySlide.JPG differ diff --git a/参考资料/NewGridAndCV/R14_MicroarrayImage_CaseStudy.m b/参考资料/NewGridAndCV/R14_MicroarrayImage_CaseStudy.m new file mode 100644 index 0000000..ff16d4e --- /dev/null +++ b/参考资料/NewGridAndCV/R14_MicroarrayImage_CaseStudy.m @@ -0,0 +1,408 @@ +%% Microarray Spot Finding Example +% This example shows a simple method for locating spots on a microarray and +% extracting the intensties of the spots. It can be downloaded from *MATLAB +% Central*. +% http://www.mathworks.com/matlabcentral +% +% Copyright 2004-2010 RBemis The MathWorks, Inc. + +%% Start with clean slate +clear %empty workspace (no variables) +close all %no figures +clc %empty command window + +%% Read image file +% MATLAB can read many standard image formats including TIFF, GIF and BMP +% using the |imread| command. In addition, the *Image Procesing Toolbox* +% provides support for working with specialized image file formats such as +% DICOM. This microarray image was stored as a J-PEG file. The image is +% much larger than the screen size, so |imshow| scales it down to fit and +% let's you know with a warning message. +% x = imread('MicroArraySlide.JPG'); +% imageSize = size(x) +% screenSize = get(0,'ScreenSize') +% +% iptsetpref('ImshowBorder','tight') +% imshow(x) +% title('original image') + +%% Crop specified region +% Next we use |imcrop| to extract a region of interest. You can repeat this +% for all print-tip blocks for a full microarray study. +% y = imcrop(x,[622 2467 220 227]); +y=imread('test.tif'); +f1 = figure('position',[40 46 285 280]); +imshow(y) + +%% Display red & green layers +% This image was stored in RGB format. We are only interested in the red +% and green planes. To extract the red plane, simply index layer 1. For the +% green plane, layer 2. Custom colormaps make visualization more intuitive. +% Notice that spot shapes are not necessarily the same in both colors. +f2 = figure('position',[265 163 647 327]); +subplot(121) +redMap = gray(256); +redMap(:,[2 3]) = 0; +subimage(y(:,:,1),redMap) +axis off +title('red (layer 1)') +subplot(122) +greenMap = gray(256); +greenMap(:,[1 3]) = 0; +subimage(y(:,:,2),greenMap) +axis off +title('green (layer 2)') + +%% Convert RGB image to grayscale for spot finding +% Initially we care more about where the spots are located than their red +% and green intensities. Converting from RGB color to grayscale allows us +% to focus first on spot locations. +z = rgb2gray(y); +figure(f1) +imshow(z) + +%% Create horizontal profile +% We are looking for a regular grid of spots so we start by looking at the +% mean intensity for each column of the image. This will help us identify +% where the centres of the spots are and where the gaps between the spots +% can be found. +xProfile = mean(z); +f2 = figure('position',[39 346 284 73]); +plot(xProfile) +title('horizontal profile') +axis tight + +%% Estimate spot spacing by autocorrelation +% Ideally the spots would be periodicaly spaced consistently printed, but +% in practice they tend to have different sizes and intensities, so the +% horizontal profile is irregular. We can use autocorrelation to enhance +% the self similarity of the profile. The smooth result promotes peak +% finding and estimation of spot spacing. The *Signal Processing Toolbox* +% allows easy computation of the autocorrelation function using the |xcov| +% command. + +ac = xcov(xProfile); %unbiased autocorrelation +f3 = figure('position',[-3 427 569 94]); +plot(ac) +s1 = diff(ac([1 1:end])); %left slopes +s2 = diff(ac([1:end end])); %right slopes +maxima = find(s1>0 & s2<0); %peaks +estPeriod = round(median(diff(maxima))) %nominal spacing +hold on +plot(maxima,ac(maxima),'r^') +hold off +title('autocorrelation of profile') +axis tight + +%% Remove background morphologically +% We can use the spacing estimate to help design a filter to remove the +% background noise from the intensity profile. We do this with the +% |imtophat| function from the *Image Processing Toolbox*. The |strel| +% command creates a simple rectangular 1D window or line shaped structuring +% element. +seLine = strel('line',estPeriod,0); +xProfile2 = imtophat(xProfile,seLine); +f4 = figure('position',[40 443 285 76]); +plot(xProfile2) +title('enhanced horizontal profile') +axis tight + +%% Segment peaks +% Now that we have clean and anchored gaps between the peaks, we can number +% each peak region with the |bwlabel| command. These regions were segmented +% by thresholding with |im2bw|. The threshold value was automatically +% determined by statistical properties of the data using |graythresh|. This +% is a good example of image processing techniques are often useful for 1D +% data analysis. +level = graythresh(xProfile2/255)*255 +bw = im2bw(xProfile2/255,level/255); +L = bwlabel(bw); +f5 = figure('position',[40 540 285 70]); +plot(L) +axis tight +title('labelled regions') + +%% Locate centers +% We can extract the centroids of the peaks. These correspond to the +% horizontal centres of the spots. This is a common blob analysis or +% feature extraction task that can be done with |regionprops|. +stats = regionprops(L); +centroids = [stats.Centroid]; +xCenters = centroids(1:2:end) +figure(f5) +hold on +plot(xCenters,1:max(L),'ro') +hold off +title('region centers') + +%% Determine divisions between spots +% The midpoints between adjacent peaks provides grid point locations. +gap = diff(xCenters)/2; +first = xCenters(1)-gap(1); +xGrid = round([first xCenters(1:end)+gap([1:end end])]) +figure(f2) +for i=1:length(xGrid) + line(xGrid(i)*[1 1],ylim,'color','m') +end +title('vertical separators') + +%% Transpose and repeat +% We just did the analysis on the vertical grid. Now we want to do the same +% for the horizontal spacing. To do this, we simply transpose the image and +% repeat all the steps used above. This time without intermediate graphics +% display commands in order to summarize the mathematical steps of this +% algorithm. + +yProfile = mean(z'); %peak profile +ac = xcov(yProfile); %cross correlation +p1 = diff(ac([1 1:end])); +p2 = diff(ac([1:end end])); +maxima = find(p1>0 & p2<0); %peak locations +estPeriod = round(median(diff(maxima))) %spacing estimate +seLine = strel('line',estPeriod,0); +yProfile2 = imtophat(yProfile,seLine); %background removed +level = graythresh(yProfile2/255); %automatic threshold level +bw = im2bw(yProfile2/255,level); %binarized peak regions +L = bwlabel(bw); %labeled regions +stats = regionprops(L); +centroids = [stats.Centroid]; %centroids +yCenters = centroids(1:2:end) %Y parts only +gap = diff(yCenters)/2; %inner region half widths +first = yCenters(1)-gap(1); + +% list defining vertical boundaries between spot regions +yGrid = round([first yCenters(1:end)+gap([1:end end])]) + +%% Put bounding boxes around each spot +% We have now found the rectangular grid. Using pairs of neighboring grid +% points we can form bounding box regions to address each spot +% individually. The position and size coordinates of each bounding box were +% tabulated for convenience into a 4-column matrix called |ROI|, which +% stands for regions of interest. +% figure(f1) +figure() +imshow(z) +line(xGrid'*[1 1],yGrid([1 end]),'color','b') +line(xGrid([1 end]),yGrid'*[1 1],'color','b') +[X,Y] = meshgrid(xGrid(1:end-1),yGrid(1:end-1)); %xGridyGridд洢˷ָӵ +[dX,dY] = meshgrid(diff(xGrid),diff(yGrid)); +ROI = [X(:) Y(:) dX(:) dY(:)]; +% first few rows of ROI table +ROI(1:5,:) + +%% Segment spots from background by thresholding +% Applying a single threshold level to the whole image so all spots are +% detected equally is generally a good idea. However, in this case is +% doesn't work so well due to large differences in spot brightness. +fSpots = figure('position',[265 163 647 327]); +subplot(121) +imshow(z) +title('gray image') +subplot(122) +bw = im2bw(z,graythresh(z)); +imshow(bw) +title('global threshold') + +%% Apply logarithmic transformation then threshold intensities +% One way to equalize large variations in magnitude is by transforming +% intensity values to logarithmic space. This works much better but some +% weak spots are still missed. +figure(fSpots) +subplot(121) +z2 = uint8(log(double(z)+1)/log(255)*255); +imshow(z2) +title('log intensity') +subplot(122) +bw = im2bw(z2,graythresh(z2)); +imshow(bw) +title('global threshold') + +%% Try local thresholding instead +% Alternatively, the bounding boxes can be used to determine local +% threshold values for each spot. The code is a little more sophisticated, +% requiring looping and indexing. Unfortunately, the results are mixed. +% Weak spots showed up well but spots with bright perimeters were as bad as +% the original global threshold before log space transformation. +figure(fSpots) +subplot(122) +bw = false(size(z)); +for i=1:length(ROI) + rows = round(ROI(i,2))+[0:(round(ROI(i,4))-1)]; + cols = round(ROI(i,1))+[0:(round(ROI(i,3))-1)]; + spot = z(rows,cols); + bw(rows,cols) = im2bw(spot,graythresh(spot)); +end +imshow(bw) +title('local threshold') + +%% Logically combine local and global thresholds +% Since both have their merits, let's combine the best of both approaches. +% This can be done using logial operation on the binary masks. These spot +% segmentation results are indeed much better. +figure(fSpots) +subplot(121) +bw = im2bw(z2,graythresh(z2)); +for i=1:length(ROI) + rows = round(ROI(i,2))+[0:(round(ROI(i,4))-1)]; + cols = round(ROI(i,1))+[0:(round(ROI(i,3))-1)]; + spot = z(rows,cols); + bw(rows,cols) = bw(rows,cols) | im2bw(spot,graythresh(spot)); +end +imshow(bw) +title('combined threshold') +subplot(122) +imshow(z) +title('linear intensity') + +%% Fill holes to solidify spots +% The silhouettes of some spots still contained pinholes. The whole image +% could be filled using a single call to |imfill| but this may not be a +% good idea. Notice that some spots run together. If four mutually adjacent +% spots (sharing a common corner) were all joined at their edges then a +% single function call would incorrectly fill in the common corner as well. +% To avoid that possibility, it's good insurance to fill each spot one +% bounding box region at a time by looping. Indeed, the spot segmentation +% now looks quite good. +figure(fSpots) +subplot(121) +warning off MATLAB:intConvertOverflow +for i=1:length(ROI) + rows = round(ROI(i,2))+[0:(round(ROI(i,4))-1)]; + cols = round(ROI(i,1))+[0:(round(ROI(i,3))-1)]; + bw(rows,cols) = imfill(bw(rows,cols),'holes'); +end +imshow(bw) %յķָbw +title('filled pinholes') + +%% Label spot masks by bounding box +% If the gridding went well, all spots should be a single color. The +% results here are pretty good. There is still room for improvement. +% +% TODO List: +% +% * Due to slightly irregular spacing, for some spots a few pixels were +% mislabeled. With additional processing, the algorithm could be extended +% to reclassify these stray pixels. +% +% * The crescent shaped spot in row 8, column 4 could be completed to be +% more circular by using the 'ConvexImage' return value from |regionprops|. +% +% * The few stray pixels that are not attached to any spots could be +% removed as well. +% +% However, in this case the spot segmentation is good enough to proceed. +L = zeros(size(bw)); +for i=1:length(ROI) + rows = ROI(i,2)+[0:(ROI(i,4)-1)]; + cols = ROI(i,1)+[0:(ROI(i,3)-1)]; + rectMask = L(rows,cols); + spotMask = bw(rows,cols); + rectMask(spotMask) = i; + L(rows,cols) = rectMask; +end +map = [0 0 0; 0.5+0.5*rand(length(ROI),3)]; +figure(f1) +imshow(L+1,map) + +%% Extract first spot for measurement +% We will now examine the first spot closely to see how we can measure its +% red and green intensities, and ultimately quantify its gene expression +% value. The measurement technique can then be repeated for all spots. +rect = ROI(1,:); %[X Y dX dY] +spot = imcrop(y,rect); %region around spot +figure(f1) +imshow(spot,'notruesize') + +%% Measure spot intensity & releative expression level +% We now simply calculate the nominal intensity over the spot for both the +% red and green layers. A measure of gene expression level can then be +% calculated from the two color intensities. Here a simple log-ratio +% measurement is shown. Other more robust measures could be used instead. +% You could also perform some analysis of the quality of the spot. +mask = imcrop(L,rect)==1; +for i=1:2 + layer = spot(:,:,i); + intensity(i) = double(median(layer(mask))); +end +intensity +expressionLevel = log(intensity(1)/intensity(2)) + +%% Remove background, calculate again and compare measurements +% If you noticed, the background intensity around the spot was not zero. +% This could bias results. To see how much difference it makes, we can +% perform background subtraction around all spots, again using |imtophat| +% but this time in 2D on the image using a disk shaped structuring element. +% Then we can calculate color intensities and relative expression level +% again to see what effect background bias had on the measurement. In this +% case the measurement shows more downregulation with background removed. +seDisk = strel('disk',round(estPeriod)); +spot2 = imtophat(spot,seDisk); +for i=1:2 + layer = spot2(:,:,i); + intensity(i) = double(median(layer(mask))); +end +intensity +expressionLevel = log(intensity(1)/intensity(2)) + +%% Set up graphical display for results +% It is helpful to see red and green intensity values overlayed onto the +% respective color images to gain confidence that measured intensities make +% sense. It is also be helpful to overlay quantitative expression levels +% onto the original image to provide additional visual assurance of +% measurement results. The rectangular grid also helps correlate measured +% values between images. The flexibility of MATLAB's powerful *Handle +% Graphics* engine allow custom graphics like this to be set up quickly and +% easily. +f7 = figure('position',[52 94 954 425]); +ax(1) = subplot(121); +subimage(y(:,:,1),redMap) +title('red intensity') +ax(2) = subplot(122); +subimage(y(:,:,2),greenMap) +title('green intensity') +f8 = figure('position',[316 34 482 497]); +ax(3) = get(imshow(y,'notruesize'),'parent'); +title('gene expression') +for i=1:3 + axes(ax(i)) + axis off + line(xGrid'*[1 1],yGrid([1 end]),'color',0.5*[1 1 1]) + line(xGrid([1 end]),yGrid'*[1 1],'color',0.5*[1 1 1]) +end + +%% Repeat measurement for all spots +% We now repeat the spot extraction and intensity calculation for all the +% spots in the grid. Here the measured values were tabulated as additional +% columns beside the ROI positions for each spot into a new matrix called +% |spotData|. +figure(f7), figure(f8) +spotData = [ROI zeros(length(ROI),5)]; +for i=1:length(ROI) + spot = imcrop(y,ROI(i,:)); %raw image + spot2 = imtophat(spot,seDisk); %background removed + mask = imcrop(L,ROI(i,:))==i; %spot mask + for j=1:2 + layer = spot2(:,:,j); %color layer + intensity(j) = double(median(layer(mask))); + text(ROI(i,1)+ROI(i,3)/2,ROI(i,2)+ROI(i,4)/2,sprintf('%.0f',intensity(j)),... + 'color','y','HorizontalAlignment','center','parent',ax(j)) + rawLayer = spot(:,:,j); + rawIntensity(j) = double(median(layer(mask))); + end + expression = log(intensity(1)/intensity(2)); + text(ROI(i,1)+ROI(i,3)/2,ROI(i,2)+ROI(i,4)/2,sprintf('%.2f',expression),... + 'color','w','HorizontalAlignment','center','parent',ax(3)) + drawnow + spotData(i,5:9) = [intensity(:)' expression rawIntensity(:)']; +end + +%% Export spot data to Excel spreadsheet +% MATLAB can write to many standard formats. We will use |xlswrite| to save +% the |spotData| to an Excel workbook. +% +% TODO list: +% +% * prepend column names first (see |xlswrite| doc example using cell arrays) +% +% * programmatically open spreadsheet in Excel (see |winopen| doc) +xlswrite('microarray.xls',spotData) diff --git a/参考资料/NewGridAndCV/cDNA.png b/参考资料/NewGridAndCV/cDNA.png new file mode 100644 index 0000000..494542a Binary files /dev/null and b/参考资料/NewGridAndCV/cDNA.png differ diff --git a/参考资料/NewGridAndCV/checkstop.m b/参考资料/NewGridAndCV/checkstop.m new file mode 100644 index 0000000..ba557e0 --- /dev/null +++ b/参考资料/NewGridAndCV/checkstop.m @@ -0,0 +1,37 @@ +function indicator = checkstop(old,new,dt) +% indicate whether we should performance further iteraions or stop + +% Copyright (c) 2009, +% Yue Wu @ ECE Department, Tufts University +% All Rights Reserved + +layer = size(new,3); + +for i = 1:layer + old_{i} = old(:,:,i); + new_{i} = new(:,:,i); +end + +if layer + ind = find(abs(new)<=.5); + M = length(ind); + Q = sum(abs(new(ind)-old(ind)))./M; + if Q<=dt*.18^2 + indicator = 1; + else + indicator = 0; + end +else + ind1 = find(abs(old_{1})<1); + ind2 = find(abs(old_{2})<1); + M1 = length(ind1); + M2 = length(ind2); + Q1 = sum(abs(new_{1}(ind1)-old_{1}(ind1)))./M1; + Q2 = sum(abs(new_{2}(ind2)-old_{2}(ind2)))./M2; + if Q1<=dt*.18^2 && Q2<=dt*.18^2 + indicator = 1; + else + indicator = 0; + end +end +return diff --git a/参考资料/NewGridAndCV/chenvese.m b/参考资料/NewGridAndCV/chenvese.m new file mode 100644 index 0000000..6203b95 --- /dev/null +++ b/参考资料/NewGridAndCV/chenvese.m @@ -0,0 +1,407 @@ +%========================================================================== +% +% Active contour with Chen-Vese Method +% for image segementation +% +% Implemented by Yue Wu (yue.wu@tufts.edu) +% Tufts University +% Feb 2009 +% http://sites.google.com/site/rexstribeofimageprocessing/ +% +% all rights reserved +% Last update 02/26/2009 +%-------------------------------------------------------------------------- +% Usage of varibles: +% input: +% I = any gray/double/RGB input image +% mask = initial mask, either customerlized or built-in +% num_iter = total number of iterations +% mu = weight of length term +% method = submethods pick from ('chen','vector','multiphase') +% +% Types of built-in mask functions +% 'small' = create a small circular mask +% 'medium' = create a medium circular mask +% 'large' = create a large circular mask +% 'whole' = create a mask with holes around +% 'whole+small' = create a two layer mask with one layer small +% circular mask and the other layer with holes around +% (only work for method 'multiphase') +% Types of methods +% 'chen' = general CV method +% 'vector' = CV method for vector image +% 'multiphase'= CV method for multiphase (2 phases applied here) +% +% output: +% phi0 = updated level set function +% +%-------------------------------------------------------------------------- +% +% Description: This code implements the paper: "Active Contours Without +% Edges" by Chan and Vese for method 'chen', the paper:"Active Contours Without +% Edges for vector image" by Chan and Vese for method 'vector', and the paper +% "A Multiphase Level Set Framework for Image Segmentation Using the +% Mumford and Shah Model" by Chan and Vese. +% +%-------------------------------------------------------------------------- +% Deomo: Please see HELP file for details +%========================================================================== + +function seg = chenvese(I,mask,num_iter,mu,method) + +%% +%-- Default settings +% length term mu = 0.2 and default method = 'chan' + if(~exist('mu','var')) + mu=0.2; + end + + if(~exist('method','var')) + method = 'chan'; + end + +%-- End default settings + +%% +%-- Initializations on input image I and mask +% resize original image + s = 200./min(size(I,1),size(I,2)); % resize scale + if s<1 + I = imresize(I,s); + end + +% auto mask settings + if ischar(mask) + switch lower (mask) + case 'small' + mask = maskcircle2(I,'small'); + case 'medium' + mask = maskcircle2(I,'medium'); + case 'large' + mask = maskcircle2(I,'large'); + case 'whole' + mask = maskcircle2(I,'whole'); + %mask = init_mask(I,30); + case 'whole+small' + m1 = maskcircle2(I,'whole'); + m2 = maskcircle2(I,'small'); + mask = zeros(size(I,1),size(I,2),2); + mask(:,:,1) = m1(:,:,1); + mask(:,:,2) = m2(:,:,2); + otherwise + error('unrecognized mask shape name (MASK).'); + end + else + if s<1 + mask = imresize(mask,s); + end + if size(mask,1)>size(I,1) || size(mask,2)>size(I,2) + error('dimensions of mask unmathch those of the image.') + end + switch lower(method) + case 'multiphase' + if (size(mask,3) == 1) + error('multiphase requires two masks but only gets one.') + end + end + + end + + +switch lower(method) + case 'chan' + if size(I,3)== 3 + P = rgb2gray(uint8(I)); + P = double(P); + elseif size(I,3) == 2 + P = 0.5.*(double(I(:,:,1))+double(I(:,:,2))); + else + P = double(I); + end + layer = 1; + + case 'vector' + s = 200./min(size(I,1),size(I,2)); % resize scale + I = imresize(I,s); + mask = imresize(mask,s); + layer = size(I,3); + if layer == 1 + display('only one image component for vector image') + end + P = double(I); + + case 'multiphase' + layer = size(I,3); + if size(I,1)*size(I,2)>200^2 + s = 200./min(size(I,1),size(I,2)); % resize scale + I = imresize(I,s); + mask = imresize(mask,s); + end + + P = double(I); %P store the original image + otherwise + error('!invalid method') +end +%-- End Initializations on input image I and mask + +%% +%-- Core function +switch lower(method) + case {'chan','vector'} + %-- SDF + % Get the distance map of the initial mask + + mask = mask(:,:,1); + phi0 = bwdist(mask)-bwdist(1-mask)+im2double(mask)-.5; + % initial force, set to eps to avoid division by zeros + force = eps; + %-- End Initialization + + %-- Display settings +% figure(); +% subplot(2,2,1); imshow(I); title('Input Image'); +% subplot(2,2,2); contour(flipud(phi0), [0 0], 'r','LineWidth',1); title('initial contour'); +% subplot(2,2,3); title('Segmentation'); +% figure(1); +% imshow(I); title('Input Image'); +% figure(2); contour(flipud(phi0), [0 0], 'r','LineWidth',1); title('initial contour'); +% figure(3); + %-- End Display original image and mask + + %-- Main loop + for n=1:num_iter + inidx = find(phi0>=0); % frontground index + outidx = find(phi0<0); % background index + force_image = 0; % initial image force for each layer + for i=1:layer + L = im2double(P(:,:,i)); % get one image component + c1 = sum(sum(L.*Heaviside(phi0)))/(length(inidx)+eps); % average inside of Phi0 + c2 = sum(sum(L.*(1-Heaviside(phi0))))/(length(outidx)+eps); % verage outside of Phi0 + force_image=-(L-c1).^2+(L-c2).^2+force_image; + % sum Image Force on all components (used for vector image) + % if 'chan' is applied, this loop become one sigle code as a + % result of layer = 1 + end + + % calculate the external force of the image + force = mu*kappa(phi0)./max(max(abs(kappa(phi0))))+1/layer.*force_image; + + % normalized the force + force = force./(max(max(abs(force)))); + + % get stepsize dt + dt=0.5; + + % get parameters for checking whether to stop + old = phi0; + phi0 = phi0+dt.*force; + new = phi0; + indicator = checkstop(old,new,dt); + + % intermediate output + if(mod(n,20) == 0) +% showphi(I,phi0,n); + end; + if indicator % decide to stop or continue +% showphi(I,phi0,n); + + %make mask from SDF + seg = phi0<=0; %-- Get mask from levelset + +% figure(4); imshow(seg); title('Global Region-Based Segmentation'); + + return; + end + end; +% showphi(I,phi0,n); + + + + %make mask from SDF + seg = phi0<=0; %-- Get mask from levelset + +% figure(4);; imshow(seg); title('Global Region-Based Segmentation'); + case 'multiphase' + %-- Initializations + % Get the distance map of the initial masks + mask1 = mask(:,:,1); + mask2 = mask(:,:,2); + phi1=bwdist(mask1)-bwdist(1-mask1)+im2double(mask1)-.5;%Get phi1 from the initial mask 1 + phi2=bwdist(mask2)-bwdist(1-mask2)+im2double(mask2)-.5;%Get phi1 from the initial mask 2 + + %-- Display settings + figure(); + subplot(2,2,1); + if layer ~= 1 + imshow(I); title('Input Image'); + else + imagesc(P); axis image; colormap(gray);title('Input Image'); + end + subplot(2,2,2); + hold on + contour(flipud(mask1),[0,0],'r','LineWidth',2.5); + contour(flipud(mask1),[0,0],'x','LineWidth',1); + contour(flipud(mask2),[0,0],'g','LineWidth',2.5); + contour(flipud(mask2),[0,0],'x','LineWidth',1); + title('initial contour'); + hold off + subplot(2,2,3); title('Segmentation'); + %-- End display settings + + %Main loop + for n=1:num_iter + %-- Narrow band for each phase + nb1 = find(phi1<1.2 & phi1>=-1.2); %narrow band of phi1 + inidx1 = find(phi1>=0); %phi1 frontground index + outidx1 = find(phi1<0); %phi1 background index + + nb2 = find(phi2<1.2 & phi2>=-1.2); %narrow band of phi2 + inidx2 = find(phi2>=0); %phi2 frontground index + outidx2 = find(phi2<0); %phi2 background index + %-- End initiliazaions on narrow band + + %-- Mean calculations for different partitions + %c11 = mean (phi1>0 & phi2>0) + %c12 = mean (phi1>0 & phi2<0) + %c21 = mean (phi1<0 & phi2>0) + %c22 = mean (phi1<0 & phi2<0) + + cc11 = intersect(inidx1,inidx2); %index belong to (phi1>0 & phi2>0) + cc12 = intersect(inidx1,outidx2); %index belong to (phi1>0 & phi2<0) + cc21 = intersect(outidx1,inidx2); %index belong to (phi1<0 & phi2>0) + cc22 = intersect(outidx1,outidx2); %index belong to (phi1<0 & phi2<0) + + f_image11 = 0; + f_image12 = 0; + f_image21 = 0; + f_image22 = 0; % initial image force for each layer + + for i=1:layer + L = im2double(P(:,:,i)); % get one image component + + if isempty(cc11) + c11 = eps; + else + c11 = mean(L(cc11)); + end + + if isempty(cc12) + c12 = eps; + else + c12 = mean(L(cc12)); + end + + if isempty(cc21) + c21 = eps; + else + c21 = mean(L(cc21)); + end + + if isempty(cc22) + c22 = eps; + else + c22 = mean(L(cc22)); + end + + %-- End mean calculation + + %-- Force calculation and normalization + % force on each partition + + f_image11=(L-c11).^2.*Heaviside(phi1).*Heaviside(phi2)+f_image11; + f_image12=(L-c12).^2.*Heaviside(phi1).*(1-Heaviside(phi2))+f_image12; + f_image21=(L-c21).^2.*(1-Heaviside(phi1)).*Heaviside(phi2)+f_image21; + f_image22=(L-c22).^2.*(1-Heaviside(phi1)).*(1-Heaviside(phi2))+f_image22; + end + + % sum Image Force on all components (used for vector image) + % if 'chan' is applied, this loop become one sigle code as a + % result of layer = 1 + + % calculate the external force of the image + + % curvature on phi1 + curvature1 = mu*kappa(phi1); + curvature1 = curvature1(nb1); + % image force on phi1 + fim1 = 1/layer.*(-f_image11(nb1)+f_image21(nb1)-f_image12(nb1)+f_image22(nb1)); + fim1 = fim1./max(abs(fim1)+eps); + + % curvature on phi2 + curvature2 = mu*kappa(phi2); + curvature2 = curvature2(nb2); + % image force on phi2 + fim2 = 1/layer.*(-f_image11(nb2)+f_image12(nb2)-f_image21(nb2)+f_image22(nb2)); + fim2 = fim2./max(abs(fim2)+eps); + + % force on phi1 and phi2 + force1 = curvature1+fim1; + force2 = curvature2+fim2; + %-- End force calculation + + % detal t + dt = 1.5; + + old(:,:,1) = phi1; + old(:,:,2) = phi2; + + %update of phi1 and phi2 + phi1(nb1) = phi1(nb1)+dt.*force1; + phi2(nb2) = phi2(nb2)+dt.*force2; + + new(:,:,1) = phi1; + new(:,:,2) = phi2; + + indicator = checkstop(old,new,dt); + + if indicator + showphi(I, new, n); + %make mask from SDF + seg11 = (phi1>0 & phi2>0); %-- Get mask from levelset + seg12 = (phi1>0 & phi2<0); + seg21 = (phi1<0 & phi2>0); + seg22 = (phi1<0 & phi2<0); + + se = strel('disk',1); + aa1 = imerode(seg11,se); + aa2 = imerode(seg12,se); + aa3 = imerode(seg21,se); + aa4 = imerode(seg22,se); + seg = aa1+2*aa2+3*aa3+4*aa4; + + subplot(2,2,4); imagesc(seg);axis image;title('Global Region-Based Segmentation'); + + return + end + % re-initializations + phi1 = reinitialization(phi1, 0.6);%sussman(phi1, 0.6);% + phi2 = reinitialization(phi2, 0.6);%sussman(phi2,0.6); + + %intermediate output + if(mod(n,20) == 0) + phi(:,:,1) = phi1; + phi(:,:,2) = phi2; + showphi(I, phi, n); + end; + end; + phi(:,:,1) = phi1; + phi(:,:,2) = phi2; + showphi(I, phi, n); + %make mask from SDF + seg11 = (phi1>0 & phi2>0); %-- Get mask from levelset + seg12 = (phi1>0 & phi2<0); + seg21 = (phi1<0 & phi2>0); + seg22 = (phi1<0 & phi2<0); + + se = strel('disk',1); + aa1 = imerode(seg11,se); + aa2 = imerode(seg12,se); + aa3 = imerode(seg21,se); + aa4 = imerode(seg22,se); + seg = aa1+2*aa2+3*aa3+4*aa4; + %seg = bwlabel(seg); + subplot(2,2,4); imagesc(seg);axis image;title('Global Region-Based Segmentation'); + + +end + diff --git a/参考资料/NewGridAndCV/choice.m b/参考资料/NewGridAndCV/choice.m new file mode 100644 index 0000000..77a7909 --- /dev/null +++ b/参考资料/NewGridAndCV/choice.m @@ -0,0 +1,51 @@ +function U1=choice(u0,Daxiao) +L = bwlabel(u0,8); +%%ȡ +STATS = regionprops(L,'Area'); + +%%ȡ +area=[STATS.Area]; +% [m_area,n_area]=size(area); +% num_big=0; +% for i=1:n_area +% if area(1,i)>100 +% num_big=num_big+1; +% end +% end + +%ѡȡĿ +area_num=area(area>Daxiao); +%Ŀ +Length_area_num=numel(area_num) ; +%趨飬洢Ŀı +array=zeros(1,Length_area_num); + +%ĿıŴ洢array +% for k=1:Length_area_num +% % array(1,k)=find(area==area_num(1,k)); +% +% end +num_big=1; +[m_area,n_area]=size(area); +for i=1:n_area + if area(1,i)>100 + array(1,num_big)=i; + num_big=num_big+1; + end +end + + +L_new=L; +[m_L_new,n_L_new]=size(L_new); +for i=1:m_L_new + for j=1:n_L_new + for k1=1:Length_area_num + if (L_new(i,j)==array(1,k1)) + L_new(i,j)=0; + end + end + end +end +U1=L_new; + + diff --git a/参考资料/NewGridAndCV/choosemaxobj.m b/参考资料/NewGridAndCV/choosemaxobj.m new file mode 100644 index 0000000..e533b84 --- /dev/null +++ b/参考资料/NewGridAndCV/choosemaxobj.m @@ -0,0 +1,33 @@ +function Lnew=choosemaxobj(L_old,night) +%%Ӷֵͼ4,8ȡ壬Ϊֵͼ + +%%ע +L = bwlabel(L_old,night); + +%%ȡ +STATS = regionprops(L,'Area'); + +%%ȡ +area=[STATS.Area]; +area_max=max(area); + +if area_max>20 %趨20ΪСֵС20ĵ㣬ȥ + %%ȷı + area_num=find(area==area_max); + + %%߼жϣ屣 + if(length(area_num)>1) + Lnew=(L==area_num(1,1)); + else + Lnew=(L==area_num); + end + + else + Lnew=0; +end + +% LL = bwlabel(Lnew,8); +% %%ȡ +% STATSLnew = regionprops(LL,'Area'); +% %%ȡ +% areanew=[STATSLnew.Area] \ No newline at end of file diff --git a/参考资料/NewGridAndCV/contour_bw.m b/参考资料/NewGridAndCV/contour_bw.m new file mode 100644 index 0000000..0ae6627 --- /dev/null +++ b/参考资料/NewGridAndCV/contour_bw.m @@ -0,0 +1,48 @@ + %%%%ֵͼ:IJһΪ1,Ϊڵ,Ϊ. + function output=contour_bw(input) + %input=imread('bw_image.bmp'); + %%ͼ:input + %%ͼ:ڱ߽紦Ϊ255,ڲΪ0 +% input=Lnew_out; + + [r,c]=size(input); + output=255*zeros(r,c); + %Ķ +% global num; + global r_avg; + r_avg=0; + global g_avg; + g_avg=0; + global b_avg; + b_avg=0; +% global rect_rgb_01;%иԭͼ + Reb_dot=0; + + + for i=2:r-1 + for j=2:c-1 + if (input(i,j)-input(i,j+1)==1 || input(i,j)-input(i,j-1)==1 || input(i,j)-input(i+1,j)==1 || input(i,j)-input(i-1,j)==1) + output(i,j)=255; +% r_avg=r_avg+double(rect_rgb_01(i,j,1)); +% g_avg=g_avg+double(rect_rgb_01(i,j,2)); +% b_avg=b_avg+double(rect_rgb_01(i,j,3)); +% num=num+1; + Reb_dot=Reb_dot+1; + else + output(i,j)=0; + end + end + Reb_dot + end + + + + + +% r_avg=uint8(r_avg/num); +% g_avg=uint8(g_avg/num); +% b_avg=uint8(b_avg/num); + + + + \ No newline at end of file diff --git a/参考资料/NewGridAndCV/cvseg.m b/参考资料/NewGridAndCV/cvseg.m new file mode 100644 index 0000000..955a613 --- /dev/null +++ b/参考资料/NewGridAndCV/cvseg.m @@ -0,0 +1,9 @@ +function seg=cvseg(imin,iter,dt) + +% Customerlized Mask +m = zeros(size(imin,1),size(imin,2)); +m(20:50,20:50) = 1; +seg = chenvese(uint8(imin),'small',iter,dt,'chan'); % ability on gray image +% seg = chenvese(uint8(imin),'whole',500,0.1,'chan'); % ability on gray image + +% imwrite(uint8(255*seg),'tvbseg.bmp'); diff --git a/参考资料/NewGridAndCV/demo_GriddingAndCV.m b/参考资料/NewGridAndCV/demo_GriddingAndCV.m new file mode 100644 index 0000000..62a4348 --- /dev/null +++ b/参考资料/NewGridAndCV/demo_GriddingAndCV.m @@ -0,0 +1,305 @@ +%% Microarray Spot Finding Example +% This example shows a simple method for locating spots on a microarray and +% extracting the intensties of the spots. It can be downloaded from *MATLAB +% Central*. +% http://www.mathworks.com/matlabcentral +% +% Copyright 2004-2010 RBemis The MathWorks, Inc. + +%% Start with clean slate +clear %empty workspace (no variables) +close all %no figures +clc %empty command window + +%% Read image file +% MATLAB can read many standard image formats including TIFF, GIF and BMP +% using the |imread| command. In addition, the *Image Procesing Toolbox* +% provides support for working with specialized image file formats such as +% DICOM. This microarray image was stored as a J-PEG file. The image is +% much larger than the screen size, so |imshow| scales it down to fit and +% let's you know with a warning message. +% x = imread('MicroArraySlide.JPG'); +% imageSize = size(x) +% screenSize = get(0,'ScreenSize') +% +% iptsetpref('ImshowBorder','tight') +% imshow(x) +% title('original image') + +%% Crop specified region +% Next we use |imcrop| to extract a region of interest. You can repeat this +% for all print-tip blocks for a full microarray study. +% y = imcrop(x,[622 2467 220 227]); + + +y=imread('test.tif'); +% y=imread('cDNA.png'); +% y=imread('test001.tif'); +% y = imread('C:\Users\Administrator\Desktop\йҽѧѧʵԱȣ\ʵ\test001.tif'); +% y = imread('C:\Users\Administrator\Desktop\йҽѧѧʵԱȣ\ʵ\test003.tif'); +f1 = figure('position',[40 46 285 280]); +imshow(y) +title('ԭʼͼ') + +%% Display red & green layers +% This image was stored in RGB format. We are only interested in the red +% and green planes. To extract the red plane, simply index layer 1. For the +% green plane, layer 2. Custom colormaps make visualization more intuitive. +% Notice that spot shapes are not necessarily the same in both colors. +% f2 = figure('position',[265 163 647 327]); +% subplot(121) +% redMap = gray(256); +% redMap(:,[2 3]) = 0; +% subimage(y(:,:,1),redMap) +% axis off +% title('red (layer 1)') +% subplot(122) +% greenMap = gray(256); +% greenMap(:,[1 3]) = 0; +% subimage(y(:,:,2),greenMap) +% axis off +% title('green (layer 2)') + +%% Convert RGB image to grayscale for spot finding +% Initially we care more about where the spots are located than their red +% and green intensities. Converting from RGB color to grayscale allows us +% to focus first on spot locations. +z = rgb2gray(y); +figure(f1) +figure,imshow(z),title('Ҷͼ') + +%% Create horizontal profile +% We are looking for a regular grid of spots so we start by looking at the +% mean intensity for each column of the image. This will help us identify +% where the centres of the spots are and where the gaps between the spots +% can be found. +xProfile = mean(z); +f2 = figure('position',[39 346 284 73]); +plot(xProfile) +title('horizontal profile') +axis tight + +%% Estimate spot spacing by autocorrelation +% Ideally the spots would be periodicaly spaced consistently printed, but +% in practice they tend to have different sizes and intensities, so the +% horizontal profile is irregular. We can use autocorrelation to enhance +% the self similarity of the profile. The smooth result promotes peak +% finding and estimation of spot spacing. The *Signal Processing Toolbox* +% allows easy computation of the autocorrelation function using the |xcov| +% command. + +ac = xcov(xProfile); %unbiased autocorrelation +f3 = figure('position',[-3 427 569 94]); +plot(ac) +s1 = diff(ac([1 1:end])); %left slopes +s2 = diff(ac([1:end end])); %right slopes +maxima = find(s1>0 & s2<0); %peaks +estPeriod = round(median(diff(maxima))) %nominal spacing +hold on +plot(maxima,ac(maxima),'r^') +hold off +title('autocorrelation of profile') +axis tight + +%% Remove background morphologically +% We can use the spacing estimate to help design a filter to remove the +% background noise from the intensity profile. We do this with the +% |imtophat| function from the *Image Processing Toolbox*. The |strel| +% command creates a simple rectangular 1D window or line shaped structuring +% element. +seLine = strel('line',estPeriod,0); +xProfile2 = imtophat(xProfile,seLine); +f4 = figure('position',[40 443 285 76]); +plot(xProfile2) +title('enhanced horizontal profile') +axis tight + +%% Segment peaks +% Now that we have clean and anchored gaps between the peaks, we can number +% each peak region with the |bwlabel| command. These regions were segmented +% by thresholding with |im2bw|. The threshold value was automatically +% determined by statistical properties of the data using |graythresh|. This +% is a good example of image processing techniques are often useful for 1D +% data analysis. +level = graythresh(xProfile2/255)*255 +bw = im2bw(xProfile2/255,level/255); +L = bwlabel(bw); +f5 = figure('position',[40 540 285 70]); +plot(L) +axis tight +title('labelled regions') + +%% Locate centers +% We can extract the centroids of the peaks. These correspond to the +% horizontal centres of the spots. This is a common blob analysis or +% feature extraction task that can be done with |regionprops|. +stats = regionprops(L); +centroids = [stats.Centroid]; +xCenters = centroids(1:2:end) +figure(f5) +hold on +plot(xCenters,1:max(L),'ro') +hold off +title('region centers') + +%% Determine divisions between spots +% The midpoints between adjacent peaks provides grid point locations. +gap = diff(xCenters)/2; +first = xCenters(1)-gap(1); +xGrid = round([first xCenters(1:end)+gap([1:end end])]) +figure(f2) +for i=1:length(xGrid) + line(xGrid(i)*[1 1],ylim,'color','m') +end +title('vertical separators') + +%% Transpose and repeat +% We just did the analysis on the vertical grid. Now we want to do the same +% for the horizontal spacing. To do this, we simply transpose the image and +% repeat all the steps used above. This time without intermediate graphics +% display commands in order to summarize the mathematical steps of this +% algorithm. + +yProfile = mean(z'); %peak profile +ac = xcov(yProfile); %cross correlation +p1 = diff(ac([1 1:end])); +p2 = diff(ac([1:end end])); +maxima = find(p1>0 & p2<0); %peak locations +estPeriod = round(median(diff(maxima))) %spacing estimate +seLine = strel('line',estPeriod,0); +yProfile2 = imtophat(yProfile,seLine); %background removed +level = graythresh(yProfile2/255); %automatic threshold level +bw = im2bw(yProfile2/255,level); %binarized peak regions +L = bwlabel(bw); %labeled regions +stats = regionprops(L); +centroids = [stats.Centroid]; %centroids +yCenters = centroids(1:2:end) %Y parts only +gap = diff(yCenters)/2; %inner region half widths +first = yCenters(1)-gap(1); + +% list defining vertical boundaries between spot regions +yGrid = round([first yCenters(1:end)+gap([1:end end])]) + +%% Put bounding boxes around each spot +% We have now found the rectangular grid. Using pairs of neighboring grid +% points we can form bounding box regions to address each spot +% individually. The position and size coordinates of each bounding box were +% tabulated for convenience into a 4-column matrix called |ROI|, which +% stands for regions of interest. +% figure(f1) +figure() +imshow(z) +line(xGrid'*[1 1],yGrid([1 end]),'color','b') +line(xGrid([1 end]),yGrid'*[1 1],'color','b') +[X,Y] = meshgrid(xGrid(1:end-1),yGrid(1:end-1)); %xGridyGridд洢˷ָӵ +[dX,dY] = meshgrid(diff(xGrid),diff(yGrid)); +ROI = [X(:) Y(:) dX(:) dY(:)]; +% first few rows of ROI table +ROI(1:5,:) +%% +%% Ӵϣȷλ +I=y; +% location_row=xGrid; +% location_col=yGrid; + +location_row=yGrid; +location_col=xGrid; + +for ii = 1:length(location_row)-1 + for jj = 1: length(location_col)-1 + subimage = I([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:);%Iȡÿgrid + subimage_double = double(subimage); + + %%%% תݵ[0,1]֮䣬Ȼת[0,255]֮ + sumimage_01 = (subimage_double-min(subimage_double(:)))./max(subimage_double(:))-min(subimage_double(:)); %ÿgridת[0,1]֮ +% subimage_norm = uint8(255*sumimage_01);%ԭ + subimage_norm = uint8(255*sumimage_01); + + %%%% TVȥ + subimage_norm(:,:,1)= tvdenoise(double(subimage_norm(:,:,1)),0.01); + subimage_norm(:,:,2)= tvdenoise(double(subimage_norm(:,:,2)),0.01); + subimage_norm(:,:,3)= tvdenoise(double(subimage_norm(:,:,3)),0.01); + + I_norm([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:)=subimage_norm;%ÿgrid + + + +%% + %CVзָ + + if mean(mean(subimage_norm(:,:,1)))<5 %%ֵС5ǿ5ͼ + subimage_norm=5*subimage_norm; + else if mean(mean(subimage_norm(:,:,1)))<30 %%ֵС30ǿ1.5ͼ + subimage_norm=1.5*subimage_norm; + end + end + + iter=500; + dt=0.1; + u0=cvseg(subimage_norm,iter,dt); %C-Vˮƽָ +% u1=im2bw(u0); + + %ֵĵΪĬֵ0 + [m_u0,n_u0]=size(u0); + for i=1:m_u0 + for j=1:n_u0 + if(isnan(u0(i,j))) + u0(i,j)=0; + end + end + end +% subimage_bw=u0; + %һеĵһһΪɫȡ + u1=u0; + if(u1(end,1)==1 || u1(end,end)) + u2=1-u1; + else + u2=u1; + end + subimage_bw=u2; + +% subimage_bwΪõĶֵͼ񣬼Ϊݵһͼ + %% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ָ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + I_bw([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:)=subimage_bw;%ÿgrid; + + + end +end +figure,imshow(I_bw) +%ʾԭͼѧȵĶֵͼ +all_bw=zeros(size(y(:,:,1))); +all_bw(1:location_row(end),1:location_col(end))=I_bw; +figure(),imshow(all_bw) +%% +I_bw1=all_bw; +I_bw_01=choice(I_bw1,100); %%޳100Ŀ꣨100Ϊ趨ֵ,ֱͼķȷֵ +figure(),imshow(I_bw_01); +%% +num_sub_I_bw_02=0; +for ii = 1:length(location_row)-1 + for jj = 1: length(location_col)-1 + sub_I_bw_01 = I_bw_01([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:); + sub_I_bw_02=choosemaxobj(sub_I_bw_01,8);%޳С20ĵ + + [m,n]=size(sub_I_bw_02); + for i=1:m + for j=1:n + if(sub_I_bw_02(1,j)==1 || sub_I_bw_02(i,1)==1|| sub_I_bw_02(m,j)==1 || sub_I_bw_02(i,n)==1) + num_sub_I_bw_02=num_sub_I_bw_02+1; + end + end + end + + if(num_sub_I_bw_02>10) + sub_I_bw_02=zeros(size(sub_I_bw_02)); + end + + num_sub_I_bw_02=0;%ÿμ + I_bw_Last_01([location_row(ii):location_row(ii+1)],[location_col(jj):location_col(jj+1)],:)=sub_I_bw_02;%ÿgrid; + end +end +all_bw_01=zeros(size(y(:,:,1))); +all_bw_01(1:location_row(end),1:location_col(end))=I_bw_Last_01; +figure(),imshow(all_bw_01); diff --git a/参考资料/NewGridAndCV/fillingholes.m b/参考资料/NewGridAndCV/fillingholes.m new file mode 100644 index 0000000..46c17f9 --- /dev/null +++ b/参考资料/NewGridAndCV/fillingholes.m @@ -0,0 +1,5 @@ +function Lout=fillingholes(Linput,neight) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%洦ȡԭͼķȻȡ壬Ȼ +Lout=choosemaxobj(1-Linput,neight); +Lout=1-Lout; \ No newline at end of file diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy.html b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy.html new file mode 100644 index 0000000..702e68a --- /dev/null +++ b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy.html @@ -0,0 +1,892 @@ + + + + + + Microarray Spot Finding Example + + + + +

Microarray Spot Finding Example

+ +

This example shows a simple method for locating spots on a microarray and extracting the intensties of the spots. It can be + downloaded from MATLAB Central. http://www.mathworks.com/matlabcentral

+ +

Contents

+
+ +
+

Start with clean slate

clear           %empty workspace (no variables)
+close all       %no figures
+clc             %empty command window
+

Read image file

+

MATLAB can read many standard image formats including TIFF, GIF and BMP using the imread command. In addition, the Image Procesing Toolbox provides support for working with specialized image file formats such as DICOM. This microarray image was stored as a J-PEG + file. The image is much larger than the screen size, so imshow scales it down to fit and let's you know with a warning message. + 

x = imread('MicroArraySlide.JPG');
+imageSize = size(x)
+screenSize = get(0,'ScreenSize')
+
+iptsetpref('ImshowBorder','tight')
+imshow(x)
+title('original image')
+
imageSize =
+        3248        1248           3
+screenSize =
+           1           1        1024         768
+Warning: Image is too big to fit on screen; displaying at 42% scale.
+

Crop specified region

+

Next we use imcrop to extract a region of interest. You can repeat this for all print-tip blocks for a full microarray study. + 

y = imcrop(x,[622 2467 220 227]);
+f1 = figure('position',[40 46 285 280]);
+imshow(y)
+

Display red & green layers

+

This image was stored in RGB format. We are only interested in the red and green planes. To extract the red plane, simply + index layer 1. For the green plane, layer 2. Custom colormaps make visualization more intuitive. Notice that spot shapes are + not necessarily the same in both colors. + 

f2 = figure('position',[265 163 647 327]);
+subplot(121)
+redMap = gray(256);
+redMap(:,[2 3]) = 0;
+subimage(y(:,:,1),redMap)
+axis off
+title('red (layer 1)')
+subplot(122)
+greenMap = gray(256);
+greenMap(:,[1 3]) = 0;
+subimage(y(:,:,2),greenMap)
+axis off
+title('green (layer 2)')
+

Convert RGB image to grayscale for spot finding

+

Initially we care more about where the spots are located than their red and green intensities. Converting from RGB color to + grayscale allows us to focus first on spot locations. + 

z = rgb2gray(y);
+figure(f1)
+imshow(z)
+

Create horizontal profile

+

We are looking for a regular grid of spots so we start by looking at the mean intensity for each column of the image. This + will help us identify where the centres of the spots are and where the gaps between the spots can be found. + 

xProfile = mean(z);
+f2 = figure('position',[39 346 284 73]);
+plot(xProfile)
+title('horizontal profile')
+axis tight
+

Estimate spot spacing by autocorrelation

+

Ideally the spots would be periodicaly spaced consistently printed, but in practice they tend to have different sizes and + intensities, so the horizontal profile is irregular. We can use autocorrelation to enhance the self similarity of the profile. + The smooth result promotes peak finding and estimation of spot spacing. The Signal Processing Toolbox allows easy computation of the autocorrelation function using the xcov command. + 

ac = xcov(xProfile);                        %unbiased autocorrelation
+f3 = figure('position',[-3 427 569 94]);
+plot(ac)
+s1 = diff(ac([1 1:end]));                   %left slopes
+s2 = diff(ac([1:end end]));                 %right slopes
+maxima = find(s1>0 & s2<0);                 %peaks
+estPeriod = round(median(diff(maxima)))     %nominal spacing
+hold on
+plot(maxima,ac(maxima),'r^')
+hold off
+title('autocorrelation of profile')
+axis tight
+
estPeriod =
+    19
+

Remove background morphologically

+

We can use the spacing estimate to help design a filter to remove the background noise from the intensity profile. We do this + with the imtophat function from the Image Processing Toolbox. The strel command creates a simple rectangular 1D window or line shaped structuring element. + 

seLine = strel('line',estPeriod,0);
+xProfile2 = imtophat(xProfile,seLine);
+f4 = figure('position',[40 443 285 76]);
+plot(xProfile2)
+title('enhanced horizontal profile')
+axis tight
+

Segment peaks

+

Now that we have clean and anchored gaps between the peaks, we can number each peak region with the bwlabel command. These regions were segmented by thresholding with im2bw. The threshold value was automatically determined by statistical properties of the data using graythresh. This is a good example of image processing techniques are often useful for 1D data analysis. + 

level = graythresh(xProfile2/255)*255
+bw = im2bw(xProfile2/255,level/255);
+L = bwlabel(bw);
+f5 = figure('position',[40 540 285 70]);
+plot(L)
+axis tight
+title('labelled regions')
+
level =
+    16
+

Locate centers

+

We can extract the centroids of the peaks. These correspond to the horizontal centres of the spots. This is a common blob + analysis or feature extraction task that can be done with regionprops. + 

stats = regionprops(L);
+centroids = [stats.Centroid];
+xCenters = centroids(1:2:end)
+figure(f5)
+hold on
+plot(xCenters,1:max(L),'ro')
+hold off
+title('region centers')
+
xCenters =
+  Columns 1 through 8 
+   23.0000   41.0000   60.0000   79.0000   98.0000  119.5000  137.0000  156.0000
+  Columns 9 through 10 
+  175.0000  193.5000
+

Determine divisions between spots

+

The midpoints between adjacent peaks provides grid point locations.

gap = diff(xCenters)/2;
+first = xCenters(1)-gap(1);
+xGrid = round([first xCenters(1:end)+gap([1:end end])])
+figure(f2)
+for i=1:length(xGrid)
+  line(xGrid(i)*[1 1],ylim,'color','m')
+end
+title('vertical separators')
+
xGrid =
+    14    32    51    70    89   109   128   147   166   184   203
+

Transpose and repeat

+

We just did the analysis on the vertical grid. Now we want to do the same for the horizontal spacing. To do this, we simply + transpose the image and repeat all the steps used above. This time without intermediate graphics display commands in order + to summarize the mathematical steps of this algorithm. + 

yProfile = mean(z');                        %peak profile
+ac = xcov(yProfile);                        %cross correlation
+p1 = diff(ac([1 1:end]));
+p2 = diff(ac([1:end end]));
+maxima = find(p1>0 & p2<0);                 %peak locations
+estPeriod = round(median(diff(maxima)))     %spacing estimate
+seLine = strel('line',estPeriod,0);
+yProfile2 = imtophat(yProfile,seLine);      %background removed
+level = graythresh(yProfile2/255);          %automatic threshold level
+bw = im2bw(yProfile2/255,level);            %binarized peak regions
+L = bwlabel(bw);                            %labeled regions
+stats = regionprops(L);
+centroids = [stats.Centroid];               %centroids
+yCenters = centroids(1:2:end)               %Y parts only
+gap = diff(yCenters)/2;                     %inner region half widths
+first = yCenters(1)-gap(1);
+
+% list defining vertical boundaries between spot regions
+yGrid = round([first yCenters(1:end)+gap([1:end end])])
+
estPeriod =
+    20
+yCenters =
+  Columns 1 through 8 
+   24.0000   43.5000   64.0000   83.5000  104.5000  123.5000  144.5000  164.0000
+  Columns 9 through 10 
+  183.0000  203.5000
+yGrid =
+    14    34    54    74    94   114   134   154   174   193   214
+

Put bounding boxes around each spot

+

We have now found the rectangular grid. Using pairs of neighboring grid points we can form bounding box regions to address + each spot individually. The position and size coordinates of each bounding box were tabulated for convenience into a 4-column + matrix called ROI, which stands for regions of interest. + 

figure(f1)
+imshow(z)
+line(xGrid'*[1 1],yGrid([1 end]),'color','b')
+line(xGrid([1 end]),yGrid'*[1 1],'color','b')
+[X,Y] = meshgrid(xGrid(1:end-1),yGrid(1:end-1));
+[dX,dY] = meshgrid(diff(xGrid),diff(yGrid));
+ROI = [X(:) Y(:) dX(:) dY(:)];
+% first few rows of ROI table
+ROI(1:5,:)
+
ans =
+    14    14    18    20
+    14    34    18    20
+    14    54    18    20
+    14    74    18    20
+    14    94    18    20
+

Segment spots from background by thresholding

+

Applying a single threshold level to the whole image so all spots are detected equally is generally a good idea. However, + in this case is doesn't work so well due to large differences in spot brightness. + 

fSpots = figure('position',[265 163 647 327]);
+subplot(121)
+imshow(z)
+title('gray image')
+subplot(122)
+bw = im2bw(z,graythresh(z));
+imshow(bw)
+title('global threshold')
+

Apply logarithmic transformation then threshold intensities

+

One way to equalize large variations in magnitude is by transforming intensity values to logarithmic space. This works much + better but some weak spots are still missed. + 

figure(fSpots)
+subplot(121)
+z2 = uint8(log(double(z)+1)/log(255)*255);
+imshow(z2)
+title('log intensity')
+subplot(122)
+bw = im2bw(z2,graythresh(z2));
+imshow(bw)
+title('global threshold')
+
Warning: Conversion rounded non-integer floating point value to nearest uint8 value.
+

Try local thresholding instead

+

Alternatively, the bounding boxes can be used to determine local threshold values for each spot. The code is a little more + sophisticated, requiring looping and indexing. Unfortunately, the results are mixed. Weak spots showed up well but spots with + bright perimeters were as bad as the original global threshold before log space transformation. + 

figure(fSpots)
+subplot(122)
+bw = false(size(z));
+for i=1:length(ROI)
+  rows = round(ROI(i,2))+[0:(round(ROI(i,4))-1)];
+  cols = round(ROI(i,1))+[0:(round(ROI(i,3))-1)];
+  spot = z(rows,cols);
+  bw(rows,cols) = im2bw(spot,graythresh(spot));
+end
+imshow(bw)
+title('local threshold')
+

Logically combine local and global thresholds

+

Since both have their merits, let's combine the best of both approaches. This can be done using logial operation on the binary + masks. These spot segmentation results are indeed much better. + 

figure(fSpots)
+subplot(121)
+bw = im2bw(z2,graythresh(z2));
+for i=1:length(ROI)
+  rows = round(ROI(i,2))+[0:(round(ROI(i,4))-1)];
+  cols = round(ROI(i,1))+[0:(round(ROI(i,3))-1)];
+  spot = z(rows,cols);
+  bw(rows,cols) = bw(rows,cols) | im2bw(spot,graythresh(spot));
+end
+imshow(bw)
+title('combined threshold')
+subplot(122)
+imshow(z)
+title('linear intensity')
+

Fill holes to solidify spots

+

The silhouettes of some spots still contained pinholes. The whole image could be filled using a single call to imfill but this may not be a good idea. Notice that some spots run together. If four mutually adjacent spots (sharing a common corner) + were all joined at their edges then a single function call would incorrectly fill in the common corner as well. To avoid that + possibility, it's good insurance to fill each spot one bounding box region at a time by looping. Indeed, the spot segmentation + now looks quite good. + 

figure(fSpots)
+subplot(121)
+warning off MATLAB:intConvertOverflow
+for i=1:length(ROI)
+  rows = round(ROI(i,2))+[0:(round(ROI(i,4))-1)];
+  cols = round(ROI(i,1))+[0:(round(ROI(i,3))-1)];
+  bw(rows,cols) = imfill(bw(rows,cols),'holes');
+end
+imshow(bw)
+title('filled pinholes')
+

Label spot masks by bounding box

+

If the gridding went well, all spots should be a single color. The results here are pretty good. There is still room for improvement.

+

TODO List:

+
+ +
+
+ +
+
+ +
+

However, in this case the spot segmentation is good enough to proceed.

L = zeros(size(bw));
+for i=1:length(ROI)
+  rows = ROI(i,2)+[0:(ROI(i,4)-1)];
+  cols = ROI(i,1)+[0:(ROI(i,3)-1)];
+  rectMask = L(rows,cols);
+  spotMask = bw(rows,cols);
+  rectMask(spotMask) = i;
+  L(rows,cols) = rectMask;
+end
+map = [0 0 0; 0.5+0.5*rand(length(ROI),3)];
+figure(f1)
+imshow(L+1,map)
+

Extract first spot for measurement

+

We will now examine the first spot closely to see how we can measure its red and green intensities, and ultimately quantify + its gene expression value. The measurement technique can then be repeated for all spots. + 

rect = ROI(1,:);                            %[X Y dX dY]
+spot = imcrop(y,rect);                      %region around spot
+figure(f1)
+imshow(spot,'notruesize')
+

Measure spot intensity & releative expression level

+

We now simply calculate the nominal intensity over the spot for both the red and green layers. A measure of gene expression + level can then be calculated from the two color intensities. Here a simple log-ratio measurement is shown. Other more robust + measures could be used instead. You could also perform some analysis of the quality of the spot. + 

mask = imcrop(L,rect)==1;
+for i=1:2
+  layer = spot(:,:,i);
+  intensity(i) = double(median(layer(mask)));
+end
+intensity
+expressionLevel = log(intensity(1)/intensity(2))
+
intensity =
+    32    81
+expressionLevel =
+   -0.9287
+

Remove background, calculate again and compare measurements

+

If you noticed, the background intensity around the spot was not zero. This could bias results. To see how much difference + it makes, we can perform background subtraction around all spots, again using imtophat but this time in 2D on the image using a disk shaped structuring element. Then we can calculate color intensities and relative + expression level again to see what effect background bias had on the measurement. In this case the measurement shows more + downregulation with background removed. + 

seDisk = strel('disk',round(estPeriod));
+spot2 = imtophat(spot,seDisk);
+for i=1:2
+  layer = spot2(:,:,i);
+  intensity(i) = double(median(layer(mask)));
+end
+intensity
+expressionLevel = log(intensity(1)/intensity(2))
+
intensity =
+    14    70
+expressionLevel =
+   -1.6094
+

Set up graphical display for results

+

It is helpful to see red and green intensity values overlayed onto the respective color images to gain confidence that measured + intensities make sense. It is also be helpful to overlay quantitative expression levels onto the original image to provide + additional visual assurance of measurement results. The rectangular grid also helps correlate measured values between images. + The flexibility of MATLAB's powerful Handle Graphics engine allow custom graphics like this to be set up quickly and easily. + 

f7 = figure('position',[52 94 954 425]);
+ax(1) = subplot(121);
+subimage(y(:,:,1),redMap)
+title('red intensity')
+ax(2) = subplot(122);
+subimage(y(:,:,2),greenMap)
+title('green intensity')
+f8 = figure('position',[316 34 482 497]);
+ax(3) = get(imshow(y,'notruesize'),'parent');
+title('gene expression')
+for i=1:3
+  axes(ax(i))
+  axis off
+  line(xGrid'*[1 1],yGrid([1 end]),'color',0.5*[1 1 1])
+  line(xGrid([1 end]),yGrid'*[1 1],'color',0.5*[1 1 1])
+end
+

Repeat measurement for all spots

+

We now repeat the spot extraction and intensity calculation for all the spots in the grid. Here the measured values were tabulated + as additional columns beside the ROI positions for each spot into a new matrix called spotData. + 

figure(f7), figure(f8)
+spotData = [ROI zeros(length(ROI),5)];
+for i=1:length(ROI)
+  spot = imcrop(y,ROI(i,:));                            %raw image
+  spot2 = imtophat(spot,seDisk);                        %background removed
+  mask = imcrop(L,ROI(i,:))==i;                         %spot mask
+  for j=1:2
+    layer = spot2(:,:,j);                               %color layer
+    intensity(j) = double(median(layer(mask)));
+    text(ROI(i,1)+ROI(i,3)/2,ROI(i,2)+ROI(i,4)/2,sprintf('%.0f',intensity(j)),...
+      'color','y','HorizontalAlignment','center','parent',ax(j))
+    rawLayer = spot(:,:,j);
+    rawIntensity(j) = double(median(layer(mask)));
+  end
+  expression = log(intensity(1)/intensity(2));
+  text(ROI(i,1)+ROI(i,3)/2,ROI(i,2)+ROI(i,4)/2,sprintf('%.2f',expression),...
+    'color','w','HorizontalAlignment','center','parent',ax(3))
+  drawnow
+  spotData(i,5:9) = [intensity(:)' expression rawIntensity(:)'];
+end
+

Export spot data to Excel spreadsheet

+

MATLAB can write to many standard formats. We will use xlswrite to save the spotData to an Excel workbook. +

+

TODO list:

+
+ +
+
+ + 
xlswrite('microarray.xls',spotData)
+


+ Published with MATLAB® 7.0

+ + + \ No newline at end of file diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy.png new file mode 100644 index 0000000..2d78f2e Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_01.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_01.png new file mode 100644 index 0000000..3567e3f Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_01.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_02.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_02.png new file mode 100644 index 0000000..0577176 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_02.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_03.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_03.png new file mode 100644 index 0000000..12b5dcb Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_03.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_04.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_04.png new file mode 100644 index 0000000..29a1768 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_04.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_05.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_05.png new file mode 100644 index 0000000..28b3b57 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_05.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_06.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_06.png new file mode 100644 index 0000000..a352f79 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_06.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_07.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_07.png new file mode 100644 index 0000000..f0ff276 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_07.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_08.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_08.png new file mode 100644 index 0000000..28fea02 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_08.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_09.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_09.png new file mode 100644 index 0000000..f9c1310 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_09.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_10.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_10.png new file mode 100644 index 0000000..7b4dd05 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_10.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_11.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_11.png new file mode 100644 index 0000000..473b5ed Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_11.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_12.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_12.png new file mode 100644 index 0000000..a39ef4d Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_12.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_13.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_13.png new file mode 100644 index 0000000..1d54291 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_13.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_14.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_14.png new file mode 100644 index 0000000..ef6eed7 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_14.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_15.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_15.png new file mode 100644 index 0000000..2fe82a0 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_15.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_16.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_16.png new file mode 100644 index 0000000..5c60c5e Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_16.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_17.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_17.png new file mode 100644 index 0000000..dffb555 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_17.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_18.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_18.png new file mode 100644 index 0000000..d16f22b Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_18.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_19.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_19.png new file mode 100644 index 0000000..9855aa9 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_19.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_20.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_20.png new file mode 100644 index 0000000..d2baa59 Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_20.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_21.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_21.png new file mode 100644 index 0000000..0dff13c Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_21.png differ diff --git a/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_22.png b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_22.png new file mode 100644 index 0000000..7caf20a Binary files /dev/null and b/参考资料/NewGridAndCV/html/R14_MicroarrayImage_CaseStudy_22.png differ diff --git a/参考资料/NewGridAndCV/kappa.m b/参考资料/NewGridAndCV/kappa.m new file mode 100644 index 0000000..8bd934e --- /dev/null +++ b/参考资料/NewGridAndCV/kappa.m @@ -0,0 +1,28 @@ +function KG = kappa(I) +% get curvature information of input image +% input: 2D image I +% output: curvature matrix KG + +% Copyright (c) 2009, +% Yue Wu @ ECE Department, Tufts University +% All Rights Reserved + +I = double(I); +[m,n] = size(I); +P = padarray(I,[1,1],1,'pre'); +P = padarray(P,[1,1],1,'post'); + +% central difference +fy = P(3:end,2:n+1)-P(1:m,2:n+1); +fx = P(2:m+1,3:end)-P(2:m+1,1:n); +fyy = P(3:end,2:n+1)+P(1:m,2:n+1)-2*I; +fxx = P(2:m+1,3:end)+P(2:m+1,1:n)-2*I; +fxy = 0.25.*(P(3:end,3:end)-P(1:m,3:end)+P(3:end,1:n)-P(1:m,1:n)); +G = (fx.^2+fy.^2).^(0.5); +K = (fxx.*fy.^2-2*fxy.*fx.*fy+fyy.*fx.^2)./((fx.^2+fy.^2+eps).^(1.5)); +KG = K.*G; +KG(1,:) = eps; +KG(end,:) = eps; +KG(:,1) = eps; +KG(:,end) = eps; +KG = KG./max(max(abs(KG))); diff --git a/参考资料/NewGridAndCV/limits.m b/参考资料/NewGridAndCV/limits.m new file mode 100644 index 0000000..69478b9 --- /dev/null +++ b/参考资料/NewGridAndCV/limits.m @@ -0,0 +1,29 @@ +function [x,y]=limits(a) +% LIMITS returns min & max values of matrix; else scalar value. +% +% [lo,hi]=LIMITS(a) returns LOw and HIgh values respectively. +% +% lim=LIMITS(a) returns 1x2 result, where lim = [lo hi] values + +% Copyright 2004-2010 RBemis The MathWorks, Inc. + +if nargin~=1 | nargout>2 %bogus syntax + error('usage: [lo,hi]=limits(a)') +end + +siz=size(a); + +if prod(siz)==1 %scalar + result=a; % value +else %matrix + result=[min(a(:)) max(a(:))]; % limits +end + +if nargout==1 %composite result + x=result; % 1x2 vector +elseif nargout==2 %separate results + x=result(1); % two scalars + y=result(2); +else %no result + ans=result % display answer +end diff --git a/参考资料/NewGridAndCV/maskcircle2.m b/参考资料/NewGridAndCV/maskcircle2.m new file mode 100644 index 0000000..07d750f --- /dev/null +++ b/参考资料/NewGridAndCV/maskcircle2.m @@ -0,0 +1,75 @@ + +function m = maskcircle2(I,type) +% auto pick a circular mask for image I +% built-in mask creation function +% Input: I : input image +% type: mask shape keywords +% Output: m : mask image + +% Copyright (c) 2009, +% Yue Wu @ ECE Department, Tufts University +% All Rights Reserved + +if size(I,3)~=3 + temp = double(I(:,:,1)); +else + temp = double(rgb2gray(I)); +end + +h = [0 1 0; 1 -4 1; 0 1 0]; +T = conv2(temp,h); +T(1,:) = 0; +T(end,:) = 0; +T(:,1) = 0; +T(:,end) = 0; + +thre = max(max(abs(T)))*.5; +idx = find(abs(T) > thre); +[cx,cy] = ind2sub(size(T),idx); +cx = round(mean(cx)); +cy = round(mean(cy)); + +[x,y] = meshgrid(1:min(size(temp,1),size(temp,2))); + +m = zeros(size(temp)); +[p,q] = size(temp); + +switch lower (type) + case 'small' + r = 10; + n = zeros(size(x)); + n((x-cx).^2+(y-cy).^21, error('USAGE: r=range(x)'), end + +if nargout>0 + r=diff(limits(x)); +else + diff(limits(x)) +end diff --git a/参考资料/NewGridAndCV/redcolorcontour.m b/参考资料/NewGridAndCV/redcolorcontour.m new file mode 100644 index 0000000..a078f84 --- /dev/null +++ b/参考资料/NewGridAndCV/redcolorcontour.m @@ -0,0 +1,19 @@ +function output=redcolorcontour(inputcolor,inputcontour) +%%ͼbwе1ʾΪɫʾڲɫͼĶӦλ +%% + +[r,c]=size(inputcontour); +% inputcolor + +for i=1:r + for j=1:c + if uint8(inputcontour(i,j))==255 + inputcolor(i,j)=255; + inputcolor(i,j,2)=0; + inputcolor(i,j,3)=0; + end + end +end +output(:,:,1)=inputcolor(:,:,1); +output(:,:,2)=inputcolor(:,:,2); +output(:,:,3)=inputcolor(:,:,3); \ No newline at end of file diff --git a/参考资料/NewGridAndCV/rgblab.m b/参考资料/NewGridAndCV/rgblab.m new file mode 100644 index 0000000..99bf4e1 --- /dev/null +++ b/参考资料/NewGridAndCV/rgblab.m @@ -0,0 +1,14 @@ +function lab=rgblab(rgb) + + +%%תΪlabռ +C = makecform('srgb2lab'); +lab = applycform(rgb,C); + +% %%ʾlabռͼͨ +% L=lab(:,:,1); +% a=lab(:,:,2); +% b=lab(:,:,3); +% figure,imshow(lab(:,:,1)),title('L') +% figure,imshow(lab(:,:,2)),title('a') +% figure,imshow(lab(:,:,3)),title('b') diff --git a/参考资料/NewGridAndCV/showphi.m b/参考资料/NewGridAndCV/showphi.m new file mode 100644 index 0000000..b800542 --- /dev/null +++ b/参考资料/NewGridAndCV/showphi.m @@ -0,0 +1,26 @@ +function showphi(I, phi, i) +% show curve evolution of phi + +% Copyright (c) 2009, +% Yue Wu @ ECE Department, Tufts University +% All Rights Reserved + +for j = 1:size(phi,3) + phi_{j} = phi(:,:,j); +end +% imshow(I,'initialmagnification','fit','displayrange',[0 255]); +% hold on; + + if size(phi,3) == 1 + contour(phi_{1}, [0 0], 'r','LineWidth',4); + contour(phi_{1}, [0 0], 'g','LineWidth',1.3); + else + contour(phi_{1}, [0 0], 'r','LineWidth',4); + contour(phi_{1}, [0 0], 'x','LineWidth',1.3); + contour(phi_{2}, [0 0], 'g','LineWidth',4); + contour(phi_{2}, [0 0], 'x','LineWidth',1.3); + end +% hold off; +% title([num2str(i) ' Iterations']); + + drawnow; \ No newline at end of file diff --git a/参考资料/NewGridAndCV/tvdenoise.m b/参考资料/NewGridAndCV/tvdenoise.m new file mode 100644 index 0000000..73cffc4 --- /dev/null +++ b/参考资料/NewGridAndCV/tvdenoise.m @@ -0,0 +1,85 @@ +function u = tvdenoise(f,lambda,Tol) +%TVDENOISE Total variation grayscale and color image denoising +% u = TVDENOISE(f,lambda) denoises the input image f. The smaller +% the parameter lambda, the stronger the denoising. +% +% The output u approximately minimizes the Rudin-Osher-Fatemi (ROF) +% denoising model +% +% Min TV(u) + lambda/2 || f - u ||^2_2, +% u +% +% where TV(u) is the total variation of u. If f is a color image (or any +% array where size(f,3) > 1), the vectorial TV model is used, +% +% Min VTV(u) + lambda/2 || f - u ||^2_2. +% u +% +% TVDENOISE(...,Tol) specifies the stopping tolerance (default 1e-2). +% +% The minimization is solved using Chambolle's method, +% A. Chambolle, "An Algorithm for Total Variation Minimization and +% Applications," J. Math. Imaging and Vision 20 (1-2): 89-97, 2004. +% When f is a color image, the minimization is solved by a generalization +% of Chambolle's method, +% X. Bresson and T.F. Chan, "Fast Minimization of the Vectorial Total +% Variation Norm and Applications to Color Image Processing", UCLA CAM +% Report 07-25. +% +% Example: +% f = double(imread('barbara-color.png'))/255; +% f = f + randn(size(f))*16/255; +% u = tvdenoise(f,12); +% subplot(1,2,1); imshow(f); title Input +% subplot(1,2,2); imshow(u); title Denoised + +% Pascal Getreuer 2007-2008 + +if nargin < 3 + Tol = 1e-2; +end + +if lambda < 0 + error('Parameter lambda must be nonnegative.'); +end + +dt = 0.25; + +N = size(f); +id = [2:N(1),N(1)]; +iu = [1,1:N(1)-1]; +ir = [2:N(2),N(2)]; +il = [1,1:N(2)-1]; +p1 = zeros(size(f)); +p2 = zeros(size(f)); +divp = zeros(size(f)); +lastdivp = ones(size(f)); + +if length(N) == 2 % TV denoising + while norm(divp(:) - lastdivp(:),inf) > Tol + lastdivp = divp; + z = divp - f*lambda; + z1 = z(:,ir) - z; + z2 = z(id,:) - z; + denom = 1 + dt*sqrt(z1.^2 + z2.^2); + p1 = (p1 + dt*z1)./denom; + p2 = (p2 + dt*z2)./denom; + divp = p1 - p1(:,il) + p2 - p2(iu,:); + end +elseif length(N) == 3 % Vectorial TV denoising + repchannel = ones(N(3),1); + + while norm(divp(:) - lastdivp(:),inf) > Tol + lastdivp = divp; + z = divp - f*lambda; + z1 = z(:,ir,:) - z; + z2 = z(id,:,:) - z; + denom = 1 + dt*sqrt(sum(z1.^2 + z2.^2,3)); + denom = denom(:,:,repchannel); + p1 = (p1 + dt*z1)./denom; + p2 = (p2 + dt*z2)./denom; + divp = p1 - p1(:,il,:) + p2 - p2(iu,:,:); + end +end + +u = f - divp/lambda; diff --git a/参考资料/NewGridAndCV/低对比度cDNA图像分割的局部水平集方法_芦碧波_刘利群_张霄宏_林忠华.pdf b/参考资料/NewGridAndCV/低对比度cDNA图像分割的局部水平集方法_芦碧波_刘利群_张霄宏_林忠华.pdf new file mode 100644 index 0000000..29255b4 Binary files /dev/null and b/参考资料/NewGridAndCV/低对比度cDNA图像分割的局部水平集方法_芦碧波_刘利群_张霄宏_林忠华.pdf differ diff --git a/参考资料/封面+低对比度cDNA图像分割的局部水平集方法_芦碧波.pdf b/参考资料/封面+低对比度cDNA图像分割的局部水平集方法_芦碧波.pdf new file mode 100644 index 0000000..599c38d Binary files /dev/null and b/参考资料/封面+低对比度cDNA图像分割的局部水平集方法_芦碧波.pdf differ diff --git a/参考资料/显微图像分割.pdf b/参考资料/显微图像分割.pdf new file mode 100644 index 0000000..40688f2 Binary files /dev/null and b/参考资料/显微图像分割.pdf differ