Deriving high contrast fluorescence microscopy images through low contrast noisy image stacks

: Contrast in fluorescence microscopy images allows for the differentiation between different structures by their difference in intensities. However, factors such as point-spread function and noise may reduce it, affecting its interpretability. We identified that fluctuation of emitters in a stack of images can be exploited to achieve increased contrast when compared to the average and Richardson-Lucy deconvolution. We tested our methods on four increasingly challenging samples including tissue, in which case results were comparable to the ones obtained by structured illumination microscopy in terms of contrast.


Introduction
Microscopy allows biologists to study life by observing at a magnified scale the objects of interest such as a cell and its inner components. In particular, fluorescent optical microscopy allows the visualization of structures of interest such as microtubules or mitochondria through special molecules called fluorophores that are attached to the sample during its preparation. However, the ability to resolve between different intensities in the image, the so-called contrast, may be affected by several factors such as blurriness due to point-spread-function of the system (PSF), Poisson noise due to photon counting, electronic noise produced by spurious charges and dark noise. While it is possible to reduce noise during acquisition by increasing the excitation laser power and/or increasing the dwell time, this can produce saturation in the sensor and photobleaching on the sample and therefore it is not always possible nor desirable. The alternative is stacking which consist on taking several images of the sample at low power and low dwell time, and then average them. However, small structures can be affected as their signal becomes closer to the noise level.
Blurriness is present in the image due to diffraction of light [1] which can also be seen as the PSF being applied to the sample by the optical system through a convolution operation. As the PSF depends also on the axial position of emitters, its effect is worse for structures in out-of-focus regions. As a result, off-focus light may introduce background and therefore worsen contrast. In this context, Koho et al. have mentioned that "a good fluorescence image should have good contrast, bright details and dark background" [2], and good efforts have been made in order to improve contrast. Different optical microscopy methods use different optical configurations to enhance the contrast by reducing the out-of-focus noise. This is obtained by performing optical sectioning such that the light arriving to the camera is predominately from the imaging plane. In the case of total internal reflection fluorescence microscopy, optical sectioning is performed by evanescent field that exponentially decays away from the surface, thereby limiting most of the illumination intensity within a few hundred nanometers. Thus, only the part of the sample within the evanescent field gets illuminated, thereby enhancing the contrast. Confocal microscopy achieves optical sectioning by focusing the illumination to a spot such that only the portion within The structure of the paper is as follows: Section 2 presents the theory behind MUSICAL with focus on the indicator function, and then the new proposed functions. The results are shown in Section 3 which summarises the results and comments of using the function on different experimental samples, in addition to a comparison with state-of-the-art denoising methods. The last part includes a discussion of the results in Section 4 and finally our conclusions and future work are presented in Section 5.

