METHOD FOR INCREASING CONTRAST

Description

BACKGROUND OF THE INVENTION

[0001] In microscopy, the deconvolution of measured image data is an important and frequent application. Commercial software is available to perform this for common microscopy types, especially confocal microscopy. However, in the special case of macroscopic light-sheet microscopy, these generalized solutions are not applicable because the required point spread function (psf) changes depending on the position within the tissue, i.e., it is not homogeneous. A sufficiently high resolution in the dimension normal to the light sheet is necessary to discretize a sufficiently fine psf or to determine iteratively. For samples on the order of several millimeters, the required amount of data increases dramatically, even though a lower resolution would suffice from a pathological perspective. The additional effort associated with data storage and analysis is a crucial economic factor.

[0002] One theoretically possible approach to deconvolution is the calculation of the inverse psf using so-called linear deconvolution methods, as described, for example, in SARDER, Pinaki; NEHORAI, Arye. Deconvolution methods for 3-D fluorescence microscopy images. IEEE Signal Processing Magazine, 2006, Vol. 23, No. 3, pp. 32-45. Approaches to the numerical calculation of the inverse psf using least squares methods are presented in PREZA, Chrysanthe, et al. Regularized linear method for reconstruction of three-dimensional microscopic objects from optical sections, JOSA A, 1992, Vol. 9, No. 2, pp. 219-228.

[0003] Classical deconvolution involves dividing the convolved image by the peak frequency function (psf) in the frequency domain, taking into account noise reduction considerations such as the application of a low-pass filter. Linear deconvolution results in image correction that improves contrast by "compensating" for the point blurring caused by the psf. In reality, this poses a challenge, particularly with three-dimensional imaging of non-homogeneous tissue, as the psf changes depending on the position within the tissue.

[0004] While a discrete but finite psf (with a finite number of entries in the convolution matrix) can be inverted in some special cases, namely when the psf is a multiplicative function, most finite discrete psfs, and all multidimensional psfs, result in inverse psfs with infinitely many entries whose magnitudes do not approach zero with increasing distance from the measurement point to be deconvolved. Since only a finite number of measurement points are available, this approach does not lead to meaningful results in practice. At microscopic resolutions, deconvolution must therefore be performed by division in the frequency domain. But even then, the problem of the position-dependent psf remains unsolved.

[0005] As mentioned above, light-sheet microscopy requires a high resolution in the dimension normal to the light sheet in order to discretize a sufficiently narrow psf, or to determine iteratively. This would require a comparatively small distance between the individual light-sheet planes, which is associated with the aforementioned problems such as a large data volume and solution instability.

[0006] MOHAMMAD SHORIF UDDIN ET AL: "Restoration of Uneven Illumination in Light Sheet Microscopy Images", MICROSCOPY AND MICROANALYSIS" Vol. 17, No. 4, August 1, 2011 (2011-08-01), pages 607-613, XP001564068, ISSN: 1431-9276, DOI: 10.1017 / S1431927611000262 discloses a method for increasing the contrast of a three-dimensional light sheet microscopy image formed from N individual images, each corresponding to a light sheet plane, spaced apart from each other in the z-direction by at least a distance d, wherein the x / y plane is the light sheet plane and the x-direction is the propagation direction of the light sheet of the light sheet plane. SUMMARY OF THE INVENTION

[0007] The present invention is defined in claim 1. The dependent claims define embodiments of the invention.

[0008] According to one embodiment, a method for increasing the contrast of a three-dimensional light-sheet microscopy image, which is formed from N individual images, each corresponding to a light-sheet plane and spaced apart from each other in the z-direction by at least a distance d, wherein the x / y plane is the light-sheet plane and the x-direction is the propagation direction of the light sheet of the light-sheet plane, comprises: Deconvolution of the three-dimensional image in the z-direction, comprising: For each intensity vector from N intensities (I x,y,1 , ... , I x,y,N ) with the same x / y value, multiplication by a tridiagonal NXN deconvolution matrix, which assigns a correction parameter f1 to each voxel (x, y, n), with which the deconvolution matrix is ​​multiplied by the intensity vector (I x,y,1 , ..., I x,y,N ) for a component I(x, y, n) of the intensity vector, the intensities I(x, y, n+1) and I(x, y, n-1) of the corresponding voxels of the adjacent image planes are multiplied before being subtracted from the intensity value I(x, y,n), where the distance between two light sheet planes is greater than the half-width of the point spread function in the z-direction.

[0009] This form of deconvolution in the z-direction exploits the macroscopic property of light-sheet microscopy that the distance between the individual light-sheet planes is so large that only the immediately adjacent light-sheet needs to be considered during deconvolution in the z-direction, because the more distant light-sheet planes do not make a significant contribution. This reduces the deconvolution matrix for the N voxels in the z-direction for a given x / y value to a tridiagonal NX N-matrix, which assigns only a correction value f1 to each of the adjacent image planes for a given voxel. This greatly simplifies deconvolution compared to conventional approaches.

