1 Introduction
Photometric stereo techniques aim at acquiring both the shape and the reflectance of a scene. To this end, multiple images are acquired under the same viewing angle but varying lighting, and a physicsbased image formation model is inverted. However, the classic way to solve this inverse problem requires lighting to be highly controlled, which restricts practical photometric stereo applications to laboratory setups where careful calibration of lighting must be carried out prior to any 3Dreconstruction.
…  
Input 
… 
Output 

Reflectance  Shape 
The objective of this research work is to simplify the overall photometric stereo pipeline, by providing an efficient solution to uncalibrated photometric stereo under general lighting, as illustrated in Figure 1.^{1}^{1}1We plan to make the code publicly available. In comparison with existing efforts in the same direction, the proposed one has the following advantages:

The joint estimation of shape, reflectance and general lighting is formulated as an endtoend, mathematically transparent variational problem;

A real 3Dsurface represented as a depth map is recovered, rather than possibly nonintegrable normals;

It is robust, due to the use of Cauchy’s robust Mestimator and HuberTV albedo regularization;

It is computationally efficient, thanks to a tailored lagged block coordinate descent scheme initialized using a simple balloonlike shape.
After reviewing related works in Section 2, we discuss in Section 3 the spherical harmonicsbased image formation model considered in this work. It can be inverted using the variational approach to the joint recovery of depth, albedo and lighting which is discussed in Section 4. A dedicated numerical solution is then introduced in Section 5 and empirically evaluated in Section 6. Section 7 eventually draws the conclusion of this research.
2 Related Work
3Dmodels of scenes are essential in many applications such as visual inspection [Farooq2005] or computeraided surgery using augmented reality [Collins2012]. A 3Dmodel consists of geometric (position, orientation, etc.) and photometric (color, texture, etc.) properties. Given a set of photographies, the aim of 3D scanning is to invert the image formation process in order to recover these geometric and photometric properties of the observed scene. This notion thus includes both those of 3Dreconstruction (geometry) and of reflectance estimation (photometry).
Many approaches to the problem of 3Dreconstruction from photographies have been studied, and they are grouped under the generic naming “shapefromX”, where X stands for the clue which is being used (shadows [Shafer1983], contours [Brady1984], texture [Witkin1981], template [Bartoli2015], structured light [Geng2011], motion [Moons2008], focus [Nayar1994], silhouettes [Hernandez2004], etc.). Geometric shapefromX techniques are based on the identification and analysis of feature point or areas in the image. In contrast, photometric techniques build upon the analysis of the quantity of light received by each photosite of the camera’s sensor. Among photometric techniques, shapefromshading
is probably the most famous one. This technique, developed in the 70s by Horn
et al. [Horn1970], consists in 3Dreconstruction from a single image of a shaded scene. It is a classic illposed inverse problem whose numerical solving usually requires the surface’s reflectance to be known [Durou2008]. In order both to limit the ambiguities of shapefromshading and to allow for automatic reflectance estimation, it has been suggested to consider not just one image of the scene, but several ones acquired from the same viewing angle but under varying lighting. This variant, which was introduced in the late 70s by Woodham [Woodham1978], is known as photometric stereo.Among the various shapefromX techniques mentioned above, photometric stereo is the only 3Dscanning technique i.e., the only one which is able to achieve both 3Dreconstruction and reflectance estimation. However, early photometric approaches strongly rely on the control of lighting. The latter is usually assumed for simplicity to be directional, although the case of nearby point light sources has recently regained some attention [logothetis2017semi, mecca2014near]. More importantly, lighting is assumed to be calibrated. Indeed, the uncalibrated problem is illposed: the underlying normal map can be estimated only up to a linear ambiguity [Hayakawa1994], which reduces to a generalized basrelief one if integrability is enforced [Belhumeur1999]. To resolve the latter ambiguity, some prior on the scene’s surface or geometry must be introduced, see [Shi2018] for a recent survey. A natural way to enforce integrability consists in following a differential approach to photometric stereo [Chand13, mecca14] i.e., directly estimate the 3Dsurface as a depth map instead of first estimating the surface normals and then integrating them. Such a differential approach to photometric stereo can be coupled with variational methods in order to iteratively refine depth, reflectance and lighting in a robust manner [CVPR2017]. In addition to the theoretical interest of enforcing integrability in order to limit ambiguities, differential approaches to photometric stereo have the advantage of bypassing the problem of integrating the estimated normal field, which is by itself a nontrivial problem [Queau_Survey]. Besides, any error in the estimated normal field might propagate during integration, and thus robustness to specularities or shadows must be enforced during normal estimation, see again [Shi2018] for some discussion.
All the research works mentioned in the previous paragraph assume that lighting is induced by a single light source. Nevertheless, many studies rather considered the case of more general illumination conditions, which finds a natural application in outdoor conditions [Sato1995]
. For instance, the apparent motion of the sun within a day induces changes in the illumination direction which, in theory, allow photometric stereobased 3Dreconstruction. However, this apparent motion is close to being planar, and thus the rank of the set of illumination vectors is equal or close to
[Shen2014] (see also [Holdgeoffroy2015] for additional discussion on the stability of singleday photometric stereo). This situation is thus similar to the twoimage case, which is known to be illposed since the early 90s [Kozera1993b, Onn1990, Yang1992], although it is still an active research area [Kozera2018]. In order to limit the instabilities due to this issue, one possibility is to consider images acquired over many seasons as in [Abrams2012, Ackermann2012], or to resort to deep neural networks
[Holdgeoffroy2018]. Another one is to consider a nondirectional illumination model to represent natural illumination, as for instance in [Jung2019]. Modeling natural illumination is a promising track, since such a model would not be restricted to sunny days, and images acquired under cloudy days are known to yield more accurate 3Dreconstructions [Holdgeoffroy2015].However, the previous approaches to photometric stereo under natural illumination assume calibrated lighting, where calibration is deduced from time and GPS coordinates or from a calibration target. The case of both general and uncalibrated lighting is much more challenging and has been fewly explored, apart from studies restricted to sparse 3Dreconstructions [Shen2009] or relying on the prior knowledge of a rough geometry [Ackermann2012b, Peng2017, Shi2014]. Uncalibrated photometric stereo under natural illumination has been revisited recently in [Mo2018], using a spatiallyvarying equivalent directional lighting model. However, results were limited to the recovery of possibly nonintegrable surface normals. Instead, the method which we propose in the present paper directly recovers the underlying surface represented as a depth map. Following the seminal work of Basri and Jacobs [Basri2007], it considers the spherical harmonics representation of general lighting in lieu of the equivalent directional approximation, as discussed in the next section.
3 Image Formation Model
In photometric stereo (PS), we are given a number of observations , each representing a multichannel image (i.e. ) over a masked pixel domain . Assuming that the object being pictured is Lambertian, the surface’s reflectance is represented by the albedo , and the general image formation model is as follows, for all , , and :
(1) 
Here is the unit sphere in , represents the channelwise intensity of the incident light, and and are the channelwise albedos and the unitlength surface normals, respectively, at the surface point conjugate to pixel . The operation in (1) encodes selfshadows. The overall integral collects elementary luminance contributions arising from all incident lighting directions . In the setup of uncalibrated PS, the quantities , , in addition to , are unknown.
Equivalent directional lighting [holdgeoffroy15] approximates (1) via
(2)  
where represents the mean lighting over the visible hemisphere at . The field is spatially variant but can be approximated by directional lighting over small local patches. Over each patch, one is thus faced with the ambiguities of directional uncalibrated PS [Hayakawa1994]. Stateoftheart patchwise methods [Mo2018] first solve this problem over each patch, then connect the patches to form a complete normal field up to rotation, and eventually estimate the rotation which best satisfies the integrability constraint. Errors may however get propagated during the sequence, resulting in a possibly nonintegrable normal field.
Instead of such an equivalent directional lighting model, we rather consider a spherical harmonic approximation (SHA) of general lighting [Basri2003, Basri2007]. By defining the halfcosine kernel as
(3) 
we can view (1) as an analog of a convolution:
(4) 
Invoking the FunkHecke theorem, we obtain the following harmonic expansion analogous to Fourier series:
(5) 
Here the spherical harmonics form an orthonormal basis on , and and are the expansion coefficients of and with respect to . Since most energy in the expansion (5) concentrates on loworder terms [Basri2003], we obtain the secondorder SHA by truncating the series up to the first nine terms (i.e., ):
(6) 
The firstorder SHA refers to the truncation up to the first four terms (i.e., ). It is shown in [Basri2003] that, for distant lighting, at least of the resulting irradiance is captured by the firstorder SHA, and by the secondorder SHA (cf. Figure 2 for a visualization).
Model (1)  Model (7) (1st order)  Model (7) (2nd order)  
Plugging (6) and specifics of spherical harmonics [Basri2003] into (4), we finalize our image formation model as:
(7)  
(8) 
Here represents the secondorder harmonic images, and represents the harmonic lighting vector whose entries have absorbed and constant factors of . A key advantage of the SHA (7) over the equivalent directional lighting model (2) lies in the spatial invariance of the lighting vectors , which yields a less illposed inverse problem [Basri2007]. The counterpart is the nonlinear dependency upon the normal components, which we will handle in Section 5 using a tailored numerical solution. In the next section, we build upon the key observations that integrability [Belhumeur1999] and perspective projection [papadhimitri2013new] both largely reduce the ambiguities of uncalibrated PS to derive a variational approach to inverting the SHA (7).
4 Variational Uncalibrated PS
In this section, we shall propose a joint variational model for uncalibrated PS. To this end, let a 3Dframe
be attached to the camera, with the optical center, the axis aligned with the optical axis such that for any 3D point in front of the camera. Further let a 2Dframe be attached to the focal plane which is parallel to the plane and contains the masked pixel domain . Under perspective projection, the surface geometry is modeled as a map given by(9) 
with the depth map and
(10) 
the calibrated camera’s intrinsics matrix. In the following we denote for convenience .
Assuming that is differentiable, the surface normal at point is the unit vector oriented towards the camera such that , which yields the following parameterization of the normal by the depth:
(11)  
(12) 
Note that the dependence of on is linear.
Based on the forward model (7) and the parameterization (11) of normals, we formulate the joint recovery of reflectance, lighting and geometry as the following variational problem:
(13) 
In the first term above, we use Cauchy’s Mestimator to penalize the datafitting discrepancy:
(14) 
It is indeed wellknown that Cauchy’s estimator, being nonconvex, is robust against outliers; see for instance
[CVPR2017] in the context of PS. The scaling parameter is used in all experiments.The second term in (13) represents a Huber totalvariation (TV) regularization on each albedo map , with the Huber loss defined by
(15) 
and being fixed in the experiments. It turns out that the Huber TV imposes desirable smoothness on the albedo maps and in turn improves the joint estimation overall. Eventually, is a weight parameter which balances the datafitting term and the Huber TV one. Its value was empirically set to (see Section 6 for some discussion).
We remark that in (13), geometry is directly optimized in terms of the depth (rather than indirectly in terms of the normal ). This is typically known as a differential approach to PS in the literature [Chand13, logothetis2017semi, mecca14, mecca2014near, CVPR2017], with the benefit of both implicitly ensuring integrability (thus limiting ambiguities) and avoiding possibly errorprone integration of normals into depths as a postprocessing step.
5 Solver and Implementation
To solve the variational problem (13) numerically, we follow a “discretizethenoptimize” approach. There, is replaced by , being the number of pixels inside , which yields discretized vectors . To alleviate notational burden, we sometimes refer to a pixel by its index and sometimes by its position . The spatial gradient is discretized using a forward difference stencil.
We shall apply a lagged block coordinate descent (LBCD) method to find a local minimum of the objective function in (23). Due to the (highly) nonconvex nature of (23), initialization of optimization variables has a strong influence on the final solution. In our implementation, we initialize for all and for all . Moreover, during the first eight iterations we freeze the secondorder spherical harmonics coefficients i.e., we reconstruct using only firstorder spherical harmonic approximation as a warm start.
Based on the observation that most realworld scenes are convex, we initialize the depth as a balloonlike surface, as discussed in the following.
5.1 Depth Initialization
It is readily seen that a trivial constant initialization of the depth yields uniform vertically aligned normals and, hence, zero entries in the initial harmonic images . This would cause nonmeaningful updates on albedos and lighting vectors ; cf. Figure 3 for an illustration.
To solve this issue, we specialize the depth initialization which undergoes two phases:
1. Following [Oswald2012], we generate a balloonlike depth map under orthographic projection.
2. We then convert the orthographic depth to a perspective depth via normal integration [Queau_Survey].
Phase 1 is pursued via the following constrained convex optimization:
(16)  
s.t. 
Intuitively, the model seeks for which has minimal surface area subject to a constant volume . A global minimizer of this model can be efficiently computed by simple projected gradient iterations:
(17)  
(18) 
where and with the spectral norm. The volume constant
is a hyperparameter which is empirically chosen, see Section
6 for discussion.Next, we convert the orthographic depth to a perspective depth . Note that complies with the orthographic projection, under which a 3Dpoint is represented by
(19) 
and the corresponding surface normal to the surface at conjugate to pixel is given by
(20) 
Since is invariant to the projection model, Eq. (11) also implies that
(21) 
where stands for the logperspective depth. This further implies the formula for :
(22) 
which can be integrated to obtain (and hence ). The overall pipeline in Phase 2 is summarized as follows:
Trivial and its result  Our and its result 
5.2 Lagged Block Coordinate Descent
Even with a reasonable initialization, the numerical resolution of Problem (23) remains challenging. Due to the appearances of the spherical harmonic approximation and the Cauchy’s Mestimator , the objective in (23) is highly nonlinear and nonconvex. To tackle these challenges, here we present a lagged block coordinate descent (LBCD) method which performs efficiently in practice.
To derive LBCD, we introduce an auxiliary variable such that . This enables us to rewrite (11) as . Then we formulate the following constrained optimization problem:
(23)  
s.t. 
where is the residual function defined by:
(24) 
Upon initialization, the proposed LBCD proceeds as follows. At iteration , we lag one iteration behind, i.e.,
(25) 
and then sequentially update each of the three blocks (namely , and ). In each resulting subproblem, we solve (lagged) weighted least squares problems as an approximation of the Cauchy loss and/or the Huber loss. This is detailed in the following:

