Microscopy module #include "diplib/microscopy.h"
Assorted tools useful in microscopy, some presumably also in astronomy and other applications.
Classes

struct dip::
ColocalizationCoefficients  Holds Colocalization Coefficients as described by Manders, see
dip::MandersColocalizationCoefficients
.
Functions

void dip::
BeerLambertMapping (dip::Image const& in, dip::Image& out, dip::Image::Pixel const& background)  Applies a logarithmic mapping to a transmittance image to obtain an absorbance image

void dip::
InverseBeerLambertMapping (dip::Image const& in, dip::Image& out, dip::Image::Pixel const& background = {255})  Applies an exponential mapping to an absorbance image to obtain a transmittance image

void dip::
UnmixStains (dip::Image const& in, dip::Image& out, std::vector<Image::Pixel> const& stains)  Unmixes stains in a brightfield absorbance image or a fluorescence emission image.

void dip::
MixStains (dip::Image const& in, dip::Image& out, std::vector<Image::Pixel> const& stains)  Composes a color image given stain densities and stain absorbance values (brightfield) or stain emission values (fluorescence)

auto dip::
MandersOverlapCoefficient (dip::Image const& channel1, dip::Image const& channel2, dip::Image const& mask = {}) > dip::dfloat  Computes the Manders Overlap Coefficient.

auto dip::
IntensityCorrelationQuotient (dip::Image const& channel1, dip::Image const& channel2, dip::Image const& mask = {}) > dip::dfloat  Computes Li’s Intensity Correlation Quotient.

auto dip::
MandersColocalizationCoefficients (dip::Image const& channel1, dip::Image const& channel2, dip::Image const& mask = {}, dip::dfloat threshold1 = 0, dip::dfloat threshold2 = 0) > dip::ColocalizationCoefficients  Computes Manders’ Colocalization Coefficients.

auto dip::
CostesColocalizationCoefficients (dip::Image const& channel1, dip::Image const& channel2, dip::Image const& mask = {}) > dip::ColocalizationCoefficients  Computes Costes’ colocalization coefficients.

auto dip::
CostesSignificanceTest (dip::Image const& channel1, dip::Image const& channel2, dip::Image const& mask, dip::Random& random, dip::UnsignedArray blockSizes = {3}, dip::uint repetitions = 200) > dip::dfloat  Computes Costes’ test of significance of true colocalization

auto dip::
CostesSignificanceTest (dip::Image const& channel1, dip::Image const& channel2, dip::Image const& mask = {}, dip::UnsignedArray blockSizes = {3}, dip::uint repetitions = 200) > dip::dfloat  Like above, using a defaultinitialized
dip::Random
object. 
void dip::
IncoherentOTF (dip::Image& out, dip::dfloat defocus = 0, dip::dfloat oversampling = 1, dip::dfloat amplitude = 1, dip::String const& method = "Stokseth")  Generates an incoherent OTF (optical transfer function)

void dip::
IncoherentPSF (dip::Image& out, dip::dfloat oversampling = 1, dip::dfloat amplitude = 1)  Generates an incoherent PSF (point spread function)

void dip::
WienerDeconvolution (dip::Image const& in, dip::Image const& psf, dip::Image const& signalPower, dip::Image const& noisePower, dip::Image& out, dip::StringSet const& options = {})  Wiener Deconvolution using estimates of signal and noise power

void dip::
WienerDeconvolution (dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 1e4, dip::StringSet const& options = {})  Wiener Deconvolution using an estimate of noisetosignal ratio

void dip::
ExponentialFitCorrection (dip::Image const& in, dip::Image const& mask, dip::Image& out, dip::dfloat percentile = 1.0, dip::String const& fromWhere = "first plane", dip::dfloat hysteresis = 1.0, dip::String const& weighting = "none")  3D fluorescence attenuation correction using an exponential fit

void dip::
AttenuationCorrection (dip::Image const& in, dip::Image& out, dip::dfloat fAttenuation = 0.01, dip::dfloat bAttenuation = 0.01, dip::dfloat background = 0.0, dip::dfloat threshold = 0.0, dip::dfloat NA = 1.4, dip::dfloat refIndex = 1.518, dip::String const& method = "DET")  3D fluorescence attenuation correction using one of three iterative algorithms

