Blue-noise sampling for human retinal cone spatial distribution modeling

This paper proposes a novel method for modeling retinal cone distribution in humans. It is based on Blue-noise sampling algorithms being strongly related with the mosaic sampling performed by cone photoreceptors in the human retina. Here we present the method together with a series of examples of various real retinal patches. The same samples have also been created with alternative algorithms and compared with plots of the center of the inner segments of cone photoreceptors from imaged retinas. Results are evaluated with different distance measure used in the field, like nearest-neighbor analysis and pair correlation function. The proposed method can effectively describe features of a human retinal cone distribution by allowing to create samples similar to the available data. For this reason, we believe that the proposed algorithm may be a promising solution when modeling local patches of retina.


Introduction
Sampling is the reduction of a continuous signal into a discrete one, or the selection of a subset from a discrete set of signals.For sampling to be effective, samples should be uniformly distributed in a way that there are no discontinuities; but at the same time, regular repeating patterns should be avoided, to prevent aliasing.In the human retina, the mosaic of the cone photoreceptor cells samples the retinal optical projection of the scene, achieving the first neural coding of the spectral information from the light that enters the eye.To solve the sampling problem, the human retina has adopted an arrangement of photoreceptors that is neither perfectly regular nor perfectly random.Local analysis of foveal mosaics shows that cones are arranged in hexagonal or triangular clusters, but extending this analysis to larger areas shows characteristics such as parallel curving and circular rows of cones associated with rotated local clusters.
There are different theories regarding the regularity and development of the cone cells mosaic.Wassle and Riemann [61] proposed two models based on mechanisms that assume the self-regulation of an original random pattern, one with a repulsive force acting between nerve cells and the other based on competition for territory for each neighboring cell.Later, Yellott [67] discovered that the photoreceptors in the human retina, especially the cones, are distributed conforming to a Poisson disk distribution.He performed spectral analysis to an array of cones treated as sampling points and observed that the spectral properties of cones mosaic are representative of a Poisson disk array, with the additional restriction of a minimum distance between the center of the cells and their nearest neighbors, because of the size of the cells.This was confirmed by Galli-Resta et al., which investigated the spatial features of the ground squirrel retinal mosaics, suggesting that a minimal-spacing rule d min in conjunction with an adequate density of receptors can adequately describe the array of rods and S cones [27].Poisson disk distribution is now regarded as one of the best sampling patterns, by virtue of its blue-noise power spectrum [38].
It is still unclear how the spatial distribution and mean density of cones can affect the sampling of a retinal image [17].An interesting evidence of this open issue is the experi-ment from Hofer [33] which tested the perception of stimuli of small spatial scale.Showing brief, monochromatic flashes of light of half the diameter of a cone size on previously characterized retinal areas of the subjects, they described the same stimuli with a large number of hue categories, including white, blue and purple, indicating that the stimulation of two different cones with the same photopigment results in different color sensations, even with no stimuli in different regions of the retina or on other wavelength-sensitive cones.
In this study, we examined the Nearest Neigbour (NN) regularity index of the population of cones in images of real human retina.We then compared the results to another measure of spatial patterning, the Pair Correlation Function.The goal of this paper is to show that the sampling properties of the cone photoreceptor mosaic can be modeled by a bluenoise algorithm, and that they can be used to generate sampling arrays with the same features of the retinal cone mosaics.More specifically, we want to identify an algorithm capable of generating sampling arrays with the same range of densities in the retina, and use specific metrics to compare the spatial and spectral properties of the cones distribution.

