|
| 1 | +import numpy as np |
| 2 | + |
| 3 | +class NumpyMultiMacenkoNormalizer: |
| 4 | + def __init__(self, norm_mode="avg-post"): |
| 5 | + self.norm_mode = norm_mode |
| 6 | + self.HERef = np.array([[0.5626, 0.2159], |
| 7 | + [0.7201, 0.8012], |
| 8 | + [0.4062, 0.5581]]) |
| 9 | + self.maxCRef = np.array([1.9705, 1.0308]) |
| 10 | + |
| 11 | + def __convert_rgb2od(self, I, Io, beta): |
| 12 | + I = np.transpose(I, (1, 2, 0)) |
| 13 | + OD = -np.log((I.reshape(-1, I.shape[-1]).astype(float) + 1) / Io) |
| 14 | + ODhat = OD[~np.any(OD < beta, axis=1)] |
| 15 | + return OD, ODhat |
| 16 | + |
| 17 | + def __find_phi_bounds(self, ODhat, eigvecs, alpha): |
| 18 | + That = np.dot(ODhat, eigvecs) |
| 19 | + phi = np.arctan2(That[:, 1], That[:, 0]) |
| 20 | + |
| 21 | + minPhi = np.percentile(phi, alpha) |
| 22 | + maxPhi = np.percentile(phi, 100 - alpha) |
| 23 | + |
| 24 | + return minPhi, maxPhi |
| 25 | + |
| 26 | + def __find_HE_from_bounds(self, eigvecs, minPhi, maxPhi): |
| 27 | + vMin = np.dot(eigvecs, [np.cos(minPhi), np.sin(minPhi)]).reshape(-1, 1) |
| 28 | + vMax = np.dot(eigvecs, [np.cos(maxPhi), np.sin(maxPhi)]).reshape(-1, 1) |
| 29 | + |
| 30 | + HE = np.concatenate([vMin, vMax], axis=1) if vMin[0] > vMax[0] else np.concatenate([vMax, vMin], axis=1) |
| 31 | + return HE |
| 32 | + |
| 33 | + def __find_HE(self, ODhat, eigvecs, alpha): |
| 34 | + minPhi, maxPhi = self.__find_phi_bounds(ODhat, eigvecs, alpha) |
| 35 | + return self.__find_HE_from_bounds(eigvecs, minPhi, maxPhi) |
| 36 | + |
| 37 | + def __find_concentration(self, OD, HE): |
| 38 | + Y = OD.T |
| 39 | + C, _, _, _ = np.linalg.lstsq(HE, Y, rcond=None) |
| 40 | + return C |
| 41 | + |
| 42 | + def __compute_matrices_single(self, I, Io, alpha, beta): |
| 43 | + OD, ODhat = self.__convert_rgb2od(I, Io, beta) |
| 44 | + |
| 45 | + cov_matrix = np.cov(ODhat.T) |
| 46 | + eigvals, eigvecs = np.linalg.eigh(cov_matrix) |
| 47 | + eigvecs = eigvecs[:, [1, 2]] |
| 48 | + |
| 49 | + HE = self.__find_HE(ODhat, eigvecs, alpha) |
| 50 | + C = self.__find_concentration(OD, HE) |
| 51 | + maxC = np.array([np.percentile(C[0, :], 99), np.percentile(C[1, :], 99)]) |
| 52 | + |
| 53 | + return HE, C, maxC |
| 54 | + |
| 55 | + def fit(self, Is, Io=240, alpha=1, beta=0.15): |
| 56 | + if self.norm_mode == "avg-post": |
| 57 | + HEs, _, maxCs = zip(*[self.__compute_matrices_single(I, Io, alpha, beta) for I in Is]) |
| 58 | + |
| 59 | + self.HERef = np.mean(HEs, axis=0) |
| 60 | + self.maxCRef = np.mean(maxCs, axis=0) |
| 61 | + elif self.norm_mode == "concat": |
| 62 | + ODs, ODhats = zip(*[self.__convert_rgb2od(I, Io, beta) for I in Is]) |
| 63 | + OD = np.vstack(ODs) |
| 64 | + ODhat = np.vstack(ODhats) |
| 65 | + |
| 66 | + cov_matrix = np.cov(ODhat.T) |
| 67 | + eigvals, eigvecs = np.linalg.eigh(cov_matrix) |
| 68 | + eigvecs = eigvecs[:, [1, 2]] |
| 69 | + |
| 70 | + HE = self.__find_HE(ODhat, eigvecs, alpha) |
| 71 | + C = self.__find_concentration(OD, HE) |
| 72 | + maxCs = np.array([np.percentile(C[0, :], 99), np.percentile(C[1, :], 99)]) |
| 73 | + |
| 74 | + self.HERef = HE |
| 75 | + self.maxCRef = maxCs |
| 76 | + elif self.norm_mode == "avg-pre": |
| 77 | + ODs, ODhats = zip(*[self.__convert_rgb2od(I, Io, beta) for I in Is]) |
| 78 | + |
| 79 | + covs = [np.cov(ODhat.T) for ODhat in ODhats] |
| 80 | + eigvecs = np.mean([np.linalg.eigh(cov)[1][:, [1, 2]] for cov in covs], axis=0) |
| 81 | + |
| 82 | + OD = np.vstack(ODs) |
| 83 | + ODhat = np.vstack(ODhats) |
| 84 | + |
| 85 | + HE = self.__find_HE(ODhat, eigvecs, alpha) |
| 86 | + C = self.__find_concentration(OD, HE) |
| 87 | + maxCs = np.array([np.percentile(C[0, :], 99), np.percentile(C[1, :], 99)]) |
| 88 | + |
| 89 | + self.HERef = HE |
| 90 | + self.maxCRef = maxCs |
| 91 | + elif self.norm_mode in ["fixed-single", "stochastic-single"]: |
| 92 | + self.HERef, _, self.maxCRef = self.__compute_matrices_single(Is[0], Io, alpha, beta) |
| 93 | + else: |
| 94 | + raise ValueError("Unknown norm mode") |
| 95 | + |
| 96 | + def normalize(self, I, Io=240, alpha=1, beta=0.15, stains=True): |
| 97 | + c, h, w = I.shape |
| 98 | + |
| 99 | + HE, C, maxC = self.__compute_matrices_single(I, Io, alpha, beta) |
| 100 | + C = (self.maxCRef / maxC).reshape(-1, 1) * C |
| 101 | + |
| 102 | + Inorm = Io * np.exp(-np.dot(self.HERef, C)) |
| 103 | + Inorm[Inorm > 255] = 255 |
| 104 | + Inorm = np.transpose(Inorm, (1, 0)).reshape(h, w, c).astype(np.int32) |
| 105 | + |
| 106 | + H, E = None, None |
| 107 | + |
| 108 | + if stains: |
| 109 | + H = Io * np.exp(-np.dot(self.HERef[:, 0].reshape(-1, 1), C[0, :].reshape(1, -1))) |
| 110 | + H[H > 255] = 255 |
| 111 | + H = np.transpose(H, (1, 0)).reshape(h, w, c).astype(np.int32) |
| 112 | + |
| 113 | + E = Io * np.exp(-np.dot(self.HERef[:, 1].reshape(-1, 1), C[1, :].reshape(1, -1))) |
| 114 | + E[E > 255] = 255 |
| 115 | + E = np.transpose(E, (1, 0)).reshape(h, w, c).astype(np.int32) |
| 116 | + |
| 117 | + return Inorm, H, E |
0 commit comments