[0010] According to one embodiment, the correction factor f1 is chosen to be larger as the penetration depth of the light sheet into the tissue increases.

[0011] This takes into account the fact that the light sheet (and thus the psf) widens in the z direction as the light sheet penetrates the tissue in the x direction.

[0012] According to one embodiment, the method further comprises: Determining the correction factor f1 at the point of penetration of the light sheet into the tissue by: Calculating the deconvolution and a quality measure of the deconvolved image based on different correction parameters f1 for a section of the total image, which, starting from the point of penetration of the light sheet, includes the first k voxels in the x-direction, where k is small compared to the extent of the image in the x-direction, and choosing the correction parameter f1 for which the calculated quality measure corresponds to the greatest quality improvement.

[0013] In this way, the correction factor at the point of penetration of the light sheet can be determined experimentally or computationally based on the recorded image.

[0014] According to one embodiment, the method further comprises: representation of the dependence of the correction value f1 on the penetration depth in the x-direction by a monotonically increasing function of the correction value f1(x); calculation, for different functions f1(x), of the respective deconvolution in the z-direction of the overall image as well as a quality measure of the deconvolved overall image,

[0015] Selection of the function f1(x) that provides the best quality measure, as a function to represent the entry depth-dependent correction value f1(x).

[0016] In this way, the penetration depth-dependent correction function and thus also the deconvolution matrix for deconvolution in the z-direction can be determined.

[0017] According to one embodiment, the method further comprises: deconvolution of the image in the x / y plane for each image plane using a two-dimensional point spreading function whose radius depends on the tissue depth of the voxel to be deconvolved in the z-direction.

[0018] In this way, lateral deconvolution can be performed, taking into account the tissue depth dependence of the point spreading function.

[0019] According to one embodiment, the method further comprises: For each image plane, perform two-dimensional deconvolution for different point spread functions that are tissue depth-dependent; calculate a quality measure for the deconvolved image; select a tissue depth-dependent point spread function that leads to an improved contrast value from the majority of point spread functions for which the deconvolution was calculated.

[0020] In this way, a suitable point spreading function for lateral deconvolution can be determined.

[0021] According to one embodiment, a function is used for the tissue depth-dependent point spreading function whose radius increases with increasing tissue depth, but preferably the increase in the radius enlargement decreases continuously.

[0022] This approach for the point spreading function takes into account the result of Monte Carlo simulation regarding the shape and tissue depth dependence of the point spreading function.

[0023] According to one embodiment, a point spread function is chosen for which the mean amplitude of the high spatial frequencies that lie beyond a spatial frequency threshold is a minimum or does not deviate from a minimum by more than a predetermined threshold.

[0024] This takes into account the fact that selecting a spread function solely based on image quality improvement might favor spread functions that lead to artifacts due to high-frequency components. By selecting a spread function where the mean amplitude of the high spatial frequencies is at a minimum or deviates from such a minimum by no more than a threshold value, artifacts can be avoided.

[0025] According to one embodiment, the distance between two light sheet planes is greater than the half-value width of the point spread function in the z-direction.

[0026] This condition ensures that it is a "macroscopic" light-sheet microscopy technique, where correction in the z-direction is possible using a tridiagonal matrix that only considers the immediate neighbors. DESCRIPTION OF THE DRAWINGS

[0027] Figs. 1 to 4The results of a simulation of a deconvolution in the z-direction according to an exemplary implementation for a set of randomly generated one-dimensional intensity vectors are shown. Fig. 5 shows the dependence of the image quality on an exponent of the basis function f_1(x) according to an exemplary implementation. Fig. 6 schematically shows the processing steps of deconvolution in the x / y plane according to an exemplary implementation. Fig. 7 shows the image contrast as a function of the correction exponent of the function r_psf(z) according to an exemplary implementation. Fig. 8 shows the mean amplitude of the high spatial frequencies depending on the correction exponent of the function r_psf(z) for two different image depths. Figs. 9 to 12 the mean intensities and the contrast along the light sheet for the process steps according to an exemplary embodiment on two different tissue samples. Figs. 13 to 19shows light sheet images of a single plane after the processing steps according to an exemplary embodiment on two different tissue samples. DETAILED DESCRIPTION

[0028] The present invention pursues a new approach to the deconvolution of light-sheet microscopy images, which is described below with reference to exemplary embodiments.

[0029] In light sheet microscopy, the illumination direction and the detection direction are orthogonal to each other. A light sheet, generated by a laser or a cylindrical optic, is aligned in the focal plane of the detection objective. By moving the sample, the light sheet "travels" through the sample in the z-direction; an image acquired at a specific z-value then corresponds to a light sheet plane. By acquiring images that each correspond to light sheet planes with different z-values, a stack of images is created, each representing a different cross-sectional plane of the sample, ultimately resulting in a three-dimensional image in which each image plane corresponds to a light sheet plane in the z-direction.

