K-factor image deshadowing for three-dimensional fluorescence microscopy

The ability to track single fluorescent particles within a three dimensional (3D) cellular environment can provide valuable insights into cellular processes. In this paper, we present a modified nonlinear image decomposition technique called K-factor that reshapes the 3D point spread function (PSF) of an XYZ image stack into a narrow Gaussian profile. The method increases localization accuracy by ~60% with compare to regular Gaussian fitting, and improves minimal resolvable distance between overlapping PSFs by ~50%. The algorithm was tested both on simulated data and experimentally.

Other factors that degrade the accuracy are out-of-focus fluorescence, scattering, background noise and shot noise (especially where the photon count is low). Overlapping PSFs are a severe problem since they disrupt the localization process. Therefore, regions of interest (ROIs) exhibiting elongated (asymmetric) intensity profiles within a frame are often discarded from the sample set, effectively reducing the sampling density and increasing the image acquisition time 13 . This problem is more severe in the case of 3D imaging, due to the nature of the PSF in the axial direction, which has long tales that penetrate near-by emitters. In addition, in focused PSFs have a high probability to overlap with out of focus PSFs 14 .
3D localization methods include analysis of the diameter of the rings appearing in the defocused images 15 . Other methods are based on breaking the axial symmetry, as is done by the biplane 16,17 and astigmatic detection 18 methods. These techniques achieve lateral and axial resolutions of ~10 nm and ~30 nm respectively. However, they require more complicated optical systems, and complex post processing algorithms. In addition, the issue of overlapping PSFs must still be addressed. 2D techniques for overlapping emitter detection 19,20 are much more complex in 3D.
In this paper we present a modified nonlinear image-factorization K-factor algorithm designed for the enhancement of different contrast levels within an image 21 . This algorithm sharpens the point source PSF, which provides two benefits. First, the 3D PSF shape is transformed into a 3D Gaussian profile, which simplifies the localization process and reaches the same high precision accuracy. Second, two overlapping PSFs become effectively separated, thus are able to be separately localized in 3D. The proposed algorithm is a computer post processing technique that can be applied to the raw data prior to emitter localization. In a previous paper, the algorithm was applied to 2D images 22 showing improvement in localization precision or data acquisition time of ~40%. Here we present a modified version, which is applied to 3D localization.

Theoretical Background
The 3D PSF of an aberration-free defocused imaging system with a finite lens aperture is given by an Airy function 23 that expands with defocus. Each optical section includes both in-focus plane and out-of-focus contributions that are determined by the corresponding pupil function 24 . The ideal PSF function has circular symmetry about the z-axis, and mirror symmetry about the central xy-plane, and can be described by the Gibson and Lanni model that assumes an optical path that includes a biological sample, a cover slip layer and an immersion layer 25 . A mathematically simpler model is a Gaussian function 26 . The measured intensity is degraded by added photon shot-noise and background noise that is created mainly by out-of-focus and extraneous fluorescence and charge coupled device (CCD) readout noise. The Gaussian model for the intensity at (x, y, z) of a fluorescent particle located at (x 0 , y 0 , z 0 ), is given by 27 : where N is the total number of collected photons during acquisition period, σ x,y and σ z are the standard deviation of the Gaussian, η B is a Poisson distributed random variable with variance N b that describes the background noise, and η shot is a Poisson distributed random variable that describes the shot noise with a mean equal to the square root of the total intensity in each pixel 28 . In the absence of aberrations, the standard deviations will be given by 24,29,30 : where N.A. is the numerical aperture of the objective, λ is the wavelength of the emitted light and n is the refractive index of the medium between the coverslip and the objective front lens element. The Gaussian model is a good approximation for confocal laser scanning microscopy (CLSM), however for wide-field and 3D PALM, a more realistic model must be used 26 .

