From 306e0322d86ad4c845241ba791c179989dd8bc2c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 17:58:19 +0100 Subject: [PATCH 01/14] added vahadane to API --- torchstain/base/normalizers/vahadane.py | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 torchstain/base/normalizers/vahadane.py diff --git a/torchstain/base/normalizers/vahadane.py b/torchstain/base/normalizers/vahadane.py new file mode 100644 index 0000000..b88de7e --- /dev/null +++ b/torchstain/base/normalizers/vahadane.py @@ -0,0 +1,10 @@ +def VahadaneNormalizer(backend='numpy'): + if backend == 'numpy': + from torchstain.numpy.normalizers import VahadaneReinhardNormalizer + return NumpyVahadaneNormalizer() + elif backend == "torch": + raise NotImplementedError + elif backend == "tensorflow": + raise NotImplementedError + else: + raise Exception(f'Unknown backend {backend}') From 82fedba9a2c0a9a98bd7dceb40b2fe5594806b41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 17:58:43 +0100 Subject: [PATCH 02/14] draft vahadane implementation --- torchstain/numpy/normalizers/vahadane.py | 27 ++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 torchstain/numpy/normalizers/vahadane.py diff --git a/torchstain/numpy/normalizers/vahadane.py b/torchstain/numpy/normalizers/vahadane.py new file mode 100644 index 0000000..da46871 --- /dev/null +++ b/torchstain/numpy/normalizers/vahadane.py @@ -0,0 +1,27 @@ +import numpy as np +from torchstain.base.normalizers import HENormalizer +from torchstain.numpy.utils.lars import lars +from torchstain.numpy.utils.stats import standardize_brightness +from torchstain.numpy.utils.extract import get_stain_matrix + +""" +Source code adapted from: https://github.com/wanghao14/Stain_Normalization/blob/master/stainNorm_Vahadane.py +""" +class NumpyVahadaneNormalizer(HENormalizer): + def __init__(self): + super().__init__() + self.stain_matrix_target = None + + def fit(self, target): + target = standardize_brightness(target) + self.stain_matrix_target = get_stain_matrix(target) + + def normalize(self, I): + I = standardize_brightness(I) + stain_matrix = get_stain_matrix(I) + concentrations = get_concentrations(I, stain_matrix) + out = np.exp(-1 * np.dot(concentrations, self.stain_matrix_target).reshape(I.shape)) + return (255 * out).astype("uint8") + + + From a35f27e2f007bf8646fffa10aa5d34147f0f9620 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 17:59:02 +0100 Subject: [PATCH 03/14] lasso numpy implementation --- torchstain/numpy/utils/lasso.py | 52 +++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 torchstain/numpy/utils/lasso.py diff --git a/torchstain/numpy/utils/lasso.py b/torchstain/numpy/utils/lasso.py new file mode 100644 index 0000000..17c090b --- /dev/null +++ b/torchstain/numpy/utils/lasso.py @@ -0,0 +1,52 @@ +import numpy as np + + +""" +LASSO implementation was adapted from: +https://www.kaggle.com/code/mcweng24/lasso-regression-using-numpy +""" +def predicted_values(X, w): + return np.matmul(X, w) + + +def rho_compute(y, X, w, j): + X_k = np.delete(X, j, 1) + w_k = np.delete(w, j) + predict_k = predicted_values(X_k, w_k) + residual = y - predict_k + return np.sum(X[:, j] * residual) + + +def z_compute(X): + z_vector = np.sum(X * X, axis=0) + return np.sum(X * X, axis = 0) + + +def coordinate_descent(y, X, w, alpha, z, tolerance): + max_step = 100 + iteration = 0 + while max_step > tolerance: + iteration += 1 + old_weights = np.copy(w) + for j in range(len(w)): + rho_j = rho_compute(y, X, w, j) + if j == 0: + w[j] = rho_j/z[j] + elif rho_j < -alpha * len(y): + w[j] = (rho_j + (alpha * len(y))) / z[j] + elif rho_j > -alpha * len(y) and rho_j < alpha * len(y): + w[j] = 0. + elif rho_j > alpha * len(y): + w[j] = (rho_j - (alpha * len(y))) / z[j] + else: + w[j] = np.NaN + step_sizes = np.abs(old_weights - w) + max_step = step_sizes.max() + return w, iteration, max_step + + +def lasso(x, y, alpha=0.1, tol=0.0001): + w = np.zeros(X.shape[1], dtype="float32") + z = z_compute(x) + w_opt, iterations, max_step = coordinate_descent(y, x, w, alpha, z, tolerance) + return w_opt From c1c860ee34a0aa0f739a5b2a8f491bb6420a8adc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 17:59:23 +0100 Subject: [PATCH 04/14] OD2RGB numpy converter --- torchstain/numpy/utils/od2rgb.py | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 torchstain/numpy/utils/od2rgb.py diff --git a/torchstain/numpy/utils/od2rgb.py b/torchstain/numpy/utils/od2rgb.py new file mode 100644 index 0000000..58bc46f --- /dev/null +++ b/torchstain/numpy/utils/od2rgb.py @@ -0,0 +1,4 @@ +import numpy as np + +def OD_to_RGB(OD): + return (255 * np.exp(-1 * OD)).astype(np.uint8) From e7fe99307fb916fe335fa494e0b651fd11e3219e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 17:59:33 +0100 Subject: [PATCH 05/14] RGB2OD numpy converter --- torchstain/numpy/utils/rgb2od.py | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 torchstain/numpy/utils/rgb2od.py diff --git a/torchstain/numpy/utils/rgb2od.py b/torchstain/numpy/utils/rgb2od.py new file mode 100644 index 0000000..9721829 --- /dev/null +++ b/torchstain/numpy/utils/rgb2od.py @@ -0,0 +1,9 @@ +import numpy as np + + +def RGB_to_OD(I): + # remove zeros + x[x == 0] = 1 + + # convert to OD and return + return -1 * np.log(I / 255) \ No newline at end of file From a0f7cf0793672b972af9f0d663b80db34aaf599b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 17:59:44 +0100 Subject: [PATCH 06/14] numpy standardize brightness --- torchstain/numpy/utils/stats.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/torchstain/numpy/utils/stats.py b/torchstain/numpy/utils/stats.py index b43ebe1..44b684c 100644 --- a/torchstain/numpy/utils/stats.py +++ b/torchstain/numpy/utils/stats.py @@ -5,3 +5,7 @@ def get_mean_std(I): def standardize(x, mu, std): return (x - mu) / std + +def standardize_brightness(x): + p = np.percentile(x, 90) + return np.clip(x * 255 / p, 0, 255).astype("uint8") From cdc64311f618d3975675b9a7094de7845cfae1bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 18:00:26 +0100 Subject: [PATCH 07/14] almost finished vahadane -DL and LarsLasso remains --- torchstain/numpy/utils/extract.py | 35 +++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 torchstain/numpy/utils/extract.py diff --git a/torchstain/numpy/utils/extract.py b/torchstain/numpy/utils/extract.py new file mode 100644 index 0000000..8d460d4 --- /dev/null +++ b/torchstain/numpy/utils/extract.py @@ -0,0 +1,35 @@ +import numpy as np +from torchstain.numpy.utils.rgb2od import rgb2od +from torchstain.numpy.utils.rgb2lab import rgb2lab +from torchstain.numpy.utils.lasso import lasso + +def extract_tissue(I, th=0.8): + LAB = rgb2lab(I) + L = I_LAB[:, :, 0] / 255.0 + return L < th + +def get_stain_matrix(I, th=0.8, lamda=0.1): + # convert RGB -> OD and flatten channel-wise + OD = rgb2od(I).reshape((-1, 3)) + + # detect glass and remove it from OD image + mask = extract_tissue(I, th).reshape((-1,)) + OD = OD[mask] + + # perform dictionary learning + # dictionary = spams.trainDL(OD.T, K=2, lambda1=lamda, mode=2, modeD=0, posAlpha=True, posD=True, verbose=False) + dictionary = None # @TODO: Implement DL train method + dictionary = dictionary.T + if dictionary[0, 0] < dictionary[1, 0]: + dictionary = dictionary[[1, 0], :] + + # normalize rows and return result + return dictionary / np.linalg.norm(dictionary, axis=1)[:, None] + +def get_concentrations(I, stain_matrix, lamda=0.01): + # convert RGB -> OD and flatten channel-wise + OD = rgb2od(I).reshape((-1, 3)) + + # perform LASSO regression + #return spams.lasso(OD.T, D=stain_matrix.T, mode=2, lambda1=lamda, pos=True).toarray().T + return lasso(OD.T, y=stain_matrix.T).T # @TODO: Implement LARS-LASSO From 603e4b0357647e39412f8186b59195e1230d36b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 18:04:42 +0100 Subject: [PATCH 08/14] updated README regarding vahadane support --- README.md | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 9047b49..6041877 100644 --- a/README.md +++ b/README.md @@ -10,6 +10,7 @@ Normalization algorithms currently implemented: - Macenko et al. [\[1\]](#reference) (ported from [numpy implementation](https://github.com/schaugf/HEnorm_python)) - Reinhard et al. [\[2\]](#reference) (only numpy & TensorFlow backend support) +- Vahadane et al. [\[3\]](#reference) (only numpy backend support) ## Installation @@ -50,6 +51,7 @@ norm, H, E = normalizer.normalize(I=t_to_transform, stains=True) |-|-|-|-| | Macenko | ✓ | ✓ | ✓ | | Reinhard | ✓ | ✗ | ✓ | +| Vahadane | ✓ | ✗ | ✗ | ## Backend comparison @@ -68,8 +70,9 @@ Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz ## Reference -- [1] Macenko, Marc, et al. "A method for normalizing histology slides for quantitative analysis." 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. IEEE, 2009. -- [2] Reinhard, Erik, et al. "Color transfer between images." IEEE Computer Graphics and Applications. IEEE, 2001. +- [1] Macenko, Marc et al. "A method for normalizing histology slides for quantitative analysis." 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. IEEE, 2009. +- [2] Reinhard, Erik et al. "Color transfer between images." IEEE Computer Graphics and Applications. IEEE, 2001. +- [3] Vahadane, Abhishek et al. "Structure-preserved color normalization for histopathological images". IEEE, 2015 ## Citing From 6bd3b44fe8c6a66b4af9bd07e6b47f267122f76c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 18:07:01 +0100 Subject: [PATCH 09/14] refactored README --- README.md | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 6041877..187f6a5 100644 --- a/README.md +++ b/README.md @@ -8,9 +8,11 @@ GPU-accelerated stain normalization tools for histopathological images. Compatible with PyTorch, TensorFlow, and Numpy. Normalization algorithms currently implemented: -- Macenko et al. [\[1\]](#reference) (ported from [numpy implementation](https://github.com/schaugf/HEnorm_python)) -- Reinhard et al. [\[2\]](#reference) (only numpy & TensorFlow backend support) -- Vahadane et al. [\[3\]](#reference) (only numpy backend support) +| Algorithm | numpy | torch | tensorflow | +|-|-|-|-| +| Macenko [\[1\]](#reference) | ✓ | ✓ | ✓ | +| Reinhard [\[2\]](#reference)| ✓ | ✗ | ✓ | +| Vahadane [\[3\]](#reference) | ✓ | ✗ | ✗ | ## Installation @@ -45,14 +47,6 @@ norm, H, E = normalizer.normalize(I=t_to_transform, stains=True) ![alt text](data/result.png) -## Implemented algorithms - -| Algorithm | numpy | torch | tensorflow | -|-|-|-|-| -| Macenko | ✓ | ✓ | ✓ | -| Reinhard | ✓ | ✗ | ✓ | -| Vahadane | ✓ | ✗ | ✗ | - ## Backend comparison Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz From 1b11f112968f4e9120db90e1c40eab5a4c25083e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 18:10:29 +0100 Subject: [PATCH 10/14] fixed typo in README [no ci] --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 187f6a5..37030a8 100644 --- a/README.md +++ b/README.md @@ -49,7 +49,7 @@ norm, H, E = normalizer.normalize(I=t_to_transform, stains=True) ## Backend comparison -Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz +Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz using Macenko. | size | numpy avg. time | torch avg. time | tf avg. time | |--------|-------------------|-------------------|------------------| @@ -62,7 +62,7 @@ Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz | 1568 | 1.1935s ± 0.0739 | 0.2590s ± 0.0088 | 0.2531s ± 0.0031 | | 1792 | 1.4523s ± 0.0207 | 0.3402s ± 0.0114 | 0.3080s ± 0.0188 | -## Reference +## References - [1] Macenko, Marc et al. "A method for normalizing histology slides for quantitative analysis." 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. IEEE, 2009. - [2] Reinhard, Erik et al. "Color transfer between images." IEEE Computer Graphics and Applications. IEEE, 2001. From 4e6d7928772dc453778d6da7060f29b07e94ec10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 18:13:38 +0100 Subject: [PATCH 11/14] minor refactor in README [no ci] --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 37030a8..e335cda 100644 --- a/README.md +++ b/README.md @@ -49,7 +49,7 @@ norm, H, E = normalizer.normalize(I=t_to_transform, stains=True) ## Backend comparison -Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz using Macenko. +Runtime results using Macenko with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz. | size | numpy avg. time | torch avg. time | tf avg. time | |--------|-------------------|-------------------|------------------| @@ -65,7 +65,7 @@ Results with 10 runs per size on a Intel(R) Core(TM) i5-8365U CPU @ 1.60GHz usin ## References - [1] Macenko, Marc et al. "A method for normalizing histology slides for quantitative analysis." 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. IEEE, 2009. -- [2] Reinhard, Erik et al. "Color transfer between images." IEEE Computer Graphics and Applications. IEEE, 2001. +- [2] Reinhard, Erik et al. "Color transfer between images." IEEE Computer Graphics and Applications. 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI), IEEE, 2001. - [3] Vahadane, Abhishek et al. "Structure-preserved color normalization for histopathological images". IEEE, 2015 ## Citing From 62d9405f41acf7623a68227ce138035e859186fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Wed, 21 Dec 2022 18:15:43 +0100 Subject: [PATCH 12/14] references typo in redirect --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index e335cda..a3c5f75 100644 --- a/README.md +++ b/README.md @@ -10,9 +10,9 @@ Normalization algorithms currently implemented: | Algorithm | numpy | torch | tensorflow | |-|-|-|-| -| Macenko [\[1\]](#reference) | ✓ | ✓ | ✓ | -| Reinhard [\[2\]](#reference)| ✓ | ✗ | ✓ | -| Vahadane [\[3\]](#reference) | ✓ | ✗ | ✗ | +| Macenko [\[1\]](#references) | ✓ | ✓ | ✓ | +| Reinhard [\[2\]](#references)| ✓ | ✗ | ✓ | +| Vahadane [\[3\]](#references) | ✓ | ✗ | ✗ | ## Installation From 2f018f7f5c9c70f4e651fe6c6bec412434227e80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Thu, 22 Dec 2022 21:15:19 +0100 Subject: [PATCH 13/14] reproduced StainTools vahadane, but dep spams --- torchstain/base/normalizers/__init__.py | 1 + torchstain/base/normalizers/vahadane.py | 2 +- torchstain/numpy/normalizers/__init__.py | 3 ++- torchstain/numpy/normalizers/vahadane.py | 26 +++++++++++++++--------- torchstain/numpy/utils/__init__.py | 3 +++ torchstain/numpy/utils/extract.py | 15 +++++++------- torchstain/numpy/utils/lasso.py | 10 ++++----- torchstain/numpy/utils/od2rgb.py | 4 +++- torchstain/numpy/utils/rgb2od.py | 8 ++++---- torchstain/numpy/utils/stats.py | 4 ++-- 10 files changed, 45 insertions(+), 31 deletions(-) diff --git a/torchstain/base/normalizers/__init__.py b/torchstain/base/normalizers/__init__.py index a8f8bb7..82f59f1 100644 --- a/torchstain/base/normalizers/__init__.py +++ b/torchstain/base/normalizers/__init__.py @@ -1,3 +1,4 @@ from .he_normalizer import HENormalizer from .macenko import MacenkoNormalizer from .reinhard import ReinhardNormalizer +from .vahadane import VahadaneNormalizer diff --git a/torchstain/base/normalizers/vahadane.py b/torchstain/base/normalizers/vahadane.py index b88de7e..4defc66 100644 --- a/torchstain/base/normalizers/vahadane.py +++ b/torchstain/base/normalizers/vahadane.py @@ -1,6 +1,6 @@ def VahadaneNormalizer(backend='numpy'): if backend == 'numpy': - from torchstain.numpy.normalizers import VahadaneReinhardNormalizer + from torchstain.numpy.normalizers import NumpyVahadaneNormalizer return NumpyVahadaneNormalizer() elif backend == "torch": raise NotImplementedError diff --git a/torchstain/numpy/normalizers/__init__.py b/torchstain/numpy/normalizers/__init__.py index d453cf1..9c9e61c 100644 --- a/torchstain/numpy/normalizers/__init__.py +++ b/torchstain/numpy/normalizers/__init__.py @@ -1,2 +1,3 @@ from .macenko import NumpyMacenkoNormalizer -from .reinhard import NumpyReinhardNormalizer \ No newline at end of file +from .reinhard import NumpyReinhardNormalizer +from .vahadane import NumpyVahadaneNormalizer diff --git a/torchstain/numpy/normalizers/vahadane.py b/torchstain/numpy/normalizers/vahadane.py index da46871..0813674 100644 --- a/torchstain/numpy/normalizers/vahadane.py +++ b/torchstain/numpy/normalizers/vahadane.py @@ -1,27 +1,33 @@ import numpy as np from torchstain.base.normalizers import HENormalizer -from torchstain.numpy.utils.lars import lars +from torchstain.numpy.utils.lasso import lasso from torchstain.numpy.utils.stats import standardize_brightness -from torchstain.numpy.utils.extract import get_stain_matrix +from torchstain.numpy.utils.extract import get_stain_matrix, get_concentrations +from torchstain.numpy.utils.od2rgb import od2rgb """ -Source code adapted from: https://github.com/wanghao14/Stain_Normalization/blob/master/stainNorm_Vahadane.py +Source code adapted from: +https://github.com/wanghao14/Stain_Normalization/blob/master/stainNorm_Vahadane.py +https://github.com/Peter554/StainTools/blob/master/staintools/stain_normalizer.py """ class NumpyVahadaneNormalizer(HENormalizer): def __init__(self): super().__init__() self.stain_matrix_target = None + self.maxC_target = None def fit(self, target): - target = standardize_brightness(target) + # target = target.astype("float32") self.stain_matrix_target = get_stain_matrix(target) + concentration_target = get_concentrations(target, self.stain_matrix_target) + self.maxC_target = np.percentile(concentration_target, 99, axis=0).reshape((1, 2)) def normalize(self, I): - I = standardize_brightness(I) + # I = I.astype("float32") + # I = standardize_brightness(I) stain_matrix = get_stain_matrix(I) concentrations = get_concentrations(I, stain_matrix) - out = np.exp(-1 * np.dot(concentrations, self.stain_matrix_target).reshape(I.shape)) - return (255 * out).astype("uint8") - - - + maxC = np.percentile(concentrations, 99, axis=0).reshape((1, 2)) + concentrations *= (self.maxC_target / maxC) + out = 255 * np.exp(-1 * np.dot(concentrations, self.stain_matrix_target)) + return out.reshape(I.shape).astype("uint8") diff --git a/torchstain/numpy/utils/__init__.py b/torchstain/numpy/utils/__init__.py index f440077..56f2aea 100644 --- a/torchstain/numpy/utils/__init__.py +++ b/torchstain/numpy/utils/__init__.py @@ -2,3 +2,6 @@ from torchstain.numpy.utils.lab2rgb import * from torchstain.numpy.utils.split import * from torchstain.numpy.utils.stats import * +from torchstain.numpy.utils.lasso import * +from torchstain.numpy.utils.od2rgb import * +from torchstain.numpy.utils.rgb2od import * diff --git a/torchstain/numpy/utils/extract.py b/torchstain/numpy/utils/extract.py index 8d460d4..f4a9809 100644 --- a/torchstain/numpy/utils/extract.py +++ b/torchstain/numpy/utils/extract.py @@ -2,10 +2,11 @@ from torchstain.numpy.utils.rgb2od import rgb2od from torchstain.numpy.utils.rgb2lab import rgb2lab from torchstain.numpy.utils.lasso import lasso +import spams def extract_tissue(I, th=0.8): - LAB = rgb2lab(I) - L = I_LAB[:, :, 0] / 255.0 + LAB = rgb2lab(I / 255) + L = LAB[:, :, 0] / 255.0 return L < th def get_stain_matrix(I, th=0.8, lamda=0.1): @@ -16,9 +17,9 @@ def get_stain_matrix(I, th=0.8, lamda=0.1): mask = extract_tissue(I, th).reshape((-1,)) OD = OD[mask] - # perform dictionary learning - # dictionary = spams.trainDL(OD.T, K=2, lambda1=lamda, mode=2, modeD=0, posAlpha=True, posD=True, verbose=False) - dictionary = None # @TODO: Implement DL train method + # perform dictionary learning @TODO: Implement DL train method + dictionary = spams.trainDL(OD.T, K=2, lambda1=lamda, mode=2, modeD=0, + posAlpha=True, posD=True, verbose=False) dictionary = dictionary.T if dictionary[0, 0] < dictionary[1, 0]: dictionary = dictionary[[1, 0], :] @@ -31,5 +32,5 @@ def get_concentrations(I, stain_matrix, lamda=0.01): OD = rgb2od(I).reshape((-1, 3)) # perform LASSO regression - #return spams.lasso(OD.T, D=stain_matrix.T, mode=2, lambda1=lamda, pos=True).toarray().T - return lasso(OD.T, y=stain_matrix.T).T # @TODO: Implement LARS-LASSO + return spams.lasso(OD.T, D=stain_matrix.T, mode=2, lambda1=lamda, pos=True).toarray().T + #return lasso(OD.T, y=stain_matrix.T).T # @TODO: Implement LARS-LASSO diff --git a/torchstain/numpy/utils/lasso.py b/torchstain/numpy/utils/lasso.py index 17c090b..ad35c89 100644 --- a/torchstain/numpy/utils/lasso.py +++ b/torchstain/numpy/utils/lasso.py @@ -22,16 +22,16 @@ def z_compute(X): return np.sum(X * X, axis = 0) -def coordinate_descent(y, X, w, alpha, z, tolerance): +def coordinate_descent(y, X, w, alpha, z, tol): max_step = 100 iteration = 0 - while max_step > tolerance: + while max_step > tol: iteration += 1 old_weights = np.copy(w) for j in range(len(w)): rho_j = rho_compute(y, X, w, j) if j == 0: - w[j] = rho_j/z[j] + w[j] = rho_j / z[j] elif rho_j < -alpha * len(y): w[j] = (rho_j + (alpha * len(y))) / z[j] elif rho_j > -alpha * len(y) and rho_j < alpha * len(y): @@ -46,7 +46,7 @@ def coordinate_descent(y, X, w, alpha, z, tolerance): def lasso(x, y, alpha=0.1, tol=0.0001): - w = np.zeros(X.shape[1], dtype="float32") + w = np.zeros(x.shape[1], dtype="float32") z = z_compute(x) - w_opt, iterations, max_step = coordinate_descent(y, x, w, alpha, z, tolerance) + w_opt, iterations, max_step = coordinate_descent(y, x, w, alpha, z, tol) return w_opt diff --git a/torchstain/numpy/utils/od2rgb.py b/torchstain/numpy/utils/od2rgb.py index 58bc46f..904a025 100644 --- a/torchstain/numpy/utils/od2rgb.py +++ b/torchstain/numpy/utils/od2rgb.py @@ -1,4 +1,6 @@ import numpy as np -def OD_to_RGB(OD): +# https://github.com/Peter554/StainTools/blob/2089900d11173ee5ea7de95d34532932afd3181a/staintools/utils/optical_density_conversion.py#L18 +def od2rgb(OD): + OD = np.maximum(OD, 1e-6) return (255 * np.exp(-1 * OD)).astype(np.uint8) diff --git a/torchstain/numpy/utils/rgb2od.py b/torchstain/numpy/utils/rgb2od.py index 9721829..bc34847 100644 --- a/torchstain/numpy/utils/rgb2od.py +++ b/torchstain/numpy/utils/rgb2od.py @@ -1,9 +1,9 @@ import numpy as np - -def RGB_to_OD(I): +# https://github.com/Peter554/StainTools/blob/2089900d11173ee5ea7de95d34532932afd3181a/staintools/utils/optical_density_conversion.py#L4 +def rgb2od(I): # remove zeros - x[x == 0] = 1 + I[I == 0] = 1 # convert to OD and return - return -1 * np.log(I / 255) \ No newline at end of file + return np.maximum(-1 * np.log(I / 255), 1e-6) diff --git a/torchstain/numpy/utils/stats.py b/torchstain/numpy/utils/stats.py index 44b684c..fcf24ed 100644 --- a/torchstain/numpy/utils/stats.py +++ b/torchstain/numpy/utils/stats.py @@ -6,6 +6,6 @@ def get_mean_std(I): def standardize(x, mu, std): return (x - mu) / std -def standardize_brightness(x): - p = np.percentile(x, 90) +def standardize_brightness(x, alpha=99): + p = np.percentile(x, alpha) return np.clip(x * 255 / p, 0, 255).astype("uint8") From e7a6746f6c7969ae77b6e07d9e8ca8a3a931f3bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pedersen?= Date: Thu, 22 Dec 2022 22:02:37 +0100 Subject: [PATCH 14/14] sklearn and spams mismatch --- torchstain/numpy/utils/extract.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/torchstain/numpy/utils/extract.py b/torchstain/numpy/utils/extract.py index f4a9809..794c562 100644 --- a/torchstain/numpy/utils/extract.py +++ b/torchstain/numpy/utils/extract.py @@ -3,13 +3,15 @@ from torchstain.numpy.utils.rgb2lab import rgb2lab from torchstain.numpy.utils.lasso import lasso import spams +from sklearn.linear_model import LassoLars +from sklearn.decomposition import DictionaryLearning def extract_tissue(I, th=0.8): LAB = rgb2lab(I / 255) L = LAB[:, :, 0] / 255.0 return L < th -def get_stain_matrix(I, th=0.8, lamda=0.1): +def get_stain_matrix(I, th=0.8, alpha=0.1): # convert RGB -> OD and flatten channel-wise OD = rgb2od(I).reshape((-1, 3)) @@ -18,8 +20,19 @@ def get_stain_matrix(I, th=0.8, lamda=0.1): OD = OD[mask] # perform dictionary learning @TODO: Implement DL train method - dictionary = spams.trainDL(OD.T, K=2, lambda1=lamda, mode=2, modeD=0, - posAlpha=True, posD=True, verbose=False) + #model = DictionaryLearning( + # fit_algorithm="lars", + # transform_algorithm="lasso_lars", n_components=2, transform_n_nonzero_coefs=0, + # transform_alpha=alpha, positive_dict=True, verbose=False, split_sign=True, + # # positive_code=True, + #) + #dictionary1 = model.fit_transform(OD) + #print(dictionary1) + dictionary = spams.trainDL(OD.T, K=2, lambda1=alpha, mode=2, modeD=0, + posAlpha=True, posD=True, verbose=False) + + #print(dictionary) + #exit() dictionary = dictionary.T if dictionary[0, 0] < dictionary[1, 0]: dictionary = dictionary[[1, 0], :] @@ -27,10 +40,14 @@ def get_stain_matrix(I, th=0.8, lamda=0.1): # normalize rows and return result return dictionary / np.linalg.norm(dictionary, axis=1)[:, None] -def get_concentrations(I, stain_matrix, lamda=0.01): +def get_concentrations(I, stain_matrix, alpha=0.01): # convert RGB -> OD and flatten channel-wise OD = rgb2od(I).reshape((-1, 3)) # perform LASSO regression - return spams.lasso(OD.T, D=stain_matrix.T, mode=2, lambda1=lamda, pos=True).toarray().T + #model = LassoLars(alpha=alpha, positive=True, fit_intercept=False) + #model.fit(X=OD.T, y=stain_matrix.T) + #print(OD.T) + #return model.predict(OD.T).T + return spams.lasso(OD.T, D=stain_matrix.T, mode=2, lambda1=alpha, pos=True).toarray().T #return lasso(OD.T, y=stain_matrix.T).T # @TODO: Implement LARS-LASSO