Retinal and Cone sampling modeling
The most recent works on retinal modeling are focused on the neural behavior [63,45,44]; for example, Virtual Retina by Wohrer and Kornprobst is a large scale simulation software that transforms a video input into spike trains, designed with a focus on nonlinearities, implementing a contrast gain control mechanism.
However, there have not been many attempts at modeling the cone sampling array.The first known sampling model for positioning cones in the retina with the same qualities as the human sampling was described by Ahumada [4].It works by placing cones, which are surrounded by circular disks representing their region of influence, starting from the center of the retina, and then applying a random jitter to each point.There is an attempt to generate a space-varying parameter model, to extend the modeling capabilities past the foveola, by varying with the eccentricity the mean radius of the cone disk, the standard deviation of the cone disk radius, and the standard deviation of postpacking jitter; but ultimately those parameters seem to be only fit for the foveola.
After their studies on human photoreceptor topography, Curcio and Sloan continued in Ahumada's direction proposing a model of cones distribution based on regular arrays subjected to a spatial compression and a jitter, to fit the actual cones mosaic [14].Their analysis was based on the distribution of distance and angles of neighboring cones, com-paring real mosaics with artificially generated ones, and evidencing anisotropies in retinal cell spacing.
Another attempt at modeling the sampling properties of the cone mosaic was proposed by Wang [60], which created a polar arranged array of cones and jittered the points according to the standard deviation of a Gaussian distribution, constrained by a minimal spacing rule.The comparison of power spectrum of the human foveal cones and the generated sampling arrays show similarities, and the generated arrays exhibit some basic features of the mosaic of foveal cones.
In Deering's [16] human eye model, cones are modeled individually as a center points surrounded by points that define a polygon constituting the boundaries of the cell, each photoreceptor is then subjected to attractive and repulsing forces to adjust its position.This retinal synthesizer is then validated by calculating the neighbor fraction ratio and by empirically measuring the cone density in cells/mm 2 and comparing it from data from Curcio et al. [15].

Blue Noise Distributions
Coined by Ulichney [58], the term blue noise refers to an even, isotropic, yet unstructured distribution of points.Blue noise was first recognized as crucial in dithering of images since it captures the intensity of an image through its local point density, without introducing artificial structures of its own.It rapidly became prevalent in various scientific fields, especially in computer graphics, where its isotropic properties lead to high-quality sampling of multidimensional signals, and its absence of structure prevents aliasing.It has even been argued that its visual efficacy (used to some extent in stippling and pointillism) is linked to the presence of a blue-noise arrangement of photoreceptors in the retina discovered by Yellott [67].Over the years, a variety of research efforts targeting both the characteristics and the generation of blue noise distributions have been conducted in computer graphics.
Arguably the oldest approach to algorithmically generate point distributions with a good balance between density control and spatial irregularity is through error diffusion [26,58], which is particularly well adapted to low-level hardware implementation in printers.Concurrently, a keen interest in uniform, regularity-free distributions appeared in computer rendering in the context of anti-aliasing [12].Cook [11] proposed the first dart-throwing algorithm to create Poisson disk distributions, for which no two points are closer together than a certain threshold.Considerable efforts followed to modify and improve this original algorithm [43,41,36,7,28].Today's best Poisson disk algorithms are very efficient and versatile [20,22], even running on GPUs [62,6,66].
Thanks to the pioneering work by Dippé and Wold [18], Mitchell [42], Cook [11], Shirley [56], the computer graph-ics community became sensitive to the fact that noise and aliasing are tightly coupled to sampling.A large variety of optimization-based approaches has been proposed since then.Two main optimization-based approaches have been developed and presented in numerous papers: (1) on-line optimization [41,21,40,5,6,23,8,54,53,25,31,68,48,32,49], and (2) off-line optimization [47,37,46,59,2,1], where the nearoptimal solution is prepared in form of lookup tables, used in runtime.The present work uses as reference the approach called Blue Noise Through Optimal Transport (BNOT), developed by de Goes et al. [31], because it allows to achieve the best Blue Noise distribution known today.
In an effort to allow fast blue noise generation, the idea of using patterns computed offline was raised in [18].To remove potential aliasing artifacts due to repeated patterns, Cohen et al. [9] recommended the use of non-periodic Wang tiles, which subsequently led to improved hierarchical sampling [37] and a series of other tile-based alternatives [47,39,46,3].Wachtel et al. [59] propose a tile-based method that incorporates spectral control over sample distributions.More recently, Ahmed et al. [1] proposed a 2-D square tilebased sampling method with one sample per tile and controllable Fourier spectra.However, all precalculated structures used in this family of approaches rely on the offline generation of high-quality blue noise.

