CytoSpectre: a tool for spectral analysis of oriented structures on cellular and subcellular levels

Background Orientation and the degree of isotropy are important in many biological systems such as the sarcomeres of cardiomyocytes and other fibrillar structures of the cytoskeleton. Image based analysis of such structures is often limited to qualitative evaluation by human experts, hampering the throughput, repeatability and reliability of the analyses. Software tools are not readily available for this purpose and the existing methods typically rely at least partly on manual operation. Results We developed CytoSpectre, an automated tool based on spectral analysis, allowing the quantification of orientation and also size distributions of structures in microscopy images. CytoSpectre utilizes the Fourier transform to estimate the power spectrum of an image and based on the spectrum, computes parameter values describing, among others, the mean orientation, isotropy and size of target structures. The analysis can be further tuned to focus on targets of particular size at cellular or subcellular scales. The software can be operated via a graphical user interface without any programming expertise. We analyzed the performance of CytoSpectre by extensive simulations using artificial images, by benchmarking against FibrilTool and by comparisons with manual measurements performed for real images by a panel of human experts. The software was found to be tolerant against noise and blurring and superior to FibrilTool when analyzing realistic targets with degraded image quality. The analysis of real images indicated general good agreement between computational and manual results while also revealing notable expert-to-expert variation. Moreover, the experiment showed that CytoSpectre can handle images obtained of different cell types using different microscopy techniques. Finally, we studied the effect of mechanical stretching on cardiomyocytes to demonstrate the software in an actual experiment and observed changes in cellular orientation in response to stretching. Conclusions CytoSpectre, a versatile, easy-to-use software tool for spectral analysis of microscopy images was developed. The tool is compatible with most 2D images and can be used to analyze targets at different scales. We expect the tool to be useful in diverse applications dealing with structures whose orientation and size distributions are of interest. While designed for the biological field, the software could also be useful in non-biological applications. Electronic supplementary material The online version of this article (doi:10.1186/s12859-015-0782-y) contains supplementary material, which is available to authorized users.


Sensitivity analysis of adjustable parameters
In order to test how sensitive CytoSpectre is to the exact values of user-adjustable parameters, we re-analyzed the simulated fibril dataset using different combinations of adjustable parameters. The value of a single parameter was varied at a time within reasonable limits while keeping all other settings at their default values. The value of the maximum iterations parameter N max was not varied, as the number of iterations performed during detail component extraction is mainly controlled by the convergence threshold parameter, as long as the maximum iterations parameter has been set at a high enough value. The parameters and their tested values are given in Table S1 with the default values made bold. The values estimated for the detail component were then compared with the true values. Relative change in RMS error was used to quantify the impact of different parameter values, allowing comparison between the different estimated parameters. The relative change was expressed as a percentage of the RMS error obtained using default settings.
The results are visualized for each parameter in Figures S1-S6. Full numerical data for the experiment are provided in Additional file 9.
Bases on the sensitivity analysis, most of the adjustable parameters may be fairly safely assumed to have relatively little impact on the results produced by CytoSpectre, as long as the adjustable values are kept within certain limits. The apparently large variations in RMS error observed for the mean orientation in the case of several parameters are mainly explained by the low absolute errors of this parameter. It is doubtful that an increase in the mean orientation error from e.g. two degrees to three degrees would be catastrophic in typical applications, even though the relative change would be 50 %. We expect that in most cases, errors of this magnitude are overshadowed by biological or technical variation. For most of the adjustable parameters and their tested values, the increases in RMS error were less than 5 % and exceeded 10 % only near the limits of the tested parameter value ranges. The prior size parameter b is the only exception, as increasing its value beyond 0.6 rapidly increased the RMS errors of most of the estimated parameters. However, there was still a relatively stable range of values for this parameter between 0.2 and 0.6, where the increase in RMS error stays below 10 % for all estimated parameters. Based on these results, we conclude that tuning of the adjustable parameters by the user is rarely necessary and the default values can be expected to represent an acceptable compromise for most images.  0, 5, 10, 15,20,25,30,35,40,45 Convergence threshold τ Controls the stopping threshold for the iterative detail component extraction procedure. 2,3,4,5,6,7,8,9,10 Prior size b Controls the initial angular width of the      algorithm. Before computing the DFT, the segment is zero-padded either to a size equaling the next power of two larger than the original segment size or to a minimum size of 1024 x 1024, whichever is larger. The power spectral density (PSD) of each segment is obtained as the squared modulus of the corresponding DFT and the final spectral estimate is formed as the arithmetic mean of the power spectra of all segments. The conventional shifting operation is then performed to place the lowest frequencies at the center of the spectrum. Finally, the resulting power spectral estimate is normalized such that its sum is forced to correspond to the variance of the input image, in accordance with Parseval's theorem.