(Update ): We evaluate the residual
(26) and then set up the (lagged) weight factors for both the Cauchy loss and the Huber loss as
(27) (28) The albedos are updated as the solution to the following linear weighted leastsquares problem:
(29) which is carried out by conjugate gradient (CG).

(Update ): The lighting subproblem is similar to the one for albedos, except for absence of the Huber TV term. Upon evaluation of the residual and the weight factor , we update by solving the following linear weighted leastsquares problem via CG:
(30) 
(Update ): The depth subproblem requires additional efforts. With and evaluated after the update, we are faced with the following weighted least squares problem:
(31) where the dependence of on is still nonlinear. Therefore, we further linearize with respect to and arrive at the following update:
(32) where is the Jacobian of the map at . The resulting linearized leastsquares problem is again solved by CG. In our experiments, we additionally incorporate backtracking line search in the update to ensure a monotonic decrease of the energy.
6 Experimental Validation
3D Shapes 


Armadillo  Joyful Yell  Lucy  Thai Statue  Exemplary environment maps  
Albedos 

Bars  Constant  Ebsd  Hippie  Lena  Pattern  Rectcircle  Voronoi  White 
This section is concerned with the evaluation of the proposed nonconvex variational approach to uncalibrated photometric stereo under general lighting.
6.1 Synthetic Experiments
To validate the impact of the initial volume in (16), the tunable hyperparameter , and the number of input images in (13), we consider 36 challenging synthetic datasets. We use four different depth maps (“Joyful Yell” [Bendansie2015], “Lucy” [Levoy2005data], “Armadillo” [Levoy2005data] and “Thai Statue” [Levoy2005data]) and nine different albedo maps and each of those 36 combinations is rendered as described in (1) using different environment maps^{2}^{2}2Environment maps are downloaded from http://www.hdrlabs.com/sibl/archive.html, cf. Figure 4. The resulting 25 RGB images per dataset are used as input, along with the intrinsic camera parameters and a binary mask . A quantitative evaluation on the triplet is carried out on four randomly chosen datasets (Armadillo & White albedo, Joyful Yell & Ebsd albedo, Lucy & Hippie albedo, and Thai Statue & Voronoi albedo), comparing the impact of each value of on the resulting mean angular error (MAE) between ground truth and estimated normals.
First, we validate the choice of the input volume using the initially fixed values of and . As the volume depends on the size of the mask, we consider a linear parametrization and evaluate a range of ratios . Figure 5 (left) indicates that the optimal value of is datasetdependent. For synthetic datasets we always selected this optimal value, yet for realworld data no such evaluation is possible and must be tuned manually. Since the ballooningbased depth initialization can be carried out in realtime (implementation is parallelized in CUDA), the user has an immediate feedback on the initial depth and thus a plausible initial shape is easily drawn. Humans excel at estimating size and shape of objects [Baldwin2016] and realworld experiments will show that a manual choice of can result in appealing geometries.
Next, we evaluate the impact of , cf. Figure 5 (right). As can be seen, the depth estimate seems to deteriorate for too small and too large values of , whereas seems to provide good depth estimates across all albedo maps. Therefore we fix for all our upcoming experimental evaluation.
Unsurprisingly, the MAE is inversely proportional to the number of input images, but runtime increases (linearly) with , cf. Figure 6. We found that represents a good tradeoff between runtime and accuracy, and fix for all our further experiments. Our Matlab implementation needs about – minutes on a computer with an Intel processor.
Having fixed the choice of , we can now evaluate our approach against other stateoftheart methods. We compare our results against those obtained by an uncalibrated photometric stereo approach assuming directional lighting [Favaro2012], and another one assuming general (firstorder spherical harmonics) illumination yet relying on an input shape prior (e.g., from an RGBD sensor) [Peng2017]. As this limiting assumption on the access to a sensorbased depth prior is not always given and to make comparison fair, we input as depth prior to this method the ballooning initialization described in Section 5.1. Furthermore, we compare against another uncalibrated photometric stereo work under natural illumination [Mo2018]^{3}^{3}3Code associated with [Favaro2012] and [Peng2017] can be found online, and the results obtained by [Mo2018] were provided by the authors., which resorts to the equivalent directional lighting instead of spherical harmonics, cf. Section 3. Table 1 shows the median and mean MAEs over all 36 datasets (a more detailed table can be found in the supplementary material). On these datasets, it can be seen that our method quantitatively outperforms the current stateoftheart by a factor of –. This gain is also evaluated qualitatively in Figure 7, which shows a selection of two results.
Approach  [Favaro2012]  [Peng2017]  [Mo2018]  Ours 

