In the late 1990s I did some work trying to understand what gives rise to the sparkly look of bright lights under various conditions. When we look at a bright light, the image we see is created from light that gets transmitted directly, plus a lot of extra light that bounces around in the eye. The end result is a point spread function that is quite broad. We don’t notice this broadness unless the light source is very bright, where indirect illumination can wash out surrounding parts of the visual image. Here is the writeup:
Light entering the eye is spread out and refracted by the optical structures of the lens, cornea and eye body. This is the source of lens flare that can be visible on bright light sources. Here I describe my simulation results for this effect.
The model for the transfer function (point spread function) comes from [1]. In this model, the point spread function consists of a sharp central peak, accounting for the directly transmitted light portion, the ciliary corona, corona rays, and lenticular halo caused by the lens, and bloom caused by light scattering from other structures in the eye. In the lenticular halo, the light from the light source is broken up into its spectral components as different wavelengths are refracted through different angles. In addition, the point spread function is dependent on pupil diameter.
The point spread function is shown above for 568nm light and a 4.9mm pupil diameter (scotopic conditions.) The side peaks are the lenticular halo, and these diminish for smaller pupil diameters.
This function has four radial components:
where L is the wavelength of the light in nm, and r is the radius in degrees.
The final point spread function is
.
The constants change with the pupil diameter and for photopic and scotopic extremes are:

w0 
w1 
w2 
w3 
Scotopic (4.9mm pupil) 
0.282 
0.478 
0.207 
0.033 
Photopic (1.9mm pupil) 
0.384 
0.478 
0.138 
0.000 
In between these pupil diameters, I linearly interpolate the constants.
In addition to the rotationally symmetric point spread function, there are also radial corona rays that are generated by mutliplying by a uniform deviate in the range 0 to 1 as a function of orientation. A different set of random numbers is used for the rays of the lenticular halo. For an equal energy white light source, the completed scotopic kernel produces the image below (with soft saturation provided by a tanh function.)
If the source image is defined in terms of a number of spectral bands then each band can be convolved with the corresponding point spread function to obtain the wavelength mix at the retina. This wavelength mix can be converted into CIE XYZ components by integration over the CIE sensitivity curves for further processing (e.g. display in RGB space.)
Here, I take a hybrid approach in which a nonlight source, i.e. lowlight scene, input image in XYZ space is used along side an additional arbitrary sized set of input images that contain the intensity of various bright light sources. The colors of the light sources are defined beforehand in CIE xy space, and their spectra are approximated by a weighted combination of equal energy white light, 440nm blue light, and another wavelength that can be freely chosen. This combination can be used to make a candidate spectrum for a light of (almost any) arbitrary chromaticity. Five convolution kernels are then calculated to transform the light source intensity map into a resulting image in XYZ space. (Three kernels are used for white light, one for 440nm and one for the additional wavelength.) The nonlight source input image is convolved with the point spread function, but without the lenticular halo. The resulting lightsource and nonlightsource images are added.
The kernels are generated at the same size and resolution as the input image, with the appropriate scaling to convert between pixels and degrees of visual angle. All convolutions are done using double precision Fourier transforms. To generate the three kernels for white light, it is necessary to integrate the spatial kernels over many wavelengths weighted by the CIE XYZ curves. This can be done in 1D as the lenticular halo kernel f3(r) is isotropic prior to the addition of its corona rays.
To display the results of this simulation work, I generate a number of views of lights of different wavelengths.
These two views are for 440nm and 500nm lights. Note how the radius of the lenticular halo increases with wavelength. Below are a number of star fields, generated with a collection of white, 580, 600 and 630nm light sources. The white light sources have rainbow colored lenticular halos due to their equal energy spectrum.
These star fields are reminiscent of real star fields from the Hubble space telescope (which of course has a different lens point spread function than the eye, and is rather “plus”shaped).
Below is a logo at four different intensities, looking like an illuminated sign. The lights, which were point sources, could be simulated easily on a GPU because it just requires a finite number of stereotyped lens flare texture stamps to be drawn, but here it is necessary to do the full convolution with the source image, and so needs to make use of the FFT to be reasonably fast.
[1] Physically based glare effects for digital images. G. Spencer, P. Shirley, K. Zimmerman, and D. P. Greenberg. Proceedings of ACM SIGGRAPH 1995:325334.