[0030] According to one embodiment, the deconvolution of such an image is first performed in the direction perpendicular to the optical plane (z-direction). The manner in which this deconvolution in the z-direction (hereinafter referred to as the direction of the lens axis perpendicular to the optical plane) is performed has some special features. First, this deconvolution in the direction perpendicular to the optical plane is based on the (justified) assumption that the macroscopic distance between two images, or rather their optical planes, is greater than the half-width of the point spread function. This means that for an nth image in the stack of a total of N optical planes corresponding to images, only the two immediately adjacent optical planes, i.e.,The (n-1)th and (n+1)th images have an influence on the nth image, and therefore only these two need to be considered in the deconvolution. This assumption significantly simplifies the deconvolution matrix for deconvolution in the z-direction.

[0031] For deconvolution, according to one embodiment, for a given voxel I(x,y,n) in the nth image plane, the intensity I(x,y,n+1) of the neighboring voxel in the adjacent image plane is multiplied by a correction factor and subtracted from the intensity I(x,y,n). The same is done with the intensity I(x,y,n-1) of the other neighboring voxel in the z-direction. In this way, the deconvolution of the N voxels in the z-direction for a given location x / y in the x / y-plane is achieved using a tridiagonal matrix that assigns the correction factor f1 to each of the intensities of the neighboring image planes as a multiplication factor.

[0032] In other words, for each intensity vector from N intensities (I x,y,1 , ... , I x,y,N ) with the same x / y value of the image, for the purpose of deconvolution, it is multiplied by a tridiagonal NXN deconvolution matrix, which assigns a correction parameter f1 to each voxel (x, y, n) of the intensity vector (I x,y,1 ,..., I x,y,N ), by which the intensities Ix, y,n+1 and Ix, y,n-1 of the corresponding voxels of the adjacent image plane are multiplied for the component I(x, y, n) of the intensity vector before being subtracted from the intensity value I(x, y,n).

[0033] In one embodiment, a constant correction factor f1 can be used as an approximation. In another embodiment, a value less than 0.5 is arbitrarily chosen as the correction factor. However, the arbitrary choice of a correction factor does not necessarily lead to an optimized result.

[0034] To select a suitable correction factor, it can be determined experimentally or by means of an optimization calculation, as described in another example. For this purpose, the deconvolution can be calculated for various correction factors, and then a contrast value (or another quality measure) for the image after deconvolution can be determined. The correction factor f1 that provides the greatest contrast improvement (or the greatest improvement in the other quality measure) can then be chosen as the correction factor for the tridiagonal deconvolution matrix.

[0035] However, such an approach does not take into account that the light sheet widens in the z-direction the more it penetrates the tissue in the x-direction (i.e., in the direction of light propagation). To account for this effect, the correction factor can be selected according to another embodiment depending on the position in the direction of light propagation (x-direction). The widening of the light sheet in the z-direction with increasing penetration of the light sheet into the tissue in the direction of light propagation (x-direction) is then taken into account, for example, by making the correction factor f1 larger as the penetration depth of the light sheet into the tissue increases.

[0036] According to one embodiment, the correction factor f1 is chosen depending on the penetration depth x, whereby f1(x) preferably increases monotonically with increasing x.

[0037] Experimentally, f1(x) can be determined for such an approach according to one embodiment as follows. First, for an image section near the point of penetration of the light sheet, i.e., for example, for an image section with small x-values, the deconvolution of the image section is calculated for various values ​​of f1. It is justifiably assumed that for small x-values, i.e., for shallow tissue depths in the x-direction, the light sheet has not yet expanded significantly, and thus the correction value f1 for these x-values ​​can be approximated by a constant factor that corresponds to the correction value f1 at the point of penetration of the light sheet into the tissue. Preferably, the number of voxels in the x-direction used for the image section is small compared to the total number of voxels of the entire image in the x-direction.For example, less than the first 20% of the voxels in the x-direction are used for the image section, alternatively less than 10%, further alternatively less than 5%, further alternatively less than 1%. Alternatively still, for example, the first two, the first three, the first four, or the first 5-10 voxels in the x-direction can be chosen as the image section to be used. In all these examples, the number of voxels in the x-direction of the selected image section for optimization is small compared to the total number of voxels in the x-direction; that is, an image section is chosen that corresponds to a small penetration depth of the light sheet in the x-direction.

[0038] For the selected image area, the deconvolved image is then calculated for various parameters f1, and subsequently a quality measure of the deconvolved image is calculated, e.g., a contrast value. The parameter f1 that yields the "best" quality measure for the deconvolved image is then chosen as the correction parameter f1(0) for the entry point of the light sheet.