1,
The power spectrum is transformed into polar coordinates for easier processing during later analysis steps. The original 2D power spectrum F (u,v), expressed in Cartesian coordinates, can be transformed into P(f,θ), expressed in polar coordinates, according to the polar coordinate transform: The transform is performed using bilinear interpolation. The number of spectral elements in the interpolated polar power spectrum is selected so as to retain or exceed the resolution of the Cartesian spectrum along both dimensions at all spatial frequencies below the Nyquist frequency f Nyquist . The numbers of bins required along the two dimensions are calculated based on the number of spatial frequency bins N uv in the Cartesian representation. Neither the zero spatial frequency bin nor frequencies above f Nyquist are included in the power spectrum after the polar coordinate transform. Moreover, only angles from 0 to 180 degrees (i.e. the upper half circle in the Cartesian power spectrum) are retained, since the lower and upper half circles of the Cartesian power spectrum are symmetric. For the angular dimension, the required number of angle bins in the polar power spectrum is calculated based on the circumference of a circle whose radius corresponds to the number of spatial frequency bins between zero frequency and f Nyquist as: For the spatial frequency dimension, the required number of frequency bins in the polar power spectrum is: Since the density of spectral elements in the Cartesian power spectrum is different at different spatial frequencies i.e. at different distances from the origin, the values of the polar power spectrum need to be adjusted by a scaling factor that is dependent on spatial frequency. An adjusted element of the polar power spectrum P(f,θ) is obtained as:

Estimation of spectral background
An estimate of spectral background is formed based on the polar power spectrum estimate to allow detection of spectral signals based on their intensity relative to the background. The background spectrum estimate represents image features, such as noise, which are not of interest when analyzing the structures present in the image. Spectral signals representing meaningful details and structures of interest in the image appear as peaks and local fluctuations on top of the smoother background spectrum. It has been widely reported, that the rotationally averaged power spectra of natural images tend to fall with increasing spatial frequency f according to a power law 1/f α with varying values of α observed in different situations [41][42] [43][44] [45]. That is, the low spatial frequency region of the power spectrum tends to dominate the high spatial frequency region. An estimate of the general shape of the spectral background is needed in order to differentiate signals from the background.
The background spectrum is estimated based on the polar power spectrum using averaging over all orientations and a smoothing procedure. The advantage of not using a parametric model for the background spectrum is greater flexibility in terms of applications and imaging systems. The downside is that since the background spectrum estimate is based on an estimate of the power spectrum of the actual image being analyzed, it is biased by spectral peaks representing signals of interest. The bias is reduced by smoothing to remove most of the spectral peaks and to retain only the general shape of the spectral background. First, the polar power spectrum is averaged over all orientations to obtain a one-dimensional spectrum. The one-dimensional spectrum is then transformed to log-log scale and linear interpolation, retaining the original number of data points, is performed to obtain equidistant data points along the spatial frequency axis. Smoothing via robust local regression using weighted linear least squares and a first degree polynomial model (robust lowess) is then performed. The smoothing procedure is implemented using MATLAB's smooth function. The method decreases the effect of outliers (i.e. elements probably belonging to spectral peaks originating from genuine signals) by assigning zero weight to data outside six mean absolute deviations. By default, the span parameter L associated with the smoothing method is 10 % of the length of the one-dimensional spectrum (i.e. the number of spatial frequency bins in the polar power spectrum). The value of the span parameter can also be adjusted by the user, but the analysis is not sensitive to the exact value. Smoothing is performed in log-log space to effectively vary the span depending on spatial frequency, allowing high frequency signal peaks to be suppressed while retaining the general spectral shape, especially the typically narrow high intensity peak near zero frequency, at low frequencies. This effect is obtained because the log-log transform followed by linear interpolation effectively corresponds to upsampling at low frequencies and downsampling at high frequencies. Finally, the smoothed one-dimensional spectrum is converted back to linear scale and tiled to form a two-dimensional polar background spectrum whose size is identical to the original polar power spectrum. The background spectrum formed this way is thus isotropic and does not introduce any orientational bias during later analysis steps.

