The dataset is composed of 492 images in total and can be divided with the following distribution: 452 images were acquired during daytime and with uncontrolled illumination, and 40 images were captured at night with a torchlight as a source of illumination. All images were registered in RGB RAW format and their size is 3936 × 2624 pixels. Table I is presenting the different types of snow composing the dataset and the dates of acquisition for images. Figure 1 is displaying some examples of images composing the dataset.

Acquisitions were performed between February 2022 and April 2022 in Norway, where our department is located. A major aspect of this dataset is that all images used were taken outside in uncontrolled conditions. It means the snow was left untouched with appropriate temperature and pressure conditions, and it was not changed due to a move in a cold room inside. Snow metamorphism [1] is a phenomenon causing snow to change based on external and mechanical properties. Thus, having kept snow in its natural state is a strength of this dataset. Acquisitions were made with a Nikon D610 DSLR camera, producing RGB images. Combined with this camera was a Sigma DG Macro lens with focal lengths varying from 28 mm to 300 mm and aperture between f/3.5 and f/6.3. The camera was mounted on a tripod and was oriented to observe the scene such as the two configurations shown in Figure 2.

Illumination and viewing conditions also vary for images within the dataset. 310 images were acquired with an elevation angle of 30∘ for the camera, and 182 images were captured with the camera facing orthogonally to the scene (i.e. an elevation angle of 0∘). Setups used could be similar to a Reflectance Transformation Imaging (RTI) setup. Here, RTI was not possible due to lack of tools to go outside in the snow to perform such acquisitions. Also, the setup used offered more flexibility to test the algorithm proposed by the work of Ferrero et al [16]. For 452 images, the sun was the light source used as images were captured during daytime. All day images were acquired under a clear sky without clouds, meaning the sun had direct illumination on the scene. Images were taken at different exposure times to cope with strong direct illumination of some setups to avoid having oversaturated images.

In the study of sparkle, azimuth angles for incidence of light sources are commonly used as parameters. For this dataset, it was possible to estimate those angles of incidence by looking at shadows cast on some of the images, or by taking pictures of the setup with the light source visible. Hence, most angles of incidence are estimated and not precisely measured, so their uncertainty could be high. We moderate this choice by conducting a statistical study on a large dataset of images to smooth these potential uncertainties. 3D lidar scans were made during the image acquisition process. And some images contain shadows cast on them, and for most of the scenes, global pictures were taken to ensure the presence of the sun in them. As a consequence, it is possible to have estimates of azimuth angles of the sun for those images. Furthermore, since those images were taken from locations close to each other, and with the knowledge of 3D point clouds, all images have been replaced in the same reference system so that azimuth angles can be used for analysis and are not tied to the camera system. Regarding the 40 images captured at night, the azimuth angle of the torchlight was controlled and is known. Those images are separated from the others for the analysis as they were not acquired in the same geographic location, and then it is not possible to replace them in the same reference system. Figure 3 shows the distribution of angles put in the same reference system across the dataset of images.

Preprocessing of images is performed according to Ferrero et al. method [16]. It consists of capturing an object’s image with known properties such as reflectance or luminance, then use that known information to compute a ratio in order to estimate the unknown reflectance or luminance of the rest of the scene or other objects. All of that is possible assuming the linearity of the camera, which is made in the case of this study as well. The most commonly used objects for these calibrations are white calibration tiles, which are considered Lambertian and with known reflective properties.

The case of snow is more challenging. In the visible range, snow is commonly white, and its reflective properties are higher than the traditional calibration tiles manufactured. The reason is mainly due to scattering and subsurface scattering phenomena occurring near the snow surface, so that the amount of light reflected is large. Since snow is the main target of this statistical study, the preprocessing needs to be refined because the difference between the sparkle (high pixel values or saturated pixels) and the snow (white, so high pixel values) is too narrow. In order to perform a white balance on the images, the Gray World assumption [18] is performed. Images from our dataset are RGB images. Then, the green channel is taken as the reference and the ratios computed are to compute the assumption. Then, those 3 channels R¯,G,B¯ are averaged following Eq. (2) to obtain an image of intensities I that is used to compute sparkle algorithms later described.