[0039] The widening of the psf can then be represented by a monotonically increasing function of the correction value f1(x). For different functions, the respective deconvolution in the z-direction and a quality measure of the deconvolved image are calculated for the overall image. The function f1(x) that yields the best quality measure (e.g., the greatest contrast improvement) is then chosen as the function to represent the depth-of-entry psf or to form the corresponding deconvolution matrix with its f1 entries.

[0040] Using the single-penetration-depth-dependent deconvolution matrix thus determined, deconvolution in the z-direction can be carried out with little effort and a significant quality-improving effect.

[0041] In a further step, which may also be carried out independently of the correction of the light sheet thickness in the z-direction, i.e. the deconvolution in the z-direction, the deconvolution then takes place in the x / y-plane.

[0042] According to one embodiment, a two-dimensional psf is used whose radius depends on the tissue depth in the z-direction, i.e., in the direction of the lens. For different approach functions for the point spread function, the two-dimensional deconvolution is performed for each light sheet as a function of tissue depth, and a quality measure for the deconvolved image is calculated, e.g., a contrast value. According to one embodiment, a tissue-depth-dependent psf is then selected that leads to an improved contrast value.

[0043] According to another embodiment, a function is used for the tissue depth-dependent psf whose radius increases with increasing tissue depth, but preferably the increase in the radius enlargement decreases steadily.

[0044] According to one embodiment, a point spread function is chosen for which the mean amplitude of the high spatial frequencies that lie beyond a spatial frequency threshold is a minimum or does not deviate from a minimum by more than a predetermined threshold.

[0045] This takes into account the fact that selecting a spread function solely based on image quality improvement might favor spread functions that lead to artifacts due to high-frequency components. Artifacts can be avoided by selecting a spread function where the mean amplitude of the high spatial frequencies is at a minimum or deviates from such a minimum by no more than a threshold value. This threshold value could, for example, be a certain percentage of the minimum value, such as 10%, 15%, or 20%. Other suitable threshold values ​​are also conceivable for a person skilled in the art, such as 30%, 40%, or 50% of the minimum. The value should be chosen so that it can still be meaningfully described as "close to the minimum." This can also be the case, for example, if the value is increased by a certain percentage.B: 20%, 30% or 40% difference between minimum and maximum deviates from the minimum value.

[0046] If a point spread function has been selected taking into account the contrast improvement and the mean amplitude of the high spatial frequencies, then two-dimensional deconvolution can be carried out using a standard procedure (such as the Richardson-Lucy algorithm).

[0047] Preferably, groups of pixels with similar tissue depths are formed from the set of local tissue depths in an image. The number of these groups is arbitrary. For each pixel group, the mean local tissue depth, and thus the local tissue depth-dependent point spread function (psf), is calculated. Deconvolution is then performed for each group on the entire image according to the local mean psf. Only then are the deconvolved images masked according to the pixel groups. When combined, they yield the fully deconvolved image.

[0048] The deconvolution of the entire image then results from the deconvolution in the z-direction (the correction of the "light sheet thickness" using a tridiagonal deconvolution matrix) as well as the two-dimensional deconvolution in the x / y-plane based on a two-dimensional psf.

[0049] A new aspect of the previously described methodology according to an exemplary implementation is the consideration of the tissue depth in two spatial directions (light sheet propagation and detection axis); the separation of the spatial directions during correction, and thus the simplification of the deconvolution to the two-dimensional case, as well as the combination of empirical and theoretical approaches for optimizing the algorithm or for determining the point spreading functions.

[0050] The correction for the light sheet thickness utilizes the macroscopic resolution of the image stack in the z-direction to calculate an exact inverse convolution matrix. This allows for improved contrast even in this resolution range, where classical algorithms would not permit correction. The algorithm is stable and fast. Additional image planes that are irrelevant to the user do not need to be acquired, resulting in time savings even during the actual measurement on the microscope. The low computational effort allows the correction to be adapted to spatial changes in the light sheet thickness resulting from scattering and deflection of the light sheet within the tissue. If the factors for spatial adaptation are determined indirectly from the image data, recursive conclusions about this spatial structure of the light sheet can be drawn.

[0051] In the application of the light-sheet microscope, the three-dimensional convolution can be decomposed, as described, into a two-dimensional convolution at the plane of the light sheet, caused by light refraction during the transmission of the fluorescent emission or scattering of the excitation by the tissue, and a one-dimensional convolution caused by the thickness of the light sheet. If the light-sheet convolution caused by the thickness of the light sheet towards the objective can be reverse-engineered, then, in principle, established deconvolution methods for the two-dimensional deconvolution at the plane of the light sheet can be applied.

[0052] According to one embodiment, the three-dimensional image is transferred into a new coordinate system by converting the distance of each voxel from the tissue surface into the detection direction. The transformed image is then deconvolved using a tissue depth-dependent psf.