void dip::
SimulatedAttenuation (dip::Image const& in, dip::Image& out, dip::dfloat fAttenuation = 0.01, dip::dfloat bAttenuation = 0.01, dip::dfloat NA = 1.4, dip::dfloat refIndex = 1.518, dip::uint oversample = 1, dip::dfloat rayStep = 1)  3D fluorescence attenuation simulation
Class documentation
struct dip::ColocalizationCoefficients
Holds Colocalization Coefficients as described by Manders, see dip::MandersColocalizationCoefficients
.
Variables  

dip::dfloat M1  Proportional to the fraction of fluorescence in channel 1 that colocalizes. 
dip::dfloat M2  Proportional to the fraction of fluorescence in channel 2 that colocalizes. 
Function documentation
void
dip::BeerLambertMapping (dip::Image const& in,
dip::Image& out,
dip::Image::Pixel const& background)
Applies a logarithmic mapping to a transmittance image to obtain an absorbance image
The BeerLambert law describes how light is attenuated as it travels through an absorbing medium. In brightfield microscopy, this law describes the relationship between the intensity of the transmitted light and the absorbance of the stains on the slide, which provide contrast. The absorbance is influenced by different factors, but for a given dye molecule, the concentration is directly proportional to the absorbance. Thus, estimating absorbance yields an estimate of the relative dye concentration at each image pixel.
This function applies the mapping to the intensities in image in
.
with the intensity of the illumination (background
), and the absorbance written to out
).
background
can be a single value or one value per tensor element (channel) in in
. out
will have the same
number of tensor elements. It should be estimated from a background region in the image, or from a
calibration image taken without a slide in the optical path.
in
must be realvalued. Values outside of the range [0,background
] will be clipped.
out
will be a floatingpoint type (do not force it to be an integer type, as the rounding will destroy
all data).
void
dip::InverseBeerLambertMapping (dip::Image const& in,
dip::Image& out,
dip::Image::Pixel const& background = {255})
Applies an exponential mapping to an absorbance image to obtain a transmittance image
Applies the inverse mapping of dip::BeerLambertMapping
, simulating the image obtained under a brightfield
microscope given the stain densities in the image in
. background
is the illumination intensity, values of
0 in the input will be mapped to the value of background
, whereas larger input values will be mapped to darker
values. Input values should be relatively small, such that background * dip::Exp10(input)
can be represented
in the output data type. Best results are obtained when the input is in the range [0,1], but larger values are
allowed.
in
must be realvalued, negative values will be clipped to 0. out
will be a floatingpoint type, unless
it was protected before calling this function (see dip::Image::Protect
).
void
dip::UnmixStains (dip::Image const& in,
dip::Image& out,
std::vector<Image::Pixel> const& stains)
Unmixes stains in a brightfield absorbance image or a fluorescence emission image.
Brightfield
A color image, obtained from a brightfield microscope, and converted to an absorbance image by
dip::BeerLambertMapping
, can be separated into individual stains as long as there are no more
stains than channels. For an RGB image, up to three stains can be separated. For a multispectral
image, this number is larger.
The stain unmixing process requires knowledge of the absorption spectrum of each of the dyes on the slide. These are usually determined using slides especially prepared with a single dye. Alternatively, find small regions in the image where each stain is on its own (not mixed with other dyes). Below is a table with values for some common dyes, which can be used as a first approximation. However, these absorbance values depend on the tissue, tissue preparation and staining protocols, and imaging equipment. Consequently, best results are always obtained with projectspecific values.
The absorption of the dyes in each channel combine linearly with the density of each of the dyes,
with the absorbance of dye in the red channel, the density (concentration) of dye , and the total absorbance in the red channel. In matrix notation this leads to
Here, is a pixel in the multichannel absorbance image (in
), is a matrix that
combines absorbance for each dye and each channel, and is a vector with the density for each
dye (a pixel in out
). To find , this linear set of equations needs to be solved. This process
is described by Ruifrok (2001). This function computes a MoorePenrose pseudoinverse of ,
and applies a perpixel matrix multiplication with in
to obtain out
.
stains
is a std::vector
that contains each of the columns of matrix . That is, each element
of stains
is the values of one column of , which we refer to as a stain vector. These
stain vectors are represented by a dip::Image::Pixel
with the same number of tensor elements as in
.
stains
cannot have more elements than channels (tensor elements) are in in
.
out
will contain one channel for each stain. For example, assuming an RGB image with 3 channels, stains
can
have one, two or three elements, each element being a dip::Image::Pixel
with exactly 3 elements (corresponding
to the 3 RGB channels).
Best results are obtained when each element of stains
is normalized (i.e. the norm of each stain vector is 1);
this function does not normalize these stain vectors. The standard brightfield stain vectors given below are
normalized.
Example:
dip::Image img = dip::ImageReadTIFF( "brightfield.tif" ); img = dip::BeerLambertMapping( img, { 255 } ); img = dip::UnmixStains( img, {{ 0.644, 0.717, 0.267 }, { 0.268, 0.570, 0.776 }} ); dip::Image nuclei = img[ 0 ]; dip::Image dab = img[ 1 ];
Fluorescence
The explanation above translates to fluorescence imaging, replacing ‘absorbance’ with ‘emission’. In the case
of fluorescence, dip::BeerLambertMapping
should not be used.
Typically, fluorescence imaging systems are set up such that each channel collects light only from a single dye,
but in practice it is not always possible to use dyes with perfectly separated emission spectra. Therefore, there
will be crosstalk, i.e. light from one dye is partially recorded in a channel set up for a different dye.
Again, it is possible to measure the emission intensity in each channel (or channel crosstalk ratios) using slides prepared for the purpose, with a single dye.
In multispectral fluorescence imaging, channels are not set up specifically for each dye. Instead, the spectrum
is divided up into a set of channels. Each dye will be visible in a subset of these channels. Measuring the
emission strength for each dye in each channel again leads to the data to be written in stains
to estimate
dye densities using this function.
Standard brightfield stain vectors
Stain name  RGB absorbance triplet 

