CC BY-NC-ND
0Contents lists available at ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing
journal homepage: www.elsevier.com/locate/isprsjprs
Effects of illumination differences on photometric stereo shape-and- q albedo-from-shading for precision lunar surface reconstruction ™
Wai Chung Liua, Bo Wua'*, Christian Wohlerb
a Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, Hung Hom, Hong Kong bImage Analysis Group, TU Dortmund, Dortmund 44227, Germany
ARTICLE INFO
Article history:
Received 10 November 2017
Received in revised form 18 December 2017
Accepted 19 December 2017
Keywords:
Surface reconstruction Shape-and-albedo-from-shading Photometric stereo Illumination differences LROC NAC
ABSTRACT
Photoclinometric surface reconstruction techniques such as Shape-from-Shading (SfS) and Shape-and-Albedo-from-Shading (SAfS) retrieve topographic information of a surface on the basis of the reflectance information embedded in the image intensity of each pixel. SfS or SAfS techniques have been utilized to generate pixel-resolution digital elevation models (DEMs) of the Moon and other planetary bodies. Photometric stereo SAfS analyzes images under multiple illumination conditions to improve the robustness of reconstruction. In this case, the directional difference in illumination between the images is likely to affect the quality of the reconstruction result. In this study, we quantitatively investigate the effects of illumination differences on photometric stereo SAfS. Firstly, an algorithm for photometric stereo SAfS is developed, and then, an error model is derived to analyze the relationships between the azimuthal and zenith angles of illumination of the images and the reconstruction qualities. The developed algorithm and error model were verified with high-resolution images collected by the Narrow Angle Camera (NAC) of the Lunar Reconnaissance Orbiter Camera (LROC). Experimental analyses reveal that (1) the resulting error in photometric stereo SAfS depends on both the azimuthal and the zenith angles of illumination as well as the general intensity of the images and (2) the predictions from the proposed error model are consistent with the actual slope errors obtained by photometric stereo SAfS using the LROC NAC images. The proposed error model enriches the theory of photometric stereo SAfS and is of significance for optimized lunar surface reconstruction based on SAfS techniques.
© 2017 The Author(s). Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). This is an open access article under the CC BY-NC-ND license (http://
creativecommons.org/licenses/by-nc-nd/4.0/).
1. Introduction
Photoclinometric surface reconstruction methods such as Shape-from-Shading (SfS) and Shape-and-Albedo-from-Shading (SAfS) estimate the three-dimensional (3D) shape of a surface on the basis of its interactions between an illumination source and a sensor through the image intensity of each pixel. Each pixel contains the combined information of the slope and the other physical properties such as albedo and roughness. SfS or SAfS aims at solving for slope information and convert it into 3D shapes. These techniques have been used for lunar and planetary surface reconstruction (Lohse and Heipke, 2004; Grumpe et al., 2014; Wu et al., 2018) to improve the quality of mapping products such as digital elevation models (DEMs) and are well-known for creating pixel-wise DEMs even in textureless areas where image matching
* Corresponding author at: Department of Land Surveying and Geo-Informatics, Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong. E-mail address: bo.wu@polyu.edu.hk (B. Wu).
is ineffective (Lohse and Heipke, 2004). Lunar and planetary DEMs are critical to most extraterrestrial scientific and engineering applications, and hence, the incorporation of SfS or SAfS methods into conventional surface reconstruction routines such as stereo pho-togrammetry and laser altimetry has gained considerable attention (Kirk et al., 2003b; Grumpe et al., 2014; Wu et al., 2018).
More robust SfS or SAfS surface reconstruction results can be achieved by incorporating images of the same location under different illumination conditions (e.g., images with sunlight coming from different angles). This is referred to as photometric stereo (Woodham, 1980). Multiple illumination directions allow for a more reliable determination of surface slopes and hence, better reconstruction. In practice, the availability of multiple illumination directions is limited depending on the target body and the mission design. For example, images of lunar equatorial areas are usually illuminated either from the east or from the west (Wohler, 2004). In this case, the quality of reconstruction is likely to depend on the azimuthal difference between the illumination sources of the images as it directly implies an abundance of excess information.
https://doi.org/10.1016/j.isprsjprs.2017.12.010
0924-2716/® 2017 The Author(s). Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
However, an explicit and systematic analysis of the relationship between the quality of the SfS or SAfS reconstruction and the illumination difference is not common, which bases quality estimation mainly on experience or the rules of thumb. In this paper, we present an algorithm for photometric stereo SAfS and an error model that describes the relationship between the illumination difference and the SAfS reconstruction quality. Doing so can allow for (1) the estimation of the SAfS slope error on the basis of known information (i.e., illumination conditions and image statistics) prior to any actual process and (2) the estimation of the relative quality between multiple pairs, which can help in choosing the best pair or adjusting the weights among pairs if multiple photometric stereo pairs are used. The rest of this paper is organized as follows: After a literature review of photoclinometric surface reconstruction in Section 2, the photometric stereo SAfS algorithm and the mathematical derivation of the error model are described in detail in Section 3. Then, the error model is verified by using real lunar images from the LROC NAC in Section 4. Finally, the findings are concluded and discussed in Section 5.
2. Related work
Single-image Shape-from-Shading (Sl-SfS) (Horn, 1990) is the most fundamental form of photoclinometric surface reconstruction. Profiles parallel to the illumination direction are referred to as characteristic strips and contain the richest shape information as they depend most strongly on the light source direction (Horn, 1977). In photoclinometry, the 3D shape of a surface is related to the reflected energy (i.e., image intensity) by reflectance models such as Lambert, lunar-Lambert (McEwen, 1991), and Hapke (Hapke, 1981, 1984, 1986, 2002); however, this relationship describes the slope of the surface, which implies an underlying shape, instead of the actual 3D shape of the surface (i.e., 3D points or heights). Kirk (1987) and Horn (1990) proposed algorithms for implementing Sl-SfS in image-based slope estimation to generate 3D shapes (e.g., DEMs) from the estimated slopes. These algorithms have laid the foundation for the use of SfS in planetary mapping and scientific applications. For example, SfS derived DEMs were used for safety assessment of the Mars Exploration Rovers (MER) landing sites (Kirk et al., 2003b) and the candidate Mars Science Laboratory (MSL) landing sites (Beyer and Kirk, 2012). SfS algorithms have also been implemented in the USGS Integrated Software for lmaging Spectrometers (lSlS) (Kirk et al., 2003a). They were also used in scientific studies such as spectral analysis of lunar surface (Wohler et al., 2014) and lunar crater analysis (Salamuniccar et al., 2014). Photoclinometric surface reconstruction performs well in local scale reconstruction but has poor large scale accuracy. Thus, methods have been developed (Grumpe et al., 2014; Wu et al., 2018) by incorporating low-resolution 3D information (e.g., low-resolution DEMs) into Sl-SfS algorithms. These algorithms have successfully improved large scale accuracy while maintaining high-resolution local scale detail of the resulting DEM. DEMs or slope maps from Sl-SfS have a resolution comparable to that of the input images, which is apparently higher than that of mapping products generated using other methods such as photogrammetry and laser altimetry, and can therefore provide better information for scientific and engineering purposes (Beyer and Kirk, 2012).
Simultaneously utilizing images acquired under different illumination conditions is referred to as photometric stereo. The idea of photometric stereo was first proposed by Woodham (1980), who defined it as images acquired under the same viewing condition but different illumination conditions. This condition generally holds for most remotely sensed images as they are usually taken in the nadir or near-nadir viewing geometry at different solar azimuth angles. Woodham (1980) demonstrated with synthetic
examples that photometric stereo can be used for determining robust surface slopes and laid the foundation of photometric stereo photoclinometric reconstruction. Heipke (1992) introduced the use of photometric stereo in topographic mapping, which is referred to as multiple-image Shape-from-Shading (Ml-SfS). This concept was further developed and investigated for planetary mapping and modeling (Lohse and Heipke, 2004; Wohler, 2004; Lohse et al., 2006). Ml-SfS gives plausible results in textureless areas and has demonstrated value in high-resolution surface modeling.
The result of photoclinometric surface reconstruction is affected by uncertainties from various sources. These error sources were critically reviewed byJankowski and Squyres (1991); they include errors from imaging devices (e.g., noise) and the relevant planetary body (e.g., albedo and photometric function). Slope error models were developed on the basis of single-image photoclinometry to relate slope errors from photoclinometry to the addressed error sources. Uncertainties in photometric functions or reflectance models have been studied rigorously, and critical improvements have been made to minimize or model these uncertainties. For example, Hapke (1984) addressed the effects of macroscopic roughness of a surface to its reflectance and proposed a modified version of the Hapke (1981) model with considerations of roughness. Labarre et al. (2017) investigated the roughness addressed by Hapke (1984) and proposed a simplified version of the model. Oren and Nayar (1994) proposed similar corrections related to roughness but on the basis of Lambert reflectance and a statistical estimation of surface roughness. These modified reflectance models have been used in photoclinometric reconstruction in later planetary mapping works such as O'Hara and Barnes (2012). The relationship between various error sources and the quality of pho-toclinometry has been studied in research conducted after Jankowski and Squyres (1991). For example, Greenberg et al. (2011) investigated the relationship between macroscopic roughness and the quality of photoclinometry, and highlighted the importance of error suppression strategies and considerations of uncertainties in photoclinometry. Liu and Wu (2017) addressed the relationship between photometric stereo SAfS and the azi-muthal illumination difference, and suggested that the best combination of a photometric stereo pair would be at 90° (e.g., one image with sunlight coming from the north and the other with sunlight coming from the east) to yield results with minimal error. However, Liu and Wu (2017) did not consider the zenithal illumination differences.
3. Effects of illumination differences on photometric stereo SAfS
3.1. Photometric stereo SAfS algorithm
As shown in Fig. 1, the photometric stereo SAfS routine begins with a photometric stereo pair under given imaging and illumination conditions. At any point (x,y), its surface slope (p,q)\x,y should satisfy the following reflectance constraints:
Gi(p, q)\Xy = G^ = A^ (1)
where for two images i = (1, 2), Gi(p,q)\xy is the reflectance model estimated using the surface slope (p,q) at point (x,y), G,-| is the expected reflectance computed by dividing the image intensity Ii by the intrinsic albedo Axy. lt is insufficient to solve for p, q and Ax,y simultaneously on the basis of these two reflectance constraints. However, Eq. (1) for i = (1, 2) can be combined as follows:
G1(p,q)\xy G1
G2(p; q)\xy G2 1 x y Axy 12 \xy h xy
Fig. 1. Workflow for photometric stereo SAfS.
This is referred to as the quotient constraint; the intrinsic albedo Axy is cancelled out in Eq. (2). This method has been used for reflectance-based reconstruction using images with coplanar light sources (Wohler, 2004) where the slope of only one direction is estimated. We adopt this idea by estimating the slope in a certain direction on the basis of the following quotient constraint:
Gl(s)|x G2(s)|x
Il L.y Ax,y Axy h^y
where s refers to the slope in the desired direction. As photometric stereo pairs are not necessarily coplanar (i.e., a is usually not equal to 180° as illustrated in Fig. 2), the direction for slope estimation (i.e., the principal direction) cannot be strictly determined. Instead, we determined it from Fig. 2, where q' bisects the azimuthal difference a between the two illuminations L1 and L2 such that (b = b2 = §) and p' is 90° clockwise rotated from q' and is treated as the principal direction. The reason of choosing p' as the principal direction instead of q' is that when L1 and L2 are projected to p' originated at point O, one of them will be on the positive side of p' and the other on the negative side, allowing a less-biased estimation. Experimental results also showed that choosing p' for the slope estimation yielded less noisy solutions.
Fig. 2. Determination of principal direction p' for quotient-based single direction slope estimation.
This method attempts to treat non-coplanar photometric stereo pairs as coplanar by estimating the slope along the direction that best describes both illumination directions; therefore, higher uncertainties are expected when the photometric stereo pair is largely non-coplanar (i.e., a strongly deviates from 0° and from 180°). However, we can initialize the pixel-wise intrinsic albedo without external information or assumptions as follows:
Mxy —
The uncertainties or the noise in the albedo can be partially suppressed by a smoothing kernel; the intrinsic albedo maps A1 and A2 are combined by using the mean or the weighted mean depending on factors such as the image quality and the imaging conditions. The images and the estimated albedo map are forwarded to an iterative photometric stereo SAfS routine that is composed of three steps: (1) albedo refinement, (2) slope from reflectance, and (3) height from slope. These three steps are similar to those of the method proposed by Wu et al. (2018) while steps (2) and (3) are common in most reflectance-based reconstruction algorithms.
''Albedo refinement" mainly re-computes the pixel-wise intrinsic albedo map Axy on the basis of the slope map of the current iteration (p,q)\x,y. It is constrained so as to not deviate largely from the albedo map of the previous iteration, leading to the following equation:
Ak — (w1)Ak-1 U + (1 - w1)
Gl(p, q)
(1 - W2)
G2(p, q).
x,y (5)
where Ak\x,y is the intrinsic albedo at the kth iteration at point (x,y); wj and w2 are the weights ranging from 0 and 1 used to adjust the contribution of various components. The determination of weights depends on the quality of each image, however it is believed that the values of each component are relative consistent, hence the value of Ak\xy will be stable within a reasonable range of weights. Ak\x,y can neither be negative nor be greater than 1. A smoothing operation is applied after computation.
''Slope from reflectance" estimates the surface slopes (p,q)\x,y on the basis of the refined intrinsic albedo Axy in the current iteration. Because the intrinsic albedo is assumed to be known in this step, two images would be sufficient to solve for p and q simultaneously by non-linear optimization.
''Height from slope" estimates the DEM height values from the slope map. The DEM was first initialized using the Frankot and Chellappa algorithm (Frankot and Chellappa, 1988). This algorithm is known for enforcing the integrability of a given surface slope map. In an integrable gradient field (i.e., slope map) the height difference between any two points is independent of the path chosen.
Hence, valid gradient fields derived by SAfS techniques must be integrable. Agrawal et al. (2006) developed a toolbox for the implementation of this algorithm, which is available at (http://www.cs. cmu.edu/~lLlM/projects/lM/aagrawal/software.html). The implementation of the Frankot and Chellappa algorithm by Agrawal et al. (2006) is incorporated in the photometric stereo SAfS implementation described in this paper. Then, the DEM is refined itera-tively by the method proposed by Wu et al. (2018) as follows: (1) Each surface slope is defined by the DEM nodes at the four corners; hence, each DEM node affects the four surrounding surface slopes. (2) Each DEM node is optimized such that the surrounding surface slopes best fit the photometric stereo SAfS estimated slopes. The DEM nodes are updated one by one according to the modified Jacobi scheme described by Horn (1990) and the iteration stops when the update values are very small or the maximum number of iterations is reached.
3.2. Illumination embedded photometric stereo SAfS
In this section, we will derive the mathematical relationship between photometric stereo SAfS and the directional difference between the involved illuminations. The mathematical analysis described here focuses on slope estimation, which is a direct and strictly forward outcome of photometric stereo SAfS. As shown in Fig. 3, the derivation focuses on a single point O on a surface and on a case wherein only one photometric stereo pair is available (i.e., double-image SAfS). It also assumes that the zenith angles of both the illumination sources are non-zero (i.e., non-overhead sun) to produce the non-zero 2D illumination vectors. The unit
vector En points towards illumination direction 1; it is also
aligned with the horizontal axis p of the Cartesian coordinate
frame. Therefore,
The unit vector
points
^11 = I"1"
towards illumination direction 2, which is rotated by an angle of a counter-clockwise from illumination direction 1. The two unit vectors can be modeled as follows:
= R(a)
pLi cos a - sin a r cos a '
.qL1. sin a cos a 0 sin a
where a is the azimuthal difference between the two illumination directions. Note that R(a) denotes a counter-clockwise rotation, which implies that a will be negative when a clockwise rotation is needed. The photometric stereo SAfS is then decomposed into two Sl-SAfS, one along each illumination direction, resulting in
slope vectors
Fig. 3. Determination of surface slope at point O through photometric stereo SAfS.
point of the perpendiculars of
PTi and PL2
Pi and P2
qi_ q2
PSAfS = mi PTi Pi" = m2 P?2 P2 "
+ +
qSAfS .q^i. .qi. q2
originating at
- sin a cos a
P?i q?
and m2 and ( m2
Pj22 to
Eqs. (7) and (8) yields
PSAfS ' Si "
qSAfS mi
scalars such that intersect. Combining where s1 is solved by
(6) SI-SAfS along illumination direction 1 and m1 can be solved by
searching m, such that projecting
" Si " to PL2
mi qL2
yield s2:
[ Si mi ]
Pi and P2 where
.qi. q2 mi
s2 - si PL2 s2 - si cos a
qL2 Therefore,
S2 cos a s2 sin a
PSAfS qSAfS
S2 -Si cos a sin a
S2 . sin a
. Si tan a
sj and s2 are scalars computed by Sl-SAfS along each illumination direction; the unit vectors will be scaled accordingly to obtain the desired slope suggested by Sl-SAfS. The two Sl-SAfS slope vectors
provide conditions for the photometric stereo SAfS
Pi and P2
qi. q2
solution, which implies that the final solution
PSAfS qSAfS
has to satisfy
and [P2 ] simultaneously when projected to q2
pTi and PL2
.où. qL2
respectively. This can be derived by computing the intersection
From Eq. (10), we inferred that pSAjS is solely dependent on sj, the Sl-SAfS slope along illumination direction 1; this finding can be attributed to that fact that it lies along the horizontal axis p of the coordinate frame. qSA[s is affected by both sj and s2 and the Sl-SAfS slope along both of the illuminations of the photometric stereo pair. Their relative contributions depend on the azimuthal illumination difference a. As a result, Eq. (10) contains the relationship between the photometric stereo SAfS solution and the azimuthal illumination difference, allowing for further analysis and propagation.
3.3. Error propagation of illumination embedded photometric stereo SAfS
3.3.1. Error model
The sources of uncertainties of the solution as in Eq. (10) stem from s1, s2, and a. Among the three, only s1 and s2 are directly related to SAfS; therefore, similar to Liu and Wu (2017), we assumed that a contains a negligible error. To propagate the error
caused by SI-SAfS along each direction, we first assume that
is the true solution for photometric stereo SAfS and that the solution where s1 and s2 contain certain errors:
pSAfS0
Qsajs'
si + Dsi
(s2 +ds2 )-(s1 +as1 ) cos a
Si + ASi
s2+as2 s1+as1 sin a tan a
pSAfS0 i< qSAfS0 ]
zero. As a result, the curve of c will be U-shaped. When (r » 1), (-4-) will dominate and overcome the effect of (-^i—), thus
sin2 a sin a tan a'7
making the curve of c U-shaped. A special case exists when r is equal to or very close to 1. At small a, sin a « tan a, and therefore,
tan2 a sin2 a sin a tan a This implies that
tan2 a sin2 a sin a tan a Hence, when (r = 1), 1 r2 2r
sin a tan a sin a tan a
tan2 a
sin2 a
sin a tan a
where Asj and As2 are positive or negative errors contained in sj and s2, leading to a distorted solution. Note that Asj and As2 are a combination of multiple errors from factors such as the reflectance modeling error and the albedo error and are implicitly represented by Eq. (11). The squared error e2 will then be the sum of the squared difference between the true solution and the distorted solution:
£2 = Apsafs + a9|ajs = (psas - psafs)2 Rearranging and substituting yields
" (qsafs' - qsafs)2
= apIajs
= (1 +
+ ^Af — As2 + . 2
As1As2
tan2 a
)As2 + )As2 - 2
sin2 a
sin a tan a tan2
As1As2 sin a tan a
Note that [(1 + ta;§^)As§ + (sm^ra)Ds2] is in fact the error model suggested by Liu and Wu (2017); the two components (1 + tan^) and (—2a implying the error curve with respect to a will always be a regular U-shape and its minimal will be at a = 90° in all cases. However the model only considered the effect of the azimuthal illumination difference a and it was concluded in Liu and Wu (2017) that the model only describes the general error pattern with considerable uncertainties. The error model presented in this paper has an additional term (-2 siAaAana). This additional term allows the involvement of additional parameters such as illumination zenith angles and image intensities of a photometric stereo pair, which results in a more realistic and robust error modeling. If one treats As2 as a scaled term of Asj such that (As2 = rAsi), the Eq. can be rewritten as follows:
tan2a 1
(sin2 ( 1
As22 - 2
As1As2
r2As21 2
sin a tan a As1rAs1
sin a tan a
sin a tan a
= c2As1
\Jc2As\ — COs1 —
sin a tan a
The coefficient c implies a slope root mean square error (RMSE) of c units with respect to a slope error of 1 because of SI-SAfS along the first illumination direction. The component + is U-
tan2 a sin2 a
shaped with respect to a for any real r, and its minimum is at 90°. However, the term (-sinf2trana) changes the shape of c depending on a and r. When (0 < r < 1), (—will dominate because (
tan2 a sin2 a
and (-sinO^na) will be scaled down by r and will be very close to
c2 = 1 +-
tan2 a
sin a tan a
This implies that the error is very small at small a when (r ffi 1) and increases with an increase in a, creating an exponential-shape curve. In general, the curve of the error coefficient c with respect to the azimuthal illumination difference a has a flat U-shape for most cases except when (r ffi 1), where smaller c can be found at smaller a.
3.3.2. Estimation of the error model parameters
In the previous section, we showed that the error model depends on the parameters r and a. The azimuthal illumination difference a can be retrieved directly from the imaging/illumination conditions, which are assumed to be given, while r denotes the ratio between the SI-SAfS errors along the two illumination directions (i.e., r — A;2) and is more difficult to retrieve. Hence, the estimation of r becomes important to understand the relationship between the azimuthal illumination difference and the SAfS quality. Assuming a Lambertian reflectance model, photometric stereo SAfS solves for a solution that simultaneously satisfies
Li(si)=si * sin bi + cosb = L
L-'l L - A
where Li(si) represents the Lambertian reflectance for illumination direction i and bi is the sun-zenith angle (i.e., 0° at overhead sun) of illumination direction i. Li is the expected reflectance value for image i computed using image intensity i and intrinsic albedo A. A is the same for all images i as it corresponds to the same point on the surface, in this case, (i — 1,2) as only two images are considered. Li(si) can be simplified by assuming s? + 1 — 1) to produce a unique solution without iterations. Hence, Eq. (13) can be simplified as follows:
Li(si) — si * sin Pi + cos Pi — Li Therefore, Li - cos Pi
sin Pi
If we assume that the sun-zenith angle (i.e., ft) contains negligible error, the SAfS solution si will solely be affected by the expected reflectance Li. If the errors in It are assumed to be negligible, which is possible for radiometrically well-calibrated images, Li will be
Fig. 4. Verification routine with real data for the proposed error model.
M1410 45774
M1363 03002
Ml 646 22650
M1122 344699
M1168 337480
Reference DEM
Fig. 5. Pixel-synchronous LROC NAC images and reference DEM for experimental image set 1.
largely affected by the intrinsic albedo A according to Eq. (19). Thus, the errors in si (i.e., As,) can be modeled as follows:
Ds(= @si M=-
( A2 sin ft
where AA corresponds to an error in A. Substituting Eq. (21) to r yields the following:
As2 _ -J2AA A2 sin ft _ J2 sin b
Ds1 A2 sin ft -J1DA J1sin b
Eq. (22) implies that the error model presented in Eq. (14) relies on not only the azimuthal illumination difference (horizontal angular difference) but also the zenith illumination difference (vertical angular difference). The image intensity ratio of the photometric stereo pair also contributes to the error model through Eq. (22). Note that this estimation does not require external information as the illumination zenith angles (i.e., ft) and the image intensities (i.e., Ii) are given at the very beginning. ln this case, r is estimated on the basis of a simplified Lambert reflectance model; the viewing conditions (e.g., viewing direction and zenith) may also contribute to r if other reflectance models such as the Lommel-Seeliger law, lunar-Lambert model, or Hapke model are utilized.
ln general, I1 and I2 are positive as they represent the energy recorded by the sensor; the illumination zenith angles (i.e., bi) are positive because (1) bi are assumed to be non-zero as stated in Section 3.1(2) negative ft implies a 180° change in a and hence, is com-
Table 1
Imaging and illumination conditions of image set 1.
lmage lD Subsolar Sun-zenith Sub-spacecraft Emission
azimuth (°) angle (°) azimuth (°) angle (°)
M141045774 5.53 69.13 292.84 19.62
M136303002 313.10 75.52 61.98 11.67
M164622650 279.50 86.23 161.92 1.15
M1122344699 59.22 79.64 4.18 1.69
M1168337480 46.46 76.77 269.97 1.14
pensated by a. Therefore, r will always be positive according to Eq. (22), which implies that the errors in Sl-SAfS along both of the directions (i.e., As1 and Ds2) will increase or decrease simultaneously. This is reasonable in the sense that when the albedo A is lower than
b) Best-fit slope error estimate and actual slope error
Fig. 6. Comparison between estimated theoretical error from the proposed error model and actual error from photometric stereo SAfS for image set 1.
Table 2
Azimuthal difference (a)and zenith difference (embedded in r) of the illuminations for photometric stereo pairs in image set 1.
^^ r a M141045774 M136303002 M164622650 M1122344699 M1168337480
M141045774 1.373 3.128 0.489 0.641
M136303002 52.4° 2.613 0.674 0.872
M164622650 86.0° 33.6° 1.384 1.535
M1122344699 53.7° 106.1° 139.7° 0.769
M1168337480 40.9° 93.4° 127.0° 12.8°
Table 3
Estimated error from the proposed error model and actual photometric stereo SAfS error for image set 1.
12.8 33.6 40.9 52.4 53.7 86.0 93.4 106.1 127.0 139.7
r 0.769 2.613 0.641 1.373 0.489 3.128 0.872 0.674 1.535 1.384
c 1.369 3.368 1.015 1.388 1.008 3.225 1.367 1.407 2.854 3.468
e (rs1 = 0.070) 0.096 0.235 0.071 0.097 0.070 0.226 0.096 0.098 0.200 0.243
Actual slope error 0.105 0.223 0.092 0.066 0.106 0.234 0.074 0.095 0.209 0.241
Actual angular error (°) 4.64 11.30 4.02 3.12 4.95 11.63 3.40 4.45 10.38 12.30
its true value, both SI-SAfS slopes will tend to incline towards the sun, leading to a positive error for SI-SAfS along both the illumination directions. When A is higher than its true value, both SI-SAfS
slopes will tend to incline away from the sun, leading to a negative error for SI-SAfS along both of the illumination directions. Therefore, the ratio of the slope error in SI-SAfS will always be positive.
f) a =86.0° g) a = 93.4° h) a =106.1° i) a =127.0° j) a =139.7°
Color scale (m)
Fig. 7. 3D views of photometric stereo SAfS DEMs at different a for dataset 1.
Fig. 8. Angular error maps from photometric stereo SAfS at different a for image set 1.
4. Experimental verifications
4.1. Overview of experimental verification routine
The proposed error model for photometric stereo SAfS described in Section 3 is verified using real lunar images. As shown in Fig. 4, the verification routine utilizes the Lunar Reconnaissance Orbiter Camera Narrow Angle Camera (LROC NAC) Calibrated Data Record (CDR) images. Multiple images covering the same area under different illumination conditions (i.e., photometric stereo image set) were obtained. They were co-registered and resampled to align and synchronize the pixels. This process was performed in ArcGlS 10.0 using manually selected tie points. The pixel-synchronous images were then geo-referenced to a reference DEM of the same area created independently by using photogram-metric methods (Burns et al., 2012). Both the images and the DEM are available in the LROC archive (http://lroc.sese.asu.edu/archive). The photometric stereo image set was decomposed into multiple photometric stereo pairs (i.e., two images per pair) with each pair containing a specific illumination azimuthal difference a. For each photometric stereo pair, the error coefficient c was estimated using the error model in Eq. (14) given a, and the parameter r was estimated using Eq. (22). The parameter r was first estimated pixel-wise and then summarized for the entire image by using the mean, median, or other appropriate statistical methods. ln this study, we used the mode values in the histogram of r to avoid outliers and retain representativeness. ln contrast, the photometric stereo pair was forwarded to photometric stereo SAfS reconstruction to produce a slope map; the photometric stereo SAfS slope map was compared with the slope map generated from the reference DEM treated as the ground truth. Thus, the actual slope error was analyzed. The verification results were obtained by comparing the estimated slope error obtained using the proposed error model and the actual slope error. The slope error was computed by first resampling the slope map from photometric stereo SAfS to the dimension of the ground truth and then applying the following equation:
1J^nY^^pf-pf^^^f-qf^ (23)
RMSEsiope =
where n is the number of samples (i.e., pixels). The angular error was also analyzed for a more intuitive presentation; it described the inclination angle in three dimensions between the photometric stereo SAfS surface normals and the normal of the ground truth as follows:
angular :
1 ^^ n
pSAfSpref + qSAfSqref '
PSAfS '
' qSAfS '
ref + ref
4.2. Experimental analysis for LROC NAC image set 1
4.2.1. Overview of image set 1
Experimental image set 1 contains five LROC NAC CDR images covering the same region. Fig. 5 shows the pixel-synchronous images (2431 x 372 pixels) processed according to Section 4.1; their imaging conditions are summarized in Table 1. The resolutions for both the pixel-synchronous images and the reference DEM are 2 m/pixel. Because of the rotation characteristics of the Moon, most image sets that covered a wide range of illumination differences are found in higher latitudes. The region is located in the southern hemisphere (i.e., 70°S latitude); therefore, the illumination mainly comes from the north of the region and the sun-zenith angle is relatively high. The five images produce 10 photometric stereo pairs with their azimuthal and zenith illumination differences summarized in Table 2.
4.2.2. Verification results for image set 1
Table 3 summarizes the essential parameters for the proposed error model, the error estimated using the proposed error model in the case of real data, and the actual error from photometric stereo SAfS. Fig. 6a shows the plot of the estimated error coefficient c with respect to a. lt is apparent that when a is 33.6°, 86.0°, 127.0°,
b) Absolute difference DEMs
Fig. 9. Vertical RMSE and the corresponding absolute difference DEMs of image set 1.
and 139.7°, c jumps significantly because of the comparatively high r. From Eq. (17), we inferred that a large r implies that the error in Sl-SAfS along illumination direction 2 (i.e., As2) is a few times larger than the Sl-SAfS error in illumination direction 1 (i.e., As3), contributing an even larger error to the final result. These photometric stereo pairs with a high error coefficient c shared image M164622650, as shown in Tables 1 and 2. This image has a sun-zenith angle of 86.23°, which implies that the sunlight is coming from almost the horizon, which makes the image darker than the others, possibly creating a considerable amount of cast shadow, hence leading to relatively large errors. Fig. 7 reveals that the photometric stereo SAfS DEM derived from these pairs (i.e., Fig. 7b, f, i, and j) have significantly fewer topographic details than the other pairs. The angular error maps in Fig. 8 also suggest similar findings; when a DEM contains fewer topographic details than the ground truth, small local topographic features will be represented by smooth surfaces, yielding a larger angular error in these regions and hence, increasing the global angular error. lt is noted that there are apparent horizontal stripes on the angular error
maps in Fig. 8. They exist because there are small but dense systematic horizontal artefacts on the reference DEM obtained from the LROC NAC archive, which is treated as ground truth. These horizontal artefacts are inherited to the angular error maps during the comparison analysis. However, as all the photometric stereo pairs used in this dataset share the same ground truth, the horizontal artefacts on the reference DEM actually do not influence the comparison between different stereo pairs.
The shape of the error coefficient curve c with respect to a (Fig. 6a) exhibits a very high correspondence with the actual slope error shown in Fig. 6b (blue1 solid line). According to Eq. (14), the error estimated using the proposed error model (i.e., e) is the product of the error coefficient c and an estimate of the Sl-SAfS slope error along illumination direction 1 (i.e., as1). Assuming that the values of as1 of all photometric stereo pairs within image set 1 are identical, the estimated error curve yields a least squares best fit at as1 = 0.070,
1 For interpretation of color in Figs. 6 and 11, the reader is referred to the web version of this article
a) b) c) d) e)
M1225787254 M1153886196 M144349709 M1188124163 Reference DEM
Fig. 10. Pixel-synchronous LROC NAC images and reference DEM for experimental image set 2.
Table 4
Imaging and illumination conditions of image set 2.
Image ID Subsolar azimuth (°) Sun-zenith angle (°) Sub-spacecraft azimuth (°) Emission angle (°)
M1225787254 3.64 69.01 282.17 2.86
M1153886196 72.43 83.11 27.21 1.69
M144349709 45.46 74.16 165.26 1.69
M1188124163 278.71 85.37 357.80 1.13
Table 5
Azimuthal difference (a)and zenith difference (embedded in r) of the illuminations for photometric stereo pairs in image set 2.
r a M1225787254 M1153886196 M144349709 M1188124163
M1225787254 0.304 0.681 3.499
M1153886196 68.8° 0.438 1.055
M144349709 41.8° 27.0° 2.217
M1188124163 84.9° 153.7° 126.8°
approximately 4°, as shown in Fig. 6b (magenta dotted line). The result shows a high consistency with the actual error curve with a fitting RMSE of 0.017; this finding implies that the proposed error model can estimate the relative error of photometric stereo SAfS between different photometric stereo pairs and estimate the
photometric stereo SAfS slope error when the estimated SI-SAfS slope error along illumination direction 1 (i.e., as1) is given.
A vertical error analysis usually gives a more strictly forward and intuitive impression of the quality of the DEMs. With reference to the vertical RMSE plots with respect to a (Fig. 9a), larger fluctuations are observed for a 6 90°. The shape of the vertical error plot with respect to a is visually different from that of the slope error curves shown in Fig. 6; this difference is reasonable because the conversion of the reconstructed slopes to DEMs involves various processes such as offsetting and integrability enforcement. As a result, the DEM is in fact the best-fit solution of the slope map instead of a perfect fit; hence, the error patterns in the SAfS slope maps are unlikely to be inherited completely after the process. However, we might notice a higher RMSE at the lowest available a; this may imply that photometric stereo pairs with a 6 90° have a higher possibility of resulting in larger vertical errors. In-depth investigation is required to confirm this observation. Fig. 9b reveals that the grayish patches (i.e., areas with a larger absolute vertical difference) are in large clusters instead of scattered small patches; they also exhibit gradual changes from bright to dark. This finding probably implies that there are slope deviations in the low-resolution components (i.e., low-frequency components) of the
Table 6
Estimated error from the proposed error model and actual photometric stereo SAfS error for image set 2.
a(°) 27.0 41.8 68.8 84.9 126.8 153.7
r 0.438 0.681 0.304 3.499 2.217 1.055
c 1.413 1.005 1.002 3.567 3.653 4.520
e (rs1 = 0.062) 0.088 0.062 0.062 0.221 0.226 0.280
Actual slope error 0.175 0.134 0.131 0.199 0.192 0.267
Actual angular error (°) 7.943 5.921 5.990 9.908 9.213 13.630
resulting SAfS DEMs. Low-frequency inconsistency is a common issue in photoclinometric surface reconstruction algorithms and DEM estimation from slope algorithms (Frankot and Chellappa, 1988; Grumpe et al., 2014). Fig. 9b also reveals obvious spatial patterns in the vertical differences of the SAfS DEMs. This is because, the SAfS DEMs produced by the proposed algorithm did not contain any vertical datum; they were vertically aligned with the reference DEM by applying a vertical offset in order to facilitate the vertical comparison. The white lines shown in Fig. 9b are the intersection lines between the SAfS DEMs and the reference DEM, while gradual changes in vertical difference imply low-resolution slope errors. Similar patterns actually imply the performance of the proposed algorithm is generally stable under a wide spectrum of a. In general, although there are fluctuations depending on different values of a, the variations are not large; the vertical RMSE is always around 15 m, which is considered relatively stable.
4.3. Experimental analysis for LROC NAC image set 2
4.3.1. Overview of image set 2
Experimental image set 2 contains four LROC NAC CDR images covering a region located around 70°S latitude. Fig. 10 shows the pixel-synchronous images (1536 x 512 pixels) used for the verification, and the imaging conditions are summarized in Table 4. The resolutions for both the pixel-synchronous images and the reference DEM are 2 m/pixel. The illumination of the images mainly comes from the north and is similar to that in the case of image set 1; the solar zenith angle is relatively high, and two images out of four have a solar zenith angle of more than 80°. The four
b) Best-fit slope error estimate and actual slope error
Fig. 11. Comparison between estimated theoretical error from the proposed error model and actual error from photometric stereo SAfS for image set 2.
images produce six photometric stereo pairs; their azimuthal and zenith illumination differences are summarized in Table 5.
4.3.2. Verification results for image set 2
Table 6 summarizes the parameter r and estimated coefficient c values of the proposed error model, the error estimated by the proposed error model according to the real data situation, and the actual error from photometric stereo SAfS, Fig. 11a shows the corresponding plot of c with respect to a. When a is 84.9°, 126.8°, and 153.7°, the estimated error coefficient c increases significantly because of high r. This increase may be attributed to the image M1188124163 shared by these three pairs. Image M1188124163 (Fig. 10d) has a large solar zenith angle according to Table 4, and although its solar zenith angle is not significantly larger than that of M1153886196 (Fig. 10b), it is apparently darker than the other images, as shown in Fig. 10. This contributes to a larger value of r and possibly larger errors because of the presence of shadows. The 3D views of the photometric stereo SAfS DEMs in Fig. 12 reveal that the aforementioned pairs (Fig. 12d, e, and f) contain few local scale topographic details; the reconstructed shape of the larger crater at the bottom-left of the image almost vanishes in Fig. 11d and e because of strong shadows. In Fig. 12f, although the larger crater is reconstructed slightly better, most of the small topographic features such as craters and boulders could not be reconstructed very well in terms of the crater depth or the boulder height. The angular error maps in Fig. 13 present a visual correspondence between the dark patches in Fig. 13d, e, and f and the shadow areas in image M1188124163 (Fig. 10d); this correspondence suggests that the errors are likely to be associated with shadows.
The error coefficient plot c with respect to a (Fig. 11a) has a significant correspondence with the actual slope error in Fig. 11b (blue solid line). Although the similarity is not as high as in image set 1 because of a considerable increase in a > 84.9°, the shape correspondence of the curves is generally maintained. Assuming that rs1 for all photometric stereo pairs are the same, we found that the estimated error plot of e with respect to a obtained a least squares best fit at rs1 = 0.062, approximately 3.5°, with a fitting RMSE of 0.057. The deviation between the estimated error curve e from the proposed error model and the actual error curve is believed to stem from the fact that the iterative photometric stereo SAfS algorithm described in Section 3.1 has implicitly enforced the inte-grability of the resulting slope map, suppressing its errors to a certain extent. Integrability enforcement is an essential component of most shape from slope applications and hence, is not within the scope of our photometric stereo SAfS error analysis with respect to the azimuthal illumination difference. This also implies that the resulting slope maps for a = 84.9°, 126.8°, and 153.7° in this image set are likely to be more strongly non-integrable.
With reference to the vertical RMSE plots with respect to a for image set 2 (Fig. 14a), the RMSE curves of image set 2 exhibit a very slight U-shape similar to the shape found by Liu and Wu (2017). As in the case of image set 1, here, the vertical error plot presents a visible deviation from its slope error plot; the reasons for this deviation have been analyzed in Section 4.2.2. Dataset 2 also presents its highest RMSE at the lowest available a; this follows from the findings and implications of the experimental analysis of image set 1. Fig. 14b reveals large grayish clusters with gradual gray level changes. They are likely to be slope errors remaining in the low-resolution components of the resulting SAfS DEMs. Similar with the previous image set, they also present obvious similar spatial patterns, which indicates the stable performances of the proposed algorithm under various a. In general, the vertical RMSE of image set 2 corresponds to 12-13 m. In conclusion, the general performance of photometric stereo SAfS is relatively stable for different values of a in terms of the vertical analysis.
d) a = 84.9C
e) a= 126.8°
f) a= 153.7°
h) Color scale (m)
Fig. 12. 3D views of photometric stereo SAfS DEMs at different a for image set 2.
5. Conclusion and discussion
In this paper, we presented the explicit relationship between the slope error in photometric stereo SAfS and the illumination difference of the image pair. A case of two images was studied, and an error model was derived and presented, which will be useful in quality control and performance improvements. We found that the slope error in photometric stereo SAfS is affected by not only the illumination azimuthal difference but also the illumination zenith angles and the general intensity (e.g., mean intensity) of the images.
The error model was verified by photometric stereo SAfS using LROC NAC images; a photometric stereo SAfS algorithm was developed and implemented for the verification. The algorithm was developed on the basis of pairwise photometric stereo but can be
used for photometric stereo sets when more images are available. The results showed a good correspondence between the slope errors estimated from the proposed error model and the actual slope errors obtained from photometric stereo SAfS. Slight inconsistencies were noted; they were attributed to the integrability constraint enforced by the photometric stereo SAfS algorithm. Inte-grability enforcement is not included in the proposed error model, which relates the azimuthal illumination difference to the photometric stereo SAfS slope error. Scalar fitting between the error model estimates and the actual error suggested that our implemented photometric stereo SAfS algorithm led to a slope error along the first illumination direction of around 0.06-0.07, approximately 3-4°. The vertical error analysis did not show a very apparent pattern that corresponded to the proposed error model; this is reasonable because the point-to-point vertical error was unable to
Fig. 13. Angular error maps from photometric stereo SAfS at different a for image set 2.
b) Absolute difference DEMs
Fig. 14. Vertical RMSE and the corresponding absolute difference DEMs of image set 2.
directly represent slope errors and the process of converting surface slopes into heights includes other constraints such as enforcement of integrability and/or vertical offsetting, which might alter the error values. The vertical error analysis revealed the possible low-resolution slope deviations in photometric stereo SAfS, which are common in photoclinometry. In general, the vertical error fluctuations with respect to the directional illumination difference are small for both datasets.
When comparing SAfS with stereo photogrammetry, it can be concluded that stereo photogrammetry can provide accurate 3D points and hence achieving high large-scale accuracy, while SAfS excels at recovering local details at finest possible scale but with relatively lower accuracy at large scale. Stereo photogrammetry is effective in well textured surface with varying albedos where image matching can provide robust solutions, while SAfS is effective in poor textured surface with relatively homogenous albedo where image matching fails to provide desirable results. Stereo photogrammetry must use two or more images with desirable imaging geometry while SAfS can work with only one image and can provide more reliable results with the use of photometric stereo.
The proposed error model enriches the theory of photometric stereo SAfS. It can serve as a component for the quantitative quality control of photometric stereo SAfS. It can be combined with the SI-SfS error models proposed by Jankowski and Squyres (1991) to provide a more comprehensive overview of uncertainties in photometric stereo photoclinometry. Further investigations can be performed using the proposed model for improving the performance of photometric stereo SAfS algorithms. The photometric stereo SAfS algorithm and the error model proposed in this paper is of significance for optimized lunar surface reconstruction based on SAfS techniques.
Acknowledgments
The work described in this paper was funded by a grant from the Research Grants Council of Hong Kong (Project No.: PolyU 152086/15E) and grants from the National Natural Science Foundation of China (Project No.: 41471345 and Project No.: 41671426).
References
Agrawal, A., Raskar, R., Chellappa, R., 2006. What is the range of surface reconstructions from a gradient field?. Computer Vision. ECCV 2006. Lect. Notes Comput. Sci. 3951, 578-591.
Beyer, R.A., Kirk, R.L., 2012. Meter-scale slopes of candidate MSL landing sites from point photoclinometry. Space Sci. Rev. 170 (1-4), 775-791.
Burns, K.N., Speyerer, E.J., Robinson, M.S., Tran, T., Rosiek, M.R., Archinal, B.A., Howington-Kraus, E., 2012. Digital elevation models and derived products from LROC NAC stereo observations. Int. Arch. Photogram. Remote Sens. Spatial Inform. Sci. XXXIX (B4), 483-488.
Frankot, R.T., Chellappa, R., 1988. A method for enforcing integrability in shape from shading algorithms. IEEE Trans. Pattern Anal. Mach. Intell. 10 (4), 439-451.
Greenberg, R., Morris, M., Foley, M.A., 2011. The effect of sub-pixel slope variations on photoclinometry. Icarus 214 (1), 348-350.
Grumpe, A., Belkhir, F., Wohler, C., 2014. Construction of lunar DEMs based on reflectance modelling. Adv. Space Res. 53,1735-1767.
Hapke, B., 1981. Bidirectional reflectance spectroscopy. 1. Theory. J. Geophys. Res. 86 (B4), 3039-3054.
Hapke, B., 1984. Bidirectional reflectance spectroscopy. 3. Correction for macroscopic roughness. Icarus 59 (1), 41-59.
Hapke, B., 1986. Bidirectional reflectance spectroscopy. 4. The extinction coefficient and the opposition effect. Icarus 67 (2), 264-280.
Hapke, B., 2002. Bidirectional reflectance spectroscopy. 5. The coherent backscatter opposition effect and anisotropic scattering. Icarus 157 (2), 523-534.
Heipke, C., 1992. Integration of digital image matching and multiple image shape from shading. ESA SP 347. 391-391.
Horn, B.P., 1977. Understanding Image Intensities. Artif. Intell. 8 (2), 201-231.
Horn, B.K.P., 1990. Height and gradient from shading. Int. J. Comput. Vision 5 (1), 37-75.
Jankowski, D.G., Squyres, S.W., 1991. Sources of error in planetary photoclinometry. J. Geophys. Res. 96 (E4), 20907-20922.
Kirk, R.L., 1987. III. A Fast Finite-element Algorithm for Two-dimensional Photoclinometry, Ph.D. Thesis, Caltech. pp. 165-258.
Kirk, R.L., Barrett, J.M., Soderblom, L.A., 2003a. Photoclinometry made simple...?. Workshop on Advances in Planetary Mapping, Houston. url: <https:// astropedia.astrogeology.usgs.gov/download/Research/ISPRS/Kirk_isprs_mar03. pdf> (last accessed: 22-0ct-2017.
Kirk, R.L., Howington-Kraus, E., Redding, B., Galuszka, D., Hare, T.M., Archinal, B.A., Soderblom, L.A., Barrett, J.M., 2003b. High-resolution topomapping of candidate MER landing sites with Mars orbiter camera narrow-angle images. J. Geophys. Res. 108 (E12), 8088.
Labarre, S., Ferrari, C., Jacquemoud, S., 2017. Surface roughness retrieval by inversion of the Hapke model: a multiscale approach. Icarus 290, 63-80.
Liu, W.C., Wu, B., 2017. Photometric stereo shape-and-albedo-from-shading for pixel-level resolution lunar surface reconstruction. Int. Arch. Photogram. Rem. Sens. Spatial Inform. Sci. XLII-3 (W1), 91-97.
Lohse, V., Heipke, C., 2004. Multi-image shape-from-shading: derivation of planetary digital terrain models using clementine images. Int. Arch. Photogram. XXXV (B4), 828-833.
Lohse, V., Heipke, C., Kirk, R.L., 2006. Derivation of planetary topography using multi-image shape-from-shading. Planet. Space Sci. 54 (7), 661-674.
McEwen, A.S., 1991. Photometric functions for photoclinometry and other applications. Icarus 92 (2), 298-311.
O'Hara, R., Barnes, D., 2012. A new shape from shading technique with application to Mars Express HRSC images. ISPRS J. Photogramm. Remote Sens. 67, 27-34.
Oren, M., Nayar, S.K., 1994. Generalization of Lambert's reflectance model.In: Proceedings of the 21st Conference on Computer Graphics and Interactive, pp. 239-246.
Salamuniccar, G., Loncaric, S., Grumpe, A., Wohler, C., 2014. Hybrid method for crater detection based on topography reconstruction from optical images and the new LU78287GT catalogue of Lunar impact craters. Adv. Space Res. 53 (12), 1783-1797.
Wohler, C., 2004. Shape from shading under coplanar light sources. Pattern Recognition. DAGM 2004. Lect. Notes Comput. Sci. 3175, 278-285.
Wohler, C., Grumpe, A., Berezhnoy, A., Bhatt, M.U., Mall, U., 2014. Integrated topographic, photometric and spectral analysis of the lunar surface: application to impact melt flows and ponds. Icarus 235 (2014), 86-122.
Woodham, R.J., 1980. Photometric method for determining surface orientation from multiple images. Opt. Eng. 19 (1), 139-144.
Wu, B., Liu, W.C., Grumpe, A., Wohler, C., 2018. Construction of pixel-level resolution DEMs from monocular images by shape and albedo from shading constrained with low-resolution DEM. ISPRS J. Photogram. Rem. Sens. https:// doi.org/10.1016/j.isprsjprs.2017.03.007.