[0053] Light sheet correction, i.e., correction of the light sheet thickness or deconvolution in the z-direction, precedes tissue depth-dependent lateral deconvolution to generate a uniformly illuminated initial image with distinct tissue boundaries. Since deconvolution requires knowledge of the tissue depth, both process steps are directly coupled. Empirical methods based on determining the information content of the deconvolved image (e.g., contrast measurement) are used to determine the parameters of the point spreading functions employed.

[0054] The following section describes in more detail a proposed solution for the algebraic correction of the light sheet thickness, i.e., the deconvolution in the z-direction, according to an exemplary embodiment.

[0055] The discretization of the convolution vector in the dimension normal to the plane of light sheet propagation forms a linear system of equations, which can be solved for any number of vertical support points. If the system is simplified by considering only the first support points above and below the center of the light sheet, the system of equations can be interpreted as a tridiagonal matrix NxN (where N is the number of recorded planes), for which computationally fast solution algorithms, such as the Thomas algorithm, are known.

[0056] Provided that the tridiagonal convolution matrix is ​​positive definite, meaning all its eigenvalues ​​are real and positive, its inverse is also a tridiagonal matrix. Then, and only then, is an algebraic correction permissible. The result is a constant factor for every combination of light sheet and measurement distance, by which the respective neighboring planes must be multiplied before being subtracted from the image plane to be corrected. This factor can be calculated once by matrix inversion, thus reducing the correction, under constant measurement conditions, to a simple, direct algebraic equation per pixel. The greater the distance between the light sheets, the smaller the magnitude of the multiplication factor, and therefore the smaller the gain in image sharpness from applying the correction.The condition of positive definiteness can be translated into a minimum distance between the measurement points, i.e., a minimum distance between the measured image planes. At smaller distances, the matrix inversion becomes unstable, and the aforementioned case of inverse psf with infinitely many entries occurs.

[0057] In this context, it's helpful to consider that at very small measurement distances, adjacent images become too similar to be corrected simply by subtracting their neighbors. The magnitude of the multiplication factor becomes so large that negative entries appear in the resulting image, which then have to be corrected by adding the next neighboring layer. This process continues across the entire image stack.

[0058] This consideration shows that the property of macroscopy—that is, the comparatively large (macroscopic) distance between the light sheet planes, which prevents the application of established deconvolution methods—is simultaneously the property that enables the computation of a finite inverse psf. Furthermore, the method can be adapted for each pixel according to the light sheet thickness, thus eliminating the problem of psf inhomogeneity. In particular, the available computational effort can be redistributed to a local adaptation of the inverted matrix because the deconvolution of the individual pixel is reduced to a simple multiplication and subtraction.

[0059] Implementing this method requires determining a correction factor f1, which will be discussed in more detail in the following section using examples. Once this factor is determined, various properties of the corrected image in the frequency domain can be examined.

[0060] The effect of the approach of correcting by subtracting only the nearest neighbor using a correction factor will now be explained in detail using a simulation. This serves to demonstrate the effect of subtracting the immediate neighbor using a correction factor for randomly generated images. For this purpose, random "one-dimensional" images, i.e., random intensity values, were generated to create a random "one-dimensional image" with 500 pixels. Based on this, the correction was then applied, in which the neighboring intensity value multiplied by the correction factor f1 is subtracted, as is the case with the tridiagonal deconvolution matrix.

[0061] For the graphs in Fig. 1 f1 was assumed to be 0.26 and the algebraic correction was applied to randomly generated one-dimensional datasets (500 "pixels" each).

[0062] This simulation was repeated one hundred thousand times. Fig. 1 This figure shows the averaged amplitudes in the frequency domain (i.e., after Fourier transformation) of the random data set before ("original") and after correction by subtracting the neighbors multiplied by the correction factor. The ratio of the original to the corrected amplitude profile is also shown. In a classical deconvolution, this would correspond to the potential difference (psf) in the frequency domain. To provide an indication of the variance within the simulations, the 99th percentiles of this ratio are also shown.

[0063] The change in the mean quotient during the simulation is represented as the residual in Fig. 2 The graph shows the residual of the simulation with standard deviation as the ratio of the change in the mean amplitude of the quotient function (when adding a simulation) to the newly calculated mean amplitude across all simulations.

[0064] The original 100,000 data points were locally averaged (along a logarithmic scale) and thus reduced to 50 data points, plus one additional data point for the first ten simulations. The running standard deviation shown was calculated in the same way. For the residual itself, the course of the averaged amplitudes of the quotient before and after adding a new simulated data point was compared. After 100,000 simulations, the quotient function changes on average by less than one ten-sixth of its own mean value.

