|
| 1 | +# Copyright 2020 - 2021 MONAI Consortium |
| 2 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 3 | +# you may not use this file except in compliance with the License. |
| 4 | +# You may obtain a copy of the License at |
| 5 | +# http://www.apache.org/licenses/LICENSE-2.0 |
| 6 | +# Unless required by applicable law or agreed to in writing, software |
| 7 | +# distributed under the License is distributed on an "AS IS" BASIS, |
| 8 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 9 | +# See the License for the specific language governing permissions and |
| 10 | +# limitations under the License. |
| 11 | + |
| 12 | +from typing import Union |
| 13 | + |
| 14 | +import numpy as np |
| 15 | + |
| 16 | +from monai.transforms.transform import Transform |
| 17 | + |
| 18 | + |
| 19 | +class ExtractHEStains(Transform): |
| 20 | + """Class to extract a target stain from an image, using stain deconvolution (see Note). |
| 21 | +
|
| 22 | + Args: |
| 23 | + tli: transmitted light intensity. Defaults to 240. |
| 24 | + alpha: tolerance in percentile for the pseudo-min (alpha percentile) |
| 25 | + and pseudo-max (100 - alpha percentile). Defaults to 1. |
| 26 | + beta: absorbance threshold for transparent pixels. Defaults to 0.15 |
| 27 | + max_cref: reference maximum stain concentrations for Hematoxylin & Eosin (H&E). |
| 28 | + Defaults to (1.9705, 1.0308). |
| 29 | +
|
| 30 | + Note: |
| 31 | + For more information refer to: |
| 32 | + - the original paper: Macenko et al., 2009 http://wwwx.cs.unc.edu/~mn/sites/default/files/macenko2009.pdf |
| 33 | + - the previous implementations: |
| 34 | +
|
| 35 | + - MATLAB: https://github.com/mitkovetta/staining-normalization |
| 36 | + - Python: https://github.com/schaugf/HEnorm_python |
| 37 | + """ |
| 38 | + |
| 39 | + def __init__( |
| 40 | + self, |
| 41 | + tli: float = 240, |
| 42 | + alpha: float = 1, |
| 43 | + beta: float = 0.15, |
| 44 | + max_cref: Union[tuple, np.ndarray] = (1.9705, 1.0308), |
| 45 | + ) -> None: |
| 46 | + self.tli = tli |
| 47 | + self.alpha = alpha |
| 48 | + self.beta = beta |
| 49 | + self.max_cref = np.array(max_cref) |
| 50 | + |
| 51 | + def _deconvolution_extract_stain(self, image: np.ndarray) -> np.ndarray: |
| 52 | + """Perform Stain Deconvolution and return stain matrix for the image. |
| 53 | +
|
| 54 | + Args: |
| 55 | + img: uint8 RGB image to perform stain deconvolution on |
| 56 | +
|
| 57 | + Return: |
| 58 | + he: H&E absorbance matrix for the image (first column is H, second column is E, rows are RGB values) |
| 59 | + """ |
| 60 | + # check image type and vlues |
| 61 | + if not isinstance(image, np.ndarray): |
| 62 | + raise TypeError("Image must be of type numpy.ndarray.") |
| 63 | + if image.min() < 0: |
| 64 | + raise ValueError("Image should not have negative values.") |
| 65 | + if image.max() > 255: |
| 66 | + raise ValueError("Image should not have values greater than 255.") |
| 67 | + |
| 68 | + # reshape image and calculate absorbance |
| 69 | + image = image.reshape((-1, 3)) |
| 70 | + image = image.astype(np.float32) + 1.0 |
| 71 | + absorbance = -np.log(image.clip(max=self.tli) / self.tli) |
| 72 | + |
| 73 | + # remove transparent pixels |
| 74 | + absorbance_hat = absorbance[np.all(absorbance > self.beta, axis=1)] |
| 75 | + if len(absorbance_hat) == 0: |
| 76 | + raise ValueError("All pixels of the input image are below the absorbance threshold.") |
| 77 | + |
| 78 | + # compute eigenvectors |
| 79 | + _, eigvecs = np.linalg.eigh(np.cov(absorbance_hat.T).astype(np.float32)) |
| 80 | + |
| 81 | + # project on the plane spanned by the eigenvectors corresponding to the two largest eigenvalues |
| 82 | + t_hat = absorbance_hat.dot(eigvecs[:, 1:3]) |
| 83 | + |
| 84 | + # find the min and max vectors and project back to absorbance space |
| 85 | + phi = np.arctan2(t_hat[:, 1], t_hat[:, 0]) |
| 86 | + min_phi = np.percentile(phi, self.alpha) |
| 87 | + max_phi = np.percentile(phi, 100 - self.alpha) |
| 88 | + v_min = eigvecs[:, 1:3].dot(np.array([(np.cos(min_phi), np.sin(min_phi))], dtype=np.float32).T) |
| 89 | + v_max = eigvecs[:, 1:3].dot(np.array([(np.cos(max_phi), np.sin(max_phi))], dtype=np.float32).T) |
| 90 | + |
| 91 | + # a heuristic to make the vector corresponding to hematoxylin first and the one corresponding to eosin second |
| 92 | + if v_min[0] > v_max[0]: |
| 93 | + he = np.array((v_min[:, 0], v_max[:, 0]), dtype=np.float32).T |
| 94 | + else: |
| 95 | + he = np.array((v_max[:, 0], v_min[:, 0]), dtype=np.float32).T |
| 96 | + |
| 97 | + return he |
| 98 | + |
| 99 | + def __call__(self, image: np.ndarray) -> np.ndarray: |
| 100 | + """Perform stain extraction. |
| 101 | +
|
| 102 | + Args: |
| 103 | + image: uint8 RGB image to extract stain from |
| 104 | +
|
| 105 | + return: |
| 106 | + target_he: H&E absorbance matrix for the image (first column is H, second column is E, rows are RGB values) |
| 107 | + """ |
| 108 | + if not isinstance(image, np.ndarray): |
| 109 | + raise TypeError("Image must be of type numpy.ndarray.") |
| 110 | + |
| 111 | + target_he = self._deconvolution_extract_stain(image) |
| 112 | + return target_he |
| 113 | + |
| 114 | + |
| 115 | +class NormalizeHEStains(Transform): |
| 116 | + """Class to normalize patches/images to a reference or target image stain (see Note). |
| 117 | +
|
| 118 | + Performs stain deconvolution of the source image using the ExtractHEStains |
| 119 | + class, to obtain the stain matrix and calculate the stain concentration matrix |
| 120 | + for the image. Then, performs the inverse Beer-Lambert transform to recreate the |
| 121 | + patch using the target H&E stain matrix provided. If no target stain provided, a default |
| 122 | + reference stain is used. Similarly, if no maximum stain concentrations are provided, a |
| 123 | + reference maximum stain concentrations matrix is used. |
| 124 | +
|
| 125 | + Args: |
| 126 | + tli: transmitted light intensity. Defaults to 240. |
| 127 | + alpha: tolerance in percentile for the pseudo-min (alpha percentile) and |
| 128 | + pseudo-max (100 - alpha percentile). Defaults to 1. |
| 129 | + beta: absorbance threshold for transparent pixels. Defaults to 0.15. |
| 130 | + target_he: target stain matrix. Defaults to ((0.5626, 0.2159), (0.7201, 0.8012), (0.4062, 0.5581)). |
| 131 | + max_cref: reference maximum stain concentrations for Hematoxylin & Eosin (H&E). |
| 132 | + Defaults to [1.9705, 1.0308]. |
| 133 | +
|
| 134 | + Note: |
| 135 | + For more information refer to: |
| 136 | + - the original paper: Macenko et al., 2009 http://wwwx.cs.unc.edu/~mn/sites/default/files/macenko2009.pdf |
| 137 | + - the previous implementations: |
| 138 | +
|
| 139 | + - MATLAB: https://github.com/mitkovetta/staining-normalization |
| 140 | + - Python: https://github.com/schaugf/HEnorm_python |
| 141 | + """ |
| 142 | + |
| 143 | + def __init__( |
| 144 | + self, |
| 145 | + tli: float = 240, |
| 146 | + alpha: float = 1, |
| 147 | + beta: float = 0.15, |
| 148 | + target_he: Union[tuple, np.ndarray] = ((0.5626, 0.2159), (0.7201, 0.8012), (0.4062, 0.5581)), |
| 149 | + max_cref: Union[tuple, np.ndarray] = (1.9705, 1.0308), |
| 150 | + ) -> None: |
| 151 | + self.tli = tli |
| 152 | + self.target_he = np.array(target_he) |
| 153 | + self.max_cref = np.array(max_cref) |
| 154 | + self.stain_extractor = ExtractHEStains(tli=self.tli, alpha=alpha, beta=beta, max_cref=self.max_cref) |
| 155 | + |
| 156 | + def __call__(self, image: np.ndarray) -> np.ndarray: |
| 157 | + """Perform stain normalization. |
| 158 | +
|
| 159 | + Args: |
| 160 | + image: uint8 RGB image/patch to be stain normalized, pixel values between 0 and 255 |
| 161 | +
|
| 162 | + Return: |
| 163 | + image_norm: stain normalized image/patch |
| 164 | + """ |
| 165 | + # check image type and vlues |
| 166 | + if not isinstance(image, np.ndarray): |
| 167 | + raise TypeError("Image must be of type numpy.ndarray.") |
| 168 | + if image.min() < 0: |
| 169 | + raise ValueError("Image should not have negative values.") |
| 170 | + if image.max() > 255: |
| 171 | + raise ValueError("Image should not have values greater than 255.") |
| 172 | + |
| 173 | + # extract stain of the image |
| 174 | + he = self.stain_extractor(image) |
| 175 | + |
| 176 | + # reshape image and calculate absorbance |
| 177 | + h, w, _ = image.shape |
| 178 | + image = image.reshape((-1, 3)) |
| 179 | + image = image.astype(np.float32) + 1.0 |
| 180 | + absorbance = -np.log(image.clip(max=self.tli) / self.tli) |
| 181 | + |
| 182 | + # rows correspond to channels (RGB), columns to absorbance values |
| 183 | + y = np.reshape(absorbance, (-1, 3)).T |
| 184 | + |
| 185 | + # determine concentrations of the individual stains |
| 186 | + conc = np.linalg.lstsq(he, y, rcond=None)[0] |
| 187 | + |
| 188 | + # normalize stain concentrations |
| 189 | + max_conc = np.array([np.percentile(conc[0, :], 99), np.percentile(conc[1, :], 99)], dtype=np.float32) |
| 190 | + tmp = np.divide(max_conc, self.max_cref, dtype=np.float32) |
| 191 | + image_c = np.divide(conc, tmp[:, np.newaxis], dtype=np.float32) |
| 192 | + |
| 193 | + image_norm: np.ndarray = np.multiply(self.tli, np.exp(-self.target_he.dot(image_c)), dtype=np.float32) |
| 194 | + image_norm[image_norm > 255] = 254 |
| 195 | + image_norm = np.reshape(image_norm.T, (h, w, 3)).astype(np.uint8) |
| 196 | + return image_norm |
0 commit comments