Median  27.16  21.14  34.06  9.17 
Mean  34.15  21.18  35.53  10.72 
[Favaro2012]  [Peng2017]  [Mo2018]  Ours  GT  

Joyful Yell & Lena 

MAE:  
Thai Statue & White 

MAE: 
6.2 RealWorld Experiments
For realworld data we use the publicly available dataset of [Haefner2018b]. It offers eight challenging realworld datasets of objects with complex geometry and albedo captured under daylight and a freely moving LED, along with intrinsics matrix and masks . Results are presented in Figure 8. Despite relying on a directional lighting model, the approach of [Favaro2012] produces reasonable results on some datasets (Face1, Ovenmitt or Shirt), but it fails on others. As [Peng2017] assumes a reliable prior on depth in order to perform a photometric refinement, this approach is biased towards its initialization and thus, only when the depth prior is very close to the objects’ rough shape (Ovenmitt, Shirt, Tabletcase, Vase) a meaningful geometry is recovered. The approach of [Mo2018] estimates a possibly nonintegrable normal field only, and it can be seen that after integration the depth map might not be satisfactory. As our approach optimizes over depth directly, such issues are not apparent and we are able to recover finescale geometric details throughout all tests.
[Favaro2012]  [Peng2017]  [Mo2018]  Ours  
Face 1 

