plot29 · 代码 10 · 图 0 · 笔记 6 · 代码库专题
ECM纤维定量
代码库整理的专题，README 见原始目录。
ECM纤维定量（自动生成占位图）
图源：（自动生成占位图）
分析代码
5 个源码文件，按目录顺序；本机绝对路径已替换为 ~。 代码可复制，但需要按自己的数据改路径与参数。
3.1_原文fork代码/analysis_pca_overlap.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# =============================================================================
# analysis_pca_overlap.py — PCA 重叠敏感性分析（非 Fig. 4H 主图复刻）
# =============================================================================
# 原文正文称 PCA 支持以下结论：Vert 与 Unwounded 的重叠由 POD14 增至
# POD30，并在 POD90 完全重叠；但 Fig. 4H 图注明确该主图为 t-SNE plots。
# 原始 MATLAB 公开脚本也执行 tsne(quantified, ...)，未提供论文中 PCA
# 图/重叠统计的确切特征子集、聚合层级或重叠定义。因此本脚本只能做 PCA
# 敏感性分析，不能据此宣称视觉或定量复刻 Fig. 4H。
#
# 本脚本报告两个探索性指标：
# 1) NN overlap: Vert 组每个样本的 k 近邻中属于 Unwounded 的比例（0~1）
# （对照: PBS 组近邻中属于 Unwounded 的比例，应始终低）
# 2) 质心距: 两组在 PCA 前3主成分空间的质心距离（越大越分离）
#
# 用法:
# python analysis_pca_overlap.py [--csv output/fiber_quantification_features.csv]
# [--dedupe] [--timepoint 1mo|2wk|3mo|all]
# [--k 10] [--low-var 0.25]
# 输出:
# 控制台表格 + output/pca_overlap_{timepoint}.png（PCA 椭圆图）
# =============================================================================
import os
import argparse
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
OUT_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "output")
def clean_X(df, feat_cols, low_var, mean_only=False):
"""inf/NaN 清理 + 低方差过滤（与 run_reproduction.py 保持一致）
mean_only=True 时仅保留每参数的 mean/全局统计列；此选项是敏感性分析，不是原文已证实的 26 特征定义。"""
if mean_only:
feat_cols = [c for c in feat_cols if c.endswith("_mean") or c in
("ori_mean", "ori_median", "ori_std", "ori_skew", "ori_kurt",
"circ_kappa", "glcm_contrast", "glcm_correlation",
"glcm_energy", "glcm_homogeneity", "n_skel", "n_branch")]
X = df[feat_cols].replace([np.inf, -np.inf], np.nan).fillna(0).values
if low_var > 0:
var = np.nanvar(X, axis=0)
keep = var >= np.percentile(var, 100 * low_var)
X = X[:, keep]
return X
def nn_overlap(emb, y, group_a, group_b, k=10):
"""A 组样本的 k 近邻中属于 B 组的比例（0=完全不混，1=完全混合）"""
from sklearn.neighbors import NearestNeighbors
nn = NearestNeighbors(n_neighbors=min(k + 1, len(emb))).fit(emb)
_, idx = nn.kneighbors(emb[y == group_a])
idx = idx[:, 1:] # 去掉自身
return float((y[idx] == group_b).mean())
def overlap_table(emb, df, timepoint_label, k):
"""对给定嵌入与子集计算 NN overlap 表（global PCA 模式复用）"""
rows = []
for a, b in (("Vert", "Unwounded"), ("PBS", "Unwounded"), ("PBS", "Vert")):
if a not in df["group"].values or b not in df["group"].values:
continue
o = nn_overlap(emb, df["group"].values, a, b, k)
c_a = emb[df["group"].values == a].mean(0)
c_b = emb[df["group"].values == b].mean(0)
d = float(np.linalg.norm(c_a - c_b))
rows.append({"timepoint": timepoint_label, "pair": f"{a}-{b}",
"NN_overlap": round(o, 3), "centroid_dist": round(d, 2),
"n_a": int((df["group"] == a).sum()), "n_b": int((df["group"] == b).sum())})
return pd.DataFrame(rows)
def analyze(df, timepoint_label, k, low_var, mean_only=False, n_pc=3):
feat_cols = [c for c in df.columns if c not in ("image", "timepoint", "group", "rep", "channel")]
X = clean_X(df, feat_cols, low_var, mean_only)
Xs = StandardScaler().fit_transform(X)
emb = PCA(n_components=min(n_pc, X.shape[1]), random_state=42).fit_transform(Xs) # 原文用 PCA 定量
rows = []
for a, b in (("Vert", "Unwounded"), ("PBS", "Unwounded"), ("PBS", "Vert")):
if a not in df["group"].values or b not in df["group"].values:
continue
o = nn_overlap(emb, df["group"].values, a, b, k)
c_a = emb[df["group"].values == a].mean(0)
c_b = emb[df["group"].values == b].mean(0)
d = float(np.linalg.norm(c_a - c_b))
rows.append({"timepoint": timepoint_label, "pair": f"{a}-{b}",
"NN_overlap": round(o, 3), "centroid_dist": round(d, 2),
"n_a": int((df["group"] == a).sum()), "n_b": int((df["group"] == b).sum())})
return pd.DataFrame(rows), emb
def plot_pca(emb, df, title, path):
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
fig, ax = plt.subplots(figsize=(6.5, 6))
colors = {"PBS": "#d62728", "Vert": "#1f77b4", "Unwounded": "#2ca02c"}
for g in sorted(df["group"].unique()):
m = (df["group"] == g).values
pts = emb[m, :2]
ax.scatter(pts[:, 0], pts[:, 1], s=10, alpha=0.6, c=colors.get(g, "gray"), label=g)
if len(pts) >= 3: # 95% 置信椭圆
cov = np.cov(pts.T)
evals, evecs = np.linalg.eigh(cov)
ang = np.degrees(np.arctan2(evecs[1, np.argmax(evals)], evecs[0, np.argmax(evals)]))
w, h = 2 * np.sqrt(5.991 * evals) # 2D 95% 卡方临界值 5.991
ax.add_patch(Ellipse(pts.mean(0), w, h, angle=ang, fill=False,
ec=colors.get(g, "gray"), lw=1.8))
ax.set_title(title)
ax.legend()
plt.tight_layout()
plt.savefig(path, dpi=150)
plt.close()
print(f"图已保存: {path}")
def main():
ap = argparse.ArgumentParser(description="复刻原文 PCA overlap 定量（Vert vs Unwounded 随时间重合趋势）")
ap.add_argument("--csv", default=os.path.join(OUT_DIR, "fiber_quantification_features.csv"))
ap.add_argument("--dedupe", action="store_true", help="全行去重（删 MD Vert 重复副本）")
ap.add_argument("--timepoint", default="all", choices=["all", "1mo", "2wk", "3mo"])
ap.add_argument("--k", type=int, default=10, help="近邻数")
ap.add_argument("--low-var", type=float, default=0.25, help="低方差过滤比例（原文25%）")
ap.add_argument("--mean-only", action="store_true",
help="仅保留每参数 mean/全局统计列（探索性敏感性分析，非原文 26 特征的已验证定义）")
ap.add_argument("--n-pc", type=int, default=3, help="PCA 主成分数（原文图 S10/S11 为 PCA 图）")
ap.add_argument("--global-pca", action="store_true",
help="全量样本拟合 PCA 后按时间点切片（跨时间点可比；独立 PCA 为默认）")
ap.add_argument("--channel", type=int, default=None, choices=[1, 2],
help="只分析指定通道（1=成熟纤维红 / 2=未成熟纤维绿；Method 要求独立分析）")
args = ap.parse_args()
df = pd.read_csv(args.csv)
if args.dedupe:
before = len(df)
df = df.drop_duplicates()
print(f"--dedupe: {before} → {len(df)}")
if args.timepoint != "all":
df = df[df["timepoint"] == args.timepoint]
if args.channel:
df = df[df["channel"] == args.channel]
print(f"--channel {args.channel}: 剩 {len(df)} 行")
# 全量 PCA 模式：所有时间点在同一 PCA 空间（跨时间点可比）
if args.global_pca:
feat_cols = [c for c in df.columns if c not in ("image", "timepoint", "group", "rep", "channel")]
Xg = clean_X(df, feat_cols, args.low_var, args.mean_only)
Xsg = StandardScaler().fit_transform(Xg)
emb_g = PCA(n_components=min(args.n_pc, Xsg.shape[1]), random_state=42).fit_transform(Xsg)
all_rows, all_emb = [], None
for tp in (["2wk", "1mo", "3mo"] if args.timepoint == "all" else [args.timepoint]):
sub = df[df["timepoint"] == tp] if "timepoint" in df.columns else df
if len(sub) < 10:
print(f"{tp}: 样本不足 ({len(sub)})，跳过")
continue
if args.global_pca:
emb = emb_g[(df["timepoint"] == tp).values]
res = overlap_table(emb, sub, tp, args.k)
else:
res, emb = analyze(sub, tp, args.k, args.low_var, args.mean_only, args.n_pc)
all_rows.append(res)
tag = "_meanonly" if args.mean_only else ""
tag += "_global" if args.global_pca else ""
tag += f"_ch{args.channel}" if args.channel else ""
plot_pca(emb, sub, f"PCA - {tp} (ECM ultrastructure{' mean-only' if args.mean_only else ''})",
os.path.join(OUT_DIR, f"pca_overlap_{tp}{tag}.png"))
table = pd.concat(all_rows, ignore_index=True)
print("\n=== PCA 探索性重叠度（NN overlap 越高 = 越重合；按 POD14→30→90 顺序）===")
print(table.to_string(index=False))
print("\n→ 原文的 PCA 结论是 POD14→30→90 重叠增强；本 NN overlap 是近似指标，不能替代原文未公开的 overlap 定义。")
if __name__ == "__main__":
main()
3.1_原文fork代码/fiberquantificationmaster.m
%% CJUN Rosa Picro Fiber Segmentation: Skeletonization
clear
clc
close all
% Get list of all TIF files in the directory
imagefiles = dir('*.tif');
nfiles = length(imagefiles);
counter = 0;
% Optional scaling factors (to convert pixels to microns)
scale_1 = 1;
scale_2 = 1;
for i = 1:2:nfiles
redfilename = imagefiles(i).name;
greenfilename = imagefiles(i+1).name;
stored_red_names{i} = redfilename;
stored_green_names{i} = greenfilename;
stored_red_im{i} = imread(redfilename);
stored_green_im{i} = imread(greenfilename);
picro_red = 255-imread(redfilename);
picro_green = 255-imread(greenfilename);
% Remove noise with adaptive filtering if needed
% picro_red = wiener2(picro_red,[3 3]);
% picro_green = wiener2(picro_green,[3 3]);
% Binarize, erode, and dilate
% picro_red_bw = im2bw(imadjust(picro_red,[0.75 1])); % For overexposed
picro_red_bw = im2bw(imadjust(picro_red));
%picro_green_bw = im2bw(imadjust(picro_green,[0.15 1])); % For overexposed
picro_green_bw = im2bw(imadjust(picro_green));
% Define diamond structuring element
seD = strel('diamond',2);
picro_red_bw = imerode(picro_red_bw,seD);
% picro_red_bw = imerode(picro_red_bw,seD); % If extra erosion needed
% picro_red_bw = imerode(picro_red_bw,seD); % If extra erosion needed
% seD = strel('diamond',1); % Erode green images if noisy
% picro_green_bw = imerode(picro_green_bw,seD);
% Close bw objects
picro_red_bw = imclose(picro_red_bw,seD);
picro_red_bw = imclose(picro_red_bw,seD);
picro_green_bw = imclose(picro_green_bw,seD);
% Skeletonize, separate at branchpoints, and delete single pixel lines
picro_red_skel = bwmorph(picro_red_bw,'skel',Inf);
picro_red_branch = bwmorph(picro_red_skel,'branchpoints');
picro_red_skel(picro_red_branch == 1) = 0;
picro_red_skel = bwmorph(picro_red_skel,'clean');
stored_red_skel{i} = picro_red_skel;
[L_red, num_red] = bwlabel(picro_red_skel);
[~, num_red_branch] = bwlabel(picro_red_branch);
props_red = regionprops(picro_red_bw,'Area','EulerNumber','Extent','ConvexArea',...
'FilledArea','Solidity','Perimeter','Eccentricity','MajorAxisLength',...
'EquivDiameter','MinorAxisLength');
props_red_skel = regionprops(picro_red_skel,'Orientation');
kappa_red = circ_kappa([props_red_skel.Orientation]*pi/180);
stored_red(i,:) = [mean(mean(picro_red)) num_red/scale_2 mean([props_red.MajorAxisLength])*scale_1 ...
mean([props_red.MinorAxisLength])*scale_1 mean([props_red.MajorAxisLength])/mean([props_red.Perimeter]) ...
kappa_red num_red_branch/scale_2 mean([props_red.EulerNumber]) mean([props_red.Extent]) ...
mean([props_red.Perimeter])*scale_1 mean([props_red.Solidity]) mean([props_red.Eccentricity]) ...
mean([props_red.EquivDiameter])*scale_1];
picro_green_skel = bwmorph(picro_green_bw,'skel',Inf);
picro_green_branch = bwmorph(picro_green_skel,'branchpoints');
picro_green_skel(picro_green_branch == 1) = 0;
picro_green_skel = bwmorph(picro_green_skel,'clean');
stored_green_skel{i} = picro_green_skel;
[L_green, num_green] = bwlabel(picro_green_skel);
[~, num_green_branch] = bwlabel(picro_green_branch);
props_green = regionprops(picro_green_bw,'Area','EulerNumber','Extent','ConvexArea',...
'FilledArea','Solidity','Perimeter','Eccentricity','MajorAxisLength',...
'EquivDiameter','MinorAxisLength');
props_green_skel = regionprops(picro_green_skel,'Orientation');
kappa_green = circ_kappa([props_green_skel.Orientation]*pi/180);
stored_green(i,:) = [mean(mean(picro_green)) num_green/scale_2 mean([props_green.MajorAxisLength])*scale_1 ...
mean([props_green.MinorAxisLength])*scale_1 mean([props_green.MajorAxisLength])/mean([props_green.Perimeter]) ...
kappa_green num_green_branch/scale_2 mean([props_green.EulerNumber]) mean([props_green.Extent]) ...
mean([props_green.Perimeter])*scale_1 mean([props_green.Solidity]) mean([props_green.Eccentricity]) ...
mean([props_green.EquivDiameter])*scale_1];
% Brightness, Number of Fibers/1000 microns^2, Length (micron), Width
% (micron), Persistence, Angle Randomness, Number of Branchpoints/1000
% micron^2, Area, Euler Number, Extent, Perimeter, Convex Area, Filled
% Area, Solidity, Eccentricity, Equivalent Diameter
counter = counter + 1;
progress = 2*100*counter/(nfiles)
end
stored_red(all(stored_red == 0,2),:) = []; % Fiber characteristics for red images
stored_red_names = stored_red_names(~cellfun(@isempty,stored_red_names))';
stored_red_skel = stored_red_skel(~cellfun(@isempty,stored_red_skel))';
stored_green(all(stored_green == 0,2),:) = []; % Fiber characteristics for green images
stored_green_names = stored_green_names(~cellfun(@isempty,stored_green_names))';
stored_green_skel = stored_green_skel(~cellfun(@isempty,stored_green_skel))';
stored_redgreen = [stored_red stored_green];
%% Define tSNE class labels
close all
% 2 week images
class_pt = {'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert'};
% 1 month images
class_pt = {'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert'};
% 3 month images
class_pt = {'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert','Vert',...
'Vert','Vert','Vert','Vert','Vert','Vert',};
% Class labels for MD Vert experiment
class_pt = {'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded','Unwounded','Unwounded',...
'Unwounded','Unwounded','Unwounded',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks','2 Doses - 2 weeks',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month','2 Doses - 1 month',...
'2 Doses - 2 weeks',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS','PBS',...
'PBS','PBS',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 2 weeks','4 Doses - 2 weeks','4 Doses - 2 weeks',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month',...
'4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month','4 Doses - 1 month'};
Y = tsne(stored_redgreen);
figure
gscatter(Y(:,1),Y(:,2),class_pt')
3.1_原文fork代码/fiberquantificationmasterv2.m
%% Fiber Quantification Algorithm v2
% Created by Shamik Mascharak
% November 28th, 2019 - Created
% December 1st, 2019 - First iteration completed
% October 28th, 2020 - Uploaded to Github
%% Read files and generate binarized images
clear
clc
close all
% Get list of all .tif files in the directory
% imagefiles = dir('*.tif');
imagefiles = dir('*.png');
nfiles = length(imagefiles);
%% Run color deconvolution
counter = 0;
% Optional scaling factors (to convert pixels to microns)
scale = 1;
% Define target image for image normalization
target_im = imread(imagefiles(1544).name); % Unwounded Extra 8
verbose = 0; % 1 to display results, 0 to not
% Define stain vector for color deconvolution
picro_vect = [0.001 0.70710677 0.70710677;...
0.70710677 0.001 0.70710677;...
0.57735026 0.57735026 0.57735026];
for i = 1:nfiles
filename = imagefiles(i).name;
stored_names{i} = filename;
stored_im{i} = imread(filename); % Store the original image
% Normalize image using RGB histogram method
stored_norm_im{i} = Norm(stored_im{i},target_im,'RGBHist',verbose);
% stored_norm_im{i} = stored_im{i}; % To bypass normalization
% Perform color deconvolution (Ruifrok method)
[Dch M] = Deconvolve(stored_norm_im{i},picro_vect,0);
[red green black] = PseudoColourStains(Dch,M);
picro_red = (255-rgb2gray(red))-(255-rgb2gray(black));
picro_green = (255-rgb2gray(green))-(255-rgb2gray(black));
% Remove noise with adaptive filtering
picro_red_noise = wiener2(picro_red,[3 3]);
picro_green_noise = wiener2(picro_green,[3 3]);
stored_picro_red_noise{i} = picro_red_noise;
stored_picro_green_noise{i} = picro_green_noise;
picro_red_bw = im2bw(imadjust(picro_red_noise));
picro_green_bw = im2bw(imadjust(picro_green_noise));
% Define diamond structuring element
seD = strel('diamond',1);
picro_red_bw_strel = imerode(picro_red_bw,seD);
picro_green_bw_strel = imerode(picro_green_bw,seD);
stored_picro_red_strel{i} = picro_red_bw_strel; % Binary map
stored_picro_green_strel{i} = picro_green_bw_strel; % Binary map
% Skeletonize, separate at branchpoints, and delete single pixel lines
picro_red_bw_skel = bwmorph(picro_red_bw_strel,'skel',Inf);
picro_red_bw_branch = bwmorph(picro_red_bw_skel,'branchpoints');
picro_red_bw_skel(picro_red_bw_branch == 1) = 0;
picro_red_bw_skel = bwmorph(picro_red_bw_skel,'clean');
stored_picro_red_skel{i} = picro_red_bw_skel; % Skeletonized image
stored_picro_red_branch{i} = picro_red_bw_branch;
picro_green_bw_skel = bwmorph(picro_green_bw_strel,'skel',Inf);
picro_green_bw_branch = bwmorph(picro_green_bw_skel,'branchpoints');
picro_green_bw_skel(picro_green_bw_branch == 1) = 0;
picro_green_bw_skel = bwmorph(picro_green_bw_skel,'clean');
stored_picro_green_skel{i} = picro_green_bw_skel; % Skeletonized image
stored_picro_green_branch{i} = picro_green_bw_branch;
counter = counter + 1;
progress = 100*counter/(nfiles)
end
stored_names = stored_names';
%% Quantify image features
clc
counter = 0;
for i = 1:nfiles
red = stored_picro_red_noise{i};
green = stored_picro_green_noise{i};
red_bw = stored_picro_red_strel{i};
green_bw = stored_picro_red_strel{i};
red_skel = stored_picro_red_skel{i};
green_skel = stored_picro_green_skel{i};
red_branch = stored_picro_red_branch{i};
green_branch = stored_picro_green_branch{i};
% Calculate values from grayscale images
props_red = regionprops(red_bw,red,'all');
props_green = regionprops(green_bw,green,'all');
% Fiber angle-related values
alpha_red = [props_red.Orientation]*pi/180;
circ_red = [mean(alpha_red) median(alpha_red) std(alpha_red) skewness(alpha_red)...
kurtosis(alpha_red) circ_kappa(alpha_red)];
alpha_green = [props_green.Orientation]*pi/180;
circ_green = [mean(alpha_green) median(alpha_green) std(alpha_green) skewness(alpha_green)...
kurtosis(alpha_green) circ_kappa(alpha_green)];
% Calculate Haralick features from grayscale images (4 different offsets)
offset = [0 1; -1 1;-1 0;-1 -1];
glcm_red = graycomatrix(red,'Offset',offset,'Symmetric',true);
glcm_green = graycomatrix(green,'Offset',offset,'Symmetric',true);
grayprops_red = graycoprops(glcm_red);
grayprops_green = graycoprops(glcm_green);
% Calculate values from binary images
props_red_bw = regionprops(red_bw,'all');
props_green_bw = regionprops(green_bw,'all');
% Fiber angle-related values
alpha_red_bw = [props_red_bw.Orientation]*pi/180;
circ_red_bw = [mean(alpha_red_bw) median(alpha_red_bw) std(alpha_red_bw) skewness(alpha_red_bw)...
kurtosis(alpha_red_bw) circ_kappa(alpha_red_bw)];
alpha_green_bw = [props_green_bw.Orientation]*pi/180;
circ_green_bw = [mean(alpha_green_bw) median(alpha_green_bw) std(alpha_green_bw) skewness(alpha_green_bw)...
kurtosis(alpha_green_bw) circ_kappa(alpha_green_bw)];
% Calculate values from skeletonized images
props_red_skel = regionprops(red_skel,'all');
props_green_skel = regionprops(green_skel,'all');
% Fiber angle-related values
alpha_red_skel = [props_red_skel.Orientation]*pi/180;
circ_red_skel = [mean(alpha_red_skel) median(alpha_red_skel) std(alpha_red_skel) skewness(alpha_red_skel)...
kurtosis(alpha_red_skel) circ_kappa(alpha_red_skel)];
alpha_green_skel = [props_green_skel.Orientation]*pi/180;
circ_green_skel = [mean(alpha_green_skel) median(alpha_green_skel) std(alpha_green_skel) skewness(alpha_green_skel)...
kurtosis(alpha_green_skel) circ_kappa(alpha_green_skel)];
% Branchpoint-related values
[L_red, num_red] = bwlabel(red_skel);
[~, num_red_branch] = bwlabel(red_branch);
[L_green, num_green] = bwlabel(green_skel);
[~, num_green_branch] = bwlabel(green_branch);
% Grayscale regionprops, grayscale circ, Haralick features, binary regionprops,
% binary circ, skeletonized regionprops, skeletonized circ, branchpoints,
quantified{i} = [mean([props_red.Area]) std([props_red.Area]) mean([props_green.Area]) std([props_green.Area])...
mean([props_red.MajorAxisLength]) std([props_red.MajorAxisLength]) mean([props_green.MajorAxisLength]) std([props_green.MajorAxisLength])...
mean([props_red.MinorAxisLength]) std([props_red.MinorAxisLength]) mean([props_green.MinorAxisLength]) std([props_green.MinorAxisLength])...
mean([props_red.Eccentricity]) std([props_red.Eccentricity]) mean([props_green.Eccentricity]) std([props_green.Eccentricity])...
mean([props_red.ConvexArea]) std([props_red.ConvexArea]) mean([props_green.ConvexArea]) std([props_green.ConvexArea])...
mean([props_red.Circularity]~=Inf) std([props_red.Circularity]~=Inf) mean([props_green.Circularity]~=Inf) std([props_green.Circularity]~=Inf)...
mean([props_red.FilledArea]) std([props_red.FilledArea]) mean([props_green.FilledArea]) std([props_green.FilledArea])...
mean([props_red.EulerNumber]) std([props_red.EulerNumber]) mean([props_green.EulerNumber]) std([props_green.EulerNumber])...
sum([props_red.EulerNumber]) sum([props_green.EulerNumber])...
mean([props_red.EquivDiameter]) std([props_red.EquivDiameter]) mean([props_green.EquivDiameter]) std([props_green.EquivDiameter])...
mean([props_red.Solidity]) std([props_red.Solidity]) mean([props_green.Solidity]) std([props_green.Solidity])...
mean([props_red.Extent]) std([props_red.Extent]) mean([props_green.Extent]) std([props_green.Extent])...
mean([props_red.Perimeter]) std([props_red.Perimeter]) mean([props_green.Perimeter]) std([props_green.Perimeter])...
mean([props_red.PerimeterOld]) std([props_red.PerimeterOld]) mean([props_green.PerimeterOld]) std([props_green.PerimeterOld])...
mean([props_red.MeanIntensity]) std([props_red.MeanIntensity]) mean([props_green.MeanIntensity]) std([props_green.MeanIntensity])...
mean([props_red.MinIntensity]) std(double([props_red.MinIntensity])) mean([props_green.MinIntensity]) std(double([props_green.MinIntensity]))...
mean([props_red.MaxIntensity]) std(double([props_red.MaxIntensity])) mean([props_green.MaxIntensity]) std(double([props_green.MaxIntensity]))...
mean([props_red.MaxFeretDiameter]) std([props_red.MaxFeretDiameter]) mean([props_green.MaxFeretDiameter]) std([props_green.MaxFeretDiameter])...
mean([props_red.MaxFeretAngle]) std([props_red.MaxFeretAngle]) mean([props_green.MaxFeretAngle]) std([props_green.MaxFeretAngle])...
mean([props_red.MinFeretDiameter]) std([props_red.MinFeretDiameter]) mean([props_green.MinFeretDiameter]) std([props_green.MinFeretDiameter])...
mean([props_red.MinFeretAngle]) std([props_red.MinFeretAngle]) mean([props_green.MinFeretAngle]) std([props_green.MinFeretAngle])...
circ_red circ_green...
[grayprops_red.Contrast] [grayprops_red.Correlation] [grayprops_red.Energy] [grayprops_red.Homogeneity]...
[grayprops_green.Contrast] [grayprops_green.Correlation] [grayprops_green.Energy] [grayprops_green.Homogeneity]...
mean([props_red_bw.Area]) std([props_red_bw.Area]) mean([props_green_bw.Area]) std([props_green_bw.Area])...
mean([props_red_bw.MajorAxisLength]) std([props_red_bw.MajorAxisLength]) mean([props_green_bw.MajorAxisLength]) std([props_green_bw.MajorAxisLength])...
mean([props_red_bw.MinorAxisLength]) std([props_red_bw.MinorAxisLength]) mean([props_green_bw.MinorAxisLength]) std([props_green_bw.MinorAxisLength])...
mean([props_red_bw.Eccentricity]) std([props_red_bw.Eccentricity]) mean([props_green_bw.Eccentricity]) std([props_green_bw.Eccentricity])...
mean([props_red_bw.ConvexArea]) std([props_red_bw.ConvexArea]) mean([props_green_bw.ConvexArea]) std([props_green_bw.ConvexArea])...
mean([props_red_bw.Circularity]~=Inf) std([props_red_bw.Circularity]~=Inf) mean([props_green_bw.Circularity]~=Inf) std([props_green_bw.Circularity]~=Inf)...
mean([props_red_bw.FilledArea]) std([props_red_bw.FilledArea]) mean([props_green_bw.FilledArea]) std([props_green_bw.FilledArea])...
mean([props_red_bw.EulerNumber]) std([props_red_bw.EulerNumber]) mean([props_green_bw.EulerNumber]) std([props_green_bw.EulerNumber])...
sum([props_red_bw.EulerNumber]) sum([props_green_bw.EulerNumber])...
mean([props_red_bw.EquivDiameter]) std([props_red_bw.EquivDiameter]) mean([props_green_bw.EquivDiameter]) std([props_green_bw.EquivDiameter])...
mean([props_red_bw.Solidity]) std([props_red_bw.Solidity]) mean([props_green_bw.Solidity]) std([props_green_bw.Solidity])...
mean([props_red_bw.Extent]) std([props_red_bw.Extent]) mean([props_green_bw.Extent]) std([props_green_bw.Extent])...
mean([props_red_bw.Perimeter]) std([props_red_bw.Perimeter]) mean([props_green_bw.Perimeter]) std([props_green_bw.Perimeter])...
mean([props_red_bw.PerimeterOld]) std([props_red_bw.PerimeterOld]) mean([props_green_bw.PerimeterOld]) std([props_green_bw.PerimeterOld])...
mean([props_red_bw.MaxFeretDiameter]) std([props_red_bw.MaxFeretDiameter]) mean([props_green_bw.MaxFeretDiameter]) std([props_green_bw.MaxFeretDiameter])...
mean([props_red_bw.MaxFeretAngle]) std([props_red_bw.MaxFeretAngle]) mean([props_green_bw.MaxFeretAngle]) std([props_green_bw.MaxFeretAngle])...
mean([props_red_bw.MinFeretDiameter]) std([props_red_bw.MinFeretDiameter]) mean([props_green_bw.MinFeretDiameter]) std([props_green_bw.MinFeretDiameter])...
mean([props_red_bw.MinFeretAngle]) std([props_red_bw.MinFeretAngle]) mean([props_green_bw.MinFeretAngle]) std([props_green_bw.MinFeretAngle])...
circ_red_bw circ_green_bw...
mean([props_red_skel.Area]) std([props_red_skel.Area]) mean([props_green_skel.Area]) std([props_green_skel.Area])...
mean([props_red_skel.MajorAxisLength]) std([props_red_skel.MajorAxisLength]) mean([props_green_skel.MajorAxisLength]) std([props_green_skel.MajorAxisLength])...
mean([props_red_skel.MinorAxisLength]) std([props_red_skel.MinorAxisLength]) mean([props_green_skel.MinorAxisLength]) std([props_green_skel.MinorAxisLength])...
mean([props_red_skel.Eccentricity]) std([props_red_skel.Eccentricity]) mean([props_green_skel.Eccentricity]) std([props_green_skel.Eccentricity])...
mean([props_red_skel.ConvexArea]) std([props_red_skel.ConvexArea]) mean([props_green_skel.ConvexArea]) std([props_green_skel.ConvexArea])...
mean([props_red_skel.Circularity]~=Inf) std([props_red_skel.Circularity]~=Inf) mean([props_green_skel.Circularity]~=Inf) std([props_green_skel.Circularity]~=Inf)...
mean([props_red_skel.FilledArea]) std([props_red_skel.FilledArea]) mean([props_green_skel.FilledArea]) std([props_green_skel.FilledArea])...
mean([props_red_skel.EulerNumber]) std([props_red_skel.EulerNumber]) mean([props_green_skel.EulerNumber]) std([props_green_skel.EulerNumber])...
sum([props_red_skel.EulerNumber]) sum([props_green_skel.EulerNumber])...
mean([props_red_skel.EquivDiameter]) std([props_red_skel.EquivDiameter]) mean([props_green_skel.EquivDiameter]) std([props_green_skel.EquivDiameter])...
mean([props_red_skel.Solidity]) std([props_red_skel.Solidity]) mean([props_green_skel.Solidity]) std([props_green_skel.Solidity])...
mean([props_red_skel.Extent]) std([props_red_skel.Extent]) mean([props_green_skel.Extent]) std([props_green_skel.Extent])...
mean([props_red_skel.Perimeter]) std([props_red_skel.Perimeter]) mean([props_green_skel.Perimeter]) std([props_green_skel.Perimeter])...
mean([props_red_skel.PerimeterOld]) std([props_red_skel.PerimeterOld]) mean([props_green_skel.PerimeterOld]) std([props_green_skel.PerimeterOld])...
mean([props_red_skel.MaxFeretDiameter]) std([props_red_skel.MaxFeretDiameter]) mean([props_green_skel.MaxFeretDiameter]) std([props_green_skel.MaxFeretDiameter])...
mean([props_red_skel.MaxFeretAngle]) std([props_red_skel.MaxFeretAngle]) mean([props_green_skel.MaxFeretAngle]) std([props_green_skel.MaxFeretAngle])...
mean([props_red_skel.MinFeretDiameter]) std([props_red_skel.MinFeretDiameter]) mean([props_green_skel.MinFeretDiameter]) std([props_green_skel.MinFeretDiameter])...
mean([props_red_skel.MinFeretAngle]) std([props_red_skel.MinFeretAngle]) mean([props_green_skel.MinFeretAngle]) std([props_green_skel.MinFeretAngle])...
circ_red_skel circ_green_skel...
num_red num_red_branch num_green num_green_branch];
counter = counter + 1;
progress = 100*counter/(nfiles)
end
quantified = cell2mat(quantified');
quantified_props = quantified(1:end,1:126);
quantified_props_bw = quantified(1:end,127:224);
quantified_props_skel = quantified(1:end,225:end);
%% Class labels - Mascharak, desJardins-Park, Januszyk et al
for i = 1:length(stored_names)
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD2')
class_pt{i} = 'PBS POD2';
color{i} = [1 0 0];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD7')
class_pt{i} = 'PBS POD7';
color{i} = [0 1 0];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD14')
class_pt{i} = 'PBS POD14';
color{i} = [0 0 1];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD30')
class_pt{i} = 'PBS POD30';
color{i} = [1 1 0];
end
if contains(stored_names{i},'- U') | contains(stored_names{i},'Unwounded')
class_pt{i} = 'Unwounded';
color{i} = [0 1 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD2')
class_pt{i} = 'Vert POD2';
color{i} = [1 0 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD7')
class_pt{i} = 'Vert POD7';
color{i} = [1 1 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD14')
class_pt{i} = 'Vert POD14';
color{i} = [0.5 0.5 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD30')
class_pt{i} = 'Vert POD30';
color{i} = [0.5 0.5 0.5];
end
end
%% t-SNE Visualization
Y = tsne(quantified,'Algorithm','exact','Distance','euclidian','Exaggeration',4);
figure
gscatter(Y(:,1),Y(:,2),class_pt')
3.1_原文fork代码/fiberquantificationmasterv3.m
%% Fiber Quantification Algorithm v3
% Created by Shamik Mascharak
% November 28th, 2019 - Created
% December 1st, 2019 - First iteration completed
% October 28th, 2020 - Uploaded to Github
% September 22nd, 2022 - Code updated
%% Read files and generate binarized images
clear
clc
close all
% Get list of all .tif files in the directory
imagefiles = dir('*.tif');
% imagefiles = dir('*.png');
nfiles = length(imagefiles);
%% Run color deconvolution
counter = 0;
% Optional scaling factors (to convert pixels to microns)
scale = 1;
% Define target image for image normalization
target_im = imread(imagefiles(1544).name); % Input index from imagefiles variable
verbose = 0; % 1 to display results, 0 to not
% Define stain vector for color deconvolution
picro_vect = [0.001 0.70710677 0.70710677;...
0.70710677 0.001 0.70710677;...
0.57735026 0.57735026 0.57735026];
% Define warning counters for images that are blank/empty after processing
warn_count_picro_green = 0;
warn_count_picro_green_bw = 0;
warn_count_picro_green_bw_branch = 0;
warn_count_picro_green_bw_skel = 0;
warn_count_picro_green_bw_strel = 0;
warn_count_picro_green_noise = 0;
warn_count_picro_red = 0;
warn_count_picro_red_bw = 0;
warn_count_picro_red_bw_branch = 0;
warn_count_picro_red_bw_skel = 0;
warn_count_picro_red_bw_strel = 0;
warn_count_picro_red_noise = 0;
for i = 1:nfiles
filename = imagefiles(i).name;
stored_names{i} = filename;
stored_im{i} = imread(filename); % Store the original image
% Normalize image using RGB histogram method
stored_norm_im{i} = Norm(stored_im{i},target_im,'RGBHist',verbose);
% stored_norm_im{i} = stored_im{i}; % To bypass normalization
% Perform color deconvolution (Ruifrok method)
[Dch M] = Deconvolve(stored_norm_im{i},picro_vect,0);
[red green black] = PseudoColourStains(Dch,M);
picro_red = (255-rgb2gray(red))-(255-rgb2gray(black));
picro_green = (255-rgb2gray(green))-(255-rgb2gray(black));
% Remove noise with adaptive filtering
picro_red_noise = wiener2(picro_red,[3 3]);
picro_green_noise = wiener2(picro_green,[3 3]);
stored_picro_red_noise{i} = picro_red_noise;
stored_picro_green_noise{i} = picro_green_noise;
picro_red_bw = im2bw(imadjust(picro_red_noise));
picro_green_bw = im2bw(imadjust(picro_green_noise));
% Define diamond structuring element
seD = strel('diamond',1);
picro_red_bw_strel = imerode(picro_red_bw,seD);
picro_green_bw_strel = imerode(picro_green_bw,seD);
stored_picro_red_strel{i} = picro_red_bw_strel; % Binary map
stored_picro_green_strel{i} = picro_green_bw_strel; % Binary map
% Skeletonize, separate at branchpoints, and delete single pixel lines
picro_red_bw_skel = bwmorph(picro_red_bw_strel,'skel',Inf);
picro_red_bw_branch = bwmorph(picro_red_bw_skel,'branchpoints');
picro_red_bw_skel(picro_red_bw_branch == 1) = 0;
picro_red_bw_skel = bwmorph(picro_red_bw_skel,'clean');
stored_picro_red_skel{i} = picro_red_bw_skel; % Skeletonized image
stored_picro_red_branch{i} = picro_red_bw_branch;
picro_green_bw_skel = bwmorph(picro_green_bw_strel,'skel',Inf);
picro_green_bw_branch = bwmorph(picro_green_bw_skel,'branchpoints');
picro_green_bw_skel(picro_green_bw_branch == 1) = 0;
picro_green_bw_skel = bwmorph(picro_green_bw_skel,'clean');
stored_picro_green_skel{i} = picro_green_bw_skel; % Skeletonized image
stored_picro_green_branch{i} = picro_green_bw_branch;
% Store indices of images that are blank/empty after processing. These
% indices can be found in the variables containing the 'warnings_' prefix.
if not(any(any(picro_green)))
warnings_picro_green(warn_count_picro_green+1) = i;
warn_count_picro_green = warn_count_picro_green + 1;
end
if not(any(any(picro_green_bw)))
warnings_picro_green_bw(warn_count_picro_green_bw+1) = i;
warn_count_picro_green_bw = warn_count_picro_green_bw + 1;
end
if not(any(any(picro_green_bw_branch)))
warnings_picro_green_bw_branch(warn_count_picro_green_bw_branch+1) = i;
warn_count_picro_green_bw_branch = warn_count_picro_green_bw_branch + 1;
end
if not(any(any(picro_green_bw_skel)))
warnings_picro_green_bw_skel(warn_count_picro_green_bw_skel+1) = i;
warn_count_picro_green_bw_skel = warn_count_picro_green_bw_skel + 1;
end
if not(any(any(picro_green_bw_strel)))
warnings_picro_green_bw_strel(warn_count_picro_green_bw_strel+1) = i;
warn_count_picro_green_bw_strel = warn_count_picro_green_bw_strel + 1;
end
if not(any(any(picro_green_noise)))
warnings_picro_green_noise(warn_count_picro_green_noise+1) = i;
warn_count_picro_green_noise = warn_count_picro_green_noise + 1;
end
if not(any(any(picro_red)))
warnings_picro_red(warn_count_picro_red+1) = i;
warn_count_picro_red = warn_count_picro_red + 1;
end
if not(any(any(picro_red_bw)))
warnings_picro_red_bw(warn_count_picro_red_bw+1) = i;
warn_count_picro_red_bw = warn_count_picro_red_bw + 1;
end
if not(any(any(picro_red_bw_branch)))
warnings_picro_red_bw_branch(warn_count_picro_red_bw_branch+1) = i;
warn_count_picro_red_bw_branch = warn_count_picro_red_bw_branch + 1;
end
if not(any(any(picro_red_bw_skel)))
warnings_picro_red_bw_skel(warn_count_picro_red_bw_skel+1) = i;
warn_count_picro_red_bw_skel = warn_count_picro_red_bw_skel + 1;
end
if not(any(any(picro_red_bw_strel)))
warnings_picro_red_bw_strel(warn_count_picro_red_bw_strel+1) = i;
warn_count_picro_red_bw_strel = warn_count_picro_red_bw_strel + 1;
end
if not(any(any(picro_red_noise)))
warnings_picro_red_noise(warn_count_picro_red_noise+1) = i;
warn_count_picro_red_noise = warn_count_picro_red_noise + 1;
end
% Display progress
counter = counter + 1;
progress = 100*counter/(nfiles)
end
stored_names = stored_names';
%% Quantify image features
clc
counter = 0;
for i = 1:nfiles
red = stored_picro_red_noise{i};
green = stored_picro_green_noise{i};
red_bw = stored_picro_red_strel{i};
green_bw = stored_picro_green_strel{i};
red_skel = stored_picro_red_skel{i};
green_skel = stored_picro_green_skel{i};
red_branch = stored_picro_red_branch{i};
green_branch = stored_picro_green_branch{i};
% Calculate values from grayscale images
props_red = regionprops(red_bw,red,'all');
props_green = regionprops(green_bw,green,'all');
% Fiber angle-related values
alpha_red = [props_red.Orientation]*pi/180;
circ_red = [mean(alpha_red) median(alpha_red) std(alpha_red) skewness(alpha_red)...
kurtosis(alpha_red) circ_kappa(alpha_red)];
alpha_green = [props_green.Orientation]*pi/180;
circ_green = [mean(alpha_green) median(alpha_green) std(alpha_green) skewness(alpha_green)...
kurtosis(alpha_green) circ_kappa(alpha_green)];
% Calculate Haralick features from grayscale images (4 different offsets)
offset = [0 1; -1 1;-1 0;-1 -1];
glcm_red = graycomatrix(red,'Offset',offset,'Symmetric',true);
glcm_green = graycomatrix(green,'Offset',offset,'Symmetric',true);
grayprops_red = graycoprops(glcm_red);
grayprops_green = graycoprops(glcm_green);
% Calculate values from binary images
props_red_bw = regionprops(red_bw,'all');
props_green_bw = regionprops(green_bw,'all');
% Fiber angle-related values
alpha_red_bw = [props_red_bw.Orientation]*pi/180;
circ_red_bw = [mean(alpha_red_bw) median(alpha_red_bw) std(alpha_red_bw) skewness(alpha_red_bw)...
kurtosis(alpha_red_bw) circ_kappa(alpha_red_bw)];
alpha_green_bw = [props_green_bw.Orientation]*pi/180;
circ_green_bw = [mean(alpha_green_bw) median(alpha_green_bw) std(alpha_green_bw) skewness(alpha_green_bw)...
kurtosis(alpha_green_bw) circ_kappa(alpha_green_bw)];
% Calculate values from skeletonized images
props_red_skel = regionprops(red_skel,'all');
props_green_skel = regionprops(green_skel,'all');
% Fiber angle-related values
alpha_red_skel = [props_red_skel.Orientation]*pi/180;
circ_red_skel = [mean(alpha_red_skel) median(alpha_red_skel) std(alpha_red_skel) skewness(alpha_red_skel)...
kurtosis(alpha_red_skel) circ_kappa(alpha_red_skel)];
alpha_green_skel = [props_green_skel.Orientation]*pi/180;
circ_green_skel = [mean(alpha_green_skel) median(alpha_green_skel) std(alpha_green_skel) skewness(alpha_green_skel)...
kurtosis(alpha_green_skel) circ_kappa(alpha_green_skel)];
% Branchpoint-related values
[L_red, num_red] = bwlabel(red_skel);
[~, num_red_branch] = bwlabel(red_branch);
[L_green, num_green] = bwlabel(green_skel);
[~, num_green_branch] = bwlabel(green_branch);
% Grayscale regionprops, grayscale circ, Haralick features, binary regionprops,
% binary circ, skeletonized regionprops, skeletonized circ, branchpoints,
quantified{i} = [mean([props_red.Area]) std([props_red.Area]) mean([props_green.Area]) std([props_green.Area])...
mean([props_red.MajorAxisLength]) std([props_red.MajorAxisLength]) mean([props_green.MajorAxisLength]) std([props_green.MajorAxisLength])...
mean([props_red.MinorAxisLength]) std([props_red.MinorAxisLength]) mean([props_green.MinorAxisLength]) std([props_green.MinorAxisLength])...
mean([props_red.Eccentricity]) std([props_red.Eccentricity]) mean([props_green.Eccentricity]) std([props_green.Eccentricity])...
mean([props_red.ConvexArea]) std([props_red.ConvexArea]) mean([props_green.ConvexArea]) std([props_green.ConvexArea])...
mean([props_red.Circularity]~=Inf) std([props_red.Circularity]~=Inf) mean([props_green.Circularity]~=Inf) std([props_green.Circularity]~=Inf)...
mean([props_red.FilledArea]) std([props_red.FilledArea]) mean([props_green.FilledArea]) std([props_green.FilledArea])...
mean([props_red.EulerNumber]) std([props_red.EulerNumber]) mean([props_green.EulerNumber]) std([props_green.EulerNumber])...
sum([props_red.EulerNumber]) sum([props_green.EulerNumber])...
mean([props_red.EquivDiameter]) std([props_red.EquivDiameter]) mean([props_green.EquivDiameter]) std([props_green.EquivDiameter])...
mean([props_red.Solidity]) std([props_red.Solidity]) mean([props_green.Solidity]) std([props_green.Solidity])...
mean([props_red.Extent]) std([props_red.Extent]) mean([props_green.Extent]) std([props_green.Extent])...
mean([props_red.Perimeter]) std([props_red.Perimeter]) mean([props_green.Perimeter]) std([props_green.Perimeter])...
mean([props_red.PerimeterOld]) std([props_red.PerimeterOld]) mean([props_green.PerimeterOld]) std([props_green.PerimeterOld])...
mean([props_red.MeanIntensity]) std([props_red.MeanIntensity]) mean([props_green.MeanIntensity]) std([props_green.MeanIntensity])...
mean([props_red.MinIntensity]) std(double([props_red.MinIntensity])) mean([props_green.MinIntensity]) std(double([props_green.MinIntensity]))...
mean([props_red.MaxIntensity]) std(double([props_red.MaxIntensity])) mean([props_green.MaxIntensity]) std(double([props_green.MaxIntensity]))...
mean([props_red.MaxFeretDiameter]) std([props_red.MaxFeretDiameter]) mean([props_green.MaxFeretDiameter]) std([props_green.MaxFeretDiameter])...
mean([props_red.MaxFeretAngle]) std([props_red.MaxFeretAngle]) mean([props_green.MaxFeretAngle]) std([props_green.MaxFeretAngle])...
mean([props_red.MinFeretDiameter]) std([props_red.MinFeretDiameter]) mean([props_green.MinFeretDiameter]) std([props_green.MinFeretDiameter])...
mean([props_red.MinFeretAngle]) std([props_red.MinFeretAngle]) mean([props_green.MinFeretAngle]) std([props_green.MinFeretAngle])...
circ_red circ_green...
[grayprops_red.Contrast] [grayprops_red.Correlation] [grayprops_red.Energy] [grayprops_red.Homogeneity]...
[grayprops_green.Contrast] [grayprops_green.Correlation] [grayprops_green.Energy] [grayprops_green.Homogeneity]...
mean([props_red_bw.Area]) std([props_red_bw.Area]) mean([props_green_bw.Area]) std([props_green_bw.Area])...
mean([props_red_bw.MajorAxisLength]) std([props_red_bw.MajorAxisLength]) mean([props_green_bw.MajorAxisLength]) std([props_green_bw.MajorAxisLength])...
mean([props_red_bw.MinorAxisLength]) std([props_red_bw.MinorAxisLength]) mean([props_green_bw.MinorAxisLength]) std([props_green_bw.MinorAxisLength])...
mean([props_red_bw.Eccentricity]) std([props_red_bw.Eccentricity]) mean([props_green_bw.Eccentricity]) std([props_green_bw.Eccentricity])...
mean([props_red_bw.ConvexArea]) std([props_red_bw.ConvexArea]) mean([props_green_bw.ConvexArea]) std([props_green_bw.ConvexArea])...
mean([props_red_bw.Circularity]~=Inf) std([props_red_bw.Circularity]~=Inf) mean([props_green_bw.Circularity]~=Inf) std([props_green_bw.Circularity]~=Inf)...
mean([props_red_bw.FilledArea]) std([props_red_bw.FilledArea]) mean([props_green_bw.FilledArea]) std([props_green_bw.FilledArea])...
mean([props_red_bw.EulerNumber]) std([props_red_bw.EulerNumber]) mean([props_green_bw.EulerNumber]) std([props_green_bw.EulerNumber])...
sum([props_red_bw.EulerNumber]) sum([props_green_bw.EulerNumber])...
mean([props_red_bw.EquivDiameter]) std([props_red_bw.EquivDiameter]) mean([props_green_bw.EquivDiameter]) std([props_green_bw.EquivDiameter])...
mean([props_red_bw.Solidity]) std([props_red_bw.Solidity]) mean([props_green_bw.Solidity]) std([props_green_bw.Solidity])...
mean([props_red_bw.Extent]) std([props_red_bw.Extent]) mean([props_green_bw.Extent]) std([props_green_bw.Extent])...
mean([props_red_bw.Perimeter]) std([props_red_bw.Perimeter]) mean([props_green_bw.Perimeter]) std([props_green_bw.Perimeter])...
mean([props_red_bw.PerimeterOld]) std([props_red_bw.PerimeterOld]) mean([props_green_bw.PerimeterOld]) std([props_green_bw.PerimeterOld])...
mean([props_red_bw.MaxFeretDiameter]) std([props_red_bw.MaxFeretDiameter]) mean([props_green_bw.MaxFeretDiameter]) std([props_green_bw.MaxFeretDiameter])...
mean([props_red_bw.MaxFeretAngle]) std([props_red_bw.MaxFeretAngle]) mean([props_green_bw.MaxFeretAngle]) std([props_green_bw.MaxFeretAngle])...
mean([props_red_bw.MinFeretDiameter]) std([props_red_bw.MinFeretDiameter]) mean([props_green_bw.MinFeretDiameter]) std([props_green_bw.MinFeretDiameter])...
mean([props_red_bw.MinFeretAngle]) std([props_red_bw.MinFeretAngle]) mean([props_green_bw.MinFeretAngle]) std([props_green_bw.MinFeretAngle])...
circ_red_bw circ_green_bw...
mean([props_red_skel.Area]) std([props_red_skel.Area]) mean([props_green_skel.Area]) std([props_green_skel.Area])...
mean([props_red_skel.MajorAxisLength]) std([props_red_skel.MajorAxisLength]) mean([props_green_skel.MajorAxisLength]) std([props_green_skel.MajorAxisLength])...
mean([props_red_skel.MinorAxisLength]) std([props_red_skel.MinorAxisLength]) mean([props_green_skel.MinorAxisLength]) std([props_green_skel.MinorAxisLength])...
mean([props_red_skel.Eccentricity]) std([props_red_skel.Eccentricity]) mean([props_green_skel.Eccentricity]) std([props_green_skel.Eccentricity])...
mean([props_red_skel.ConvexArea]) std([props_red_skel.ConvexArea]) mean([props_green_skel.ConvexArea]) std([props_green_skel.ConvexArea])...
mean([props_red_skel.Circularity]~=Inf) std([props_red_skel.Circularity]~=Inf) mean([props_green_skel.Circularity]~=Inf) std([props_green_skel.Circularity]~=Inf)...
mean([props_red_skel.FilledArea]) std([props_red_skel.FilledArea]) mean([props_green_skel.FilledArea]) std([props_green_skel.FilledArea])...
mean([props_red_skel.EulerNumber]) std([props_red_skel.EulerNumber]) mean([props_green_skel.EulerNumber]) std([props_green_skel.EulerNumber])...
sum([props_red_skel.EulerNumber]) sum([props_green_skel.EulerNumber])...
mean([props_red_skel.EquivDiameter]) std([props_red_skel.EquivDiameter]) mean([props_green_skel.EquivDiameter]) std([props_green_skel.EquivDiameter])...
mean([props_red_skel.Solidity]) std([props_red_skel.Solidity]) mean([props_green_skel.Solidity]) std([props_green_skel.Solidity])...
mean([props_red_skel.Extent]) std([props_red_skel.Extent]) mean([props_green_skel.Extent]) std([props_green_skel.Extent])...
mean([props_red_skel.Perimeter]) std([props_red_skel.Perimeter]) mean([props_green_skel.Perimeter]) std([props_green_skel.Perimeter])...
mean([props_red_skel.PerimeterOld]) std([props_red_skel.PerimeterOld]) mean([props_green_skel.PerimeterOld]) std([props_green_skel.PerimeterOld])...
mean([props_red_skel.MaxFeretDiameter]) std([props_red_skel.MaxFeretDiameter]) mean([props_green_skel.MaxFeretDiameter]) std([props_green_skel.MaxFeretDiameter])...
mean([props_red_skel.MaxFeretAngle]) std([props_red_skel.MaxFeretAngle]) mean([props_green_skel.MaxFeretAngle]) std([props_green_skel.MaxFeretAngle])...
mean([props_red_skel.MinFeretDiameter]) std([props_red_skel.MinFeretDiameter]) mean([props_green_skel.MinFeretDiameter]) std([props_green_skel.MinFeretDiameter])...
mean([props_red_skel.MinFeretAngle]) std([props_red_skel.MinFeretAngle]) mean([props_green_skel.MinFeretAngle]) std([props_green_skel.MinFeretAngle])...
circ_red_skel circ_green_skel...
num_red num_red_branch num_green num_green_branch];
% Display progress
counter = counter + 1;
progress = 100*counter/(nfiles)
end
quantified = cell2mat(quantified');
quantified_props = quantified(1:end,1:126);
quantified_props_bw = quantified(1:end,127:224);
quantified_props_skel = quantified(1:end,225:end);
%% Class labels - Mascharak, desJardins-Park, Januszyk et al
for i = 1:length(stored_names)
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD2')
class_pt{i} = 'PBS POD2';
color{i} = [1 0 0];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD7')
class_pt{i} = 'PBS POD7';
color{i} = [0 1 0];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD14')
class_pt{i} = 'PBS POD14';
color{i} = [0 0 1];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD30')
class_pt{i} = 'PBS POD30';
color{i} = [1 1 0];
end
if contains(stored_names{i},'- U') | contains(stored_names{i},'Unwounded')
class_pt{i} = 'Unwounded';
color{i} = [0 1 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD2')
class_pt{i} = 'Vert POD2';
color{i} = [1 0 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD7')
class_pt{i} = 'Vert POD7';
color{i} = [1 1 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD14')
class_pt{i} = 'Vert POD14';
color{i} = [0.5 0.5 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD30')
class_pt{i} = 'Vert POD30';
color{i} = [0.5 0.5 0.5];
end
end
%% t-SNE Visualization
Y = tsne(quantified,'Algorithm','exact','Distance','euclidian','Exaggeration',4);
figure
gscatter(Y(:,1),Y(:,2),class_pt')
3.1_原文fork代码/fiberquantificationmasterv3_trichrome.m
%% Fiber Quantification Algorithm v3
% Created by Shamik Mascharak
% November 28th, 2019 - Created
% December 1st, 2019 - First iteration completed
% October 28th, 2020 - Uploaded to Github
% September 22nd, 2022 - Code updated
%% Read files and generate binarized images
clear
clc
close all
% Get list of all .tif files in the directory
imagefiles = dir('*.tif');
% imagefiles = dir('*.png');
nfiles = length(imagefiles);
%% Run color deconvolution
counter = 0;
% Optional scaling factors (to convert pixels to microns)
scale = 1;
% Define target image for image normalization
target_im = imread(imagefiles(1544).name); % Input index from imagefiles variable
verbose = 0; % 1 to display results, 0 to not
% Define stain vector for color deconvolution
trichrome_vect = [0.7995107 0.5913521 0.10528667;...
0.099971585 0.73738605 0.6680326;...
0.59227383 0.3264422 0.7366459];
% Define warning counters for images that are blank/empty after processing
warn_count_blue = 0;
warn_count_blue_bw = 0;
warn_count_blue_bw_branch = 0;
warn_count_blue_bw_skel = 0;
warn_count_blue_bw_strel = 0;
warn_count_blue_noise = 0;
for i = 1:nfiles
filename = imagefiles(i).name;
stored_names{i} = filename;
stored_im{i} = imread(filename); % Store the original image
% Normalize image using RGB histogram method
stored_norm_im{i} = Norm(stored_im{i},target_im,'RGBHist',verbose);
% stored_norm_im{i} = stored_im{i}; % To bypass normalization
% Perform color deconvolution (Ruifrok method)
[Dch M] = Deconvolve(stored_norm_im{i},trichrome_vect,0);
[blue red green] = PseudoColourStains(Dch,M);
trichrome_blue = (255-rgb2gray(blue))-(255-rgb2gray(green));
% Remove noise with adaptive filtering
trichrome_blue_noise = wiener2(trichrome_blue,[3 3]);
stored_trichrome_blue_noise{i} = trichrome_blue_noise;
trichrome_blue_bw = im2bw(imadjust(trichrome_blue_noise));
% Define diamond structuring element
seD = strel('diamond',1);
trichrome_blue_bw_strel = imerode(trichrome_blue_bw,seD);
stored_trichrome_blue_strel{i} = trichrome_blue_bw_strel; % Binary map
% Skeletonize, separate at branchpoints, and delete single pixel lines
trichrome_blue_bw_skel = bwmorph(trichrome_blue_bw_strel,'skel',Inf);
trichrome_blue_bw_branch = bwmorph(trichrome_blue_bw_skel,'branchpoints');
trichrome_blue_bw_skel(trichrome_blue_bw_branch == 1) = 0;
trichrome_blue_bw_skel = bwmorph(trichrome_blue_bw_skel,'clean');
stored_trichrome_blue_skel{i} = trichrome_blue_bw_skel; % Skeletonized image
stored_trichrome_blue_branch{i} = trichrome_blue_bw_branch;
% Store indices of images that are blank/empty after processing. These
% indices can be found in the variables containing the 'warnings_' prefix.
if not(any(any(trichrome_blue)))
warnings_trichrome_blue(warn_count_blue+1) = i;
warn_count_blue = warn_count_blue + 1;
end
if not(any(any(trichrome_blue_bw)))
warnings_trichrome_blue_bw(warn_count_blue_bw+1) = i;
warn_count_blue_bw = warn_count_blue_bw + 1;
end
if not(any(any(trichrome_blue_bw_branch)))
warnings_trichrome_blue_bw_branch(warn_count_blue_bw_branch+1) = i;
warn_count_blue_bw_branch = warn_count_blue_bw_branch + 1;
end
if not(any(any(trichrome_blue_bw_skel)))
warnings_trichrome_blue_bw_skel(warn_count_blue_bw_skel+1) = i;
warn_count_blue_bw_skel = warn_count_blue_bw_skel + 1;
end
if not(any(any(trichrome_blue_bw_strel)))
warnings_trichrome_blue_bw_strel(warn_count_blue_bw_strel+1) = i;
warn_count_blue_bw_strel = warn_count_blue_bw_strel + 1;
end
if not(any(any(trichrome_blue_noise)))
warnings_trichrome_blue_noise(warn_count_blue_noise+1) = i;
warn_count_blue_noise = warn_count_blue_noise + 1;
end
% Display progress
counter = counter + 1;
progress = 100*counter/(nfiles)
end
stored_names = stored_names';
%% Quantify image features
clc
counter = 0;
for i = 1:nfiles
blue = stored_trichrome_blue_noise{i};
blue_bw = stored_trichrome_blue_strel{i};
blue_skel = stored_trichrome_blue_skel{i};
blue_branch = stored_trichrome_blue_branch{i};
% Calculate values from grayscale images
props_blue = regionprops(blue_bw,blue,'all');
% Fiber angle-related values
alpha_blue = [props_blue.Orientation]*pi/180;
circ_blue = [mean(alpha_blue) median(alpha_blue) std(alpha_blue) skewness(alpha_blue)...
kurtosis(alpha_blue) circ_kappa(alpha_blue)];
% Calculate Haralick features from grayscale images (4 different offsets)
offset = [0 1; -1 1;-1 0;-1 -1];
glcm_blue = graycomatrix(blue,'Offset',offset,'Symmetric',true);
grayprops_blue = graycoprops(glcm_blue);
% Calculate values from binary images
props_blue_bw = regionprops(blue_bw,'all');
% Fiber angle-related values
alpha_blue_bw = [props_blue_bw.Orientation]*pi/180;
circ_blue_bw = [mean(alpha_blue_bw) median(alpha_blue_bw) std(alpha_blue_bw) skewness(alpha_blue_bw)...
kurtosis(alpha_blue_bw) circ_kappa(alpha_blue_bw)];
% Calculate values from skeletonized images
props_blue_skel = regionprops(blue_skel,'all');
% Fiber angle-related values
alpha_blue_skel = [props_blue_skel.Orientation]*pi/180;
circ_blue_skel = [mean(alpha_blue_skel) median(alpha_blue_skel) std(alpha_blue_skel) skewness(alpha_blue_skel)...
kurtosis(alpha_blue_skel) circ_kappa(alpha_blue_skel)];
% Branchpoint-related values
[L_blue, num_blue] = bwlabel(blue_skel);
[~, num_blue_branch] = bwlabel(blue_branch);
% Grayscale regionprops, grayscale circ, Haralick features, binary regionprops,
% binary circ, skeletonized regionprops, skeletonized circ, branchpoints,
quantified{i} = [mean([props_blue.Area]) std([props_blue.Area])...
mean([props_blue.MajorAxisLength]) std([props_blue.MajorAxisLength])...
mean([props_blue.MinorAxisLength]) std([props_blue.MinorAxisLength])...
mean([props_blue.Eccentricity]) std([props_blue.Eccentricity])...
mean([props_blue.ConvexArea]) std([props_blue.ConvexArea])...
mean([props_blue.Circularity]~=Inf) std([props_blue.Circularity]~=Inf)...
mean([props_blue.FilledArea]) std([props_blue.FilledArea])...
mean([props_blue.EulerNumber]) std([props_blue.EulerNumber])...
sum([props_blue.EulerNumber])...
mean([props_blue.EquivDiameter]) std([props_blue.EquivDiameter])...
mean([props_blue.Solidity]) std([props_blue.Solidity])...
mean([props_blue.Extent]) std([props_blue.Extent])...
mean([props_blue.Perimeter]) std([props_blue.Perimeter])...
mean([props_blue.PerimeterOld]) std([props_blue.PerimeterOld])...
mean([props_blue.MeanIntensity]) std([props_blue.MeanIntensity])...
mean([props_blue.MinIntensity]) std(double([props_blue.MinIntensity]))...
mean([props_blue.MaxIntensity]) std(double([props_blue.MaxIntensity]))...
mean([props_blue.MaxFeretDiameter]) std([props_blue.MaxFeretDiameter])...
mean([props_blue.MaxFeretAngle]) std([props_blue.MaxFeretAngle])...
mean([props_blue.MinFeretDiameter]) std([props_blue.MinFeretDiameter])...
mean([props_blue.MinFeretAngle]) std([props_blue.MinFeretAngle])...
circ_blue...
[grayprops_blue.Contrast] [grayprops_blue.Correlation] [grayprops_blue.Energy] [grayprops_blue.Homogeneity]...
mean([props_blue_bw.Area]) std([props_blue_bw.Area])...
mean([props_blue_bw.MajorAxisLength]) std([props_blue_bw.MajorAxisLength])...
mean([props_blue_bw.MinorAxisLength]) std([props_blue_bw.MinorAxisLength])...
mean([props_blue_bw.Eccentricity]) std([props_blue_bw.Eccentricity])...
mean([props_blue_bw.ConvexArea]) std([props_blue_bw.ConvexArea])...
mean([props_blue_bw.Circularity]~=Inf) std([props_blue_bw.Circularity]~=Inf)...
mean([props_blue_bw.FilledArea]) std([props_blue_bw.FilledArea])...
mean([props_blue_bw.EulerNumber]) std([props_blue_bw.EulerNumber])...
sum([props_blue_bw.EulerNumber])...
mean([props_blue_bw.EquivDiameter]) std([props_blue_bw.EquivDiameter])...
mean([props_blue_bw.Solidity]) std([props_blue_bw.Solidity])...
mean([props_blue_bw.Extent]) std([props_blue_bw.Extent])...
mean([props_blue_bw.Perimeter]) std([props_blue_bw.Perimeter])...
mean([props_blue_bw.PerimeterOld]) std([props_blue_bw.PerimeterOld])...
mean([props_blue_bw.MaxFeretDiameter]) std([props_blue_bw.MaxFeretDiameter])...
mean([props_blue_bw.MaxFeretAngle]) std([props_blue_bw.MaxFeretAngle])...
mean([props_blue_bw.MinFeretDiameter]) std([props_blue_bw.MinFeretDiameter])...
mean([props_blue_bw.MinFeretAngle]) std([props_blue_bw.MinFeretAngle])...
circ_blue_bw...
mean([props_blue_skel.Area]) std([props_blue_skel.Area])...
mean([props_blue_skel.MajorAxisLength]) std([props_blue_skel.MajorAxisLength])...
mean([props_blue_skel.MinorAxisLength]) std([props_blue_skel.MinorAxisLength])...
mean([props_blue_skel.Eccentricity]) std([props_blue_skel.Eccentricity])...
mean([props_blue_skel.ConvexArea]) std([props_blue_skel.ConvexArea])...
mean([props_blue_skel.Circularity]~=Inf) std([props_blue_skel.Circularity]~=Inf)...
mean([props_blue_skel.FilledArea]) std([props_blue_skel.FilledArea])...
mean([props_blue_skel.EulerNumber]) std([props_blue_skel.EulerNumber])...
sum([props_blue_skel.EulerNumber])...
mean([props_blue_skel.EquivDiameter]) std([props_blue_skel.EquivDiameter])...
mean([props_blue_skel.Solidity]) std([props_blue_skel.Solidity])...
mean([props_blue_skel.Extent]) std([props_blue_skel.Extent])...
mean([props_blue_skel.Perimeter]) std([props_blue_skel.Perimeter])...
mean([props_blue_skel.PerimeterOld]) std([props_blue_skel.PerimeterOld])...
mean([props_blue_skel.MaxFeretDiameter]) std([props_blue_skel.MaxFeretDiameter])...
mean([props_blue_skel.MaxFeretAngle]) std([props_blue_skel.MaxFeretAngle])...
mean([props_blue_skel.MinFeretDiameter]) std([props_blue_skel.MinFeretDiameter])...
mean([props_blue_skel.MinFeretAngle]) std([props_blue_skel.MinFeretAngle])...
circ_blue_skel...
num_blue num_blue_branch];
counter = counter + 1;
progress = 100*counter/(nfiles)
end
quantified = cell2mat(quantified');
quantified_props = quantified(1:end,1:126);
quantified_props_bw = quantified(1:end,127:224);
quantified_props_skel = quantified(1:end,225:end);
%% Class labels - Mascharak, desJardins-Park, Januszyk et al
for i = 1:length(stored_names)
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD2')
class_pt{i} = 'PBS POD2';
color{i} = [1 0 0];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD7')
class_pt{i} = 'PBS POD7';
color{i} = [0 1 0];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD14')
class_pt{i} = 'PBS POD14';
color{i} = [0 0 1];
end
if contains(stored_names{i},'- P') & contains(stored_names{i},'POD30')
class_pt{i} = 'PBS POD30';
color{i} = [1 1 0];
end
if contains(stored_names{i},'- U') | contains(stored_names{i},'Unwounded')
class_pt{i} = 'Unwounded';
color{i} = [0 1 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD2')
class_pt{i} = 'Vert POD2';
color{i} = [1 0 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD7')
class_pt{i} = 'Vert POD7';
color{i} = [1 1 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD14')
class_pt{i} = 'Vert POD14';
color{i} = [0.5 0.5 1];
end
if contains(stored_names{i},'- V') & contains(stored_names{i},'POD30')
class_pt{i} = 'Vert POD30';
color{i} = [0.5 0.5 0.5];
end
end
%% t-SNE Visualization
Y = tsne(quantified,'Algorithm','exact','Distance','euclidian','Exaggeration',4);
figure
gscatter(Y(:,1),Y(:,2),class_pt')