Probabilistic representation of the power spectrum
Based on the estimated background spectrum, the values of the polar power spectrum are converted to probabilities. After the conversion, each element ( ,θ) of the spectrum is assigned a value representing the probability that the element belongs to a signal. In this context, 'a signal' should be understood loosely as some feature of the image represented by power spectral values larger than the corresponding values of the estimated spectral background. That is, elements of the power spectrum having a high magnitude relative to the corresponding element of the background spectrum should be assigned a high probability value, while spectral elements with values close to the background should be assigned low probabilities. To convert the values of the polar power spectrum to probabilities, we first utilize a normalized WOSA spectrum, defined as [46]: where ( ,θ) is the background spectrum estimate. In the assumption of vanishing correlations caused by the windowing and overlapping process, the normalized WOSA spectrum is gamma distributed with σ = 1/N S and h = N S [46] [47]. The probability ( ,θ) that the element ( ,θ) represents a signal is then obtained as the value of the gamma cumulative density function evaluated at the corresponding value of ( , ): where Γ(•) is the Gamma function.

Selection of spatial frequency cutoffs
Some of the lowest spatial frequencies are usually excluded from analysis because the definition of angle is very coarse near the zero spatial frequency i.e. close to the origin of the Cartesian power Based on this criterion and the Cartesian representation of the power spectrum, a suitable lower spatial frequency cutoff f L can be selected by excluding a circular element of the power spectrum, centered on the zero frequency bin. By requiring the number of bins along the circumference of this circle to be such that the angular resolution criterion is fulfilled, the corresponding number of spatial frequency bins to exclude, starting from the origin or zero frequency, can be obtained as: The actual spatial frequency value corresponding to the 12 th bin then depends on the physical dimensions of the image. Provided that random noise is not an issue, specifying a value for the upper cutoff frequency f H is somewhat arbitrary [9]. By default, the highest spatial frequency included in the analysis in CytoSpectre is limited only by the Nyquist frequency. Both of the cutoff frequencies can also be freely modified by the user by selecting corresponding cutoff wavelengths of choice. This may be desirable e.g. to exclude high-frequency noise or to focus the analysis on a particular wavelength range of interest.

Mixed component extraction
The mixed component represents all features of the image present in the range of spatial frequencies dictated by the upper and lower cutoff frequencies.
To extract the parts of the spectrum associated with the mixed component, all spectral elements outside the cutoff frequencies are set to zero. The mixed component can thus be expressed as:

Estimation of initial parameter values for detail component detection
The detail component represents some particular group of structures with a limited range of spatial frequencies and/or orientations, such as intracellular fibrils of given size. Such structures are reflected in the power spectrum as more or less contiguous regions of high intensity relative to background. The aim of the detail component extraction procedure is to separate these regions from the power spectrum. CytoSpectre uses an iterative method, bearing some resemblance to a Bayesian classifier and segmentation algorithms based on region growing, to separate the spectral region corresponding to the detail component. Before the iterative procedure, a number of initial parameter values are estimated.
First, the dominant orientation of the low frequency range of the power spectrum is estimated using a sliding window approach. Mean orientation is first computed for each spatial frequency bin of the polar power spectrum using the corresponding function of the CircStat toolbox [22], excluding bins below and above the lower and higher cutoff frequencies, respectively. Angle doubling is performed to account for the axial nature of the data [22] and the mean orientation angle is afterwards divided by two. The mean orientation values are then transformed to absolute differences relative to the mean orientation of the first bin to avoid problems near zero or 180 degrees. A sliding window, whose length is equal to 10 % of the total number of spatial frequency bins by default, is then applied to detect the spatial frequency at which a transition in the mean direction takes place. The window is incrementally moved by one bin at a time, starting from the low frequency end, towards higher frequencies, and a Mann-Whitney test is performed at each location to test if the mean orientations within the bins overlapping the window differ significantly from all previous bins. By default, a significance level of 0.05 is used to detect the transition. The window length L and the significance level α can also be adjusted by the user, but the analysis is not sensitive to the exact values of these parameters. Mean orientation θ within the range of spatial frequencies defined by the lower cutoff frequency and the detected transition frequency is then calculated as described above. If no transition is detected, the whole spatial frequency range between the cutoff frequencies is used for calculating the mean direction. The obtained mean direction represents the dominant orientation at low frequencies, which is mostly dominated by gross properties of the image, such as the overall orientation of a complete cell imaged at high magnification.
The spatial frequency range of the detail component is estimated by first obtaining a slice of the polar power spectrum along the expected mean orientation of the detail component, followed by a peak fitting procedure. By default, the slice contains the elements of the polar power spectrum whose orientations are within ten degrees from the dominant orientation estimated during the previous step. This is the case for typical structures, which exhibit intensity variations along orientations that are transverse to the longitudinal axis of the structures. Because orientations in the power spectrum reflect the orientation of intensity variation in the image, spectral signals originating from such structures are offset by 90 degrees relative to the orientation of the structures themselves. By default, this offset is corrected during the analysis. For structures exhibiting intensity variations along their longitudinal axis, such as myofibrils and their characteristic striated intensity patterns, the correction should not be made, as the orientations seen in the power spectrum directly correspond to the orientation of the structures of interest. The selection between these two options can be accomplished by adjusting the analysis settings accordingly. In the latter case, the spectral slice is obtained along orientations within ten degrees from the orientation transverse to the dominant low frequency orientation. The default tolerance of ten degrees can also be adjusted by the user by modifying the value of the angular tolerance parameter ε, but the exact value is not critical for the analysis. Averaging over the selected orientations is then performed to obtain a one-dimensional vector of spectral values.
Next, the values of corresponding elements of the background spectrum are subtracted from the one-dimensional spectral slice. A Gaussian peak is then fitted to this background subtracted spectrum using a non-linear least squares fitting routine. If the semi-automatic detail component detection setting is enabled and the user supplies a guess for the wavelength range of the detail component, spatial frequencies corresponding to the user-specified limits are used to constrain the peak fitting procedure. In this case, the mean of the Gaussian peak is constrained between the user-specified limits and the standard deviation of the Gaussian is constrained between ¼ and ½ of the spatial frequency range defined by the limits. Peak height is constrained between zero and the actual maximum value of the background subtracted spectrum between the user-specified limits. If the fully automatic setting is used and no guesses are available, the highest peak between lower and upper cutoff frequencies is searched. The peak is then treated as if it were a perfect Gaussian to find locations at distances corresponding to two standard deviations from the peak location. In other words, nearest locations at both sides of the peak with values less than e 2 times the peak height are searched. The spatial frequencies corresponding to these locations are then used to constrain the peak fitting procedure as in the case of user-specified limits. In both cases, the initial values of all three parameters are set halfway between the lower and upper limits. The mean ̅ and standard deviation produced by the fitting procedure, reflecting the estimated spatial frequency range of the detail component, are then stored.
After the expected mean orientation θ and spatial frequency range of the detail component have been estimated, an initial estimate of the spread of the orientation distribution is computed. This is accomplished by calculating the angular standard deviation of the polar power spectrum at spatial frequencies within two standard deviations from the mean of the Gaussian peak. The values obtained for ̅ and during the previous peak fitting step are used here. The angular standard deviation is computed with the corresponding function from the CircStat toolbox [22].

Detail component extraction
The where prior size b is an adjustable parameter with a default value of 0.5. The reason for incorporating the parameter b is to allow iterative inflation of the detail component spectral region starting from an initial smaller region, in a manner resembling segmentation algorithms based on region-growing strategies. This avoids the inclusion of noncontiguous spectral regions which are likely to originate from multiple different sources in the image. Initially, the Gaussian selector function is centered at the location ̅ ,θ , as these values represent the expected mean orientation and mean spatial frequency of the detail component. Moreover, the Gaussian function is truncated to assign zero prior probability to spectral regions outside the specified spatial frequency cutoffs.
The truncated Gaussian selector function with the initial parameter values can thus be expressed as: After the initial selector function has been formed, the initial spectral region of interest (ROI) P D,0 representing the detail component is obtained by first computing the 'posterior' probabilities of each spectral element as: ( ,θ) = ( ,θ) ( , θ) (15) and then accepting all elements of the power spectrum for which X 0 > 0.5: After P D,0 has been obtained, the process proceeds in iterative manner. By default, the maximum number of iterations N max is 10 and the convergence threshold τ is one degree, but these values can also be adjusted by the user. The iterations i = 1…N max proceed in the following steps: 1. Calculate mean orientation θ and angular standard deviation , of the current spectral ROI P D,i . Use P D,0 during the first iteration.
2. Calculate the change in mean orientation Δθ and change in angular standard deviation Δ relative to the previous values. If Δθ < τ and Δ < τ or i ≥ N max , stop and return P D,i .
3. Get updated selector function K i according to Eq. 14 by substituting θ with θ and , with , . 4. Get updated 'posterior' probabilities X i according to Eq. 15 by substituting K 0 with K i . 5. Get updated spectral ROI P D,i+1 . according to Eq. 16 by substituting X 0 with X i . 6. Go to step 1.