Face 2 

Rucksack 

Backpack 

Ovenmitt 

Shirt 

Tabletcase 

Vase 
7 Conclusion
We proposed a variational approach to uncalibrated photometric stereo (PS) under general lighting. Assuming a perspective camera setup, our method jointly estimates shape, reflectance and lighting in a robust manner. The possible nonintegrability of normals is bypassed by the direct estimation of the underlying depth map, and robustness is ensured by resorting to Cauchy’s Mestimator and HuberTV albedo regularization. Although the problem is nonconvex and thus numerically challenging and initializationdependent, we tackled it efficiently through a tailored lagged block coordinate descent algorithm and ballooningbased depth initialization. Over a series of evaluations on synthetic and real data, we demonstrated that our method outperforms existing methods in terms of MAE by a factor of – and provides highly detailed reconstructions even in challenging realworld settings.
In future research, a more automated balloonlike depth initialization is desirable. Exploring the theoretical foundations (uniqueness of a solution) of differential perspective uncalibrated PS under spherical harmonic lighting and analyzing the convergence properties of the proposed numerical scheme constitute two other promising perspectives.
Appendix A Further Details on Synthetic Experiments
To provide further insights on the synthetic experiments (in Section 6.1), we visualize the environment lighting , , used to render each image. Figure 9 shows all environment maps^{4}^{4}4All environment maps were downloaded from http://www.hdrlabs.com/sibl/archive.html.
The impact of each incident lighting , , is illustrated in Figure 10 showing the Joyful Yell with a White () albedo. Thus, color changes in the images are caused by lighting only, as depicted in model (1) and (7) in the main paper.
Table 2 shows the mean angular error (MAE) of each dataset on the stateoftheart approaches [Favaro2012, Mo2018, Peng2017] and our proposed methodology. It can be seen that our approach consistently overcomes [Favaro2012, Mo2018, Peng2017] by a factor of –. Only the Pattern albedo seems to bias the resulting depth negatively, yet even in this case our approach estimates the geometry more faithfully than the current stateoftheart.
Dataset  [Favaro2012]  [Peng2017]  [Mo2018]  Ours  