Methods
MUSICAL is an image processing technique based on the decomposition of a sequence of images that exhibit the intensity fluctuation of the independent fluorophores in fluorescence microscopy into a sequence of orthogonal vectors, also called eigenimages, that explain the observed dataset. This can be considered similar to the Fourier Analysis where signals can be represented as a sum of weighted exponentials, with the difference that the eigenimages represent the most prominent spatial structures in the temporal stack. Every eigenimage u i has an eigenvalue σ i associated that can be seen as a relevance factor in explaining the observed data with higher values representing the more relevant ones. This value also allows to order the set of eigenimages and we refer to this as the order of the eigenimage. In this definition, higher order means higher index but lower eigenvalue as they are sorted decreasingly. By doing this, we can expand on the relation with the Fourier domain as an interesting feature of eigenimages is that high order ones can be linked to high spatial frequencies [11,12].
By using the eigenvalue value the eigenimages are classified into two subsets that then form two mutually orthogonal subspaces: the signal subspace S and the noise subspace N . The assumption is that eigenimages with eigenvalues above certain threshold correspond to actual signal, while the ones with values below correspond to noise (therefore the name of the subspaces).
Regarding threshold selection, which is in principle chosen manually by the user, two rules have been suggested [11], both based on the second eigenvalue distribution across patches. From these rule, the simplest one is to use the lowest second eigenvalue (the so-called rule A) which corresponds to the first guideline published with MUSICAL and stated as the knee criterion. Unless the contrary is indicated, this is rule we followed during this work.
The workflow of MUSICAL is based on processing small overlapping sections or three dimensional patches (the third dimension being the time index) of the image of the size of the main lobe of the PSF of the system. This is estimated using four parameters given by the user: the emission wavelength of the fluorophore, the numerical aperture, magnification of the microscope and the pixel size of the microscope camera. Each 3D patch is processed to obtain a two dimensional image, that is then aggregated to the results of all the patches in the set. In the conventional MUSICAL used for super-resolution, the super-resolution is achieved during this patch-processing step by the use of an indicator function f whose value depends on whether an emitter is present or not at a particular test location. MUSICAL's original publication presented a single function, which later was generalized by Acuña et. al. [11] to: The function f is evaluated at an arbitrary test point r test and its definition is made in such a way that it produces high values when tested at emitters locations, and low otherwise. When the aim of MUSICAL is super-resolution, the test points are chosen as belonging to a grid finer than the original image grid by a factor of a parameter called sub-pixelation in the original MUSICAL paper. However, since the aim of this paper is contrast enhancement only, the r test in our case correspond to the image grid itself i.e the center of each of the pixels. For each point, the estimated PSF g (considered as a vector) is generated and then projected into the eigenimage space of the patch being processed. Therefore, in this space, the PSF can be represented by its projections elements as g = [g 1 , . . . , g N ] where g i represents the projection of the PSF in the eigenimage u i . The variables a i and b i are part of the generalized framework [11] as they can be defined flexibly allowing the definition of different indicator functions, from where the original function is just a particular case of the general function. The last parameter, α, while it plays a role in contrast, it can be considered as fix as no further studies have been made after MUSICAL publication and there its recommended value was 4 .
We provide MUSICAL's coefficients definition in the context of the generalized indicator function next.
MUSICAL: in its original form, a i and b i are binary and complementary (i.e. a i + b i = 1) and their value depends on whether the eigenimage u i is classified as signal (a i = 1 and b i = 0) or noise (a i = 0 and b i = 1). Therefore, the coefficients can be defined as follows: In this definition, the projections of the PSF in the signal subspace are allocated in the numerator while the projections in the noise subspace are placed in the denominator. This allows to exploit the fact that every frame in the dataset is the result of a linear combination of the actual emitter's PSF and therefore belongs to the span of the collection of emitters' PSF. In MUSICAL, this space is also estimated using the eigenimages in S. Hence, as N is by definition orthogonal to S it is also orthogonal to the PSFs. As a result, for actual emitters, the components in the denominator will be low as the PSF should be orthogonal to the noise subspace and therefore, the indicator function is high. It is importance to remark that in MUSICAL design, it is the denominator the component that enables super-resolution. The numerator purpose is to provide a local context and it is used exclusively for stitching between nearby patches.
Based on the generalized version of MUSICAL, we designed two indicator functions whose purpose is to enhance the contrast of fluorescence image. Their main difference with MUSICAL is that they do not rely on the noise subspace but only the signal subspace, which reflects their main purpose of contrast enhancement and no super-resolution. The new methods are defined next.
Contrast enhancement 1 (CE1): while in MUSICAL's indicator function the numerator contains only the projections in the signal subspace and the denominator only the projections in the noise subspace, in this new design shown in Eq. (3) only the signal part is used as the denominator contains all the projections regardless how they have been classified.
As a result, we have removed high spatial frequencies associated to eigenimages with low eigenvalues as they are filtered out from the numerator which is now the main component of the indicator function. On the other hand, it provides a normalization component that is not present in MUSICAL: eigenimages with high eigenvalues are weighted the same in the numerator regardless of their relevance and also without being affected by their projections on the noise subspace. This can help to enhance edges or dim structures with good fluctuation levels.
Contrast enhancement 2 (CE2): instead of weighting every eigenimage equally as in CE1, we assign a weight based on their corresponding eigenvalue σ i . As a result, while all the eigenimages are present in the denominator equally, the subspace associated with the signal will be weighted, with higher order eigenimages (the ones with lowest eigenvalues) being enhanced compared to the more prominent ones. By doing this, we expect that projectiona associated with lower eigenvalues but still containing signal to be over enhanced in comparison to CE1.
We have also expanded the two methods described previously (CE1 and CE2) by relaxing their conditions on the threshold by including a smooth transition between signal and noise subspaces. We followed the soft-threshold scheme suggested in [11] and defined a linear decay for coefficients in CE1 and CE2. We call to these methods CE1s, and CE2s respectively and their definitions are shown are shown in Eq. (5) and Eq. (6). CE1s:

CE2s
: The values for σ max and σ min can be picked programmatically by selecting two patches based on the mean image of the stack and computing their second eigenvalues. These patches correspond to the ones centered in the pixels with highest and lowest values, and their computed eigenvalue is assigned as σ max and σ min respectively.
As result, this methods are hybrids in terms of weighting. In the case of CE1s, in comparison to CE1, it adds a weight to the eigenimages proportional to their relevance. In the case of CE2, its the opposite as the extra term competes with the eigenvalue value. Therefore it alleviates the enhancement of lower order eigenimages.
Further insights into the presented contrast enhancement functions: It has been shown earlier in [13] that the the eigenimages and spatial frequencies are indeed related. The authors of [13] showed how the decomposition of the temporal auto-correlation matrix encodes not only temporal variance but also spatial distributions. The same article shows that there are multiple Fourier components of the sample present in each eigenimage, but the maximum frequency in the eigenimages increases with the order. In this sense, eigendecomposition is pseudospectral, and not precisely a spectral decomposition. Also, the presented contrast enhancement functions involve a pseudospectral low pass eigen-filtering approach, which is roughly analogous to a spectral filtering approach.
One may argue that high pass filter on the averaged image is analogous to the approach presented here. However, it is notable that the averaging operator first removes the information encoded in the temporal fluctuations arising from fluorescence photokinetics, causing a loss of information. Further, performing a high pass filter inadvertently suppresses or removes low frequency information without any consideration of statistical significance or presence of features with low frequency support, causing further loss of information. Similar arguments may be considered for bilateral filters as well. Further, the question of designing such spectral filters is also demanding. In comparison, the contrast enhancements proposed in this paper do exclude eigenimages of higher order, but the exclusion is based on the statistical evidence of their diminished significance. Therefore, similar to dimensionality reduction performed using principal component analysis, our approach chooses the statistically most significant information, which is replete with the frequency support in the eigenimages that do not get excluded.
Indeed MUSICAL is a non-quantitative and non-linear approach in its original form. However, due to the design choice b i = 1, the contrast enhancement indicator functions presented here are non-linear only in the sense of alpha, vaguely similar to the gamma correction. Nonetheless, MUSICAL's indicator functions, the original as well as the ones proposed in this paper indicate a confidence in the presence of an emitter better than the diffraction limited image. In this sense, the proposed approaches provide opportunities of better suppressing the out-of-focus light and providing a better contrast when the densities of intensities of the fluorophores vary across the sample. Examples to this effect are included in Supplement 1 (Note S3 and S4), where simulated datasets with well-known ground truths are used. It is noted from the example in note S3 that CE2 and CE2s suffer from the non-linear weighing terms (of the form 1/σ i in a i ) which gives increasing weightage to the terms of lesser statistical significance. This property indeed indicates that the presented approach can be considered as eigen-filtering where the choice of the coefficients a i can be considered as the eigen-filter design. At the same time, the importance of g i , i.e. the projection of the PSF at a candidate location upon the eigenimages is essentially analogous to identifying the coefficients of principal component analysis, and thereby indicating the significance of that location contributing to the measurements.