Orientation analysis
For each spectral component, the distribution of orientations, normalized to assume the form of a probability density, is obtained as: where P C (f,θ) is substituted either with the mixed component P M (f,θ) or the detail component

Wavelength analysis
For each spectral component, the distribution of wavelengths is obtained by first converting the modified polar power spectrum into wavelength form. As the power spectrum is a density distribution function, expressing it in terms of wavelength λ rather than spatial frequency f is not simply a matter of making the substitution f = 1/λ [48]. In addition to this substitution, the values in each spatial frequency bin need to be multiplied with a correction factor of 1/λ 2 (or equivalently f 2 ) during the transformation. This correction is simply based on the requirement that total energy be conserved. The power spectrum as a function of wavelength can thus be expressed as: (λ,θ) = ,θ where P C (f,θ) is substituted either with the mixed component P M (f,θ) or the detail component P D (f,θ) spectral ROI. The wavelength distribution, normalized to assume the form of a probability density, is then obtained as: Linear interpolation is then performed to obtain equidistant spacing of data points. The following summary statistics are then calculated: Mean wavelength: The mean of the wavelength distribution.

Median wavelength:
The median of the wavelength distribution.

Mode:
The mode of the wavelength distribution.

Standard deviation:
The standard deviation of the wavelength distribution.

Generation of synthetic phase contrast microscopy images of cell clusters
Synthetic images with bright targets on a dark background, resembling clusters of cells, were generated using a custom MATLAB script. An image size of 1024 x 1024 pixels, 10x magnification and an image pixel size of 0.68 μm x 0.68 μm were used as the basis of the images. The targets were randomly located two-dimensional Gaussians with varying sizes and orientations. First, several parameter values were obtained for each image by random sampling from uniform distributions. The maximum number of targets per image was sampled from the interval [10, 100].
The mean width of the targets i.e. the full width at half maximum (FWHM) of the Gaussians, was sampled from the interval (10, 50) μm. The standard deviation of the target width was sampled from the interval (0, 20) μm. The aspect ratio of the targets, i.e. the ratio of the target length to the target width, was sampled from the interval (1, 5). Parameters for the von Mises orientation distribution i.e. mean direction and concentration parameter kappa were sampled from the intervals (0, π) and (0, 100), respectively.
For each image, targets were located one by one, first sampling the center of each target from a uniform distribution and the direction of the target from a von Mises distribution with the given mean direction and concentration parameters using the corresponding function from the CircStat toolbox [22]. Directions were converted to orientations by subtracting π from all direction values larger than or equal to π. The width of each target was then sampled from a normal distribution having the previously obtained mean and standard deviation parameters. Widths less than 5 μm were rejected and the sampling was repeated to ensure that each target had an acceptable number of pixels. Target length was then calculated on the basis of the target width and aspect ratio values. Points belonging to the target were defined as the set of pixels located within an elliptical region at most two standard deviations from the center of the Gaussian. If any of these pixels were located in a forbidden region, the target was rejected and the random placement was repeated. Otherwise, the target was accepted. Initially, the forbidden region comprised all points located within ten pixels from the edges of the image. Whenever a target was accepted, all of the pixels belonging to the target were added to the forbidden region, preventing placement of new targets on existing ones. The target placement process was continued until the maximum number of targets had been placed or more than 1000 consecutive iterations were rejected. The synthetic images were output in TIF format and the orientation and size distributions of the targets were saved for performance calculations. In total, 1000 images were generated.

