-
Notifications
You must be signed in to change notification settings - Fork 17
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Overhaul the tool and add support for single-image processing
- Loading branch information
Showing
9 changed files
with
1,059 additions
and
560 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1 @@ | ||
| ../../macros/creators.xml |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,88 +1,122 @@ | ||
| """ | ||
| Copyright 2021-2022 Biomedical Computer Vision Group, Heidelberg University. | ||
| Author: Qi Gao ([email protected]) | ||
| Authors: | ||
| - Qi Gao ([email protected]) | ||
| - Leonid Kostrykin ([email protected]) | ||
| Distributed under the MIT license. | ||
| See file LICENSE for detail or copy at https://opensource.org/licenses/MIT | ||
| """ | ||
|
|
||
| import argparse | ||
|
|
||
| import imageio | ||
| import giatools.io | ||
| import numpy as np | ||
| import pandas as pd | ||
| from skimage.feature import peak_local_max | ||
| from skimage.filters import gaussian, laplace | ||
|
|
||
|
|
||
| def getbr(xy, img, nb, firstn): | ||
| ndata = xy.shape[0] | ||
| br = np.empty((ndata, 1)) | ||
| for j in range(ndata): | ||
| br[j] = np.NaN | ||
| if not np.isnan(xy[j, 0]): | ||
| timg = img[xy[j, 1] - nb - 1:xy[j, 1] + nb, xy[j, 0] - nb - 1:xy[j, 0] + nb] | ||
| br[j] = np.mean(np.sort(timg, axis=None)[-firstn:]) | ||
| return br | ||
|
|
||
|
|
||
| def spot_detection(fn_in, fn_out, frame_1st=1, frame_end=0, | ||
| typ_filter='Gauss', ssig=1, th=10, | ||
| typ_br='smoothed', bd=10): | ||
| ims_ori = imageio.mimread(fn_in, format='TIFF') | ||
| ims_smd = np.zeros((len(ims_ori), ims_ori[0].shape[0], ims_ori[0].shape[1]), dtype='float64') | ||
| if frame_end == 0 or frame_end > len(ims_ori): | ||
| frame_end = len(ims_ori) | ||
|
|
||
| for i in range(frame_1st - 1, frame_end): | ||
| ims_smd[i, :, :] = gaussian(ims_ori[i].astype('float64'), sigma=ssig) | ||
| ims_smd_max = np.max(ims_smd) | ||
|
|
||
| txyb_all = np.array([]).reshape(0, 4) | ||
| for i in range(frame_1st - 1, frame_end): | ||
| tmp = np.copy(ims_smd[i, :, :]) | ||
| if typ_filter == 'LoG': | ||
| tmp = laplace(tmp) | ||
|
|
||
| tmp[tmp < th * ims_smd_max / 100] = 0 | ||
| coords = peak_local_max(tmp, min_distance=1) | ||
| idx_to_del = np.where((coords[:, 0] <= bd) | (coords[:, 0] >= tmp.shape[0] - bd) | | ||
| (coords[:, 1] <= bd) | (coords[:, 1] >= tmp.shape[1] - bd)) | ||
| coords = np.delete(coords, idx_to_del[0], axis=0) | ||
| xys = coords[:, ::-1] | ||
|
|
||
| if typ_br == 'smoothed': | ||
| intens = getbr(xys, ims_smd[i, :, :], 0, 1) | ||
| elif typ_br == 'robust': | ||
| intens = getbr(xys, ims_ori[i], 1, 4) | ||
| else: | ||
| intens = getbr(xys, ims_ori[i], 0, 1) | ||
|
|
||
| txyb = np.concatenate(((i + 1) * np.ones((xys.shape[0], 1)), xys, intens), axis=1) | ||
| txyb_all = np.concatenate((txyb_all, txyb), axis=0) | ||
|
|
||
| df = pd.DataFrame() | ||
| df['FRAME'] = txyb_all[:, 0].astype(int) | ||
| df['POS_X'] = txyb_all[:, 1].astype(int) | ||
| df['POS_Y'] = txyb_all[:, 2].astype(int) | ||
| df['INTENSITY'] = txyb_all[:, 3] | ||
| import scipy.ndimage as ndi | ||
| from numpy.typing import NDArray | ||
| from skimage.feature import blob_dog, blob_doh, blob_log | ||
|
|
||
| blob_filters = { | ||
| 'dog': blob_dog, | ||
| 'doh': blob_doh, | ||
| 'log': blob_log, | ||
| } | ||
|
|
||
|
|
||
| def mean_intensity(img: NDArray, y: int, x: int, radius: int) -> float: | ||
| assert img.ndim == 2 | ||
| assert radius >= 0 | ||
| if radius == 0: | ||
| return float(img[y, x]) | ||
| else: | ||
| mask = np.ones(img.shape, bool) | ||
| mask[y, x] = False | ||
| mask = (ndi.distance_transform_edt(mask) <= radius) | ||
| return img[mask].mean() | ||
|
|
||
|
|
||
| def spot_detection( | ||
| fn_in: str, | ||
| fn_out: str, | ||
| frame_1st: int, | ||
| frame_end: int, | ||
| filter_type: str, | ||
| min_scale: float, | ||
| max_scale: float, | ||
| abs_threshold: float, | ||
| rel_threshold: float, | ||
| boundary: int, | ||
| ) -> None: | ||
|
|
||
| # Load the single-channel 2-D input image (or stack thereof) | ||
| stack = giatools.io.imread(fn_in) | ||
|
|
||
| # Normalize input image so that it is a stack of images (possibly a stack of a single image) | ||
| assert stack.ndim in (2, 3) | ||
| if stack.ndim == 2: | ||
| stack = stack.reshape(1, *stack.shape) | ||
|
|
||
| # Slice the stack | ||
| assert frame_1st >= 1 | ||
| assert frame_end >= 0 | ||
| stack = stack[frame_1st - 1:] | ||
| if frame_end > 0: | ||
| stack = stack[:-frame_end] | ||
|
|
||
| # Select the blob detection filter | ||
| assert filter_type.lower() in blob_filters.keys() | ||
| blob_filter = blob_filters[filter_type.lower()] | ||
|
|
||
| # Perform blob detection on each image of the stack | ||
| detections = list() | ||
| for img_idx, img in enumerate(stack): | ||
| blobs = blob_filter(img, threshold=abs_threshold, threshold_rel=rel_threshold, min_sigma=min_scale, max_sigma=max_scale) | ||
| for blob in blobs: | ||
| y, x, scale = blob | ||
| radius = scale * np.sqrt(2) * 2 | ||
| intensity = mean_intensity(img, round(y), round(x), round(radius)) | ||
| detections.append( | ||
| { | ||
| 'frame': img_idx + 1, | ||
| 'pos_x': round(x), | ||
| 'pos_y': round(y), | ||
| 'scale': scale, | ||
| 'radius': radius, | ||
| 'intensity': intensity, | ||
| } | ||
| ) | ||
|
|
||
| # Build and save dataframe | ||
| df = pd.DataFrame.from_dict(detections) | ||
| df.to_csv(fn_out, index=False, float_format='%.2f', sep="\t") | ||
|
|
||
|
|
||
| if __name__ == "__main__": | ||
|
|
||
| parser = argparse.ArgumentParser(description="Spot detection") | ||
| parser.add_argument("fn_in", help="Name of input image sequence (stack)") | ||
| parser.add_argument("fn_out", help="Name of output file to save the coordinates and intensities of detected spots") | ||
| parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack)") | ||
| parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack)") | ||
| parser.add_argument("filter", help="Detection filter") | ||
| parser.add_argument("ssig", type=float, help="Sigma of the Gaussian for noise suppression") | ||
| parser.add_argument("thres", type=float, help="Percentage of the global maximal for thresholding candidate spots") | ||
| parser.add_argument("typ_intens", help="smoothed or robust (for measuring the intensities of spots)") | ||
| parser.add_argument("bndy", type=int, help="Number of pixels (Spots close to image boundaries will be ignored)") | ||
|
|
||
| parser.add_argument("fn_in", help="Name of input image or image sequence (stack).") | ||
| parser.add_argument("fn_out", help="Name of output file to write the detections into.") | ||
| parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack).") | ||
| parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack).") | ||
| parser.add_argument("filter_type", help="Detection filter") | ||
| parser.add_argument("min_scale", type=float, help="The minimum scale to consider for multi-scale detection.") | ||
| parser.add_argument("max_scale", type=float, help="The maximum scale to consider for multi-scale detection.") | ||
| parser.add_argument("abs_threshold", type=float, help=( | ||
| "Filter responses below this threshold will be ignored. Only filter responses above this thresholding will be considered as blobs. " | ||
| "This threshold is ignored if the relative threshold (below) corresponds to a higher response.") | ||
| ) | ||
| parser.add_argument("rel_threshold", type=float, help=( | ||
| "Same as the absolute threshold (above), but as a fraction of the overall maximal filter response of an image. " | ||
| "This threshold is ignored if it corresponds to a response below the absolute threshold.") | ||
| ) | ||
| parser.add_argument("boundary", type=int, help="Width of image boundaries (in pixel) where spots will be ignored.") | ||
|
|
||
| args = parser.parse_args() | ||
| spot_detection(args.fn_in, args.fn_out, | ||
| frame_1st=args.frame_1st, frame_end=args.frame_end, | ||
| typ_filter=args.filter, ssig=args.ssig, th=args.thres, | ||
| typ_br=args.typ_intens, bd=args.bndy) | ||
| filter_type=args.filter_type, | ||
| min_scale=args.min_scale, max_scale=args.max_scale, | ||
| abs_threshold=args.abs_threshold, rel_threshold=args.rel_threshold, | ||
| boundary=args.boundary) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
Binary file not shown.
Oops, something went wrong.