Results
The methods were tested on experimental data available from previous works [9,10,12]. Four samples were used corresponding to in-vitro F-actin filaments [10], microtubules in U2OS cells [12], placental tissue and myocardial tissue [9]. We have added an extensive description of the acquisition methods used for tissue in Note S1 of our Supplement 1. Details for the remaining samples can be found in their respective publications.
As a means of comparison, we included the results of applying Richardson-Lucy [5,6] (referred to as RL from now on) using the restoration package of scikit-image [14] and ACsN using the ImageJ/Matlab implementation. We have also included MUSICAL result with subpixelaiton of 1 for the different samples in order to analyze the difference between the proposed indicator functions (CE1, CE2, CE1s and CE2s) and the original one. The results of MUSICAL and our new methods were obtained using Python, by modifying just the window processing section of the MUSICAL implementation. MUSICAL results are included as the starting point as any implementation available could be used to obtain these results by setting the subpixelation to 1. The RL and ACsN results were obtained by applying both methods to the mean image which was obtained by taking the average for each pixel across time which as mentioned in the introduction, it already reduces noise. The RL was applied with the iteration parameter set to 30 and using the main lobe of the theoretical PSF for the corresponding optical system. These parameters were also the input for ACsN. For CE1, CE2 and original MUSICAL we set the threshold following the first rule of selection as described in Section 2. A special case is the threshold for placenta's MUSICAL result for which we generated an additional result using a higher value as it provided a better result. We have also included a close view of two regions of interest in most cases. These are labeled as R1 and R2 and are only marked in the mean image to provide more clarity in the results.
In-vitro F-actin filaments: the sample shows actin tagged with Phalloidin-Atto 565 imaged in a TIRF microscope, with numerical aperture of 1.49 and a pixel size of 65 nm (more details can be found in [10]). As the sample has already been shown suitable for MUSICAL and therefore shown to have good intensity fluctuations, it provides a good example for testing our CE methods both in terms of contrast and also in terms of the original MUSICAL's indicator function. It is also a sample with low background as TIRF illuminates a small volume of the sample and therefore it a good assessment of how CE methods deal with images with high difference between foreground and background. From a total of 10000 frames we extracted the first 500 frames and processed them to obtain the results shown in Fig. 1. The diameter of the strands is approximately 6 nm, therefore it is not possible to measure them through optical microscopy as the point-spread function is 2 orders of magnitude larger. However, it is still valuable as single strands can be looked at separately and also globally as a network or distribution. Our CE methods provide good definition of the strands and a good delimitation between foreground and background. It is also interesting to note how dim structures that do not appear as part of strands but as single points are also enhanced. An example of this is shown in region 1 (marked in Fig. 1(a) as R1). In both MUSICAL and our methods, two strands are clearly visible, in contrast with the mean image. It is interest to note that MUSICAL not only emphasizes the two strands but make them thinner. This is expected as MUSICAL provides super-resolution and strands are supposed to be thin. We note that since the mean image do not show this double strand, it is expected that such structure do not appear in results of RL and ACsN. We also observed the expected local normalization properties of CE methods as dim and small structures are also enhanced when compared to the mean. This can be appreciated in the structures placed on above of the two strands shown in R1. The second region (Fig. 1 R2) shows some strand cluttered and crossing each other. While this network can be partially seen in the mean image ( Fig. 1(a) R2), its intricate structure is easily revealed after improving its contrast with CE1 and CE1s. In the case of CE2, we observed that the network becomes denser with no clear definition of single strands and with a structure that resembles a mesh. This is an example of over-enhancement where normalization of dim structure do not provide better insight into the sample. For this case, we consider this result as an artifact as the image provides an almost homogeneous region instead of lines. Such artifact is produced by the over-representation of high-frequency content derived from higher order eigenimages in CE2 due to the use of 1/σ i in the definition of a i . This is not the case for CE2s whose result is similar to the one obtained by CE1. Therefore, CE methods in general, manage to enhance the local contrast and details of bright single strands in the foreground while also enhancing the visualization of dimmer structures that are hard to observe in the mean image. In addition, the distinction between foreground and background is clear for every case.
Microtubules: this sample corresponds to U2OS cells in which microtubules were tagged with Alexa 647 and acquired in widefield mode at numerical aperture of 1.49 and pixel size of 108 nm. It provides a more challenging case than the actin sample presented before as the cell is a three-dimensional structure and the illumination can activate fluorescence in a larger volume than for TIRF. This contributes to poorer contrast and provides an example where contrast-enhancement can provide a different insight into the sample as several as the brightness can change across the sample due to both density and axial plane of emitters. The whole acquisition comprises 10000 frames (details can be obtained from [12]) from which we extracted and processed 500 frames to obtain the results shown in Fig. 2. From the mean it is possible to observe how the brightness of the sample is concentrated along a circular structure where the density of strands is larger than in other regions. From region 1 (labeled as R1 in Fig. 2(a)) it is possible to observe how single strand are hard to observe without saturating brighter regions of the image. The same happens in MUSICAL ( Fig. 2(b)), RL (Fig. 2(c)) and ACsN (Fig. 2(d)).
On the other hand, our CE methods allow their visualization without saturating other regions , which provides a better assessment of the whole image as no fine tuning of particular regions is needed. The best results are obtained with CE1 and CE2s as the strands are enhanced when compared to the background. The results are not that clear in region 2 where the density of strands is higher. In this situation even though the methods contribute to improve contrast between strands and background, they look very similar to the mean and also to the RL results. Regarding the results of MUSICAL, there is subtle suppression of the structures in the middle and also in the periphery of the whole structure. This is alleviated by CE1 and CE1s without enhancing the already bright structures as much as the other methods. These results provide an example of intensity normalization where CE methods generate better local distribution of intensities than the mean. Background is clear and the dim structures are enhanced without affecting surrounding bright structures.
Placenta tissue: this sample corresponds to chorionic tissue from placenta and constitutes a different type of example when compared to the ones with strand-like structures presented before. Tissue samples are generally quite dense and therefore, differentiation between signal and background can be a problem. They may also have autofluorescence from a variety of proteins in the extra-cellular matrix. Therefore, contrast in tissue samples is a bigger challenge than in cells. For that reason, this sample provides a good study case where contrast between foreground and background is poor and enhancement may provide better insight. The results are shown on a single channel corresponding to actin labeled with Phalloidin-Atto647N and acquired in widefield illumination, with numerical aperture of 1.42 and pixel size of 80 nm. The sample thickness is 1 µm and a stack of 500 frames was processed. As the sample is thick, contrast is compromised as the out-of-focus light and background compete with the signal of interest. The structure of interest corresponds to the microvilli brush border (MVBB), which are cellular membrane protrusions in the form of small hairs located at the periphery of the chorionic tissue, on the fetal side of the placenta. These are more apparent after processing the stack as can be seen on the results labeled with the R1 suffix at the middle rows of Fig. 3. Our methods produce a range of results, with the best one being CE1 and CE1s. These show a clear definition of the edges of the MVBB against the background, reducing blurriness in a similar way to the RL method across the entire structure. An example is shown in region 1 (R1). However, we note that background is not suppressed to a zero value. However, in the interface between MVBB and the background it is possible to observe darker areas than the background which effectively helps to differentiate foreground from background. On the contrary, CE2 and CE2s completely suppress the MVBB, leaving only the brightest regions visible and several artifacts in form of artificial double layers. From the methods used for comparison, RL shows a great improvement over the mean image as the MVBB becomes visible and the background is suppressed. In the case of ACsN, this method produces a smoother results but not noticeable improvement.
Region 2 (marked as R2 in Fig. 3) shows a situation where CE2, CE1s and CE2s fail in reconstructing the right structure and is also an example where fine tuning may be necessary. On the other hand, CE1 provides the best result with some artifacts appearing for low threshold which can be improved upon increasing this value.
Finally, this sample also shows a case where threshold required tuning as the best results were obtained using a higher value than the default one (rule A). In this case, a higher value allows better results both for MUSICAL and CE1 as can be observed in Fig. 3(b) and 3(e) respectively. This also provides an example of the different sensitivity that CE1 and MUSICAL may have, as the former provided a good visualization of the MVBB even with rule A. Therefore, we observed that CE2 and CE2s failed in obtaining a valid result. On the other hand, CE1 and CE1s allow to effectively obtain better contrast on placenta tissue when compared to the mean, even if background is not removed.
In addition to the visual inspection, we have estimated the signal-to-background ratio for every image. In order to do so, we picked 5 pairs of signal and background locations, and took a region of 5 by 5 pixels around each. Then, we obtained the mean for each region and obtained the SBR for each pair of regions as a division of the mean intensity in the foreground region over the background region. The points are used for the estimation are shown in Fig. 3(a) in blue and orange color which represent foreground and background respectively. For the mean image, this results in a SBR of 2.01. In the case of CE2, the value was not possible to compute as the background was precisely zero in the background and therefore the division is undefined. Something similar occurs for CE2s which produces a rather large SBR due to small background values. However, this should not be considered as an indicator of good results as the visual inspection reveals easily artifacts on the image. A similar situation occurs with MUSICAL. On the other hand, CE1s seems to do a pretty good job at increasing the SBR slightly without introducing artifacts. This result correlates with the observation described before.
Pig heart tissue: this sample corresponds to a 100 nm thick sample, where the lipid membrane was tagged with the CellMask Orange fluorescent dye, acquired in widefield illumination, with numerical aperture of 1.42 and pixel size of 80 nm. While this sample is thinner than the placenta sample shown before, it is also denser in structures. As the lipid membrane is present across different organelles of the cell, the density of emitters can also be high which affects the contrast directly as PSF of large number of emitters overlap. In addition to the image stack, we have also used the same region taken using Structured Illumination Microscopy (SIM) [15]. SIM allows to obtain an image with improved spatial resolution compared to conventional fluorescence microscopy by a factor of two by combining several images through image processing taken under special illumination patterns (more details can be found at Note S2 of our Supplement 1). While SIM is a super-resolution technique, by producing in practice a smaller PSF, edges and also contrast are also better than conventional microscopy. Therefore, as fluctuation are not exploited, it provides an independent reference for both super-resolution and high-contrast.
Results are shown in Fig. 4 for all the methods, and the first we can observe is how the mean image ( Fig. 4(a)) displays an almost uniform structure which require fine-tuning of the histogram for a nice visualization. We have marked a region of interest on it (R1) where the striations are visible as a ladder-like structure formed by the sarcomeres. The SIM image (Fig. 4(i)) makes a better job at showing such structure due to the improved resolution and contrast as mentioned before. However, we note that that this is not a problem of resolution as the strides separation is not below the resolution limit but contrast. From the tested methods, we obtained the best results with CE1s ( Fig. 3(f)) which produced clear definition of details through contrast enhancement. As in the placenta case, background was not suppressed. A close result was obtained by deconvolution, where striation is partially visible. On the other hand, MUSICAL which performs super-resolution, reconstructs a mesh with all the details lost. This may indicate that fluorescence characteristics required by the method are not met for this dataset.
The remaining methods (CE1, CE2 and CE2s) do not show an improvement on the mean and several artifacts are visible (CE1 and CE2) or even missing information (CE2s). In the case of ACsN the image is smoother but in the process, details as striation are completely lost. Finally, as a means of comparison, we include a profile line (normalized) across the striation in Fig. 4(j). The results are shown only for the methods that worked for this sample corresponding to the mean, CE1s, RL and ACsN. This provides an example of how contrast is improved which is reflected on the clear apparition of the ladder structure of sarcomeres.