Algorithms for the detection of sparkle and its study were developed by Ferrero et al. along several articles [14–17]. These methods were originally designed for the use of goniometric measuring tools to ensure precise values for azimuth angles and more control over illumination and viewing conditions, but can be applied to the study case we built (with higher uncertainties). The requirements are: digital images with apparent sparkle coming from the material studied, information on azimuth angles of incident light and reference for the calibration. The following description of steps of the algorithm is coming from Ferrero et al. [16].

The first step is to calibrate images obtained to get luminance factors β for each pixel of the images. As mentioned previously, the case of snow is slightly more delicate than using traditional white calibration tiles. Therefore, after applying the Gray World assumption, luminance factors are computed for each image of the dataset (as referred in [16]). Each image of the collection is large with size of 3936 × 2624 pixels. Then, in an attempt to reduce computational time, smaller areas are selected from original images. Patches were selected from the center of the image to ensure a maximum of focused snow on patches and no artefacts visible that are not snow related (presence of a colour checker as seen in Fig. 1). A parameter shw controls the half-width size of the patches and can be modified to reduce or expand the area for the sparkle algorithm. Once the patches are selected, a procedure to detect sparkles spots is applied and follows those steps:

The choice of the size of elementary area lhw is to be discussed as it cannot be chosen small for physical meaning, but not too large as well. Otherwise, there is a risk of avoiding some relevant sparkle spots if the distance between two consecutive spots is smaller than lhw. An impact study on the choice of shw and lhw has been conducted and is presented in a later part. Results of this procedure gives luminance factors of sparkle spots and luminance factors for backgrounds of elementary areas.

Ferrero et al. introduced three indicators to quantify sparkle distribution of metallic samples [15, 17], and these same quantities can be applied in the case of snow. The first indicator is an illumination contrast Csp which can be computed by illuminance or luminance factors. The illumination contrast of a sparkle spot can be defined for an elementary area by with βsp and βbg luminance factors previously computed. The second indicator is called ensemble contrast of sparkle spots CE and is defined as the median value of all contrasts of sparkle spots that are above a threshold value Cth. For this study case of snow, this quantity Cth is evaluated for each patch as the middle of the total range of contrast value computed in the image and follows the definition given by Eq. (4)

The third and last indicator introduced is the density of sparkle spots dsp which represents the number of sparkle spots in the area considered with contrast values Csp higher than the threshold value Cth. This value is given in mm−2. As one could expect, this value is linked to the choice of shw so it is important to choose an area of study large enough to provide a stable statistic for this metric. Among those three indicators presented, only two of them, being the ensemble contrast of sparkle spots CE and the density of sparkle spots dsp, are used as metrics in this study. The contrast of sparkle spots Csp is indirectly used in the two metrics cited previously.

To adapt the algorithm to the case of this snow dataset, two parameters are introduced: shw and lhw. In their articles, Ferrero et al. suggested choosing a patch size large enough to ensure materials for the statistics, and the choice of the elementary area lhw was decided to reproduce the circular area visible by a human observer at a distance of 40 cm from the scene in their case. However, they do not provide the values they actually used. Therefore, several values of shw and lhw are chosen in a given range, and the algorithm is run on part of the dataset to test the influence of those parameters on the resulted contrast and density values.

The patch half-width size is chosen among the following values: 80, 100, 120, 150, and 200 pixels. The range for the elementary area is: 11, 13, 15, 17, and 19 pixels. The algorithm was run 25 times to cover all scenarios on a small part of the dataset to avoid long computation time. Only two scenes from the dataset are considered for showing the results displayed in Figure 4 for the first scene and in Figure 5 for the second scene. Units of the parameters are given in pixels to have a common unit for all images of the dataset. Although the pixel unit could be linked to the resolution of the camera, acquisition setups for all images of the dataset are different. As mentioned, the dataset is composed of images taken with different setups of acquisition. Then, the first images (chronologically speaking) were taken with a setup which has evolved. We do not have the required information for those images to compute the ratio that would link the pixel size to an international unit size such as the meter (or centimeter). To avoid having inhomogeneous data, we therefore chose to not include it here, but it is considered for future works.