Generation of synthetic fluorescence microscopy images of cells with intracellular structures
Synthetic images of cells with intracellular fibrils were generated using the MATLAB implementation (version 1.0) of SimuCell [24] and custom MATLAB scripts. An image size of 1024 x 1024 pixels, 40x magnification and an image pixel size of 170 nm x 170 nm were used as the basis of the images. For each image, SimuCell was first used to create a single cell with the cytoplasm stained red and the nucleus blue. Random nuclear placement was constrained within two pixels from the midpoint of the image. If any parts of the cell were located less than ten pixels away from the image edge, the cell generation procedure was repeated. For the nuclei, the default nuclear model was used with the following parameters: radius 10 μm, eccentricity 0.5, randomness 0.1, mean intensity 0.5 and intensity standard deviation 0. For the cytoplasms, a nuclear-centric model was used with the following parameters: radius 50 μm, eccentricity 0.7, randomness 0.7,  Directions were converted to orientations by subtracting π from all direction values larger than or equal to π. Fibril length was then sampled from a normal distribution with an initial mean of 50 μm and an initial standard deviation of 10 μm. End points of the fibril were then calculated based on the known midpoint, length and orientation. Points along the entire length of the fibril were subsequently computed using linear interpolation with tenfold upsampling. The resulting region was then thickened by dilating using a square structuring element of ones with size 3x3 pixels. If any of the points of this region were located in a region of the image that is forbidden for fibrils, the fibril was rejected. Otherwise, the fibril was accepted. Initially, the forbidden region consisted of all pixels outside the cell or within the nucleus, dilated by a square structuring element of ones with size equal to twice the subunit length. The dilation was performed to avoid parts of the fibrils extending out from the cell or into the nucleus. When a fibril was accepted, the region consisting of the interpolated points of the fibril was first dilated using a square structuring element of ones with size equal to the given fibril spacing. This region was then added to the existing forbidden region, preventing the placement of new fibrils less than the given fibril spacing away from existing ones. If more than 10 000 consecutive iterations were rejected, the mean and the standard deviation of the fibril length were decreased by 10 % to fit more fibrils into the cell. The fibril placement process was continued until either a maximum number of 20 fibrils were accepted, the mean fibril length became shorter than the mean subunit-to-subunit distance times 15 or more than 100 000 consecutive iterations were rejected.
After the locations of the fibrils had been chosen, the image was formed by locating the fibril subunits. The subunits were located fibril by fibril, starting from one end of each fibril and proceeding towards the other end in steps whose sizes were sampled from a normal distribution with the given mean and standard deviation for subunit-to-subunit distance. The orientation of each Gaussian was sampled from a von Mises distribution having a mean orientation equal to the orientation of the fibril in question, and a concentration parameter equal to the concentration parameter of the fibrils, again using the corresponding function from the CircStat toolbox [22]. As in the case of the complete fibrils, π was subtracted from directions larger than or equal to π. The resulting image, consisting of Gaussians organized into fibrils, was used to construct the green channel of the synthetic image.
To add some interference to the green channel, resulting from the nucleus and the cytoplasm in the blue and red channels, respectively, an RGB image was first formed directly from the three color channels. The final green channel image was formed by converting the RGB image to grayscale format using MATLAB's rgb2gray function i.e. by eliminating the hue and saturation information while retaining only the luminance. As the relative intensities of the red, green and blue channels in the RGB image were 0.3, 1 and 0.5, respectively, the high intensity features of the resulting image represent the target fibrils, while weaker background signals originate from the cytoplasm and the nucleus. The motivation for this mixing operation is that in real images, the objects of interest are often accompanied by other interfering features, which arise e.g. in the case of fluorescence microscopy from unspecific binding of antibodies to different parts of the cell. The final RGB image was then constructed from the green channel image with the added interference, and the original red and blue channel images created by SimuCell. Finally, the red, green and blue channels of the image were rescaled to have relative intensities of 0.8, 1 and 0.9, respectively, in accordance with observations made from our test set images. The synthetic images were output in TIF format. In total, 1000 images were generated.

Degradation of synthetic images
Synthetic images were blurred by filtering with with MATLAB's imfilter function and a Gaussian filter having a prescribed standard deviation. The standard deviation value was adjusted to obtain images with varying degrees of blurring. The dimensions of the square filter were set equal to six times the standard deviation. Image boundaries were handled by assuming values outside the image to equal the value of the nearest pixel on the image border. Blurring was applied separately to each color channel in the case of RGB images.
Additive zero mean white Gaussian noise and Poisson noise were added using MATLAB's imnoise function. The variance of the white Gaussian noise was adjusted to obtain images with different amounts of white Gaussian noise. The amount of Poisson noise was adjusted by first scaling the image intensity values down by division with the desired noise level and then scaling the values up with the same noise level value after the noise generation step. Noise was separately generated for each color channel in the case of RGB images.

Patient-specific human iPSC lines and culturing
The collection of biopsies for generating patient specific human induced pluripotent stem cell