Microscopy module

#include "diplib/microscopy.h"

struct dip::ColocalizationCoefficients

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

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 Beer-Lambert 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 real-valued. Values outside of the range [0,background] will be clipped. out will be a floating-point 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 real-valued, negative values will be clipped to 0. out will be a floating-point 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 multi-spectral 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 project-specific 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 multi-channel 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 Moore-Penrose pseudo-inverse of , and applies a per-pixel 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 cross-talk, 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 cross-talk ratios) using slides prepared for the purpose, with a single dye.

In multi-spectral 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

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 per-pixel 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 real-valued.

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 multi-channel 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 real-valued. If mask is not forged, all input pixels are considered.

To compute Li’s Intensity Correlation Quotient between two channels in a multi-channel 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 real-valued. 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 multi-channel 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 two-dimensional 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 real-valued. 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 multi-channel 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 P-value) 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 point-spread 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 real-valued. 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 = S::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 compile-time 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 in-focus 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::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 non-masked 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 real-valued 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 real-valued 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 real-valued image. For images with fewer than 3 dimensions, the input is returned unchanged.