[0065] For Fig. 1 and Fig. 2Only the magnitudes of the frequency domain were used. If the complex components, including the phase, are obtained when calculating the quotient, this new complex frequency domain function can be inverted to a real function in the intensity domain. Since both the original function and the corrected function are real in the intensity domain, their corresponding values ​​in the frequency domain have complex conjugate entries. Their quotient thus remains complex conjugate, and the inverse (in the intensity domain) becomes real again (apart from numerical rounding errors). This function is in Fig. 3 The figure shows the Fourier inverse of the simulated mean quotient function. Fig. 1 (f1=0.26) and for further correction parameters in the frequency domain.

[0066] A spacing of 10 µm between two images in the image stack was assumed, and this corresponds to the resolution in the representation. The distribution shown can be assigned a full width at half maximum (FWHM) of 13.8 µm through linear interpolation between the data points. If a different correction value is used instead of 0.26, a different FWHM results, as also shown.

[0067] The respective half-widths were calculated for various factors f1; these are plotted in Fig. 4 As expected, the full width at half maximum (FWHM) increases with increasing correction factor f1 until a local maximum of 18.2 µm FWHM is reached at f1 = 0.43. After this point, the function no longer changes around 0 µm; instead, the offset value increases in the periphery (see...). Fig. 3). Approximately at f1 = 0.49, the half-widths increase again, and further local maxima form in the quotient function. The mathematical limit of stable invertibility has been reached.

[0068] However, the simulation shows that the approach of subtracting the nearest neighbor using a correction factor is fundamentally suitable for stable deconvolution, as long as the distance between the planes is sufficiently large compared to the half-width of the point spread function, for example, at least as large as the half-width of the point spread function. This is the case in macroscopic light-sheet microscopy, which is why this method can be applied to deconvolution in the z-direction, i.e., light-sheet thickness correction.

[0069] However, the selection or determination of the correction factor f1 is required for the application. This empirical determination of the local (tissue depth-dependent) correction factor is described in more detail below using an example.

[0070] The thickness of the light sheet before entering and after exiting the tissue can be determined directly through experimentation. This experimental determination can be used to restrict the approximate function f_1 both above and below, thereby improving the empirical determination of the function. However, the local progression of the light sheet's expansion within the tissue is not trivial and cannot be readily determined. The same applies to any potential deflection of the light sheet from the original exposure plane.

[0071] Based on an existing macroscopic light sheet measurement, the local expansion of the light sheet is empirically determined according to one embodiment. The initial assumption is that expansion occurs only within the tissue. A global thresholding method (e.g., Otsu) is then used to distinguish between tissue and background. This makes it possible to assign a local tissue depth to each pixel within the tissue in the direction of light sheet propagation.

[0072] To model the expansion of the light sheet in the tissue, a monotonically increasing function is used according to one embodiment, which provides positive values ​​for all tissue depths.

[0073] Furthermore, it is advantageous to use no more than one or two empirically determined parameters in order to limit the computational effort. For example, a power function of the following form can be used: f 1 x T = 1 2 − 1 x T n

[0074] Here, f1 is the required correction factor, i.e., the value by which the adjacent layers are multiplied before being subtracted from the layer to be corrected. XT is the tissue depth in pixels, and n is the previously unknown power. The described function approaches the theoretical maximum correction value of 0.5 as the tissue depth XT increases (for n > 0). To avoid negative f1, the function is set to the correction value of the undisturbed light sheet in tissue, f1.0, for small XT values. This can be determined empirically by performing the correction with different values ​​for f1.0 at the tissue edge facing the light sheet. The value f1.0 that yields the highest image contrast (or subjectively best image impression) is selected as the parameter. The lower limit of the tissue depth XT0,0 for a known undisturbed light sheet thickness f1.0 (0.26 in this example) is calculated from the equation above as follows: x T , 0 = 1 1 2 − f 1,0 n

[0075] After each correction, the contrast in the foreground is calculated as the standard deviation of the corrected intensities of the foreground, relative to the mean intensity (see PELI, Eli. Contrast in complex images. JOSA A, 1990, Vol. 7, No. 10, pp. 2032-2040). For demonstration purposes, the method was applied to a photogrammetric measurement of a human appendix, which had previously been optically cleared using a method as described in DE 102016123458 B3.

[0076] This results in the in Fig. 5 The contrast value is shown as a function of the exponent of the correction function. The optimal overall contrast is achieved here at an exponent of -0.65.

[0077] The described method can, in principle, be varied and extended as desired by using a different criterion for image quality (instead of contrast, or with a different definition of image contrast) or a different function for the tissue depth-dependent correction factor. Furthermore, the number of free variables can be increased. By assuming different correction factors above and below the image, a deflection of the light sheet away from the exposure plane can be modeled. According to one embodiment, a numerical iterative method for solving nonlinear minimization problems, such as the Gauss-Newton method, can then be used, albeit with the associated increased computational effort.

[0078] Once the correction function has been determined, the same, once determined correction function can always be used for multiple light sheet measurements with the same tissue types and clearing protocols.