K-Factor Algorithm
The K-factor algorithm is an nonlinear image decomposition method 21,31 , that divides the image into a set of contrast based factors. The reconstructed image is the joint product of all harmonics in the decomposition. Here we demonstrate a modification of the original algorithm. where k n is the contrast depth parameter and its value varies between 0 and 1. In the original decomposition we used, the parameter k was constant for all the different factors. However, varying k yields improved results, as will be discussed further on. g n is a binary image computed as: This is an iterative process. For each n, g n is calculated as a binary image with pixel values that are either 1 or 0 according the threshold set by k n . Once g n is calculated, f n can be found. The next iteration will use the original image divided by the factor f n as the threshold for the next iteration of g n . Due to the nature of the decomposition, when a fixed parameter k is used, the first few factors contain the details with the highest contrast in the image. These details are the desired image information along with some spatial noise. Higher order factors contain mostly spatial noise, together with some fine spatial information associated with low contrast levels. In a previous paper, we reported that the joint product of the original image and the first few factors of the decomposition, produces a sharper image, with data that is de-emphasized and reduced noise 22 . The choice of k controls the contrast depth of the factors. A high value close to one produces an image version with fine geometrical details of the image, but requires many factors. A value close to zero requires relatively few factors but produces large contrast steps and a spatially coarse version of the image. Since the desired image information is usually concentrated mostly in the first few factors, the factors were assigned a high k value. The k-value was gradually reduced for the higher factors. In order to find an optimal set of parameters for the decomposition, an object consisting of two closely spaced Gaussian point sources was tested. The parameter k for each factor was allowed to vary, together with the number of factors M and the additive noise level of the image. The search was designed to maximize the height of the maxima of the PSFs, while minimizing the saddle between them, therefore maximizing the PSF peak to minimum saddle ratio. The optimal choice was found to be given by: for the previous algorithm. In addition, the σ of the PSF was reduced by 50%. This reduction enables the detection of overlapping emitters with reduced separation. Furthermore, here, only 9 factors are needed for the decomposition, compared to 48 in the previous work. The proposed technique can be added to the conventional PALM workflow. The 3D image acquisition technique begins with the capturing of a set of z-stack images of the sample. The images are then processed with the K-factor algorithm. Next, a localization routine using methods of single emitter fitting produces the 3D emitter positions. The final super-resolved image is the linear combination of all emitters coordinates.