Discussion
The original MUSICAL indicator function is defined in a way that the denominator and numerator components have different purposes. While the denominator is the part that enables superresolution as its value approaches zero when tested at emitters location, the numerator provides local context which enables stitching of the results of different patches. Therefore, as a superresolution technique, it is a method that exploits predominantly the noise subspace of the imaging signal. On the other hand, our CE methods exploit only the signal part with the purpose of enabling only contrast enhancement. After testing the new methods on a variety of experimental samples, we have proven this concept for two of our methods. CE1 and CE1s present the most promising results as improvement in the quality of the images was consistent across the four samples. On the other hand, CE2 and CE2s appear to over-enhance features. This is noticeable for both tissue samples where several artifacts were obtained across the entire resulting image ( Fig. 3(f,h), and 4(f,h)). However, CE1 and CE1s worked for every situation with almost no artifacts and great separation of foreground and background. In the case of CE1, we notice that even if the threshold is chosen roughly using rule A, which can be implemented programmatic, results were good. This is completely different to what was observed for the same sample in MUSICAL image where the result varied dramatically for different thresholds. An example of this is shown in Fig. 3(b) and Fig. 3(e) where results for two different thresholds are displayed. In addition, for CE1 the threshold appears to affect the local saturation but not structures recovered. This can be translated into a simple rule of thumb where the threshold can be increased until the desired saturation is obtained.
Regarding contrast, we noticed that CE1 and CE1s equalize the brightness of different structures for every sample. This is noticeable specially where structures that appear dim in the mean, are clearer and crisper after processing the stack. This can observed in Fig. 4 where structures and details are enhanced regardless of their average brightness, with results very similar to the ones obtained by SIM. This can be explained by how MUSICAL's family of indicator function are sensitive not only to the mean but also to the second moment [12]. By decomposing the image into eigenimages and re-weighting the basis (in CE1 all the eigenimages weigh the same while in CE2 they are scaled according to their eigenvalues). As a result, CE methods can enhance dimmer structures provided they show intensity variations.
It was also observed that depending on the sample, background is not entirely suppressed by the CE methods. This is particularly visible in the tissue sample, which were precisely challenging samples as the background can be indeed a source of light. This is expected as this samples are dense or have several sources from off-focus planes. Therefore, as these method work locally, they also enhance the low signal found in the background. However, edges were not compromised and therefore, the effect on the background is not entirely negative at least for CE1 and CE1s. In order to provide the user a darker background if needed, CE methods and the mean could be combined using image multiplication which would weight down the background and increase the foreground. However, the user should be careful as the mean could still have dim structures in the background.
As the result for CE1 and CE1s, it might be interesting to try on a larger number of samples or even revisit previous results where methods as MUSICAL may have failed. We include a small study about the sensitivity of CE1s to the accurate knowledge of the PSF in Note S5 of Supplement 1. It shows an added advantage of less sensitivity of CE1s in particular and the contrast enhancement indicator functions in general over MUSICAL.
In order to avoid subjectivity, we computed some metrics on the image with the purpose of quantifying contrast enhancement. In particular, for each image we calculated the BRISQUE index, standard deviation of the pixel intensities, and the entropy of the pixels. However, the results were inconclusive as no index was coherent with the observations or had a definitive trend that could be insightful. The numerical results are included in Note S2 of Supplement 1 for the sake of completeness.

Conclusion
We have identified that it is possible to process a stack in order to obtained higher definition and contrast in fluorescence images. This work presents a novel approach for enhancing the local contrast of fluorescence microscopy images as alternatives to the mean. These methods, referred to as CE methods, are based on the indicator function of MUSICAL, originally designed for super-resolution. Such function was modified by dropping the super-resolution component and leaving only the local contrast part. We have tested the different designed methods in four samples with different degree of complexity and compare these results with Richardson-Lucy deconvolution and ACsN. Particularly, the results of our methods can be compared to the ones obtained with deconvolution but with an improved contrast of dimmer structures as they appear enhanced in the final image. From the four CE methods designed in this work, we checked that for tissue samples, CE1 and CE1s provide very detailed results. We concluded this after comparing the results of the same sample under SIM as it provided a good reference in terms of contrast and also resolution.