AEC  0.274, 0.679, 0.680 
Alcian blue  0.875, 0.458, 0.158 
Aniline blue  0.853, 0.509, 0.113 
Azocarmine  0.071, 0.977, 0.198 
DAB  0.268, 0.570, 0.776 
Eosin  0.093, 0.954, 0.283 
Fast blue  0.749, 0.606, 0.267 
Fast red  0.214, 0.851, 0.478 
Feulgen  0.464, 0.830, 0.308 
Hematoxylin  0.644, 0.717, 0.267 
Hematoxylin + PAS  0.553, 0.754, 0.354 
Methyl blue  0.799, 0.591, 0.105 
Methyl green  0.980, 0.144, 0.133 
Methylene blue  0.553, 0.754, 0.354 
OrangeG  0.107, 0.368, 0.923 
PAS  0.175, 0.972, 0.155 
PonceauFuchsin  0.100, 0.737, 0.668 
void
dip::MixStains (dip::Image const& in,
dip::Image& out,
std::vector<Image::Pixel> const& stains)
Composes a color image given stain densities and stain absorbance values (brightfield) or stain emission values (fluorescence)
This function does the opposite of what dip::UnmixStains
does: it applies the perpixel
matrix multiplication to obtain (out
) from (in
) and
(composed from the values in stains
).
stains
is a vector with these absorbance/emission values, and should have the same number of
elements as channels (tensor elements) in the image in
. Each element of the vector should have the
same number of channels, and these dictate the number of channels in the output image out
. If out
has three channels, it will be tagged as an RGB image. Call dip::InverseBeerLambertMapping
with out
to create an image as seen through a brightfield microscope:
Example:
dip::Image img( { 1024, 1024 }, 2 ); dip::DrawBandlimitedBall( img, 300, { 400, 500 }, { 1.0, 0.0 } ); dip::DrawBandlimitedBall( img, 200, { 500, 600 }, { 0.0, 1.0 } ); img = dip::MixStains( img, {{ 0.644, 0.717, 0.267 }, { 0.268, 0.570, 0.776 }} ); img = dip::InverseBeerLambertMapping( img, { 255 } );
If there are more stains than channels, this process is irreversible (that is, it will not be possible to unmix the stains again).
dip::dfloat
dip::MandersOverlapCoefficient (dip::Image const& channel1,
dip::Image const& channel2,
dip::Image const& mask = {})
Computes the Manders Overlap Coefficient.
The Manders Overlap Coefficient is defined similarly to the Pearson Correlation Coefficient, but without subtracting the means from each of the variables,
with and the two channels. Thus, it returns a value proportional to the fraction of pixels where both channels have a large value. Do note the arguments against this method by Adler and Parmryd (2010).
The images must be scalar and realvalued.
If mask
is not forged, all input pixels are considered. For this measure, it is very important to select
only relevant pixels, and exclude any background pixels with background staining. Furthermore, the zero
level is important, any offset should be subtracted first.
To compute the Manders Overlap Coefficient between two channels in a multichannel image (a tensor image):
MandersOverlapCoefficient( in[ 0 ], in[ 1 ], mask );
dip::dfloat
dip::IntensityCorrelationQuotient (dip::Image const& channel1,
dip::Image const& channel2,
dip::Image const& mask = {})
Computes Li’s Intensity Correlation Quotient.
Li’s Intensity Correlation Quotient is proportional to the fraction of pixels where the two channels vary in a dependent manner. For each pixel, is computed. Then the ICQ is , with the count operator.
The images must be scalar and realvalued. If mask
is not forged, all input pixels are considered.
To compute Li’s Intensity Correlation Quotient between two channels in a multichannel image (a tensor image):
IntensityCorrelationQuotient( in[ 0 ], in[ 1 ], mask );
dip::ColocalizationCoefficients
dip::MandersColocalizationCoefficients (dip::Image const& channel1,
dip::Image const& channel2,
dip::Image const& mask = {},
dip::dfloat threshold1 = 0,
dip::dfloat threshold2 = 0)
Computes Manders’ Colocalization Coefficients.
Manders’ Colocalization Coefficients separate out the contributions in Manders Overlap Coefficient of the two channels and , defining two coefficients as the fraction of staining in one channel that appears where the other channel has some staining,
Note that if the two input images are binary, this is equivalent to computing the precision and sensitivity using
dip::Precision
and dip::Sensitivity
.
Here, instead of thresholding at 0, we apply threshold1
for channel1
, and threshold2
for channel2
. These
thresholds default to 0 to match the method proposed by Manders.
The images must be scalar and realvalued. Any negative values in the input images will cause wrong output,
make sure to clamp the input to 0 before calling this function.
If mask
is not forged, all input pixels are considered.
To compute Manders’ Colocalization Coefficients between two channels in a multichannel image (a tensor image):
MandersColocalizationCoefficients( in[ 0 ], in[ 1 ], mask );
dip::ColocalizationCoefficients
dip::CostesColocalizationCoefficients (dip::Image const& channel1,
dip::Image const& channel2,
dip::Image const& mask = {})
Computes Costes’ colocalization coefficients.
Costes’ Colocalization Coefficients are similar to Manders’ colocalization coefficients, but uses a threshold for each channel under which the correlation is zero. This threshold cuts out the background signal. Staining is colocalized at those pixels where both channels are above their respective threshold. The two coefficients are defined as the fraction of total staining that is colocalized,
with the two input channels. are the thresholds for each channel, with , and and the slope and intercept of the regression line of the twodimensional histogram. The thresholds are successively lowered until the pixels that are not in exhibit no correlation between the two channels.
The images must be scalar and realvalued. Any negative values in the input images will cause wrong output,
make sure to clamp the input to 0 before calling this function.
If mask
is not forged, all input pixels are considered.
To compute Costes’ Colocalization Coefficients between two channels in a multichannel image (a tensor image):
CostesColocalizationCoefficients( in[ 0 ], in[ 1 ], mask );
dip::dfloat
dip::CostesSignificanceTest (dip::Image const& channel1,
dip::Image const& channel2,
dip::Image const& mask,
dip::Random& random,
dip::UnsignedArray blockSizes = {3},
dip::uint repetitions = 200)
Computes Costes’ test of significance of true colocalization
This test verifies whether there is colocalization in the image pair by comparing the correlation between the two channels to that of a randomly shuffled version of the channels. When randomly shuffling one of the channels, there no longer exists correlation. This test gives a significance value to colocalization estimates, but does not quantify the amount of colocalization. Use one of the methods listed below to quantify colocalization.
The algorithm shuffles one of the channels by dividing it into blocks of blockSizes
pixels,
and randomly permuting these blocks. This is repeated repetitions
times. The correlation
between the shuffled channel and the other channel is computed for each of these repetitions,
and a normal distribution is fitted to the obtained values. From this distribution, we compute
the probability that a correlation not larger than the correlation between the two channels is
obtained randomly. This probability (a Pvalue) is returned, and can be compared to, for example,
0.95 to determine with a 5% confidence whether there exists true colocalization in the image
pair.
blockSizes
should be set to the size of the pointspread function or the size of the texture
in the image. An appropriate value can be estimated as the smaller of the widths of the
autocorrelation functions for the two channels. If the block size is too small, the method
will overestimate the significance of the colocalization.
The images must be scalar and realvalued. If mask
is not forged, all input pixels are considered.
If mask
is forged, only blocks that overlap the masked area by at least 3/4 are used.
However, the full block is used, including the portion that falls outside the mask.
void
dip::IncoherentOTF (dip::Image& out,
dip::dfloat defocus = 0,
dip::dfloat oversampling = 1,
dip::dfloat amplitude = 1,
dip::String const& method = "Stokseth")
Generates an incoherent OTF (optical transfer function)
This function implements the formulae for a (defocused) incoherent OTF as described by Castleman.
The defocus
is defined as the maximum defocus path length error divided by the wave length
(See Castleman for details).
When defocus
is unequal to zero, either the Stokseth approximation or the Hopkins approximation
is used, depending on the value of method
(which can be either "Stokseth"
or "Hopkins"
).
The summation over the Bessel functions in the Hopkins formulation is stopped when the change is
smaller than 0.0001 (this is a compiletime constant).
oversampling
is the oversampling rate. If set to 1, the OTF is sampled at the Nyquist rate. Increase
the value to sample more densely.
amplitude
is the value of the OTF at the origin, and thus equivalent to the integral over the PSF.
out
will be scalar and of type dip::DT_SFLOAT
. It should have 1 or 2 dimensions, its sizes will be preserved.
If out
has no sizes, a 256x256 image will be generated.
void
dip::IncoherentPSF (dip::Image& out,
dip::dfloat oversampling = 1,
dip::dfloat amplitude = 1)
Generates an incoherent PSF (point spread function)
This function generates an incoherent infocus point spread function of a diffraction limited objective.
oversampling
is the oversampling rate. If set to 1, the OTF is sampled at the Nyquist rate. Increase
the value to sample more densely.
amplitude
is the integral over the PSF.
out
will be scalar and of type dip::DT_SFLOAT
. It should have 1 or 2 dimensions, its sizes will be preserved.
For 1D images, the PSF returned is a single line through the middle of a 2D PSF.
If out
has no sizes, a square image of size ceil(19*oversampling)
will be generated.
void
dip::WienerDeconvolution (dip::Image const& in,
dip::Image const& psf,
dip::Image const& signalPower,
dip::Image const& noisePower,
dip::Image& out,
dip::StringSet const& options = {})
Wiener Deconvolution using estimates of signal and noise power
If is the Fourier transform of in
, is the Fourier transform of psf
,
and is the Fourier transform of out
, then this function estimates the that optimally
(in the least squares sense) satisfies (that is, in
is the result of the convolution of
out
with psf
).
Finding out
requires knowledge of the power spectrum of the signal and the noise. The Wiener deconvolution
filter is defined in the frequency domain as
where is signalPower
, and is noisePower
. These functions are typically not known, but:

signalPower
can be estimated as the Fourier transform of the autocorrelation ofin
. If a raw image is passed for this argument (dip::Image{}
), then it will be computed as such. 
noisePower
can be estimated as a flat function. A 0D image can be given here, it will be expanded to the size of the other images.noisePower
should not be zero anywhere, as that might lead to division by zero and consequently meaningless results.
The other syntax for dip::WienerDeconvolution
takes an estimate of the noisetosignal
ratio instead of the signal and noise power spectra. Note that can be rewritten as
where is the noisetosignal ratio.
psf
is given in the spatial domain, and will be zeropadded to the size of in
and Fourier transformed.
The PSF (point spread function) should sum to one in order to preserve the mean image intensity.
If the OTF (optical transfer function, the Fourier transform of the PSF) is known, it is possible to pass
that as psf
; add the string "OTF"
to options
.
All input images must be realvalued and scalar, except if the OFT is given instead of the PSF, in which
case psf
could be complexvalued.
void
dip::WienerDeconvolution (dip::Image const& in,
dip::Image const& psf,
dip::Image& out,
dip::dfloat regularization = 1e4,
dip::StringSet const& options = {})
Wiener Deconvolution using an estimate of noisetosignal ratio
See the description of the function with the same name above. The difference here is that a single number,
regularization
, is given instead of the signal and noise power spectra. We then set (the
noisetosignal ratio) to regularization * dip::Maximum(P)
, with P
equal to .
void
dip::ExponentialFitCorrection (dip::Image const& in,
dip::Image const& mask,
dip::Image& out,
dip::dfloat percentile = 1.0,
dip::String const& fromWhere = "first plane",
dip::dfloat hysteresis = 1.0,
dip::String const& weighting = "none")
3D fluorescence attenuation correction using an exponential fit
This routine implements a simple correction of absorption, reflection and bleaching in 3D fluorescence
imaging based upon the assumption that the sum of these effects result in a exponential extinction of
the signal as a function of depth. Only pixels within mask
, if given, are taken into account to determine
this attenuation, but the whole image is corrected.
The attenuation is estimated based on the mean of the nonmasked pixels as a function of depth. If
percentile
is in the valid range [0, 100], the corresponding percentile is used instead of the mean.
An exponential function is fitted to these values. The starting point of the fit is determined by fromWhere
,
which can be "first plane"
, "first max"
, or "global max"
. In the case of "first max"
, the first maximum
is found with point[z+1] > hysteresis * point[z]
.
If the mean variant is chosen, one can choose to apply a variance weighting to the fit by setting weighting
to "variance"
.
in
must be a 3D, scalar and realvalued image. For images with fewer than 3 dimensions, the input is returned
unchanged.
void
dip::AttenuationCorrection (dip::Image const& in,
dip::Image& out,
dip::dfloat fAttenuation = 0.01,
dip::dfloat bAttenuation = 0.01,
dip::dfloat background = 0.0,
dip::dfloat threshold = 0.0,
dip::dfloat NA = 1.4,
dip::dfloat refIndex = 1.518,
dip::String const& method = "DET")
3D fluorescence attenuation correction using one of three iterative algorithms
This function implements an attenuation correction using three different recursive attenuation correction
algorithms. The method is selected with the method
parameter, which must be one of "DET"
, "LT2"
or "LT1"
.
DET stands for Directional Extinction Tracking. LT2 is the Two Light Cone convolutions method, and LT1 is the
One Light Cone convolution.
The DET algorithm is the most accurate one, since it takes both forward and backward attenuation
into account (specified through the fAttenuation
and bAttenuation
parameters). It is however considerably
slower that the LT2 and LT1 algorithms, which take only forward attenuation into account (fAttenuation
).
These last two algorithms assume a constant attenuation (background
) for pixels with an intensity lower
than the threshold
.
The system is characterized by parameters NA
(numerical aperture) and refIndex
(refractive index of the
medium), as well as the pixel size information in in
(the x and y pixel size must be the same, the z size
must be have the same units, but can be different).
in
must be a 3D, scalar and realvalued image. For images with fewer than 3 dimensions, the input is returned
unchanged.
void
dip::SimulatedAttenuation (dip::Image const& in,
dip::Image& out,
dip::dfloat fAttenuation = 0.01,
dip::dfloat bAttenuation = 0.01,
dip::dfloat NA = 1.4,
dip::dfloat refIndex = 1.518,
dip::uint oversample = 1,
dip::dfloat rayStep = 1)
3D fluorescence attenuation simulation
Simulates an attenuation based on the model of a CSLM, using a ray tracing method.
The system is characterized by parameters NA
(numerical aperture) and refIndex
(refractive index of the
medium), as well as the pixel size information in in
(the x and y pixel size must be the same, the z size
must be have the same units, but can be different).
fAttenuation
and bAttenuation
are the forward and backward attenuation factors, respectively.
The ray tracing method uses a step size of rayStep
, and a ray casting oversampling of oversample
.
in
must be a 3D, scalar and realvalued image. For images with fewer than 3 dimensions, the input is returned
unchanged.