[0079] Following the description of the deconvolution in the z-direction, the tissue depth-dependent deconvolution perpendicular to the lens axis is described in more detail below according to an exemplary embodiment.

[0080] After the light sheet thickness correction has been performed for all planes in the image stack, the "actual deconvolution" takes place, i.e., deconvolution in the lateral direction (in the x / y plane). Due to the separation between axial and lateral psf (light sheet thickness and scattering during transmission through the tissue), a two-dimensional psf remains after the light sheet thickness correction.

[0081] The extent of this psf depends on the light path length between the emission source and the tissue boundary in the direction of the objective lens; in the case of a light-sheet microscope, this is perpendicular to the illumination plane. Experimental determination of this psf using fluorescent beads allows the determination of its shape, but not its size in actual tissue.

[0082] Similar to the case of the penetration-depth-dependent correction factor for the light sheet correction, one embodiment determines the extent of the psf as a function of tissue depth. Analogous to the procedure described above, each voxel within the tissue is assigned a tissue depth, now in the axial (z) direction instead of in the direction of the light sheet's propagation. According to one embodiment, a power function again serves as the basis function for the psf radius rpsf.

[0083] Monte Carlo simulations of scattering within tissues showed that light beams diverge monotonically with increasing tissue depth, but the rate of this diverging decreases steadily (see FLOCK, Stephen T., et al. Monte Carlo modeling of light propagation in highly scattering tissues. I. Model predictions and comparison with diffusion theory. IEEE Transactions on Biomedical Engineering, 1989, Vol. 36, No. 12, pp. 1162-1168). A power function with an exponent 0 < n < 1 satisfies this condition. r psf = k z T Δ z n

[0084] In the equation above, ZT denotes the local tissue depth and Δz the image resolution in the axial direction. Knowing the exponent n, this equation allows the calculation of the psf expansion. The constant normalization factor k serves to scale the psf radius and can initially be chosen arbitrarily. In the example shown here, k was set to 0.5. It can be corrected based on the empirical determination of the exponent n, as described below.

[0085] To implement this local deconvolution process, the image is divided into three to five regions of similar tissue depth at each layer. The local psf is then calculated as an Airy function with a depth-dependent radius for each region. The radius rpsf corresponds to the position of the first zero of the Airy function. Deconvolution is then performed for each psf for the entire image using the Richardson-Lucy algorithm. With five regions, this yields five deconvolved images. These are each masked (using the defined regions of similar depth as masks) and combined to form the deconvolved overall image. This process is performed for each image in the image stack. The entire described procedure is illustrated in [reference to be added]. Fig. 6 depicted.

[0086] The value of the correction exponent n depends on the optical properties of the tissue under investigation at the given excitation or emission wavelength. How Fig. 7 This shows that the image contrast, in terms of the standard deviation of the intensities in the tissue, increases with n. However, an excessively "large" psf leads to ring artifacts in the image, visible as an increase in the amplitudes of high spatial frequencies (cf. Fig. 8 ).

[0087] An exponent n = 0.37 shows no ring artifacts in the two image planes considered for cleared lymphatic tissue and is therefore used for deconvolution in the present embodiment.

[0088] If a deconvolved image without ring artifacts cannot be generated for any exponent n, the normalization factor k is preferably adjusted downwards. If the frequency contrast shows a broad plateau of low amplitudes over large ranges of n, k can be increased. In both cases, the empirical determination of n is then repeated.

[0089] According to one exemplary embodiment, the complete algorithm for correcting and increasing the contrast of light sheet measurements consists of six individual steps: determination of a global threshold intensity, low-pass filtering, local light sheet correction, tissue depth-dependent deconvolution, and adjustment of the intensity distribution along the light sheet.

[0090] In the Figs. 9 to 12The mean intensity and contrast (as standard deviation of intensity) along the light sheet for the most important process steps are shown using two exemplary tissues. Each data point represents the set of all points within the tissue that are equidistant from the tissue edge in the direction of light sheet propagation. The intensity profile of the appendix in Fig. 9 The image shows a tissue indentation between 4mm and 6mm deep after the application of light-based correction. The cause is the intestinal lumen of an appendix located in the middle of the tissue.

[0091] The calculation of the global threshold intensity is performed in two steps, according to one example implementation. First, a mean-value z-projection is calculated along the lens axis. Then, a global Otsu threshold method is applied to this resulting image. The calculated threshold value is then corrected downwards by a safety factor of 0.6.

[0092] The low-pass filter is an ideal FFT filter with a circular mask that sets the 20% highest spatial frequencies in the image to zero. This is followed by light sheet correction and tissue depth-dependent deconvolution, as described above.