From Figs. 4(b) and 5(b), one can notice a variation of density values with the choice of lhw as the linear regression plots are parallel but not superimposed. The smaller this lhw value gets, the larger the density, as it is possible to detect more sparkle spots and have fewer overlaps. However, if lhw gets too small, it does not hold a physical meaning as pointed by Ferrero et al. when they chose it. Similarly, by observing results from Figs. 4(c) and 5(c), contrast values are impacted by the choice of the size of the patch observed shw. Expanding the size of the patch opens the possibility to detect other sparkle spots and so increases the contrast values as estimated by the algorithm. However, large values of shw influence the computational time of the algorithm. As a consequence, a trade-off is made in the choice of the size of the elementary area. Ultimately, the values of shw = 150 pixels and lhw = 15 pixels are chosen and fixed for the rest of the analysis for the dataset.

This part tackles the variation in the estimated contrast and density distributions depending on the illumination conditions. One aspect to remind is that day illumination is symbolizing the sunlight, i.e. there was no control over the quantity of light nor the incident angle on the scene. In this study, all images were taken under direct illumination. As mentioned by Kirchner et al. [19], being under direct (no clouds) or diffuse (cloudy overcast sky) illumination highly influences the type of sparkle seen on materials. Azimuth angles for incident light are estimated as explained previously. Regarding night illumination, images were acquired under a torchlight powerful enough to produce sparkle on the snow surface. The algorithm is run on both subsets (day and night) of the dataset to provide contrast and density distributions of sparkle spots.

Figure 6(a) and (b) represent the scatter between the density of sparkle spots (x-axis) and the contrast of sparkle spots (y-axis) respectively for day illumination and night illumination. One thing to note is that the scales for density values are different. The maximum density under daylight is 2.5 mm−2 while the maximum density under torchlight at night is close to 0.02 mm−2 so a factor of 100 between them. Even though the gap is important, it is relevant to notice that the shape of the scatter clouds looks similar in both cases, and could either be an inverse proportionality function or a decreasing exponential. Our goal is not to estimate this correlation. However, one can note that high contrast values are mostly achieved with small values of densities, while low contrast values are spread in the range of density values. The impact of illumination on density values is massive. High density means there are more sparkle events happening on the surface considered. Under sunlight, it is less surprising to see those high density values. Even though the daylight setups are uncontrolled, the behaviour between contrast and density values remains stable between both illumination modes considered. For the rest of the analysis, the night subset is excluded, as it does not hold a lot of images and contains only one type of snow.

In the description of the acquisition of data, two configurations for the position of the camera are mentioned, such as shown in Fig. 2. Therefore, a small study is conducted by separating images taken in the two configurations to check the influence of the elevation angle on the result of contrast and density for sparkle spots.

Figure 7 shows the results obtained for those subsets. For both, the trend for a decreasing exponential or an inverse proportionality function can be observed, even though the scales for contrast and density values are not similar. Since results between the two subsets are quite similar, the distinction on the elevation of the camera is not maintained for the rest of the study. However, the elevation information is to be considered if an experiment is designed to try to identify snow grain shapes.

Table I is already presenting the different types of snow composing the dataset and the repartition between images. The reason why the type of snow is important in this study is the link between the type of snow, the age of snow and the grain structure. The age of snow is directly linked to the morphology of the snow grain [1]. When a snow grain is freshly fallen, its size is rather small. Over time and with increasing temperature, it expands and therefore both his size and shape are evolving. Therefore, the older the snow, the bigger the grains with less complex shapes. The shape and size of snow grains should have an interaction with how sparkle events are emitted. Therefore, the aim of this part is to see whether the sparkle spots seen on images are varying with the type of snow.

First, all images used to obtain the results in this part are images taken from the day illumination subset. As stated in Table I, there is a total of 452 images used for the results displayed in Figures 8 and 9. Even though the fresh snow label is in Fig. 8, we chose to add Fig. 9 for a better readability due to smaller scales. Then, as they are now, an accurate description can be made related to the type of snow.