Simulations
Monte-Carlo (MC) simulations were used in order to generate simulated data sets with emitters at different lateral and axial positions in each set. In these simulations, the model was a fluorescence point source emitting at λ = 540 nm that was imaged through an objective lens onto a CCD camera. The model parameters were chosen to correspond to an optical system consisting of a 63x/NA = 1.4 oil immersion objective lens, imaged though a 1.0x relay lens onto a CCD sensor array with 6.45μm × 6.45 μm pixels, which translates to 102 nm × 102 nm in the object plane with the 63x objective. The z-step size was Δz i = 200 nm. Background noise was introduced by adding a Poisson distributed random value with the parameter N b , that can be measured experimentally and therefore is considered to be known 32 . In our simulations we choose an average sized background noise of N b = 5 photons. Shot noise was added as a Poisson process with an expected value which corresponds to the noiseless pixel values and a standard deviation (STD) that equals the square root of the value of each pixel. The value for σ x,y and σ z for the given imaging parameters was calculated from Equation (2) to be 167 nm and 413 nm respectively. Figure 2 is the simulated lateral and axial PSF profiles with added shot noise (N = 5000) and N b = 5. Figure 2(a) is a simulation of the lateral profile of the PSF and Figure 2(d) is a simulation of the axial profile of the PSF. When the K-factor algorithm is applied to the PSF, the areas with highest contrast in the PSF are enhanced and low level regions are discarded. The resulting PSF shape is well approximated by a Gaussian, both in the lateral (Fig. 2b) and the axial (Fig. 2e) directions. In order to compare the original and K-factor processed images, the line profiles of the yellow dashed lines in Fig. 2(a,b) is presented in Fig. 2(c), and the line profiles of the yellow dashed lines in Fig. 2(d,e) is presented in Fig. 2(f).
The transform of the PSF into a Gaussian shape improves the localization precision of each emitter. In order to test this assumption, the raw simulated data set of z-stack images of a single emitter with added shot noise and background noise (N b = 5 photons), and the K-factor filtered data, were processed with a least square localization routine that was implemented in MATLAB. The position of the fluorophore and the number of arriving photons N were allowed to vary. The initial z position of the emitter was chosen Where N, σ x,y , σ z , x 0 , y 0 and z 0 are all fit parameters. The algorithm's initial coordinates were the local maximum. The fit was carried out by taking 9 × 9 × 9 pixels around the initial location 25 . The fit produces the best estimate of the position of the fluorophore and was repeated for 1000 MC iterations yielding a set of L = 1000 localization positions ( , ,ˆx y z i i i ) that can be compared with the known original position (x i , y i , z i ). The root mean square (RMS) localization error was computed for both the raw and K-factor filtered simulated data: Figure 3 compares the RMS localization error of the fitting process for the lateral direction (Fig. 3a) and for the axial direction (Fig. 3b) as a function of detected photons with added shot noise and background noise (N b = 5). Gaussian fitting applied to raw data is presented in the black line. Gibson & Lanni fitting applied to the raw data is presented in the blue line and Gaussian fitting applied to the raw data that underwent processing with the K-factor algorithm is the red line. As was expected, for the raw data, the Gibson and Lanni model yielded good results for both the lateral and axial directions, whereas for the processed data, the Gaussian model yielded similar results. These results are more accurate by a factor of two with compare to a Gaussian fitting of the raw data. This is due to the K-factor algorithm's ability to transform the 3D PSF shape to fit that of a Gaussian. With regard to the photon number, as can be seen in Fig. 3, an increase in N lowers the RMS error. This is due to the shot noise that is a Poisson random process with a rate that depends on the total number of photons detected, and is proportional to √N and therefore is reduced with higher N 33 .
Another advantage of the proposed technique is that it effectively narrows the PSF width. When the width is smaller, overlapping PSF can be identified and localized separately. Since frames that contain overlapping PSFs are usually discarded, the ability to detect them will decrease image acquisition time. This problem is more pronounced in 3D imaging, since the contribution from each emitter is spread over a larger volume. In order to determine the minimal distance at which two emitters can be resolved, the distance between two emitters for both the lateral and axial axes and the added shot noise and background noise parameters were varied. The results are summarized in Table 1.
It is seen that for the lateral axis, the K-factor reduces the minimum distance required to resolve two point emitters by more than a factor of two, and reaches 250 nm.The improvement is even greater in the axial axis, due to the widening of the PSF in the z-direction. The axial PSFs profile for the distances that are presented in Table 1 is shown in Fig. 4 with N = 5000, N b = 5. Figure 4(a) is the original axial profile. Figure 4(b) is the K-factor algorithm applied to Fig. 4(a,c) is line profiles of the dashed lines in  . (a,b) are the lateral and axial RMS localization error, respectively. Gaussian fitting applied to raw data is presented in the black line.
Gibson & Lanni fitting applied to the raw data is presented in the blue line and Gaussian fitting applied to the raw data that underwent processing with the K-factor algorithm is the red line. Fig. 4(a,b). In the original image, the two PSFs were undistinguishable, whereas with the K-factor processed image, there was a distinct separation between the two PSFs, as well as the reshaping of the PSF into a Gaussian shape.

