Repository: Project-MONAI/MONAI
File: monai/apps/pathology/transforms/stain/array.py (Python)

Code before:
        absorbance = -np.log(image.clip(max=self.tli) / self.tli)

        # remove transparent pixels
        absorbance_hat = absorbance[np.all(absorbance > self.beta, axis=1)]
        if len(absorbance_hat) == 0:
            raise ValueError("All pixels of the input image are below the absorbance threshold.")

        # compute eigenvectors
        _, eigvecs = np.linalg.eigh(np.cov(absorbance_hat.T).astype(np.float32, copy=False))

        # project on the plane spanned by the eigenvectors corresponding to the two largest eigenvalues
        t_hat = absorbance_hat.dot(eigvecs[:, 1:3])

        # find the min and max vectors and project back to absorbance space
        phi = np.arctan2(t_hat[:, 1], t_hat[:, 0])
        min_phi = np.percentile(phi, self.alpha)
        max_phi = np.percentile(phi, 100 - self.alpha)
        v_min = eigvecs[:, 1:3].dot(np.array([(np.cos(min_phi), np.sin(min_phi))], dtype=np.float32).T)
        v_max = eigvecs[:, 1:3].dot(np.array([(np.cos(max_phi), np.sin(max_phi))], dtype=np.float32).T)

        # a heuristic to make the vector corresponding to hematoxylin first and the one corresponding to eosin second
        # Hematoxylin: high blue, lower red (low R/B ratio)
        # Eosin: high red, lower blue (high R/B ratio)
        ε = np.finfo(np.float32).eps
        v_min_rb_ratio = v_min[0, 0] / (v_min[2, 0] + ε)
        v_max_rb_ratio = v_max[0, 0] / (v_max[2, 0] + ε)
        if v_min_rb_ratio < v_max_rb_ratio:
            he = np.array((v_min[:, 0], v_max[:, 0]), dtype=np.float32).T
        else:
            he = np.array((v_max[:, 0], v_min[:, 0]), dtype=np.float32).T

        return he

    def __call__(self, image: np.ndarray) -> np.ndarray:
        """Perform stain extraction.

        Args:
            image: uint8 RGB image to extract stain from

        return:
            target_he: H&E absorbance matrix for the image (first column is H, second column is E, rows are RGB values)
        """
        if not isinstance(image, np.ndarray):
            raise TypeError("Image must be of type numpy.ndarray.")

        target_he = self._deconvolution_extract_stain(image)
        return target_he


class NormalizeHEStains(Transform):
    """Class to normalize patches/images to a reference or target image stain (see Note).

Review comment (performance):
```suggestion
        eps = np.finfo(np.float32).eps
        v_min_rb_ratio = v_min[0] / (v_min[2] + eps)
        v_max_rb_ratio = v_max[0] / (v_max[2] + eps)
```

Code after:
        absorbance = -np.log(image.clip(max=self.tli) / self.tli)

        # remove transparent pixels
        absorbance_hat = absorbance[np.all(absorbance > self.beta, axis=1)]
        if len(absorbance_hat) == 0:
            raise ValueError("All pixels of the input image are below the absorbance threshold.")

        # compute eigenvectors
        _, eigvecs = np.linalg.eigh(np.cov(absorbance_hat.T).astype(np.float32, copy=False))

        # project on the plane spanned by the eigenvectors corresponding to the two largest eigenvalues
        t_hat = absorbance_hat.dot(eigvecs[:, 1:3])

        # find the min and max vectors and project back to absorbance space
        phi = np.arctan2(t_hat[:, 1], t_hat[:, 0])
        min_phi = np.percentile(phi, self.alpha)
        max_phi = np.percentile(phi, 100 - self.alpha)
        v_min = eigvecs[:, 1:3].dot(np.array([(np.cos(min_phi), np.sin(min_phi))], dtype=np.float32).T)
        v_max = eigvecs[:, 1:3].dot(np.array([(np.cos(max_phi), np.sin(max_phi))], dtype=np.float32).T)

        # a heuristic to make the vector corresponding to hematoxylin first and the one corresponding to eosin second
        # Hematoxylin: high blue, lower red (low R/B ratio)
        # Eosin: high red, lower blue (high R/B ratio)
        eps = np.finfo(np.float32).eps
        v_min_rb_ratio = v_min[0] / (v_min[2] + eps)
        v_max_rb_ratio = v_max[0] / (v_max[2] + eps)
        if v_min_rb_ratio < v_max_rb_ratio:
            he = np.array((v_min[:, 0], v_max[:, 0]), dtype=np.float32).T
        else:
            he = np.array((v_max[:, 0], v_min[:, 0]), dtype=np.float32).T

        return he

    def __call__(self, image: np.ndarray) -> np.ndarray:
        """Perform stain extraction.

        Args:
            image: uint8 RGB image to extract stain from

        return:
            target_he: H&E absorbance matrix for the image (first column is H, second column is E, rows are RGB values)
        """
        if not isinstance(image, np.ndarray):
            raise TypeError("Image must be of type numpy.ndarray.")

        target_he = self._deconvolution_extract_stain(image)
        return target_he


class NormalizeHEStains(Transform):
    """Class to normalize patches/images to a reference or target image stain (see Note).