[0093] Due to absorption and scattering, the excitation intensity decreases with increasing distance from the light sheet origin. Therefore, after light sheet correction, the intensities are adjusted along the light sheet. For each image in the stack, the average intensity is calculated as a function of tissue depth. The assignment of voxels to tissue depth is analogous to the description above for light sheet correction and deconvolution, in this case in the direction of light sheet propagation. The individual voxels in the tissue are then divided by the average intensity in the tissue just determined. Subsequently, the entire individual image is normalized back to the order of magnitude of the desired image depth (e.g., 8-bit or 16-bit).

[0094] For the Figs. 13 to 17Light sheet images of a single plane were presented as examples after the respective processing steps, i.e., original image, Fourier filtering, light sheet thickness correction, lateral deconvolution, and intensity adjustment. The curves in the tissue shown are also derived from this tissue. Fig. 9 and 10 The effects of each processing step become clearly evident. Light sheet correction adds new information to the individual image, making new structures visible. Tissue depth-dependent deconvolution renders these new structures more sharply. Adjusting the intensity distribution emphasizes image areas located deep within the tissue.

[0095] As another example, compare Fig. 18 and Fig. 19 A light-sheet microscopy image of a human thymus before and after the entire described procedure. The curves for this tissue image are also derived from this image. Fig. 11 and 12The 16-bit grayscale images were converted into contour diagrams, each with ten iso-lines of equal intensity.

[0096] The person skilled in the art recognizes that the described methods for contrast enhancement can be implemented using a computer program running on a commercially available computer. Based on the preceding description, it is readily possible for the person skilled in the art to design a corresponding computer program. For this purpose, the image data, i.e., the intensities of the voxels of the light-sheet microscope image, are subjected to the described mathematical operations, which are implemented using a suitable computer program.

[0097] The image processed by the described methods is a technical entity; the contrast improvement of an image, in particular a microscope image, is a technical effect or the solution of a technical problem that is solved in a new, non-obvious way by the methods described in the exemplary embodiments and defined in the claims.

Claims

1. A method of contrast enhancement of a three-dimensional light sheet microscopy image formed from N individual images each corresponding to a light sheet plane and spaced apart from each other in the z-direction by at least a distance d, the x / y-plane being the light sheet plane and the z-direction being the propagation direction of the light sheet of the light sheet plane, comprising: Deconvolution of the three-dimensional image in the z-direction, comprising: For each intensity vector of N intensities (Ix,y,1, ... , Ix,y,N) having the same x / y value, performing a multiplication with a tridiagonal N X N deconvolution matrix, which assigns to each voxel (x, y, n) a correction parameter f1, with which, by the multiplication of the deconvolution matrix with the intensity vector (Ix,y,1, ... , Ix,y,N), for a component I(x, y, n) of the intensity vector the intensities I(x, y,n+1) and I(x, y,n-1) of the corresponding voxels of the neighboring image plane are multiplied, before they are subtracted from the intensity value I(x, y,n), wherein the distance between two light sheet planes is larger than the half-width of the point spread function in the z-direction.

2. The method of claim 1, wherein the correction factor f1 is selected to be larger with increasing penetration depth of the light sheet into the tissue in the x-direction.

3. The method of claim 1 or 2, further comprising: Determining the correction factor f1 at the point of entry of the light sheet into the tissue by: calculating the deconvolution as well as a quality measure of the deconvolved image based on different correction parameters f1 for a section of the whole image comprising the first k voxels in x-direction starting from the penetration location of the light sheet, where k is small compared to the extent of the image in x-direction, and choosing the correction parameter f1 for which the calculated quality measure corresponds to the largest quality improvement.

4. The method of any one of the preceding claims, further comprising: Representing the dependence of the correction value f1 on the penetration depth in the x-direction by a monotonically increasing function of the correction value f1(x); calculating, for different functions f1(x), the respective deconvolution in z-direction of the whole image as well as a quality measure of the deconvolved whole image, selection of that function f1(x), which provides the best quality measure, as a function to represent the entry depth dependent correction value f1(x).

5. The method according to any one of the preceding claims, further comprising: Deconvolving the image in the x / y plane for each image plane using a two-dimensional point spread function whose radius depends on the tissue depth of the voxel to be deconvolved in the z-direction.

6. The method of claim 5, comprising: For each image plane, performing the two-dimensional deconvolution for different point spread functions that are tissue depth dependent; computing a quality measure for the deconvolved image, selecting a tissue depth dependent point spread function that leads to an improved contrast value from the plurality of point spread functions for which the deconvolution has been computed.

7. The method of claim 6, wherein for the tissue depth-dependent point spread function, a function is used the radius of which increases with increasing tissue depth, whereby, however, preferably the increase of the radius increase decreases steadily.

8. The method of one of claims 5 to 7, wherein as the point spread function a function is chosen for which the average amplitude of the high spatial frequencies lying beyond a spatial frequency threshold is a minimum or does not deviate from a minimum by more than a predetermined threshold value.

9. A computer program comprising computer program code which, when executed by a computer, causes the computer to perform a method according to any one of claims 1 to 8.