Materials and Methods
Peripheral-blood lymphocytes (previously extracted from whole blood) were placed on a glass slide, and the telomere labeling was carried out using a Dako K5326 Telomere PNA FISH Kit/Cy3 (Dako, Carpinteria, Canada). Briefly, the cells were fixed and stripped on the slide as per the protocol, and hybridization was done with the HYBrite system (Abbot Laboratories, Des Plaines, IL, USA). The slides were then counterstained with DAPI and mounted using DAPI II counterstain medium (Abbot Laboratories, Des Plaines, IL, USA). The slides were then covered with a #1.5 glass coverslip and sealed with nail-polish. All experiments which involved human blood samples were carried out in accordance with the approved guidelines and regulations, and the experimental protocols were approved by the Meir Hospital Institutional Committee for Ethics. Informed consent was obtained from all study participants.
Z-stacks of the FISH labeled telomeres were acquired using a fully automated Nikon TE2000E inverted fluorescence widefield microscope, through a 60x/NA = 1.4 oil immersion objective. The excitation source was an xCite 120 W metal halide source (Excelitas Technologies Corp, Waltham, MA, USA), and a narrowband Cy3 filter set (#41003, Chroma Technology, Bellows Falls, Vermont USA) was used. Images were acquired with a Retiga 2000 R cooled CCD camera (QImaging, Surrey, BC, Canada) with 7.4 μm × 7.4 μm pixel size. The system was controlled with Nikon's NIS elements software. Multiple fields were acquired automatically, and in each field, a Z-stack of 17 slices at 200 nm spacing was acquired. The multifield ND2 images were preprocessed using the Fiji distribution of ImageJ 34 . The raw ND2 files were read using the Bioformats plugin 35 . The desired field and channel was extracted, and the image stack was cropped to an area of 135 × 165 pixels, selected to include one cell nucleus, in order to reduce computation time, with an effective pixel size of 117.5 nm. The cropped image stack was saved as a TIFF image sequence for further downstream processing in MATLAB.
The K-factor algorithm was implemented in MATLAB (version 2012b, MathWorks, Natick, MA, USA). The program was run on a HP Compaq Elite 8300 Microtower PC with Windows 7 Professional 64 bit operation system, Intel ® Core ™ i5-3470 processor, 3.20 GHz, 12 GB RAM.

Experimental Results
The raw z-stack FISH labeled telomeres images (described in the methods section) were filtered with the K-factor algorithm, after which a LS localization procedure was applied in MATLAB. The raw data was fitted to a Gibson and Lanni PSF profile, while the K-factor processed images were fitted to a Gaussian profile. The fluorescence data typically included dozens emitting of telomeres simultaneously. Figure 5(a-e) presents a sequence of XY images at two planes apart, with Δz = 400 nm between frames. Figure 5(c) is the most focused image, as determined by eye. Figure 5(f-j) are the frames after the application of the K-factor algorithm. The obtained profile of the raw data PSF, matched the one expected from Fig. 2(d), and is widening as the defocus grows. For the K-factor processed images, the PSF profile matches Fig. 2(e) and is narrowing as the defocus grows. The localization routine was done both on the raw data images and on the K-factor processed images. The focused images were chosen as the initial z coordinates for the localization routine, and the 3D neighborhood around each local maxima were fitted using the LS algorithm with the corresponding PSF shape. The unprocessed focused frame from the z-stack measurements is shown in Figure 6(a). Overlapping emitters can be seen as spots with a larger radius. One such region of closely spaced molecules is marked by a yellow frame, which is magnified in Figure 6(b). By applying the K-factor algorithm on each individual frame (Fig. 6d), areas with overlapping molecules became distinguishable, as can be seen in the magnified area (Fig. 6e). Figure 6(c,f) show the cross section of the dashed line that passes through the center of the overlapping two emitters.
The PSFs recorded in Fig. 6(b) are longitudinally extended and lacks any form of symmetry, therefore in the absence of serious optical aberrations or distortions, it clearly originates from multiple emitters, and therefore the localization process is not applicable in this case. In the K-factor processed image (Fig. 6e), these PSFs were localized as two distinct emitters. This is due to the reduction in the PSF width, which reduced the saddle between the PSFs and enabled their localization. The same effect occurs in the axial direction, as can be seen in the axial direction as can be seen in Fig. 7. Each pixel in the axial direction corresponds to 200 nm, whereas in the lateral direction it is 117 nm.  Table 1. Minimal distance at which two emitters can be resolved by both the standard LS method and the proposed K-factor.
The Gibson and Lanni PSF profile is clearly seen in Fig. 7(a), whereas it becomes a Gaussian shape in the K-factor processed image Fig. 7(b,c) is the line profiles of the red dashed lines in Fig. 7(a,b). In the original image there are overlapping PSFs that appear as large spots that resolves into two separate PSFs in Fig. 7(b). Furthermore, in the processed image, background noise is reduced which makes the background appear totally black. An important remark is the existence of Spherical aberration when performing 3D imaging of thick samples. This aberration is due to the refractive index mismatch between the imaging medium and the objective, which effectively increasing the apparent z position of each emitter. Such aberration is clearly present in the data in the defocussed images of Fig. 5(a) compared to Fig. 5(e). There are numerous ways to overcome this aberration 36 . One is to use an imaging medium with a higher refractive index. The other is to apply the localization procedure only on emitters that are below the focal plane. Alternatively, there are objective lenses available with correction ring for spherical aberration 37 . Here, when applying the K-factor algorithm to the PSFs image, the result is a Gaussian shape. Simulations show that the effect of the spherical aberrations with the proposed technique may cause a shift of the Gaussian shape which causes a constant drift error in the z position throughout the entire image. Since this is a constant drift, it does not affect the obtained 3D reconstruction and can be corrected using a rescaling factor, without any decrease in the localization precision.
In order to calculate the localization precision at all directions, the Gaussian localization routine was applied to 100 different emitters in each of two sets, and the standard deviation of particles coordinates was calculated for 20 repeated experiments. The standard deviation σ x,y of the raw data was 5.2 nm, which  corresponds to a resolution of (FWHM = 2.35σ) 12.22 nm. The K-factor algorithm applied on the row data yielded a standard deviation of 2.3 nm, which corresponds to resolution of 5.4 nm, an improvement of ~60%. With the same analysis, the axial resolution was 86.15 nm with compare to 27.2 nm in the processed data which is an improvement of ~70%. The Gibson and Lanni localization routine applied to the raw data yielded resolution of 28.6 nm. This relatively low resolution stems from the sampling in the z-direction that was done with steps of 200 nm. A smaller Z-step would increase the resolution. In addition, for each set, the number of detected PSFs was measured and averaged. The results show an increase of 50% in the density of the PSFs with the K-factor algorithm, thus enables an increase in the activated fluorophore density inside the sample. As a last step to prove the viability of the proposed technique to a 3D image of a biological sample, samples of human epidermoid carcinoma cell line, A431 38 , that were loaded with 20 nm spheres GNPs immobilized on a coverslip, using a known protocol were used 39,40 . Each sample was illuminated using a green laser at 532 nm (Photop DPGL-2100 F). The Z-stack set of images  of the scattered light from the sample was acquired using a fully automated Nikon TE2000E inverted fluorescence widefield microscope, through a 40x/NA = 0.6 long working distance air objective. Images were acquired with a Retiga 2000 R cooled CCD camera (QImaging, Surrey, BC, Canada). Multiple fields were acquired automatically, and in each field, a Z-stack of 40 slices at 100 nm spacing was acquired. The multifield ND2 images were preprocessed using the Fiji distribution of ImageJ 34 . The raw ND2 files were read using the Bioformats plugin 35 , where they were cropped to contain a single cell. The obtained 2D z-stack set is merged in order to create a 3D image of the cell using post processing and is presented in Fig. 8. Figure 8(a) is a phase image of the single cell. Figure 8(b) is the obtained 3D image of the cell using the original z-stack. The K-factor algorithm was applied to each z-stack frame individually and the entire K-factor processed set was used to create the 3D cell's image (Fig. 8c). The scale bar is identical for both Fig. 8(b,c). The original 3D image contains sharp peaks due to the PSF shape in the axial direction, which broadens with defocus, whereas in the K-factor 3D processed image the result is smoother without sharp edges, due to the reshaping of the PSF into a Gaussian shape. Thus, the technique is applicable to 3D data.

Summary and Conclusions
This paper presents the use of the K-factor algorithm for enhanced performance of 3D localization microscopy techniques. There are three major advantages of the K-factor algorithm over the used Gibson-Lanni model. For the case of a single PSF, the Gibson-Lanni model is very noise sensitive, which results in an increased localization error, whereas the K-factor algorithm is a noise reduction technique, thus is much less sensitive to noise. In addition, the computational complexity involving the localization using the Gibson-Lanni model is very high and time consuming. The fit of the PSF to a Gaussian profile is fast and simple, however results in an enlarged localization error. The K-factor algorithm reshapes the PSF into a Gaussian profile, enabling its localization using a regular Gaussian shape with the same precision as the Gibson-Lanni model and at a shorter period. For the case of closely spaced PSFs, the Gibson-Lanni model fails to perform their localization, whereas the K-factor algorithm separates the PSFs, which enables their localization. The proposed approach is generic and can be applied to any 3D microscope configuration, where the ability to use Gaussian fitting with high accuracy on 3D data can facilitate the computational complexity, hence reduce the processing time required for the generation of the 3D superresolved image.