From 4eaec6dc665e85a4111f90c309a464a6f467916f Mon Sep 17 00:00:00 2001 From: David Manthey Date: Fri, 16 Feb 2024 11:22:01 -0500 Subject: [PATCH] Nuclei detection should be more stable across regions This does a better job of detecting the average color of the area being analyzed to reduce individual tile effects. It also fixes a potential infinite loop in a degenerate case. --- .../cli/NucleiDetection/NucleiDetection.py | 28 +- .../color_normalization/__init__.py | 3 +- .../color_normalization/reinhard_stats.py | 9 + .../label/_trace_object_boundaries_cython.pyx | 300 +++++++++--------- .../nuclear/detect_tile_nuclei.py | 4 +- histomicstk/utils/sample_pixels.py | 9 +- 6 files changed, 197 insertions(+), 156 deletions(-) diff --git a/histomicstk/cli/NucleiDetection/NucleiDetection.py b/histomicstk/cli/NucleiDetection/NucleiDetection.py index 906f9bb51..933b3cfef 100644 --- a/histomicstk/cli/NucleiDetection/NucleiDetection.py +++ b/histomicstk/cli/NucleiDetection/NucleiDetection.py @@ -187,7 +187,7 @@ def detect_nuclei_with_dask(ts, tile_fgnd_frac_list, it_kwargs, args, return nuclei_list -def main(args): +def main(args): # noqa # Flags invert_image = False @@ -220,7 +220,7 @@ def main(args): 'tile_size': {'width': args.analysis_tile_size}, 'scale': {'magnification': args.analysis_mag}, 'tile_overlap': {'x': tile_overlap, 'y': tile_overlap}, - 'style': {args.style}, + 'style': args.style, } # retrieve frame @@ -230,7 +230,7 @@ def main(args): msg = 'The given frame value is not an integer' raise Exception(msg) else: - it_kwargs['frame'] = args.frame + it_kwargs['frame'] = int(args.frame) # # Initiate Dask client @@ -305,6 +305,28 @@ def main(args): src_mu_lab, src_sigma_lab = compute_reinhard_norm( args, invert_image=invert_image, default_img_inversion=default_img_inversion) + if src_mu_lab is None: + smallImage, _ = ts.getRegion( + output=dict(maxWidth=4096, maxHeight=4096), resample=False, + region=it_kwargs.get('region', {}), + format=large_image.tilesource.TILE_FORMAT_NUMPY, + frame=args.frame) + if len(smallImage.shape) == 2: + smallImage = np.resize(smallImage, (smallImage.shape[0], smallImage.shape[1], 1)) + if smallImage.shape[2] < 3: + if default_img_inversion or invert_image: + smallImage = (np.iinfo(smallImage.dtype).max + if smallImage.dtype.kind == 'u' + else np.max(smallImage)) - smallImage + smallImage = np.repeat(smallImage[:, :, :1], 3, 2) + elif not default_img_inversion and invert_image: + smallImage = (np.iinfo(smallImage.dtype).max + if smallImage.dtype.kind == 'u' + else np.max(smallImage)) - smallImage + smallImage = smallImage[:, :, :3] + src_mu_lab, src_sigma_lab = htk_cnorm.reinhard_stats_rgb(smallImage) + print(src_mu_lab, src_sigma_lab) + # # Detect nuclei in parallel using Dask # diff --git a/histomicstk/preprocessing/color_normalization/__init__.py b/histomicstk/preprocessing/color_normalization/__init__.py index 6b5301703..7168b0605 100644 --- a/histomicstk/preprocessing/color_normalization/__init__.py +++ b/histomicstk/preprocessing/color_normalization/__init__.py @@ -6,7 +6,7 @@ from .deconvolution_based_normalization import \ deconvolution_based_normalization from .reinhard import reinhard -from .reinhard_stats import reinhard_stats +from .reinhard_stats import reinhard_stats, reinhard_stats_rgb # list out things that are available for public use __all__ = ( @@ -15,5 +15,6 @@ 'background_intensity', 'reinhard', 'reinhard_stats', + 'reinhard_stats_rgb', 'deconvolution_based_normalization', ) diff --git a/histomicstk/preprocessing/color_normalization/reinhard_stats.py b/histomicstk/preprocessing/color_normalization/reinhard_stats.py index fba427aae..1175d26b9 100644 --- a/histomicstk/preprocessing/color_normalization/reinhard_stats.py +++ b/histomicstk/preprocessing/color_normalization/reinhard_stats.py @@ -82,3 +82,12 @@ def reinhard_stats( stats = ReinhardStats(Mu, Sigma) return stats + + +def reinhard_stats_rgb(rgb_image): + rgb_pixels = np.reshape(rgb_image, (1, rgb_image.shape[0] * rgb_image.shape[1], 3)) + + Mu, Sigma = color_conversion.lab_mean_std(rgb_pixels) + + # build named tuple for output + return collections.namedtuple('ReinhardStats', ['Mu', 'Sigma'])(Mu, Sigma) diff --git a/histomicstk/segmentation/label/_trace_object_boundaries_cython.pyx b/histomicstk/segmentation/label/_trace_object_boundaries_cython.pyx index 347df41d4..236a599fe 100644 --- a/histomicstk/segmentation/label/_trace_object_boundaries_cython.pyx +++ b/histomicstk/segmentation/label/_trace_object_boundaries_cython.pyx @@ -275,179 +275,183 @@ def _isbf(long[:, :] mask, long[:, :] mask_90, long[:, :] mask_180, long[:, :] m cdef long i, j, t - with nogil: - # push the first x and y points - list_bx.push_back(x_start); - list_by.push_back(y_start); - - while(True): - - h = np.zeros((row_isbf, col_isbf), dtype=int) + # check degenerate case where mask contains 1 pixel + cdef long sum = np.sum(mask) + if sum > 1: with nogil: - with cython.boundscheck(False): - # initialize a and b which are indices of ISBF - a = 0 - b = 0 - - # get length of the current linked list - size_boundary = list_bx.size() + # push the first x and y points + list_bx.push_back(x_start); + list_by.push_back(y_start); - x = list_bx[size_boundary-1] - y = list_by[size_boundary-1] + while(True): - if (DX == 1) & (DY == 0): - for i in range(ncols-x-2, ncols-x+1): - for j in range(y-1, y+1): - h[a, b] = mask_90[i, j] - b = b + 1 - b = 0 - a = a + 1 - angle = M_PI/2 + h = np.zeros((row_isbf, col_isbf), dtype=int) - elif (DX == 0) & (DY == -1): - for i in range(y-1, y+2): - for j in range(x-1, x+1): - h[a, b] = mask[i, j] - b = b + 1 - b = 0 - a = a + 1 - angle = 0 + with nogil: + with cython.boundscheck(False): + # initialize a and b which are indices of ISBF + a = 0 + b = 0 - elif (DX == -1) & (DY == 0): - for i in range(x-1, x+2): - for j in range(nrows-y-2, nrows-y): - h[a, b] = mask_270[i, j] - b = b + 1 - b = 0 - a = a + 1 - angle = 3*M_PI/2 + # get length of the current linked list + size_boundary = list_bx.size() + + x = list_bx[size_boundary-1] + y = list_by[size_boundary-1] + + if (DX == 1) & (DY == 0): + for i in range(ncols-x-2, ncols-x+1): + for j in range(y-1, y+1): + h[a, b] = mask_90[i, j] + b = b + 1 + b = 0 + a = a + 1 + angle = M_PI/2 + + elif (DX == 0) & (DY == -1): + for i in range(y-1, y+2): + for j in range(x-1, x+1): + h[a, b] = mask[i, j] + b = b + 1 + b = 0 + a = a + 1 + angle = 0 + + elif (DX == -1) & (DY == 0): + for i in range(x-1, x+2): + for j in range(nrows-y-2, nrows-y): + h[a, b] = mask_270[i, j] + b = b + 1 + b = 0 + a = a + 1 + angle = 3*M_PI/2 - else: - for i in range(nrows-y-2, nrows-y+1): - for j in range(ncols-x-2, ncols-x): - h[a, b] = mask_180[i, j] - b = b + 1 - b = 0 - a = a + 1 - angle = M_PI - - cX = vector[int](1) - cY = vector[int](1) - - if h[1, 0] == 1: - # 'left' neighbor - cX[0] = -1 - cY[0] = 0 - DX = -1 - DY = 0 - else: - if (h[2][0] == 1) & (h[2][1] != 1): - # inner-outer corner at left-rear + else: + for i in range(nrows-y-2, nrows-y+1): + for j in range(ncols-x-2, ncols-x): + h[a, b] = mask_180[i, j] + b = b + 1 + b = 0 + a = a + 1 + angle = M_PI + + cX = vector[int](1) + cY = vector[int](1) + + if h[1, 0] == 1: + # 'left' neighbor cX[0] = -1 - cY[0] = 1 - DX = 0 - DY = 1 + cY[0] = 0 + DX = -1 + DY = 0 else: - if h[0, 0] == 1: - if h[0, 1] == 1: - # inner corner at front + if (h[2][0] == 1) & (h[2][1] != 1): + # inner-outer corner at left-rear + cX[0] = -1 + cY[0] = 1 + DX = 0 + DY = 1 + else: + if h[0, 0] == 1: + if h[0, 1] == 1: + # inner corner at front + cX[0] = 0 + cY[0] = -1 + cX.push_back(-1) + cY.push_back(0) + DX = 0 + DY = -1 + else: + # inner-outer corner at front-left + cX[0] = -1 + cY[0] = -1 + DX = 0 + DY = -1 + elif h[0, 1] == 1: + # front neighbor cX[0] = 0 cY[0] = -1 - cX.push_back(-1) - cY.push_back(0) - DX = 0 - DY = -1 + DX = 1 + DY = 0 else: - # inner-outer corner at front-left - cX[0] = -1 - cY[0] = -1 + # outer corner DX = 0 - DY = -1 - elif h[0, 1] == 1: - # front neighbor - cX[0] = 0 - cY[0] = -1 - DX = 1 - DY = 0 - else: - # outer corner - DX = 0 - DY = 1 + DY = 1 - # transform points by incoming directions and add to contours - s = sin(angle) - c = cos(angle) + # transform points by incoming directions and add to contours + s = sin(angle) + c = cos(angle) - if (cX[0]!=0) | (cY[0]!=0): + if (cX[0]!=0) | (cY[0]!=0): - for t in range(cX.size()): + for t in range(cX.size()): - cx = c*cX[t] - s*cY[t] - cy = s*cX[t] + c*cY[t] + cx = c*cX[t] - s*cY[t] + cy = s*cX[t] + c*cY[t] - with gil: - list_bx.push_back(list_bx.back()+int(round(cx))) - list_by.push_back(list_by.back()+int(round(cy))) + with gil: + list_bx.push_back(list_bx.back()+int(round(cx))) + list_by.push_back(list_by.back()+int(round(cy))) - ci = c*DX - s*DY - cj = s*DX + c*DY + ci = c*DX - s*DY + cj = s*DX + c*DY - DX = int(round(ci)) - DY = int(round(cj)) + DX = int(round(ci)) + DY = int(round(cj)) - # get length of the current linked list - size_boundary = list_bx.size() + # get length of the current linked list + size_boundary = list_bx.size() + + if size_boundary > 3: + + fx1 = list_bx[0] + fx2 = list_bx[1] + fy1 = list_by[0] + fy2 = list_by[1] + lx1 = list_bx[size_boundary-1] + ly1 = list_by[size_boundary-1] + lx2 = list_bx[size_boundary-2] + ly2 = list_by[size_boundary-2] + lx3 = list_bx[size_boundary-3] + ly3 = list_by[size_boundary-3] + lx4 = list_bx[size_boundary-4] + ly4 = list_by[size_boundary-4] - if size_boundary > 3: - - fx1 = list_bx[0] - fx2 = list_bx[1] - fy1 = list_by[0] - fy2 = list_by[1] - lx1 = list_bx[size_boundary-1] - ly1 = list_by[size_boundary-1] - lx2 = list_bx[size_boundary-2] - ly2 = list_by[size_boundary-2] - lx3 = list_bx[size_boundary-3] - ly3 = list_by[size_boundary-3] - lx4 = list_bx[size_boundary-4] - ly4 = list_by[size_boundary-4] - - # check if the first and the last x and y are equal - if (size_boundary > max_length) | \ - ((lx1 == fx2)&(lx2 == fx1)&(ly1 == fy2)&(ly2 == fy1)): - # remove the last element - list_bx.pop_back() - list_by.pop_back() - break - if cX.size() == 2: - if (lx2 == fx2)&(lx3 == fx1)&(ly2 == fy2)&(ly3 == fy1): + # check if the first and the last x and y are equal + if (size_boundary > max_length) | \ + ((lx1 == fx2)&(lx2 == fx1)&(ly1 == fy2)&(ly2 == fy1)): + # remove the last element list_bx.pop_back() list_by.pop_back() + break + if cX.size() == 2: + if (lx2 == fx2)&(lx3 == fx1)&(ly2 == fy2)&(ly3 == fy1): + list_bx.pop_back() + list_by.pop_back() + list_bx.pop_back() + list_by.pop_back() + break + # detect cycle + if (lx1 == lx3)&(ly1 == ly3)&(lx2 == lx4)&(ly2 == ly4): list_bx.pop_back() list_by.pop_back() - break - # detect cycle - if (lx1 == lx3)&(ly1 == ly3)&(lx2 == lx4)&(ly2 == ly4): - list_bx.pop_back() - list_by.pop_back() - list_bx.pop_back() - list_by.pop_back() - # change direction from M_PI to 3*M_PI/2 - if (DX == 0) & (DY == 1): - DX = -1 - DY = 0 - # from M_PI/2 to M_PI - elif (DX == 1) & (DY == 0): - DX = 0 - DY = 1 - # from 0 to M_PI/2 - elif (DX == 0) & (DY == -1): - DX = 1 - DY = 0 - else: - DX = 0 - DY = -1 + list_bx.pop_back() + list_by.pop_back() + # change direction from M_PI to 3*M_PI/2 + if (DX == 0) & (DY == 1): + DX = -1 + DY = 0 + # from M_PI/2 to M_PI + elif (DX == 1) & (DY == 0): + DX = 0 + DY = 1 + # from 0 to M_PI/2 + elif (DX == 0) & (DY == -1): + DX = 1 + DY = 0 + else: + DX = 0 + DY = -1 return list_bx, list_by diff --git a/histomicstk/segmentation/nuclear/detect_tile_nuclei.py b/histomicstk/segmentation/nuclear/detect_tile_nuclei.py index acffb4724..3e2606d60 100644 --- a/histomicstk/segmentation/nuclear/detect_tile_nuclei.py +++ b/histomicstk/segmentation/nuclear/detect_tile_nuclei.py @@ -44,7 +44,9 @@ def detect_tile_nuclei(tile_info, args, src_mu_lab=None, # perform image inversion if invert_image: - im_tile = np.max(im_tile) - im_tile + im_tile = (np.iinfo(im_tile.dtype).max + if im_tile.dtype.kind == 'u' + else np.max(im_tile)) - im_tile # perform color normalization im_nmzd = htk_cnorm.reinhard(im_tile, diff --git a/histomicstk/utils/sample_pixels.py b/histomicstk/utils/sample_pixels.py index 31970c6a5..2737cbe0c 100644 --- a/histomicstk/utils/sample_pixels.py +++ b/histomicstk/utils/sample_pixels.py @@ -83,7 +83,9 @@ def sample_pixels(slide_path, sample_fraction=None, magnification=None, # perform image inversion if invert_image: - im_lres = np.max(im_lres) - im_lres + im_lres = (np.iinfo(im_lres.dtype).max + if im_lres.dtype.kind == 'u' + else np.max(im_lres)) - im_lres # compute foreground mask of whole-slide image at low-res. # it will actually be a background mask if background is set. @@ -174,8 +176,9 @@ def _sample_pixels_tile(slide_path, iter_args, positions, sample_fraction, # perform image inversion if invert_image: - im_tile = np.max(im_tile) - im_tile - + im_tile = (np.iinfo(im_tile.dtype).max + if im_tile.dtype.kind == 'u' + else np.max(im_tile)) - im_tile # get tile foreground mask at resolution of current tile tile_fgnd_mask = np.array(PIL.Image.fromarray(tile_fgnd_mask_lres).resize( im_tile.shape[:2],