MATLAB光密度,H&；E染色切片patch的标准化算法

最近读了一片论文A method for normalizing histology slides for quantitative analysis论文下载地址http://wwwx.cs.unc.edu/~mn/sites/default/files/macenko2009.pdf

博客上先记录一下matlab和python语法的区别

论文里给出详细的解答步骤

Input: RGB Slide

1. Convert RGB to OD

2. Remove data with OD intensity less thanβ

3. Calculate SVD on the OD tuples

4. Create plane from the SVD directions corresponding to the two largest singular values

5.Project data onto the plane, and normalize to unit length

6. Calculate angle of each point wrt the first SVD direction

7. Find robust extremes (αth and(100−α)th percentiles) of the angle

8. Convert extreme values back to OD space

Output

: Optimal Stain Vectors

首先就需要将所有的颜色值转换成对应的光密度值(optical density，记为OD)，对应转换的公式

OD = -log(I).其中这里的I是rgb颜色空间向量，并且归一化到(0-1)

matlab代码首先将I三个通道拉直，也就是将rgb三个颜色空间拉成一个(w*h，3)的矩阵matlab代码如下

I = double(I);

I = reshape(I, [], 3);

python运用numpy的reshape函数

I = I.reshape((w*h,3))

下面一步将图像从颜色空间向量转换成对应的光密度空间向量，matlab代码如下

OD = -log((I+1)/Io);

python里利用numpy的log函数实现这一步操作

OD = -np.log((I+1)/Io)

第二步、就需要将OD敏感度小于β的删除，这个值默认设置为0.15

matlab代码如下

ODhat = OD(~any(OD < beta, 2), :);

主要是为了OD列里面第2 dim中小于beta的去掉

python代码和matlab相似

ODhat = OD[~(OD[:,2] < beta), :]

第三步、对这个ODhat进行SVD分解，也就是奇异值分解，matlab代码如下

[V, ~] = eig(cov(ODhat));

python中的cov函数主要是利用nump中的np.cov函数，这里用法和matlab有很大的区别，matlab中cov(A)中的A是以列向量为单元，最终解的是(列*列的矩阵)，

而python是以行为单元，最终解的是(行*行)的矩阵，比如我输入(65525*3)的矩阵，matlab会得出(3*3)的协方差，而python则会得出(65525*65525)的值，

这里同样有和matlab不同之处，matlab得到的特征向量是按照特征值从小到大排布的，而python是按照从大到小排布的，下面是python的代码

_ , V = np.linalg.eig(np.cov(ODhat.T))

第四步、矩阵相乘，matlab直接

That = ODhat*V(:,2:3);

而python的*是点乘，不是矩阵的乘，这里我调用了numpy中的np.dot函数

That = np.dot(ODhat,V[:,0:2])

得到的That是一个n*2的矩阵，然后计算方向，目前也不是很清楚这个方向是什么意思，matlab里直接利用

phi = atan2(That(:,2), That(:,1));

python里我利用math.atan2这个包，math.atan2(x,y)就是求的点(x，y)的角，但是这个函数只能接受单个值，所一我采用了最简单的方法放这个函数

m ,_ = That.shape

phi = np.zeros((m,1))

for i in range(int(m)):

phi[i] = math.atan2(That[i,1], That[i,0])

python 是利用np.percentile(phi, alpha).

A\B这个说法，在python里我没有发现，无奈只能用其他等式代替，A\B = A逆*B，但是这里A不是非奇异矩阵，所以不能用inv函数，于是我求解一个广义逆矩阵

利用

np.linalg.pinv(HE)

可以求得HE的广义逆矩阵，总之，复习了一下矩阵论，以及相关的函数