# @ String (visibility=MESSAGE, value="The SMO provides a robust estimation of the background dsitribution.
Consult the documentation at https://github.com/maurosilber/SMO", required=false) msg # noqa: E501
# @ImagePlus (label="Input") image
# @String (label="Output", choices={"Background-corrected image", "SMO image", "Background mask"}, style="radioButtonVertical") output_choice # noqa: E501
# @Float (label="Smoothing (sigma)", description="Size of the gaussian smoothing kernel.", style="slider,format:0.0", min=0, max=10, stepSize=0.1, value=0, persist=true) sigma # noqa: E501
# @Integer (label="Averaging window", description="Size of the gradient direction averaging kernel
Must be smaller than the foreground features..", style="slider", min=3, max=21, stepSize=2, value=7, persist=true) size # noqa: E501
# @Float (label="SMO threshold", description="Percentage of background pixels to include in the background distribution estimation.
A higher threshold leads to more foreground pixels being (incorrectly) included.", style="slider,format:0.0", min=0, max=100, stepSize=1, value=5, persist=true) threshold # noqa: E501
# @Float (label="Background percentile", description="Subtract this percentile of the background distribution to the input image.
Note that the resulting image will contain negative values in some background regions,
but the background will be centered at 0.", style="slider,format:0.0", min=0, max=100, stepSize=1, value=50, persist=true) percentile # noqa: E501
# @Boolean (label="Add info to log", value=false) log
# @Boolean (label="Show histogram", description="Show histogram of the background distribution.", value=false) show_histogram # noqa: E501
# @output ImagePlus output
from ij import IJ, ImagePlus
from ij.gui import HistogramWindow
from ij.measure import Measurements
from ij.plugin import ImageCalculator
from ij.process import FloatProcessor
from java.lang import Double
from java.util import Random
def create_mask(image, low_threshold, high_threshold):
"""
Parameters
----------
image: FloatProcessor
low_threshold, high_threshold: float
Returns
-------
BinaryProcessor
"""
image.setThreshold(low_threshold, high_threshold, image.NO_LUT_UPDATE)
return image.createMask()
def apply_mask(image, mask):
"""Applies mask in-place to image.
Parameters
----------
image: ImagePlus
mask: BinaryProcessor
Returns
-------
None
"""
ip = image.getProcessor()
ip.setValue(Double.NaN)
ip.fill(mask)
def smo(image, sigma, size):
"""
image: FloatProcessor
sigma: float
size: int (odd)
"""
if sigma > 0:
image.blurGaussian(sigma)
# Gradient: y
grad_y = image.duplicate()
grad_y.convolve3x3([0, 1, 0, 0, 0, 0, 0, -1, 0])
# Gradient: x
grad_x = image
grad_x.convolve3x3([0, 0, 0, -1, 0, 1, 0, 0, 0])
# Normalize gradient
for x in range(image.getWidth()):
for y in range(image.getHeight()):
xv = grad_x.getf(x, y)
yv = grad_y.getf(x, y)
norm = (xv**2 + yv**2) ** 0.5
if norm == 0:
continue
grad_x.setf(x, y, xv / norm)
grad_y.setf(x, y, yv / norm)
# Average gradient direction
kernel = [1.0 / size**2 for _ in range(size**2)]
grad_x.convolve(kernel, size, size)
grad_y.convolve(kernel, size, size)
# Gradient norm
grad_x.sqr()
grad_y.sqr()
ImageCalculator().run(
"add", ImagePlus("grad_x", grad_x), ImagePlus("grad_y", grad_y)
)
image = grad_x
image.sqrt()
return image
def create_random_image():
"""Creates a image of random uniform noise.
Returns
-------
FloatProcessor
"""
pixels = Random(42).doubles(1024 * 1024).toArray()
return FloatProcessor(1024, 1024, pixels)
def get_quantile(image, quantile):
"""Compute quantile from image.
Parameters
----------
image: ImagePlus
quantile: float
Returns
-------
float
"""
if quantile == 0.5:
stats = image.getStatistics(Measurements.MEDIAN)
return stats.median
else:
stats = image.getStatistics()
fraction = quantile * stats.pixelCount
cdf = 0
for i, x in enumerate(stats.histogram()):
if cdf >= fraction:
break
cdf += x
return i * stats.binSize + stats.histMin
def main(image, sigma, size, threshold, percentile, output_choice, log, show_histogram):
image = image.getProcessor().convertToFloat().duplicate()
# Calculate SMO image
smo_image = image.duplicate()
smo_image = smo(smo_image, sigma, size)
if output_choice == "SMO image":
if log:
IJ.log(
", ".join(
(
"Output: %s" % output_choice,
"Smoothing (sigma): %r" % sigma,
"Averaging window: %r" % size,
)
)
)
smo_image.setMinAndMax(0.0, 1.0)
return ImagePlus(output_choice, smo_image)
# Calculate SMO threshold
random_image = create_random_image()
random_smo_image = smo(random_image, sigma, size)
random_smo_image = ImagePlus("random_smo_image", random_smo_image)
smo_threshold = get_quantile(random_smo_image, threshold / 100)
# Threshold SMO image and apply mask
smo_mask = create_mask(smo_image, smo_threshold, 1.0)
masked_image = ImagePlus("masked_image", image.duplicate())
apply_mask(masked_image, smo_mask)
if output_choice == "Background mask":
if log:
IJ.log(
", ".join(
(
"Output: %s" % output_choice,
"Smoothing (sigma): %r" % sigma,
"Averaging window: %r" % size,
"SMO threshold: %r" % threshold,
)
)
)
return masked_image
# Compute background
background = get_quantile(masked_image, percentile / 100)
if show_histogram:
HistogramWindow(masked_image)
# Correct image
image.add(-background)
if output_choice == "Background-corrected image":
if log:
IJ.log(
", ".join(
(
"Output: %s" % output_choice,
"Smoothing (sigma): %r" % sigma,
"Averaging window: %r" % size,
"SMO threshold: %r" % threshold,
"Background percentile: %r" % percentile,
"Background: %r" % background,
)
)
)
return ImagePlus(output_choice, image)
if __name__ == "__main__":
output = main(
image, # noqa: F821
sigma, # noqa: F821
size, # noqa: F821
threshold, # noqa: F821
percentile, # noqa: F821
output_choice, # noqa: F821
log, # noqa: F821
show_histogram, # noqa: F821
)