Background Correction#
Open in Colab
# /// script
# requires-python = ">=3.10"
# dependencies = [
# "matplotlib",
# "ndv[jupyter,vispy]",
# "numpy",
# "scikit-image",
# "scipy",
# "tifffile",
# "imagecodecs",
# ]
# ///
Overview#
In this notebook, we will explore different approaches to background correction in fluorescence microscopy images. Background correction is a crucial pre-processing step that helps remove unwanted background signal and improves the quality of quantitative analysis. We will use the scikit-image library to perform the background correction.
Background subtraction is useful when the background is uniform and the signal to noise ratio is high. We will demonstrate a simple background subtraction method using a sample fluorescence image. The main approaches we’ll cover are:
Subtracting a constant background value (e.g. mode or median of the image)
Selecting and averaging background regions to determine background level
The choice of method depends on your specific imaging conditions and the nature of the background in your images. Here we’ll demonstrate a basic approach that works well for images with relatively uniform background and distinct fluorescent signals.
Note: Background correction should be done on raw images before any other processing steps. The corrected images can then be used for further analysis like segmentation and quantification.
The images we will use for this section can be downloaded from the Measurements and Quantification Dataset.
Importing libraries#
import matplotlib.pyplot as plt
import ndv
import numpy as np
import scipy
import skimage
import tifffile
Background subtraction: mode subtraction#
Background subtraction can be done in different ways. If the background dominates the image as in the example we will use, the most common pixel value (the mode value) can serve as a rough background estimate to subtract from the image.
Let’s first load the image, then display it with ndv to explore the pixel values of the 07_bg_corr_nuclei.tif image.
# raw image and labeled mask
image = tifffile.imread("../../_static/images/quant/07_bg_corr_nuclei.tif")
ndv.imshow(image)
As you can notice, most of the pixels in the image belong to the background (everything but the nuclei). Therefore, we can try to use the mode value of the image as a background estimate. We can use the scipy.stats.mode function to compute the mode of the image.
# Flatten image and get the mode
mode_val = scipy.stats.mode(image.ravel(), keepdims=False).mode
print(f"Estimated background (mode): {mode_val:.3f}")
Estimated background (mode): 33055.000
Then, we can subtract the mode value from the image and print the minimum and maximum pixel values of the resulting image.
Important: Before performing subtraction, convert the image to a floating-point (e.g., image.astype(np.float32). This prevents unsigned integer underflow, which occurs when subtracting the mode value from pixels with intensities lower than the mode. In unsigned integer formats (like uint16), negative results wrap around to very large positive values (e.g., -1 becomes 65535), leading to incorrect results.
# Subtract mode from the image (after converting to float32)
image_mode_sub = image.astype(np.float32) - mode_val
print(f"Min: {image_mode_sub.min():.2f}, Max: {image_mode_sub.max():.2f}")
Min: -71.00, Max: 3808.00
As you can see there are some negative values in the resulting image. This is because some pixels in the original image had values lower than the mode value, and when we subtract the mode from these pixels, we get negative values.
To keep working with the image, we need to handle these negative values. One common approach is to clip the negative values to zero, effectively setting any negative pixel values to zero. This is appropriate since we want background-corrected intensities to be greater than or equal to zero.
For that, we can use the numpy np.clip function to set any negative values to zero and then print the minimum and maximum pixel values of the resulting image.
image_mode_sub_to_zero = np.clip(image_mode_sub, 0, None)
print(
f"Min: {image_mode_sub_to_zero.min():.2f}, Max: {image_mode_sub_to_zero.max():.2f}"
)
Min: 0.00, Max: 3808.00
Finally, we can visualize the background-corrected image (either with matplotlib or ndv):
plt.figure(figsize=(10, 8))
plt.subplot(121)
plt.imshow(image)
plt.title("Original")
plt.axis("off")
plt.subplot(122)
plt.imshow(image_mode_sub_to_zero)
plt.title("Background subtracted (mode)")
plt.axis("off")
plt.tight_layout()
plt.show()
Background subtraction: selected regions#
Sometimes the background isn’t uniform, or the mode isn’t representative. In these cases, we can manually choose one (or more) region we believe contains only background, estimate the average intensity in that region, and then ubtract this average value from the image.
To select regions, it might be helpful to first plot the images with the axes turned on, so we can estimate the pixel coordinates of the regions we want to select. Using matplotlib, we can even visualize and draw on the image the pixel values in the selected regions with the plt.gca().add_patch() function.
# plot the image with axis turned on to select a region
plt.figure(figsize=(8, 8))
plt.imshow(image)
# draw a rectangle around the region we want to select
x = 600
y = 200
width = 100
height = 100
plt.gca().add_patch(
plt.Rectangle(
(x, y), # top-left corner
width, # width
height, # height
edgecolor="yellow",
facecolor="none",
linewidth=2,
)
)
plt.show()
Now we can calculate the mean within the selected region.
# Choose a top-left corner patch assumed to be background
# The region selected above is (600, 200) with size 100x100 therefore we can extract the
# patch as follows: [y:y+height, x:x+width]
bg_patch = image[y : y + height, x : x + width] # image[200:300, 600:700]
bg_mean = np.mean(bg_patch)
print(f"Estimated background (mean of selected region): {bg_mean:.3f}")
Estimated background (mean of selected region): 33055.170
As we did before, we can subtract this value from the image, and clip the result to ensure no negative values remain.
image_mean_sub = image.astype(np.float32) - mode_val
print(f"Min: {image_mean_sub.min():.2f}, Max: {image_mean_sub.max():.2f}")
image_mean_sub_to_zero = np.clip(image_mean_sub, 0, None)
print(
f"Min: {image_mean_sub_to_zero.min():.2f}, Max: {image_mean_sub_to_zero.max():.2f}"
)
Min: -71.00, Max: 3808.00
Min: 0.00, Max: 3808.00
Finally, we can visualize the background-corrected image (either with matplotlib or ndv):
plt.figure(figsize=(8, 8))
plt.subplot(121)
plt.imshow(image)
plt.title("Original")
plt.axis("off")
plt.subplot(122)
plt.imshow(image_mean_sub_to_zero)
plt.title("Background subtracted (regions)")
plt.axis("off")
plt.tight_layout()
plt.show()
Background subtraction: rolling ball algorithm#
Another way of performing background subtraction is the rolling ball algorithm. This is a method that uses a rolling ball to estimate the background. It is a good method to use when the background is not uniform.
The radius parameter configures how distant pixels should be taken into account for determining the background intensity and should be a bit bigger than the size of the structures you want to keep.
Let’s first load another image of a Drosophila embryo that has a non-uniform background:
# raw image
image = tifffile.imread("../../_static/images/quant/07_bg_corr_WF_drosophila.tif")
Let’s explore the pixel values of the image with ndv:
ndv.imshow(image)
We can now estimate the background by using the rolling ball algorithm. Since this function returns an image (which is the background we want to subtract), we can also visualize it with matplotlib:
background_residue = skimage.restoration.rolling_ball(image, radius=100)
plt.figure(figsize=(8, 4))
plt.imshow(background_residue, cmap="gray")
plt.title("Background (rolling ball)")
plt.axis("off")
plt.show()
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[17], line 7
5 plt.title("Background (rolling ball)")
6 plt.axis("off")
----> 7 plt.show()
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/pyplot.py:614, in show(*args, **kwargs)
570 """
571 Display all open figures.
572
(...) 611 explicitly there.
612 """
613 _warn_if_gui_out_of_main_thread()
--> 614 return _get_backend_mod().show(*args, **kwargs)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib_inline/backend_inline.py:90, in show(close, block)
88 try:
89 for figure_manager in Gcf.get_all_fig_managers():
---> 90 display(
91 figure_manager.canvas.figure,
92 metadata=_fetch_figure_metadata(figure_manager.canvas.figure)
93 )
94 finally:
95 show._to_draw = []
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/IPython/core/display_functions.py:278, in display(include, exclude, metadata, transient, display_id, raw, clear, *objs, **kwargs)
276 publish_display_data(data=obj, metadata=metadata, **kwargs)
277 else:
--> 278 format_dict, md_dict = format(obj, include=include, exclude=exclude)
279 if not format_dict:
280 # nothing to display (e.g. _ipython_display_ took over)
281 continue
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/IPython/core/formatters.py:238, in DisplayFormatter.format(self, obj, include, exclude)
236 md = None
237 try:
--> 238 data = formatter(obj)
239 except:
240 # FIXME: log the exception
241 raise
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/decorator.py:235, in decorate.<locals>.fun(*args, **kw)
233 if not kwsyntax:
234 args, kw = fix(args, kw, sig)
--> 235 return caller(func, *(extras + args), **kw)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/IPython/core/formatters.py:282, in catch_format_error(method, self, *args, **kwargs)
280 """show traceback on failed format call"""
281 try:
--> 282 r = method(self, *args, **kwargs)
283 except NotImplementedError:
284 # don't warn on NotImplementedErrors
285 return self._check_return(None, args[0])
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/IPython/core/formatters.py:402, in BaseFormatter.__call__(self, obj)
400 pass
401 else:
--> 402 return printer(obj)
403 # Finally look for special method names
404 method = get_real_method(obj, self.print_method)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/IPython/core/pylabtools.py:170, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
167 from matplotlib.backend_bases import FigureCanvasBase
168 FigureCanvasBase(fig)
--> 170 fig.canvas.print_figure(bytes_io, **kw)
171 data = bytes_io.getvalue()
172 if fmt == 'svg':
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
3254 # ValueError can occur when resizing a window.
3256 self.patch.draw(renderer)
-> 3257 mimage._draw_list_compositing_images(
3258 renderer, self, artists, self.suppressComposite)
3260 renderer.close_group('figure')
3261 finally:
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axes/_base.py:3190, in _AxesBase.draw(self, renderer)
3187 for spine in self.spines.values():
3188 artists.remove(spine)
-> 3190 self._update_title_position(renderer)
3192 if not self.axison:
3193 for _axis in self._axis_map.values():
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axes/_base.py:3134, in _AxesBase._update_title_position(self, renderer)
3132 if title.get_text():
3133 for ax in axs:
-> 3134 ax.yaxis.get_tightbbox(renderer) # update offsetText
3135 if ax.yaxis.offsetText.get_text():
3136 bb = ax.yaxis.offsetText.get_tightbbox(renderer)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axis.py:1351, in Axis.get_tightbbox(self, renderer, for_layout_only)
1349 if renderer is None:
1350 renderer = self.get_figure(root=True)._get_renderer()
-> 1351 ticks_to_draw = self._update_ticks()
1353 self._update_label_position(renderer)
1355 # go back to just this axis's tick labels
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axis.py:1283, in Axis._update_ticks(self)
1281 major_locs = self.get_majorticklocs()
1282 major_labels = self.major.formatter.format_ticks(major_locs)
-> 1283 major_ticks = self.get_major_ticks(len(major_locs))
1284 for tick, loc, label in zip(major_ticks, major_locs, major_labels):
1285 tick.update_position(loc)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axis.py:1664, in Axis.get_major_ticks(self, numticks)
1660 numticks = len(self.get_majorticklocs())
1662 while len(self.majorTicks) < numticks:
1663 # Update the new tick label properties from the old.
-> 1664 tick = self._get_tick(major=True)
1665 self.majorTicks.append(tick)
1666 self._copy_tick_props(self.majorTicks[0], tick)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axis.py:1592, in Axis._get_tick(self, major)
1588 raise NotImplementedError(
1589 f"The Axis subclass {self.__class__.__name__} must define "
1590 "_tick_class or reimplement _get_tick()")
1591 tick_kw = self._major_tick_kw if major else self._minor_tick_kw
-> 1592 return self._tick_class(self.axes, 0, major=major, **tick_kw)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axis.py:429, in YTick.__init__(self, *args, **kwargs)
428 def __init__(self, *args, **kwargs):
--> 429 super().__init__(*args, **kwargs)
430 # x in axes coords, y in data coords
431 ax = self.axes
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/axis.py:176, in Tick.__init__(self, axes, loc, size, width, color, tickdir, pad, labelsize, labelcolor, labelfontfamily, zorder, gridOn, tick1On, tick2On, label1On, label2On, major, labelrotation, grid_color, grid_linestyle, grid_linewidth, grid_alpha, **kwargs)
170 self.gridline.get_path()._interpolation_steps = \
171 GRIDLINE_INTERPOLATION_STEPS
172 self.label1 = mtext.Text(
173 np.nan, np.nan,
174 fontsize=labelsize, color=labelcolor, visible=label1On,
175 fontfamily=labelfontfamily, rotation=self._labelrotation[1])
--> 176 self.label2 = mtext.Text(
177 np.nan, np.nan,
178 fontsize=labelsize, color=labelcolor, visible=label2On,
179 fontfamily=labelfontfamily, rotation=self._labelrotation[1])
181 self._apply_tickdir(tickdir)
183 for artist in [self.tick1line, self.tick2line, self.gridline,
184 self.label1, self.label2]:
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/text.py:139, in Text.__init__(self, x, y, text, color, verticalalignment, horizontalalignment, multialignment, fontproperties, rotation, linespacing, rotation_mode, usetex, wrap, transform_rotates_text, parse_math, antialiased, **kwargs)
137 self._x, self._y = x, y
138 self._text = ''
--> 139 self._reset_visual_defaults(
140 text=text,
141 color=color,
142 fontproperties=fontproperties,
143 usetex=usetex,
144 parse_math=parse_math,
145 wrap=wrap,
146 verticalalignment=verticalalignment,
147 horizontalalignment=horizontalalignment,
148 multialignment=multialignment,
149 rotation=rotation,
150 transform_rotates_text=transform_rotates_text,
151 linespacing=linespacing,
152 rotation_mode=rotation_mode,
153 antialiased=antialiased
154 )
155 self.update(kwargs)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/text.py:174, in Text._reset_visual_defaults(self, text, color, fontproperties, usetex, parse_math, wrap, verticalalignment, horizontalalignment, multialignment, rotation, transform_rotates_text, linespacing, rotation_mode, antialiased)
157 def _reset_visual_defaults(
158 self,
159 text='',
(...) 172 antialiased=None
173 ):
--> 174 self.set_text(text)
175 self.set_color(mpl._val_or_rc(color, "text.color"))
176 self.set_fontproperties(fontproperties)
File ~/work/bobiac-book/bobiac-book/.venv/lib/python3.12/site-packages/matplotlib/text.py:1287, in Text.set_text(self, s)
1275 def set_text(self, s):
1276 r"""
1277 Set the text string *s*.
1278
(...) 1285 ``None`` which is converted to an empty string.
1286 """
-> 1287 s = '' if s is None else str(s)
1288 if s != self._text:
1289 self._text = s
KeyboardInterrupt:
We can now subtract the background residue from the image and plot with matplotlib the raw image, the background image, and the background-corrected image:
image_rb_sub = image - background_residue
plt.figure(figsize=(10, 10))
plt.subplot(131)
plt.imshow(image, cmap="gray")
plt.title("Original")
plt.axis("off")
plt.subplot(132)
plt.imshow(background_residue, cmap="gray")
plt.title("Background (rolling ball)")
plt.axis("off")
plt.subplot(133)
plt.imshow(image_rb_sub, cmap="gray")
plt.title("Background subtracted (rolling ball)")
plt.axis("off")
plt.tight_layout()
plt.show()
Other Background Correction Techniques#
These are more advanced or specialized techniques you can explore:
Morphological opening: Removes small foreground objects to approximate the background.
Gaussian/median filtering: Smooths out the image to isolate large-scale variations.
Polynomial surface fitting: Useful when background varies gradually across the field.
Tiled/local background subtraction: Estimate and subtract background patch-by-patch.
Your method choice should depend on image modality, signal-to-noise, and application.
A good reference for background correction is the scikit-image documentation.