From Fig. 8, the difference made between dense snow (in black) and old snow (in blue) is related to the melting process. Dense snow is an accumulation of fallen snow grains that could have been accumulated for a long time, hence they could be qualified as old. However, ambient temperature is cold enough to maintain the current state of snow grains. Most of the time, those grains get refrozen due to colder temperature and the wind. For the case of old snow, temperatures are getting above zero degree and then grains are expanding and the transformation from solid to liquid starts. As it can be observed in Fig. 8 with black dots, dense snow tends to produce sparkle with high densities while maintaining contrast values which could be qualified as average. For dense snow as defined here, there are more sparkle spots visible and observable on snow surface even though they remain homogeneous. In the context of old snow, the observation seems to be opposite, as shown in Fig. 8 with blue stars. As opposed to dense snow, older snow whose melting process has started tends to produce sparkle spots with high contrast values and low densities. So, sparkle events may be more visible and shinier on old snow, but they are less likely to happen, or they would be less packed together. Regarding the fresh snow, as illustrated by Fig. 9, it does not involve high contrast values nor density values. Actually, fresh snow is described as snow grains falling on the ground and accumulating quickly. The structure of fresh layer of snow on the surface is not really well-formed, and so fresh snow is less likely to generate sparkles on its surface. However, this conclusion needs to be slightly balanced by the small number of images available, only 47 according to Table I.

As a consequence, it is more difficult to see sparkle spots which correlates with low values both for contrast and density: sparkle events are less likely to happen on fresh snow as the snow grains are too small (under 50 μm) and not well-formed. As opposed to fresh snow, dense snow (with snow grain of size around 500 μm and 1 mm) produces a lot of visible sparkles spots and old snow (snow grains above 4–5 mm) generates less sparkle events, but they are more intense.

Some images from the dataset capture the same scene but observed from two different positions. An example of such a scene is displayed in Figure 10.

Then, all images from this particular scene are gathered and studied for their contrast and density of sparkle spots. Since all images are coming from the dense snow subset, performing the analysis on the variation of illumination and viewing conditions on a similar scene may provide some information related to the geometry and the shape of snow grains.

Figures 11 and 12 display the results obtained for the contrast and density of sparkle spots for this study case. These results are presented differently from the previous ones. Here, it is more interesting to focus on the variations of contrast and density due to different viewing conditions. The camera is looking at the same scene under two perspectives. Hence, snow grains remain constant in the scene, but their effect vary with the viewing conditions. Both radial plots in Figs. 11 and 12 are results correlated to the estimated azimuth angle of the sun on images. On the left image of Fig. 10 (chosen as a reference for the reference system), one can notice shadows cast from the tripod of the camera. Using these shadows, the angle is estimated to 30∘ and then 210∘ for its counterpart. The angles displayed on the radial plots are the relative azimuth angle between the camera and the light source.

Focusing on the radial scatter plots, the interesting aspect to notice is that contrast and density values do not vary much even though sparkle spots are observed from a different viewing angle. This statement is also supported by looking to the shapes of contrast and density distributions. Even though the number of images used (66) is not large enough to make statistical estimations, the shapes of unimodal distributions can be distinguished. More accurately, distributions seem to have a centered value and other values dispatched around it. It can be interpreted as it follows: contrast and density values accumulate towards one average value in a range controlled by a small standard deviation.

A reason why such results can be observed could be correlated to the randomness of the distribution of snow grains onto the snow sample considered. Due to this randomness, potential effects due to geometry of snow grains would be averaged and therefore have less impact. However, if we assume the distribution could be estimated, another reason would then be the shape of the snow grains. In fact, it could even give information on the shape. If one assumes a spherical shape for a grain, then such a grain would reflect light the same way in all possible directions. Then, regardless of the viewing conditions, sparkle spots would have a small variation and could remain stable while the camera rotates around the scene. Due to the lack of number of viewing positions for this scene, it is impossible to conclude accurately on the type of the shape. However, it opens possibilities for designing an acquisition protocol to do so. By selecting one single snow scene and fixing the illumination conditions, several captures can be performed at various positions around the scene. Then, similar radial scatter plots can be obtained, and by observing potential symmetries, one might conclude on the nature of the grain shape.