Methods
The cone mosaics used for this work are from previously published images of patches of real human retinas, as shown in the leftmost boxes of Figures 2 through 5; they were acquired from the pdf versions of the papers or html, if available, and saved as png images.The pictures are from different subjects of various ages and were obtained with different techniques, from histological tissue prepared for electronic microscopic imaging in [15,35,13,30], to the most recent in vivo imaging techniques, using adaptive optics like deformable mirrors coupled with a wavefront sensor to compensate for the ocular aberrations of the eye [51,57,55,64].The x and y coordinates of the cells inner segments were manually plotted using WebPlotDigitizer [50].This preliminary work has been based on a relatively small dataset sue to the difficulty of finding wide collections of retinal images.We understand these difficulties related also with problem of the use of different imaging techniques and tissue preparation and we hope to have larger datasets in the future.When analyzing the points distribution, the distance between the cone centers was converted in real µm on the retina by multiplying them with the appropriate scale factor of the image, determined by the size of the sample window's side.Conversion from degrees was performed according to the model from [19], with one degree of visual angle equal to 288 µm on the retina.Cone spacing values are compatible with Wyszecki and Styles [65], with the exception of the data from [30] exhibiting lower values, probably due to post mortem shrinkage.Retinas 6, 7-A and 7-B have been cropped during analysis because they didn't fully fill the sampling window, and would have included uncharacterized areas.

