-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathsmo.py
209 lines (181 loc) · 7.05 KB
/
smo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
# @ String (visibility=MESSAGE, value="<html>The SMO provides a robust estimation of the background dsitribution.<br>Consult the documentation at <a>https://github.com/maurosilber/SMO</a></html>", 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="<html>Size of the gradient direction averaging kernel<br>Must be smaller than the foreground features.</html>.", style="slider", min=3, max=21, stepSize=2, value=7, persist=true) size # noqa: E501
# @Float (label="SMO threshold", description="<html>Percentage of background pixels to include in the background distribution estimation.<br>A higher threshold leads to more foreground pixels being (incorrectly) included.</html>", style="slider,format:0.0", min=0, max=100, stepSize=1, value=5, persist=true) threshold # noqa: E501
# @Float (label="Background percentile", description="<html>Subtract this percentile of the background distribution to the input image.<br>Note that the resulting image will contain negative values in some background regions,<br>but the background will be centered at 0.</html>", 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
)