forked from bnsreenu/python_for_microscopists
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path122_normalizing_HnE_images.py
169 lines (124 loc) · 6.04 KB
/
122_normalizing_HnE_images.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#!/usr/bin/env python
__author__ = "Sreenivas Bhattiprolu"
__license__ = "Feel free to copy, I appreciate if you acknowledge Python for Microscopists"
# https://youtu.be/yUrwEYgZUsA
"""
This code normalizes staining appearance of H&E stained images.
It also separates the hematoxylin and eosing stains in to different images.
Workflow based on the following papers:
A method for normalizing histology slides for quantitative analysis.
M. Macenko et al., ISBI 2009
http://wwwx.cs.unc.edu/~mn/sites/default/files/macenko2009.pdf
Efficient nucleus detector in histopathology images. J.P. Vink et al., J Microscopy, 2013
Original MATLAB code:
https://github.com/mitkovetta/staining-normalization/blob/master/normalizeStaining.m
Other useful references:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5226799/
https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0169875
PROPOSED WORKFLOW:
Input: RGB image
Step 1: Convert RGB to OD
Step 2: Remove data with OD intensity less than β
Step 3: Calculate singular value decomposition (SVD) on the OD tuples
Step 4: Create plane from the SVD directions corresponding to the
two largest singular values
Step 5: Project data onto the plane, and normalize to unit length
Step 6: Calculate angle of each point wrt the first SVD direction
Step 7: Find robust extremes (αth and (100−α)th 7 percentiles) of the
angle
Step 8: Convert extreme values back to OD space
Output: Optimal Stain Vectors
"""
import numpy as np
import cv2
from matplotlib import pyplot as plt
############### INPUT RGB IMAGE #######################
#Using opencv to read images may bemore robust compared to using skimage
#but need to remember to convert BGR to RGB.
#Also, convert to float later on and normalize to between 0 and 1.
#Image downloaded from:
#https://pbs.twimg.com/media/C1MkrgQWQAASbdz.jpg
img=cv2.imread('images/HnE_Image.jpg', 1)
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
Io = 240 # Transmitted light intensity, Normalizing factor for image intensities
alpha = 1 #As recommend in the paper. tolerance for the pseudo-min and pseudo-max (default: 1)
beta = 0.15 #As recommended in the paper. OD threshold for transparent pixels (default: 0.15)
######## Step 1: Convert RGB to OD ###################
## reference H&E OD matrix.
#Can be updated if you know the best values for your image.
#Otherwise use the following default values.
#Read the above referenced papers on this topic.
HERef = np.array([[0.5626, 0.2159],
[0.7201, 0.8012],
[0.4062, 0.5581]])
### reference maximum stain concentrations for H&E
maxCRef = np.array([1.9705, 1.0308])
# extract the height, width and num of channels of image
h, w, c = img.shape
# reshape image to multiple rows and 3 columns.
#Num of rows depends on the image size (wxh)
img = img.reshape((-1,3))
# calculate optical density
# OD = −log10(I)
#OD = -np.log10(img+0.004) #Use this when reading images with skimage
#Adding 0.004 just to avoid log of zero.
OD = -np.log10((img.astype(np.float)+1)/Io) #Use this for opencv imread
#Add 1 in case any pixels in the image have a value of 0 (log 0 is indeterminate)
"""
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(16, 8))
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(img[:,0],img[:,1],img[:,2])
ax2 = fig.add_subplot(122, projection='3d')
ax2.scatter(OD[:,0],OD[:,1],OD[:,2])
plt.show()
"""
############ Step 2: Remove data with OD intensity less than β ############
# remove transparent pixels (clear region with no tissue)
ODhat = OD[~np.any(OD < beta, axis=1)] #Returns an array where OD values are above beta
#Check by printing ODhat.min()
############# Step 3: Calculate SVD on the OD tuples ######################
#Estimate covariance matrix of ODhat (transposed)
# and then compute eigen values & eigenvectors.
eigvals, eigvecs = np.linalg.eigh(np.cov(ODhat.T))
######## Step 4: Create plane from the SVD directions with two largest values ######
#project on the plane spanned by the eigenvectors corresponding to the two
# largest eigenvalues
That = ODhat.dot(eigvecs[:,1:3]) #Dot product
############### Step 5: Project data onto the plane, and normalize to unit length ###########
############## Step 6: Calculate angle of each point wrt the first SVD direction ########
#find the min and max vectors and project back to OD space
phi = np.arctan2(That[:,1],That[:,0])
minPhi = np.percentile(phi, alpha)
maxPhi = np.percentile(phi, 100-alpha)
vMin = eigvecs[:,1:3].dot(np.array([(np.cos(minPhi), np.sin(minPhi))]).T)
vMax = eigvecs[:,1:3].dot(np.array([(np.cos(maxPhi), np.sin(maxPhi))]).T)
# a heuristic to make the vector corresponding to hematoxylin first and the
# one corresponding to eosin second
if vMin[0] > vMax[0]:
HE = np.array((vMin[:,0], vMax[:,0])).T
else:
HE = np.array((vMax[:,0], vMin[:,0])).T
# rows correspond to channels (RGB), columns to OD values
Y = np.reshape(OD, (-1, 3)).T
# determine concentrations of the individual stains
C = np.linalg.lstsq(HE,Y, rcond=None)[0]
# normalize stain concentrations
maxC = np.array([np.percentile(C[0,:], 99), np.percentile(C[1,:],99)])
tmp = np.divide(maxC,maxCRef)
C2 = np.divide(C,tmp[:, np.newaxis])
###### Step 8: Convert extreme values back to OD space
# recreate the normalized image using reference mixing matrix
Inorm = np.multiply(Io, np.exp(-HERef.dot(C2)))
Inorm[Inorm>255] = 254
Inorm = np.reshape(Inorm.T, (h, w, 3)).astype(np.uint8)
# Separating H and E components
H = np.multiply(Io, np.exp(np.expand_dims(-HERef[:,0], axis=1).dot(np.expand_dims(C2[0,:], axis=0))))
H[H>255] = 254
H = np.reshape(H.T, (h, w, 3)).astype(np.uint8)
E = np.multiply(Io, np.exp(np.expand_dims(-HERef[:,1], axis=1).dot(np.expand_dims(C2[1,:], axis=0))))
E[E>255] = 254
E = np.reshape(E.T, (h, w, 3)).astype(np.uint8)
plt.imsave("images/HnE_normalized.jpg", Inorm)
plt.imsave("images/HnE_separated_H.jpg", H)
plt.imsave("images/HnE_separated_E.jpg", E)