Analysis of point process
In this section, we briefly introduce basic notions from Stochastic Point Processes [34].A point process S is a stochastic generating point in a given domain Ω (here, [0, 1) s ).We denote by P n := {x (1) , x (2) and isotropic if any rotation or translation of S have the same statistical properties.We also define the density of a point set as the average number of samples inside a region B of volume V B around a sample x.
This density is constant for isotropic and stationary point processes.A sampler generating sets with a non constant density is sometimes called a non-uniform sampler.To characterize isotropic stationary point processes, the Pair Correlation Function (PCF) is a widely used tool.Such function is a characterization of the distribution of pair distances of a point process.Oztireli [48] devised a simplified estimator for this measure in the particular case of isotropic and stationary point processes.The PCF of a pointset P n in the unit domain [0, 1) s is given by where d(x (i) , x ( j) ) is a distance measure between x (i) and x ( j) .The factor k σ is used to smooth out the function.Oztireli relies on this smoothing to assume ergodicity for all sets.He uses the Gaussian function as a smoothing kernel, but one could use a box kernel or a triangle kernel instead.To compute a PCF, we use this estimator with 3 parameters, the minimal r, r min , the maximal r, r max and the smoothing value σ .Those values are usually chosen empirically.Note that as the number of samples increase, the distances between samples will be very different for similar distributions.To alleviate this, we normalize the distances in our estimations using the maximal possible radius for n samples ( [29], Eq (5)).
In Figure 1, we illustrate how the PCF of several point processes captures the spectral content of the point distribution: a pure uniform sampling, Greean-Noise and Pink-Noise samplers obtained using [3], a jittered sampler (for N samples, subdivision of the domain into regular √ N × √ N square tile and a uniform random sample is drawn in each tile), a Poisson-Disk sampler [7] and a Blue-noise sampler (BNOT) [31].

Results and discussion
Regularity index, or conformity ratio is a quantitative method used for assessing spatial regularity of photoreceptor distributions [61,24,10].A k-d tree structure has been used to find the nearest neighbor for each point, the euclidean distance was calculated for each pair found this way and all the results are classified in histograms.Each distribution of nearest neighbours can be described by a normal Gaussian distribution described by the equation where µ is the mean of the distribution and σ the standard deviation of the measurements.The regularity index is expressed by the ratio of the mean µ by the standard deviation σ .This index is reported to be 1.9 for a full random sampling and the more regular the arrangement, the higher the value, usually 3-8 for retinal mosaics.Regularity indexes for retinal data are shown in Table 1.In contrast with previous claims, our calculated indexes range from 8 to 12.In the lower bound there is data obtained from [13], which instead of a retinal image shows the marked locations of the inner segments of photoreceptors; meanwhile in the upper bound, close to 12, most of the data is from foveal centers in [30], with the exception of retina 8-G, where the different sizes of the photoreceptor profiles reflect different levels of sectioning through the inner segments.
The indexes for data generated with Green noise, Pink noise and BNOT samplers are presented in the same table.As expected, the indexes for Green and Pink noise are assimilable to those of a full random sampling, in fact they are even lower, averaging 1.3 and 1.4 respectively; meanwhile, for the BNOT data, the indexes values are much higher, more than the double of the highest values for retinal RIs.It is not very surprising that, thanks to the the uniformity optimization of BNOT, the indexes are this high; but still very far from the infinite RI of regular lattices.Given the fact that fully regular hexagonal or square patterns are proven to have poor sampling properties and therefore not suitable for simulating cones distribution, in the scope of this paper a higher RI indicates that BNOT is better at generating point processes than the other analyzed point processes.A more recent and reliable method for assessing the goodness of these processes is the previously mentioned Pair Correlation Function.In Table 2, we present the l ∞ distance, between our generated point sets and the measured PCF.From two PCFs ρ and ρ 2 , we denote their l ∞ distance as the maximal distance between the two functions: where r is a given radius.We rely on this measure as it was already used in [48] to compare PCFs, two distributions can be considered the same if this difference is under 0.1.The closest results are from comparison with BNOT and Dart Throwing samplers, moreover, the higher the measured RI for the retinal distribution of photoreceptors, the lower the distance from BNOT PCF.The opposite happens when comparing with Dart throwing algorithm, the closer to the reported RI of 8, the lower the l ∞ distance.This evidences that not only the indexes are actually higher than the ones previously measured, but also that the most effective method to simulate these distributions comes from Blue-noise samplers.

Conclusions
Blue noise sampling can describe features of a human retinal cone distribution with a certain degree of similarity to the available data and can be efficiently used for modeling local patches of retina.We hope this work can be useful to understand how spatial distribution affects the sampling of a retinal image, or the mechanisms underlying the develop- ment of this singular distribution of neuron cells and the implications it has on human vision.Given the nature of bluenoise algorithms, it should be possible to develop an adaptive sampling model that spans the whole retina.However, there would be issues in validating the cone sampling, since imaging of the whole retina is difficult to obtain and analyze.All validation in fact should also be local.Future works will explore the possibility of applying a smooth sampling across the retina to obtain an adaptive sampling, given the PCF and spectra of local patches, the patches can be reproduced [69] and correlated with a heat map that represents interpolation in space [52].

Figure 1
Figure 1 PCF of various 2-D samplers.Fist row from left to right: Realizations of 1024 samples from a uniform (a), a Green-Noise sampler (b), a Pink-Noise sampler (c), a jittered (d), a Poisson-disk (e) and a Blue-Noise sampler ( f ).The second row shows the Fourier specrtum (power spectrum) of each sampler ((g) − (l), spectrum computed on 4096 samples).The PCFs capture the spectral content of each sampler as shown in (m).

Figure 2
Figure 2 From left to right: The picture of the patch of retina, the point samples extracted from the cones' location, Nearest neighbor analysis with mean and standard deviation, Pair Correlation Function.a-d.Images from Scoles et al. [55] e. Image from Roorda & Williams [51].

Figure 3
Figure3From left to right: The picture of the patch of retina, the point samples extracted from the cones' location, Nearest neighbor analysis with mean and standard deviation, Pair Correlation Function.Images from Curcio et al.[15]

Figure 4 Figure 5
Figure 4 From left to right: The picture of the patch of retina, the point samples extracted from the cones' location, Nearest neighbor analysis with mean and standard deviation, Pair Correlation Function.a-d.Images from Jonas et al. [35] e. Image from Curcio et al.[13]

Figure 6
Figure6From left to right: The picture of the patch of retina, the point samples extracted from the cones' location, Nearest neighbor analysis with mean and standard deviation, Pair Correlation Function.Images from Gao & Hollyfield[30] realization of a point process with n samples.A point process S is stationary if it is invariant by translation, and isotropic if it is invariant by rotation.More formally, if we assume P a probability measure, S is stationary if ∀x ∈ R s