Shape  Albedo  
Armadillo 
Bars  26.22  27.84  36.91  16.78 
Constant  25.84  26.64  36.87  13.97  
Ebsd  25.34  26.88  27.80  14.26  
Hippie  28.21  27.30  25.82  14.52  
Lena  27.07  27.33  28.36  14.78  
Pattern  45.87  26.82  24.01  19.06  
Rectcircle  26.97  26.71  36.23  14.06  
Voronoi  25.62  26.91  50.70  14.07  
White  26.19  26.64  52.04  14.13  
Joyful Yell 
Bars  21.84  16.26  31.80  8.69 
Constant  23.95  14.93  33.47  5.96  
Ebsd  26.08  15.63  15.91  7.28  
Hippie  28.67  16.23  22.96  7.49  
Lena  21.33  16.33  19.70  9.21  
Pattern  26.07  18.76  26.67  16.97  
Rectcircle  35.27  15.19  52.41  7.34  
Voronoi  22.27  16.42  45.74  6.57  
White  27.12  14.32  33.06  6.20  
Lucy 
Bars  49.13  21.90  36.51  8.16 
Constant  54.98  19.89  36.57  8.71  
Ebsd  62.33  20.81  23.56  9.61  
Hippie  58.61  21.29  32.38  7.87  
Lena  64.01  22.24  30.93  9.56  
Pattern  48.83  22.25  32.68  17.78  
Rectcircle  24.68  20.99  43.13  8.98  
Voronoi  61.53  22.10  48.14  7.59  
White  64.43  19.33  44.76  8.76  
Thai Statue 
Bars  25.53  21.91  66.17  8.55 
Constant  27.20  18.91  38.47  9.58  
Ebsd  27.85  20.22  34.11  9.47  
Hippie  21.91  21.86  30.62  8.83  
Lena  33.53  19.66  34.00  9.19  
Pattern  26.77  22.06  28.81  15.27  
Rectcircle  29.36  19.92  43.86  8.84  
Voronoi  30.65  21.56  36.58  8.69  
White  28.02  18.64  37.31  9.16  
Median  27.16  21.14  34.06  9.17  
Mean  34.15  21.18  35.53  10.72 
Two more qualitative results on synthetic data are shown in Figure 11. While [Favaro2012] gives more meaningful results on Armadillo with Constant albedo, depth deteriorates strongly on Lucy with Hippie albedo. Methods of [Mo2018, Peng2017] both result in rather flattened shapes (cf. Lucy). Most accurate results are achieved using the proposed method where fine geometric details, as well as non flattened depth estimates are shown.
[Favaro2012]  [Peng2017]  [Mo2018]  Ours  GT  

Armadillo & Constant 

MAE:  
Lucy & Hippie 

MAE: 
Appendix B Further Details on RealWorld Results
Supplementary to the realworld experiments (in Section 6.2), Figure 12 shows alternative viewpoints of the realworld results. The estimated albedos, which are mapped onto the surfaces, appear satisfactory.
Face 1  Face 2  Rucksack  Backpack 
Ovenmitt  Shirt  Tabletcase  Vase 
Comments
There are no comments yet.