Shape reconstruction of symmetric surfaces using Photometric Stereo Roberto Mecca and Silvia Tozza

Abstract The reconstruction of a 3D surface through one gray scale digital image does not admit a unique solution in the orthographic Shape from Shading (SfS) framework. With the aim to make this type of problem well-posed it is possible to use the Photometric Stereo (PS) technique. It allows to add information about the surface introducing other images of the object taken from the same point of view but modifying, for each photo, the direction of the light source. The methods that use the PS technique with the orthographic model of SfS need of, at least, three images. However, even if three images are used, there is the possibility that the SfS-PS problem continues to be ill-posed. This is the case when the three images are taken using three coplanar light vectors. This work analyses this kind of illposedness in order to understand how it is possible to establish a connection among the images that do not guarantee uniqueness. A further result in this paper is given by a geometrical characterization of the surfaces for which it is possible to solve the classic SfS problem.

1 Introduction to the Shape from Shading - Photometric Stereo model and symmetric surfaces The Shape from Shading (SfS) problem consists on the 3D reconstruction of an object starting from only one gray scale photo. In the framework of the shape reconstruction the classical orthographic SfS problem fails its aims using only one digital

Roberto Mecca Dipartimento di Matematica “G. Castelnuovo”, Sapienza - University of Rome e-mail: [email protected] Silvia Tozza Dipartimento di Matematica “G. Castelnuovo”, Sapienza - University of Rome e-mail: [email protected]

1

2

Roberto Mecca and Silvia Tozza

image for the reconstruction of the corresponding surface. Our contribution is given by the SfS model where we are considering the following hypotheses: 1. the light sources are at infinity; 2. the surface is Lambertian; 3. the optical point is sufficiently far from the surface so that perspective deformations can be neglected. In this case the SfS model is given by the following irradiance equation: ρ(x, y)(n(x, y) · ω) = I(x, y)

(1)

where ρ(x, y) is the albedo (supposed constant and equal to one), n(x, y) = (n1 (x, y), n2 (x, y), n3 (x, y)) is the unit outgoing normal to the unknown surface e ω3 ) is the direction of the light source and z = u(x, y), ω = (ω1 , ω2 , ω3 ) = (ω, I : Ω → [0, 1] is the image function. Under the previous assumptions this inverse problem is ill-posed due to the concave/convex ambiguity. Many papers show the impossibility to have a unique solution for this SfS model (see [2, 6]). In order to avoid this ambiguity we can add information of the surface considering other data of the same type. The Shape from Shading problem with the Photometric-Stereo techniques (SfSPS) uses more than one image of the same surface taken from the same point of view, but using a different light source for each image, that is: n(x, y) · ω i = Ii (x, y).

(2)

Let us clarify by now the meaning of coplanar light sources: since we consider each light source as a unit vector in R3 , we will talk about coplanar light sources with the same meaning of coplanar unit vectors. That is, we will say that three or more light sources are coplanar if the representing vectors are coplanar. There are many possible approaches to model the SfS-PS problem (e.g. [11, 12, 15]). The main two ways that allow to reconstruct the surface are the normal vector field approach [4] and the differential approach [13] which we briefly describe later. The main difference between these approaches is that the first one is a local approach that works pixel by pixel in order to find the normal vector field of the unknown surface and then it reconstructs the surface z = u(x, y) all over the domain using the gradient field [4]. The second approach instead works globally: it is based on a differential problem that, using only two images, can recognize the shape of the object but, as we will see, this differential problem is well-posed only if the Dirichlet boundary condition (i.e. the height of the surface on the boundary) is known. In order to use it we need to approximate the boundary condition and another image must be added making the problem well-posed with three images. We emphasize one more time that both approaches need three image functions and our goal is to characterize this kind of functions with the aim to reduce the number of images needed considering the symmetry of the surface. There are several previous studies [7, 3] where this geometric constraint is used together with

Shape reconstruction of symmetric surfaces using Photometric Stereo

3

the photometric features. In particular there are works that focus on this field detecting the symmetries in the images [14, 8]. Other papers in which a single image of a symmetric object can be regarded as two images of another object taken from two different points of view and two different illumination directions [20] or other similar approaches as [9, 10]. In Section 2 we recall the local and the differential approach. We explain the impossibility to solve the SfS-PS problem with both approaches using three images obtained with three coplanar light source vectors. In this section we give a definition of linear dependent images used in all the subsequent sections. In Section 3 we introduce an algorithm for the “artificial” reconstruction of linear dependent images. Finally in Section 4 we study how it is possible to reduce the number of necessary images for the shape reconstruction of a symmetric surface. The main result in this section is the theorem that explains how to solve the SfS problem (with only one image) through the SfS-PS formulation when the surface under observation has at least four axes of symmetry. In the last section we present some numerical tests related to the algorithms presented.

2 Condition of linear independent images for the SfS-PS reconstruction We recall briefly the mathematical formulation of the approaches introduced in the previous section.

2.1 Normal vector approach As we have anticipated, the local approach solves the problem finding the outgoing normal vector to the surface locally, that is pixel by pixel [4, 5]. This means that, once we discretize the image domain Ω with an uniform lattice Ω d = Ωd ∪ ∂ Ωd , the following system is solved for every (xi , y j ) ∈ Ω d : 0 0 0 ω1 n1 (xi , y j ) + ω2 n2 (xi , y j ) + ω3 n3 (xi , y j ) = I1 (xi , y j ) 00 00 00 ω1 n1 (xi , y j ) + ω2 n2 (xi , y j ) + ω3 n3 (xi , y j ) = I2 (xi , y j ) (3) 000 000 000 ω1 n1 (xi , y j ) + ω2 n2 (xi , y j ) + ω3 n3 (xi , y j ) = I3 (xi , y j ). As explained before, next step for this method is to approximate the values of the height of the surface starting from the knowledge of its gradient field (i.e. knowing ∇u(xi , y j ) ∀(xi , y j ) ∈ Ω d ) [4]. In fact, using the definition of the outgoing normal vector: ∂u ∂u − (x, y), − (x, y), 1 ∂x ∂y (−∇u(x, y), 1) p n(x, y) = p = , (4) 1 + ||∇u(x, y)||2 1 + ||∇u(x, y)||2

4

Roberto Mecca and Silvia Tozza

it follows: ∇u(xi , y j ) =

! n1 (xi , y j ) n2 (xi , y j ) ,− − . n3 (xi , y j ) n3 (xi , y j )

(5)

In the local method it is very easy to see that the uniqueness of the normal field can be reached using the non-singularity of the matrix with respect to the linear system (3), namely: 0 ω1 ω20 ω30 A = ω100 ω200 ω300 , (6) ω1000 ω2000 ω3000 that is supposing the non-coplanarity of the three light sources. An important example that explains why the case of coplanar light sources is not trivial, is represented by the natural phenomenon of the sun illumination of the solar system planets (like Earth). In fact, since these planets are moving on a elliptic trajectory around the sun, we can see the sun as a light source moving on a plane (see Fig. 1).

Fig. 1 Schematic representation of the Earth lighted up in three different times by the Sun.

Shape reconstruction of symmetric surfaces using Photometric Stereo

5

2.2 PDE approach The differential approach is based on a global method to solve the problem given by the following non-linear PDEs system of the Hamilton-Jacobi type: e 0 + ω30 −∇u(x, y) · ω p = I1 (x, y), ∀(x, y) ∈ Ω 1 + ||∇u(x, y)||2 e 00 + ω300 y) · ω −∇u(x, p = I2 (x, y), ∀(x, y) ∈ Ω 1 + ||∇u(x, y)||2

(7)

with the Dirichlet boundary condition u(x, y) = g(x, y) known for all (x, y) ∈ ∂ Ω . This system can also be studied eliminating the non-linearity [13]. We can then consider the following hyperbolic problem: ( b(x, y) · ∇u(x, y) = f (x, y), a.e. (x, y) ∈ Ω (8) u(x, y) = g(x, y), ∀(x, y) ∈ ∂ Ω where and

b(x, y) = (I2 (x, y)ω10 − I1 (x, y)ω100 , I2 (x, y)ω20 − I1 (x, y)ω200 )

(9)

f (x, y) = I2 (x, y)ω30 − I1 (x, y)ω300 .

(10)

For the differential method some other passages have to be done in order to see the necessity to have non-coplanar light sources. It is clear that only taking pictures is not possible to get the height of the surface on the boundary points (i.e. the boundary condition g(x, y)). This supplementary information is necessary to make (8) well-posed and that is why we need to compute it. In fact the approximation of the Dirichlet boundary condition g(x, y) binds us to use another image (we call it I3 , as in the previous approach). This step allows us to start solving the problem with the same starting data set (i.e. ω 0 , ω 00 , ω 000 and the respective images I1 , I2 , I3 ) for both approaches. It is clear that we can approximate the boundary condition g(x, y) using the previous local method that is using (3) for all pixels (xi , y j ) ∈ ∂ Ωd . With the aim to emphasize that also using the differential approach the condition of non-coplanarity is required, we use the combinations of the three images in the hyperbolic model in order to approximate the boundary condition. We consider the following system of hyperbolic equations: (1,2) (x, y) · ∇u(x, y) = f (1,2) (x, y), b b(1,3) (x, y) · ∇u(x, y) = f (1,3) (x, y), (2,3) b (x, y) · ∇u(x, y) = f (2,3) (x, y),

a.e. (x, y) ∈ Ω a.e. (x, y) ∈ Ω a.e. (x, y) ∈ Ω

where b(h,k) (x, y) = (Ik (x, y)ω1h − Ih (x, y)ω1k , Ik (x, y)ω2h − Ik (x, y)ω2k ) and

(11)

6

Roberto Mecca and Silvia Tozza

f (x, y)(h,k) = Ik (x, y)ω3h − Ih (x, y)ω3k where (h, k) is the combination of two of the first three natural integer without repetition. In other words (h, k) contains the indices (and the respective light sources) of the images we are using. Then, using these equations, a way to approximate the boundary condition is to choose two of the previous three equations and to consider them in a generic boundary point (xi , y j ) ∈ ∂ Ωd obtaining the following linear system: ∂u ∂u (1,2) (1,2) b1 (xi , y j ) (xi , y j ) + b2 (xi , y j ) (xi , y j ) = f (1,2) (xi , y j ) ∂x ∂y b(1,3) (xi , y j ) ∂ u (xi , y j ) + b(1,3) (xi , y j ) ∂ u (xi , y j ) = f (1,3) (xi , y j ). 2 1 ∂x ∂y

(12)

It permits to compute the gradient field of u on the boundary points and then we can approximate the value of the height of u using some well know algorithm [4, 5]. Considering the coefficient matrix of the latter linear system it is possible to find a relation between the non-singularity of the matrix A (from (6)) and the matrix of the linear system (12) proving the following proposition: Proposition 1. Let (1,2)

L(xi , y j ) =

b1 (xi , y j ) (1,3) b1 (xi , y j )

! (1,2) b2 (xi , y j ) . (1,3) b2 (xi , y j )

(13)

Then det(L(xi , y j )) = 0 ∀(xi , y j ) ∈ Ω d if and only if ω 0 , ω 00 , ω 000 are coplanar. Proof. (⇒) Let us prove that: det(L(xi , y j )) = 0 implies ω 0 , ω 00 , ω 000 coplanar. We generalize the proof by simplifying the notation instead of considering the dependence on a specific point (xi , y j ). We explicit the functions that appear in the matrix L in terms of the images and of the light sources and we obtain: det(L) =ω10 ω20 I2 I3 − ω10 ω2000 I1 I2 − ω100 ω20 I1 I3 + ω100 ω2000 (I1 )2 −

− [ω10 ω20 I2 I3 − ω1000 ω20 I1 I2 − ω10 ω200 I1 I3 + ω1000 ω200 (I1 )2 ]

=I1 [−ω10 ω2000 I2 − ω100 ω20 I3 + ω100 ω2000 I1 + ω1000 ω20 I2 + ω10 ω200 I3 − ω1000 ω200 I1 ]

=I1 [(−ω10 ω2000 + ω1000 ω20 )I2 + (−ω100 ω20 + ω10 ω200 )I3 + (−ω1000 ω200 + ω100 ω2000 )I1 ]. (14) Let us consider the case where the three images I1 , I2 and I3 are defined by the surface u and the respective light sources as follows:

Shape reconstruction of symmetric surfaces using Photometric Stereo

I1 (x, y) = I2 (x, y) = I3 (x, y) =

7

−ω10 ∂∂ ux (x, y) − ω20 ∂∂ uy (x, y) + ω30 p , 1 + ||∇u(x, y)||2

−ω100 ∂∂ ux (x, y) − ω200 ∂∂ uy (x, y) + ω300 p , 1 + ||∇u(x, y)||2

(15)

−ω1000 ∂∂ ux (x, y) − ω2000 ∂∂ uy (x, y) + ω3000 p . 1 + ||∇u(x, y)||2

Replacing the previous equalities in (14) we have: det(L) = ∂ u 0 00 000 1 ∂u I1 p (ω1 ω1 ω2 − ω20 ω100 ω1000 ) + (ω10 ω200 ω2000 − ω20 ω200 ω1000 )− 2 ∂x ∂y 1 + ||∇u(x, y)|| ∂u ∂u ω10 ω300 ω2000 + ω20 ω300 ω1000 + (ω20 ω100 ω1000 − ω10 ω200 ω1000 ) + (ω20 ω100 ω2000 − ω10 ω200 ω2000 )− ∂x ∂y ∂ u ∂ u ω20 ω100 ω3000 + ω10 ω200 ω3000 + (ω10 ω200 ω1000 − ω10 ω100 ω2000 ) + (ω20 ω200 ω1000 − ω20 ω100 ω2000 )− ∂x ∂y ω30 ω200 ω1000 + ω30 ω100 ω2000 . (16) Since the common denominator (1 + ||∇u(x, y)||2 ) is always different from zero, we can continue the computation considering only the numerator of (16). Explicit I1 using (15) we can write as follows: ∂u − ω10 − ω10 ω300 ω2000 + ω20 ω300 ω1000 − ω20 ω100 ω3000 + ω10 ω200 ω3000 − ω30 ω200 ω1000 + ω30 ω100 ω2000 + ∂x ∂u − ω20 − ω10 ω300 ω2000 + ω20 ω300 ω1000 − ω20 ω100 ω3000 + ω10 ω200 ω3000 − ω30 ω200 ω1000 + ω30 ω100 ω2000 + ∂y ω30 − ω10 ω300 ω2000 + ω20 ω300 ω1000 − ω20 ω100 ω3000 + ω10 ω200 ω3000 − ω30 ω200 ω1000 + ω30 ω100 ω2000 . (17) With the aim to prove the coplanarity of ω 000 with respect to ω 0 and ω 00 , we define the direction parameters of the plane generated by these last two vectors: s = ω20 ω300 − ω30 ω200 , l = ω30 ω100 − ω10 ω300 ,

(18)

m = ω10 ω200 − ω20 ω100 . The direction parameters (s, l, m) of a plane represent a vector orthogonal to it, therefore, orthogonal to every vector v = (v1 , v2 , v3 ) belonging to the plane. This means that: sv1 + lv2 + mv3 = 0. (19)

8

Roberto Mecca and Silvia Tozza

Now, considering det(L)(x, y) (without denominator) as a polynomial in the variables ξ = ∂∂ ux (x, y) and η = ∂∂ uy (x, y), we have: det(L)(ξ , η) = 0 000 0 00 0 00 000 0 00 0 00 000 0 00 0 00 + ξ − ω1 ω1 ω2 ω3 − ω3 ω2 + ω2 ω3 ω1 − ω1 ω3 + ω3 ω1 ω2 − ω2 ω1 + η − ω20 ω1000 ω20 ω300 − ω30 ω200 + ω2000 ω30 ω100 − ω10 ω300 + ω3000 ω10 ω200 − ω20 ω100 ω30 ω1000 ω20 ω300 − ω30 ω200 + ω2000 ω30 ω100 − ω10 ω300 + ω3000 ω10 ω200 − ω20 ω100 . (20) Since det(L) = 0 ∀(ξ , η) then this means that all the coefficients are zero: ω1000 ω20 ω300 − ω30 ω200 + ω2000 ω30 ω100 − ω10 ω300 + ω3000 ω10 ω200 − ω20 ω100 = 0 (21) which implies the coplanarity. (⇐) Let us prove now that: ω 0 , ω 00 , ω 000 coplanar implies det(L) = 0. We start seeing how, considering two images I1 and I2 obtained by the light vectors ω 0 and ω 00 , it is possible to get all the images obtainable with a light vector that is coplanar to the first two (i.e. considering the vectorial equality ω 000 = αω 0 +β ω 00 ), that is I3 = αI1 + β I2 as α and β change. Considering the image definition of the orthographic SfS model we can write: I3 = =

− ∂∂ ux ω1000 − ∂∂ uy ω2000 + ω3000 p 1 + ||∇u||2

− ∂∂ ux (αω10 + β ω100 ) − ∂∂ uy (αω20 + β ω200 ) + (αω30 + β ω300 ) p 1 + ||∇u||2

=α

(22)

− ∂∂ ux ω10 − ∂∂ uy ω20 + ω30 − ∂∂ ux ω100 − ∂∂ uy ω200 + ω300 p p +β 1 + ||∇u||2 1 + ||∇u||2

= αI1 + β I2 .

We continue replacing in (14) the coplanarity equality ω 000 = αω 0 + β ω 00 and, in consequence, I3 with αI1 + β I2 . We obtain: det(L) = − ω10 (αω20 + β ω200 )I1 I2 − ω100 ω20 I1 (αI1 + β I2 ) + ω100 (αω20 + β ω200 )(I1 )2 +

(αω10 + β ω100 )ω20 I1 I2 + ω10 ω200 I1 (αI1 + β I2 ) − (αω10 + β ω100 )ω200 (I1 )2 = − αω10 ω20 I1 I2 − β ω10 ω200 I1 I2 − αω100 ω20 (I1 )2 − β ω100 ω20 I1 I2 +

αω100 ω20 (I1 )2 + β ω100 ω200 (I1 )2 + αω10 ω20 I1 I2 + β ω100 ω20 I1 I2 +

αω10 ω200 (I1 )2 + β ω10 ω200 I1 I2 − αω10 ω200 (I1 )2 − β ω100 ω200 (I1 )2 = 0.

Shape reconstruction of symmetric surfaces using Photometric Stereo

9

Now we want to remark why, from the point of view of the image data, this is a required condition in order to solve this PS formulation of the SfS problem. In particular we proved that, if we consider a third image I3 which is obtained by a light source ω 000 coplanar with ω 0 and ω 00 (that we use to obtain the first two images) we do not add any information. Rather, we will see that it is possible to compute such I3 “artificially”. Remark 1. The modelization of the orthographic SfS problem allows us to define the linear dependence of the images through their associated light vectors. In fact the way we consider the vector ω 000 coplanar with the other two allows us to define, in the same way, the concept of linearly dependent image functions. In order to give a better explanation of the linearly independent images and their importance for the resolution of the problem, we start giving the following definition: Definition 1. Three image functions I1 , I2 and I3 , obtained from the orthographic SfS model, are linearly independent if they are generated by three non coplanar light vectors. The idea behind this kind of definition uses the orthographic SfS model in order to compare the linear dependence/independence of vectors in R3 with respect to the image functions.

3 Linear dependent image reconstruction An interesting application regarding the orthographic SfS model is related to the possibility of generating the images of an object by varying the light source between two points ω 0 and ω 00 for which we already have the respective images. We suppose to have two images I1 and I2 obtained with two light sources ω 0 and ω 00 respectively, it is possible to deduce all the linearly dependent images. In the geometrical sense, if we consider these two light vectors like points of the upper hemisphere B((0, 0, 0), 1), then it is possible to calculate all the obtainable images with a light source in the geodesic identified by ω 0 and ω 00 . The main idea is related to the approximation of the points of a geodesic defined between the two light source vectors on the unit sphere. The method that we use to calculate the vector defined from the geodesic is based on the Remark 1. In fact, once we fix the extremal points, we discretize the geodesic through an affine transformation that allows us to see this curve as an arc of a circumference on the xy plane. Let us start by observing that the angle between the two unit light vectors θˆ is computable using the dot product between them. In fact: θˆ = arccos(ω 0 · ω 00 )

(23)

y

10

Roberto Mecca and Silvia Tozza

z

z ω!

ω !!

ω! θˆ

i=n

ω !!

ω !!

ω !! θˆ

(x1 , y1 )

x

j = 1, . . . , n

x

y

y

y

Fig. 2 Examples of geodesics (in blue) on the unit sphere between ω 0 and ω 00 . z

Now, we determine the affine transformation which allows the use of one angular parameter for the discretization of all the curve. For this purpose, we specify the linear mapping T like a matrix in R3x3 . In order to compute the coefficients of this matrix (i.e. the degree of freedom) we fix the following three conditions: T (1, 0, 0)T = ω 0 T (cos θˆ , sin θˆ , 0)T = ω 00

(24)

T

T (0, 0, 1) = q. y

z

z

ω!

ω

ω !!

!

ω !!

ω! θˆ θˆ

θˆ θˆ

x

x

x

θˆ

x

y

y

Fig. 3 Examples of affine transformations T between the area defined by the arc of the circumferz ence (in green) and the one z relates to the geodesic (in blue). !

We denote q asωthe unit vector (opportunely oriented) orthogonal to the plane where the vectors ω 0 and ω 00 lie (such that they are outgoing from the plane, i.e. q3 > 0). We consider the plane, generated by the vectors ω 0 and ω 00 , to have the following direction parameters:

Shape reconstruction of symmetric surfaces using Photometric Stereo

s = ω20 ω300 − ω30 ω200 l = ω30 ω100 − ω10 ω300

11

(25)

m = ω10 ω200 − ω20 ω100 . We consider the normalization of the vector (s, l, m) like a possible vector q depending on the sign of m. If m < 0 we take the opposite direction to the normalization. Let us make explicit the condition imposed by (24) writing the coefficients of the matrix that represents the transformation:

ω10

0 T = ω2 ω30

ω100 −ω10 cos θˆ sin θˆ ω200 −ω20 cos θˆ sin θˆ ω300 −ω30 cos θˆ sin θˆ

q1

q2 .

(26)

q3

If we consider an angle θ such that 0 < θ < θˆ , we obtain a third light source ω¯ using just the previous application: ¯ T (cos θ , sin θ , 0)T = ω.

(27)

Now, if we have built ω¯ coplanar to the previous two light vectors, then we have ω¯ = αω 0 + β ω 00 . The next step is to determine the coefficients α, β ∈ R. We write the previous equality and we obtain the following overdetermined linear system: 0 00 ω¯ 1 = αω1 + β ω1 ω¯ 2 = αω20 + β ω200 (28) ¯ 0 00 ω3 = αω3 + β ω3 With the aim to calculate the coefficients that determine the third light source, we solve the following linear system of two equations in two unknowns: α C =b (29) β where C ∈ R2×2 and clearly b ∈ R2 . Selecting two of the three equations of (28) in order to avoid singularity, the matrix C depends on the first two light vectors that we can express in spherical coordinates as follows: ω 0 = (cos θ1 sin ϕ1 , sin θ1 sin ϕ1 , cos ϕ1 ) and

ω 00 = (cos θ2 sin ϕ2 , sin θ2 sin ϕ2 , cos ϕ2 )

(with θ1 , θ2 ∈ [0, 2π] and ϕ1 , ϕ2 ∈ [0, π2 )). Our aim is to solve this overdetermined linear system considering only two of the three equations in order to avoid that one of the first two equations has all coefficients close or equal to zero (for example it

12

Roberto Mecca and Silvia Tozza

can happen if we have θ1 = θ2 = 0 with respect to the first equation or if we have θ1 = θ2 = π2 with respect to the second equation). Then, in order to avoid a bad conditioned linear system, if we consider the following set: π 3π 5π 7π ¯ A = 0, ∪ , ∪ , 2π (30) 4 4 4 4 we take

while

0 00 ω2 ω2 α ω¯ 2 = ω30 ω300 β ω¯ 3 0 00 ω1 ω1 α ω¯ 1 = ω30 ω300 β ω¯ 3

if

θ1 , θ2 ∈ A¯ otherwise.

(31)

(32)

Once we determine the coefficients α and β , we can compute the linear dependent image as follows: I¯ = αI1 + β I2 . (33) It is important to note that there is no approximation when we compute the linear ¯ In fact, the only computation that can introduce an error in the dependent image I. computation of the third image is related to the approximation of the linear system (29) for which we consider two different formulations ((31) and (32)) in order to make unimportant the error committed.

4 Reduction of the number of the images using symmetries Supposing that we are able to determine univocally the surface u(x, y) of the SfSPS problem using three images obtained with three non coplanar light sources, now we want to consider some classes of particular surfaces for which it is possible to resolve the SfS-PS problem using less than three images.

4.1 Symmetric surfaces We will prove a theorem that gives us the possibility to produce an image (of a particular surface) by using another image of the same surface, obtained from a light source appropriately chosen. Let us start to define the type of surface we can use in order to apply that. Definition 2. Let z = u(x, y) be a function defined in Ω ⊂ R2 . Let πs (symmetry plane) be a plane passing through the z axis and let rs (symmetry straight line) be its intersection with the xy plane. We say that u is symmetric with respect to the plane πs if, for every point (x0 , y0 ) ∈ Ω , it results:

Shape reconstruction of symmetric surfaces using Photometric Stereo

13

u(x0 , y0 ) = u(x1 , y1 )

(34)

where (x1 , y1 ) ∈ Ω is the symmetric point i =1, . . .of , n(x0 , y0 ) with respect to the straight line rs (Fig. 4). j=0

y rs

ω !!

(x0 , y0 ) α Fig. 4 Orthogonal view of the plane (x, y): example of a positioning of the points (x0 , y0 ) and (x1 , y1 ) with respect to the symmetry straight line rs .

(x1 , y1 )

x

With regard to the orthographic SfS problem, since the model has not a direct dependence on the value of the function u in a point, but only on its gradient, it is essential to determine a symmetry relation with respect to the partial derivatives of the surface u. We identify rs like a particular straight line passing through the origin of the reference system xyz, that is x sin α = y cos α, with α ∈ [0, π].

(35)

We consider this straight line like a subspace generated by the vector (cos α, sin α). It is possible to identify the orthogonal straight line to rs from its base vector. Then, let π π v = cos α + , sin α + = (− sin α, cos α) (36) 2 2 be this orthogonal straight line. The main relation about the symmetry that we use in the SfS model is the following: v · ∇u(x0 , y0 ) = −v · ∇u(x1 , y1 ).

(37)

4.2 Uniqueness theorem for the symmetric surfaces We aim at obtaining further information relative to the symmetric surfaces by deducing the grey scale values of an image corresponding to a particular light source

y

14

Roberto Mecca and Silvia Tozza

whose position depends on the symmetry straight line rs . These information, together with a third image, allow us to resolve the SfS-PS problem using three images linearly independent. Theorem 1. Let u(x, y) be a symmetric surface with respect to any straight line rs 00 be two light vectors such that: such that (0, 0) ∈ rs . Let ω 0 and ωded 1. 2.

e 0 is orthogonal with respect to the symmetry straight line rs ; ω 00 =−ω 00 = (−ω e e 0 (that is ωded e 0 , ω30 ) since the light vectors are constrained to ωded be taken in the upper hemisphere).

00 only by Then it is possible to deduce the image associated to the light source ωded I1 as follows: (38) I2ded (x0 , y0 ) = I1 (x1 , y1 ).

Proof. Considering the image definition of the SfS problem, we have: 00 · ∇u(x , y ) + ω 00 eded e 0 · ∇u(x0 , y0 ) + ω30 ∗ −ω ω 0 0 ded 3 ∗∗ p p = = 1 + ||∇u(x0 , y0 )||2 1 + ||∇u(x0 , y0 )||2 e 0 · ∇u(x1 , y1 ) + ω30 e 0 · ∇u(x1 , y1 ) + ω30 −ω −ω = p = p = I1 (x1 , y1 ) (39) 1 + ||∇u(x0 , y0 )||2 1 + ||∇u(x1 , y1 )||2

I2ded (x0 , y0 ) =

where in ∗∗ we use the assumption 2 of the theorem, whereas in ∗ , as the surface is symmetric, we exploit (37). In the last equality we consider furthermore that the normal vector in the point (x0 , y0 ) has the same norm as in (x1 , y1 ). t u Let us explain with an example the previous theorem considering the symmetric surface of Fig. 5 that we will use again in the numerical test section. n3

ω�

SR ω

Pˆ n2

Pˆ � χ

R χ

n2

n1

ω

n(x, y) (x, y, u(x, y))

(a)

(b)

y

Fig. 5 On the left the 3D shape of the symmetric surface taken into account. On the right the orthographic view of the surface is shown in order to make clear the symmetry (the red line). The color is proportional to surface height, starting from the cold color (for the lowest values) to the hotter color.

Shape reconstruction of symmetric surfaces using Photometric Stereo

15

It is clear that its symmetry straight line rs drown in Fig. 5 (b) defined by the angle α = 3π 4 , permits us to apply the previous theorem starting with a first light e 0 orthogonal to rs . For example we source ω 0 having its projection on the xy plane ω π can choose θ1 = 4 and, with the only aim to avoid shadows, we choose ϕ1 = 0.1 obtaining the image I1 shown in Fig. 6 (on the left). Theorem 1 says that we can deduce (and not approximate) a second image I2ded shown in Fig. 6 (on the right) 00 = (−ω 0 , ω 0 ). In this example the angles that associated to the light source ωded 3 ded = ϕ . 00 define the polar coordinates of ωded are θ2ded = θ1 + π = 5π 1 4 and ϕ2 As we will remark the deduced image will help to have the sufficient information (about the three linear independent images) in order to solve the SfS-PS problem using less than three images.

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

00 : ϕ ded = 0.1, θ ded = ωded 2 2

5π 4

→ (I2ded )

Fig. 6 On the left the starting image obtained by the light source ω 0 and on the right the deduced 00 . image computed using the light source ωded

16

Roberto Mecca and Silvia Tozza

Once we choose (x0 , y0 ) ∈ Ω , the problem consists in determining the coordinates of the symmetric point (x1 , y1 ) with respect to rs . The point (x1 , y1 ) belongs to the straight line rs⊥ , orthogonal to rs and passing through the point (x0 , y0 ), whose equation is (x − x0 ) cos α = −(y − y0 ) sin α. (40) Finally, observing that the symmetric straight line rs passes through the origin, it results necessarily that the distances with respect to (0, 0) of both points are the same.

i =1, . . . , n j=0

y

y rs

(x0 , y0 )

(x0 , y0 ) α

(x1 , y1 )

x Fig. 7 In this schematic example is shown the circumference containing the two symmetric points (x0 , y0 ) and (x1 , y1 ) with respect to the straight line rs .

That is, both points belong to the same circumference centered in the origin of the axis as shown in Fig. 7. This formulation of the problem allows that the following non-linear system ( (x − x0 ) cos α = −(y − y0 ) sin α (41) x2 + y2 = x02 + y20 admits only two solutions, that is (x0 , y0 ) and (x1 , y1 ). The procedure we use to solve the system (41) is related to an algorithm that takes into account the numerical instability. To this purpose, some particularities appear. For first in the computation of the solutions of (41) we have to take into account the angle α. In order to avoid numerical instability, we consider two possible sets to which this angle could belong: π 3π π 3π A = 0, ∪ , π and B = , . (42) 4 4 4 4 Let us suppose that α ∈ A; we use the first equation of (41), and we write the first coordinate of the point (since the division by cos α does not involve numerical

Shape reconstruction of symmetric surfaces using Photometric Stereo

17

instability). We have then: x = (y0 − y)

sin α + x0 cos α

which, replaced in the second equation of (41), gives a second order equation (in the general case written as ay2 + by + c = 0) where the coefficients are: a = 1 + (tan α)2 b = −2(y0 (tan α)2 + x0 tan α)

(43)

c = y20 (tan α)2 + 2y0 x0 tan α − y20 .

The two solutions of this equation are y0 and y1 . It is possible to determine the exact solution since y0 is known. In the same way, if α ∈ B it is possible to divide by sin α. We write then the second coordinate of the point: y = (x0 − x)

cos α + y0 sin α

and, by substitution, we solve the second order equation in the variable x with the following coefficients: a = 1 + (cot α)2 b = −2(x0 (cot α)2 + y0 cot α) c

(44)

= x02 (cot α)2 + 2x0 y0 cot α − x02 .

The choice of the right solution can be done exactly as explained before, that is excluding the value x0 already known. In the end, if we consider an image I1 obtained with a light source that respects the hypotheses of the Theorem 1, it is possible to deduce another image I2 in the way we explained before. Corollary 1. Let u(x, y) be a symmetric surface with respect to a straight line rs passing through the point (0, 0). e 0 , ω30 ) such that ω e 0 is 1. Let I1 be an image obtained with a light vector ω 0 = (ω orthogonal to the symmetry straight line of u(x, y), e 000 , ω3000 ) such that 2. let I3 be an image obtained with a light source ω 000 = (ω 0 000 0 000 0 000 e ·ω e 6= ±||ω e ||||ω e || (that is ω e and ω e have not the same direction). ω Then it is possible to univocally determine the surface u(x, y) using the SfS-PS model. The use of the image I2 defined in the Theorem 1 is necessary to the application of the reconstruction technique using three images obtained with non coplanar light vectors.

18

Roberto Mecca and Silvia Tozza

Other simplifications in the calculus of the point (x1 , y1 ) are possible when the surface u(x, y) has a particular kind of symmetry. As a first particular case we consider now a symmetry with respect to the orthogonal straight line: let rs0 and rs00 be two symmetry straight lines of u orthogonal between them. Thus we have the following relation: ∇u(x0 , y0 ) = −∇u(−x0 , −y0 ) (45) for which we note the facility of calculus in the point previously called (x1 , y1 ) = (−x0 , −y0 ).

4.3 Surfaces with four symmetry straight lines A very interesting particular case is related to the surfaces with four symmetry straight lines. In relation to these surfaces, in fact, it is possible to solve the classical SfS problem considering only one image generated from any light source. For this case the constraint concerning the light sources direction vanishes. As before, we follow the procedure of reconstructing the missing information (that is the three linear independent images). Now, because of the particular surface geometry, we are able to obtain the necessary three images only from one. That is we can reconstruct the two missing images (generated with non coplanar light vectors) starting from only one image of the surface lighted up by any light source. Let us consider the image I1 and let ω 0 be its relative light source. The other two images can be calculated in the following way: e 0 , ω30 ) I2 (x, y) = I1 (−x, −y) with ω 00 = (−ω

I3 (x, y) = I1 (−y, x) with ω 000 = (ω20 , −ω10 , ω30 ).

(46)

Relating to the third image, we use (v1 , v2 ) · ∇u(x, y) = −(v2 , −v1 ) · ∇u(−y, x), where (v1 , v2 ) is any vector of the plane (x, y). The choice of the light placing, together with the way of determining the images, is not unique. The aim of this construction is to obtain linearly independent reconstructed images for which the uniqueness of solution for the SfS problem stands. We can adapt to this particular situation based on symmetries the uniqueness theorem in [13] that proves the uniqueness of a Lipschitz solution (i.e. a solution in the Sobolev space W 1,∞ (Ω )) relatively to the differential approach: Theorem 2. Let I(x, y) be an image of the surface u(x, y) ∈ W 1,∞ (Ω ), with four straight lines of symmetry obtained by the SfS model such that 0 < I(x, y) ≤ 1 ∀(x, y) ∈ Ω with the light source ω. Then it is possible to determine this surface univocally solving the SfS classic problem. The procedure followed until now is a constructive proof of the previous theorem and allows to solve, for a particular class of surfaces, the classical SfS problem passing through the SfS-PS one. Furthermore, our analysis of the problem allows us to give also a numerical method for the approximation of this unique weak solution.

Shape reconstruction of symmetric surfaces using Photometric Stereo

19

Remark 2. Following this procedure it is possible to prove the uniqueness of a weak solution of the classic SfS problem even the surface is radial by symmetries. Remark 3. For both cases (one or four axes) the construction of the other images can be done only if the light source of the first image is not vertical. This means that in this case, even using this image deduction technique, the concave/convex ambiguity cannot be solved.

5 Numerical Tests This section is divided in two parts. The first one is related to the computation of the linear dependent images explained in Section 3. In the second part we present some shape recovery tests relative to the local and differential method for the SfSPS using symmetric surfaces in order to compute the error of both methods using less than three images.

5.1 Numerical computation of linear dependent images The tests on the linear independent images concept is carried out on synthetic and real surfaces. For the synthetic case we start using the surface (with the respective initial images) shown in Fig. 8. In Fig. 9 it is possible to see the images of the same surface for which the light sources are equidistant in the geodesic defined by ω 0 and ω 00 .

3D surface

ω 0 : ϕ1 = 0.26, θ1 = 0.0 → (I1 ) ω 00 : ϕ2 = 0.19, θ3 = π → (I2 )

Fig. 8 From left to right: the 3D shape of the surface and the starting images with the respective light sources.

For the real case we start from the images shown in Fig. 10, namely I1real and In Fig. 11 it is possible to see three linear dependent images obtained with the

I2real .

20

Roberto Mecca and Silvia Tozza

I1

I2

Fig. 9 Between I1 and I2 there are three reconstructed images with the intermediate light sources taken equidistant in the geodesic.

algorithm explained in the previous section. Also in this case we are considering 0 00 . three light sources equidistant in the geodesic defined by ωreal and ωreal

00 : ϕ = 0.281, θ = 3.502 → (I real ) 0 : ϕ1 = 0.263, θ1 = −0.305 → (I1real ) ωreal ωreal 2 2 2

Fig. 10 Real starting data.

I1real

I2real

Fig. 11 The results of the computation using our algorithm. The three images on the middle are lighted up from equidistant light sources in the geodesic defined by the starting images.

Shape reconstruction of symmetric surfaces using Photometric Stereo

21

5.2 Shape reconstruction for symmetric surfaces Let us start the numerical tests by studying several surfaces. Every single one is of a typology that allows us to apply each variant of the previous results explained. These surfaces have symmetries that allows us to apply all the previous techniques that solve the SfS-PS problem.

vsymm (x, y)

vlip (x, y)

vrad (x, y)

Fig. 12 Set of surfaces used for the numerical tests.

The different typologies of the three surfaces represented in Fig. 12 are such that their geometrical characteristics differentiate the number of necessary images for their tridimensional reconstruction. In fact: • vsymm (x, y), used also in the previous numerical test, is a surface with one symmetry axis (the straight line y + x = 0) for which we need two images; • vlip (x, y) is a Lipschitz surface with four axis of symmetry (that is y = ±x, y = 0, x = 0). For its reconstruction we need only one image; • vrad (x, y) is a particular case of the previous type of surfaces. In this radial case the curves where it is not differentiable are concentric circles. For these starting data, summarized in Fig. 13, we make different tests modifying the size of the grid ∆ , uniform for all the domain Ω . In order to solve the differential problem (8), once boundary condition g(xi , y j ) = u(xi , y j ) ∀(xi , y j ) ∈ ∂ Ωd is computed, we use the two following semi-lagrangian fix point schemes [13]: n un+1 i, j = u (xi − hγ1 (xi , y j ), y j − hγ2 (xi , y j )) +

fi, j h, ||bi, j ||

∀(xi , y j ) ∈ Ωd

(47)

n un+1 i, j = u (xi + hγ1 (xi , y j ), y j + hγ2 (xi , y j )) −

fi, j h, ||bi, j ||

∀(xi , y j ) ∈ Ωd

(48)

and

22

Roberto Mecca and Silvia Tozza

i =1, . . . , n j=1

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

j=0

ω 00 : ϕ2 = 0.1, θ2 = 34 π → (I2 )

(x0 , y0 )

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

Fig. 13 Set of images used with the respective light sources described by their spherical coordinates. For each surface it is possible to see the straight line of symmetry (in red) and the curve where the surface is not differentiable (in blue). In the first line we have the only two images used to compute the shape of vsym . On the second line we have the only image needed to compute the surfaces vlip and vrad (respectively from left to right).

b(x,y) , fi, j is the discrete version of (10) and where γ(x, y) = (γ1 (x, y), γ2 (x, y)) = ||b(x,y)|| h is a small parameter that discretizes the directional derivative. Using a bilinear interpolation with respect to the points of the lattice Ωd at point (xi ±hγ1 (xi , y j ), y j ± hγ2 (xi , y j )) ∈ / Ωd , we can reach the first order of consistency choosing h = ∆ . The previous two semi-lagrangian numerical schemes will be respectively called in the sequel as forward (47) and backward (48). For these iterative schemes we use the same initial guess, that is u0i, j = 0 in Ωd and u0i, j = gi, j in ∂ Ωd . The numerical approximation is carried out with both semi-lagrangian schemes and we calculate the errors in L∞ (Ωd ) norm by varying the uniform discretization step ∆ . In particular we consider a domain Ω d = [−1, 1]2 using synthetic images of size starting from 100 × 100 pixels (∆ = 0.02) to 800 × 800 pixels (∆ = 0.0025). In all the following tables we compute also the convergence order with the aim to verify the convergence of the methods. In the last column we show the error committed using the approximation by the local normal approach [4] always in the same L∞ (Ωd ) norm. In all the tables it is possible to see how both the numerical schemes work well for all the cases considered. In particular it is possible to see how the semi-lagrangian iterative schemes have an order of about one for all the tests. This means that the

Shape reconstruction of symmetric surfaces using Photometric Stereo

23

∆ SL forward error SL forward order SL backward error SL backward order 0.02 5.254 × 10−2 6.856 × 10−2 1.0055 1.0021 0.01 2.618 × 10−2 3.428 × 10−2 1.0033 1.0008 0.005 1.307 × 10−2 1.713 × 10−2 1.0011 0.9992 −3 −3 0.0025 6.530 × 10 8.570 × 10

Local Normal 1.63 × 10−2 1.74 × 10−2 6.75 × 10−2 1.289 × 10−1

∆ SL forward error SL forward order SL backward error SL backward order 0.02 2.731 × 10−2 2.731 × 10−2 0.8975 0.8975 −2 0.01 1.466 × 10 1.466 × 10−2 0.8630 0.8630 −3 −3 0.005 8.060 × 10 8.060 × 10 0.9165 0.9165 −3 −3 0.0025 4.270 × 10 4.270 × 10

Local Normal 1.01 × 10−2 5.1 × 10−3 3.0 × 10−3 4.32 × 10−2

SL forward error SL forward order SL backward error SL backward order 2.5491 × 10−1 3.7057 × 10−1 1.0313 1.0423 1.2472 × 10−1 1.7993 × 10−1 1 0.9833 −2 −2 6.2360 × 10 9.1010 × 10 0.9963 0.9878 3.1280 × 10−2 4.5890 × 10−2

Local Normal 1.136 × 10−1 4.53 × 10−2 2.0 × 10−2 8.3 × 10−3

Table 1 Error in L∞ (Ωd ) norm and order of convergence for vsym with angles: θ1 = ϕ1 = ϕ2 = 0.1.

Table 2 Error in L∞ (Ωd ) norm and order of convergence for vlip with angles: θ1 =

∆ 0.02 0.01 0.005 0.0025

Table 3 Error in L∞ (Ωd ) norm and order of convergence for vrad with angles: θ1 =

π 4

π 4

π 4

and θ2 = 34 π;

and ϕ1 = 0.1.

and ϕ1 = 0.1.

approximation of the symmetric surfaces can be done without loss of precision considering less than three images.

6 Conclusion and perspectives We have seen how it is possible to resolve the orthographic SfS problem with the PS technique uniquely, getting over the concave/convex ambiguity. We have shown how the linear independence of images is fundamental to obtain a unique reconstruction of the surface. Also an interesting application is explained: the possibility to generate linear dependent images of an object changing the light source between two points ω 0 and ω 00 for which we already have the relative images. The numerical tests give good results for both the two approaches: normal vector field approach and differential one. The main perspective is to apply these two results to a weaker model that take into account a less number of assumptions. For example, the well-know perspective SfS model [18, 19, 17], represented by the following equation e · ∇u(x)) − ω3 (x · ∇u(x) + u(x)) f (ω p = I(x) f 2 |∇u(x)|2 + (x · ∇u(x) + u(x))2

24

Roberto Mecca and Silvia Tozza

(where f is the focal distance), contains some degree of ill-posedness. A further step forward can be done solving the perspective problem with the same technique explained before instead of the orthographic projection assumption. With the aim to eliminate the Lambertian surface assumption the specular model can be considered (see [16, 21, 1, 22, 23] ).

References 1. J. F. Blinn. Models of light reflection for computer synthesized pictures. Computer Graphics, 11(2):192–198, 1977. 2. M. J. Brooks, W. Chojnacki, and R. Kozera. Impossible and ambiguous shading patterns. International Journal of Computer Vision, 7(2):119–126, 1992. 3. T.-J. Cham and R. Cipolla. Symmetry detection through local skewed symmetries. Image and Vision Computing, 13(5):439 – 450, 1995. 4. J.-D. Durou, J.-F. Aujol, and F. Courteille. Integrating the Normal Field of a Surface in the Presence of Discontinuities. In Proceeding of Energy Minimization Method in Computer Vision and Pattern Recognition 2009, volume 5681 of Lecture Notes in Computer Science, pages 261–273. 5. J.-D. Durou and F. Courteille. Integration of a normal field without boundary condition. In Proceedings of the First International Workshop on Photometric Analysis For Computer Vision, 2007. 6. J.-D. Durou, M. Falcone, and M. Sagona. Numerical methods for shape from shading: a new survey with benchmarks. Computer Vision and Image Understanding, 109(1):22–43, 2008. 7. S. A. Friedberg. Finding axes of skewed symmetry. Computer Vision, Graphics, and Image Processing, 34(2):138 – 155, 1986. 8. Y. Gofman and N. Kiryati. Detecting symmetries in grey level images: The global optimization approach. In Proceedings of the International Conference on Pattern Recognition, 1996. 9. G.G. Gordon. Shape from symmetry. In Proceedings of SPIE Intellent Robots and Computer Vision VIII, 1989. 10. G.G. Gordon. 3d reconstruction from a single view of an object ant its image in a plane mirror. In Proceedings of International Conference on Pattern Recognition, 1998. 11. R. Kozera. Existence and uniqueness in photometric stereo. Appl. Math. Comput., 44(1):103, 1991. 12. R. Kozera. Uniqueness in shape from shading revisited. Journal of Mathematical Imaging and Vision, 7(2):123–138, 1997. 13. R. Mecca. Uniqueness for shape from shading via photometric stereo technique. In Proceedings of the 18th IEEE International Conference on Image Processing, pages 2994–2997, 2011. 14. D.P. Mukherjee, A. Zisserman, and J.M. Brady. Shape from symmetry - detecting and exploiting symmetry in affine images. Technical Report, Report No. OUEL 1988/93. University of Oxford, Department of Engineering Science. 15. R. Onn and A. M. Bruckstein. Integrability disambiguates surface recovery in two-image photometric stereo. International Journal of Computer Vision, 5(1):105–113, 1990. 16. B. T. Phong. Illumination for computer generated pictures. Communications of the ACM, 18(6):311–317, 1975. 17. E. Prados, F. Camilli, and O. Faugeras. A viscosity solution method for shape-from-shading without image boundary data. Math. Model. Numer. Anal. 40, pages 393–412, 2006. 18. E. Prados and O. Faugeras. Perspective shape from shading and viscosity solutions. In Proceedings of the 9th IEEE International Conference on Computer Vision (vol. 2), pages 826– 831, 2003.

Shape reconstruction of symmetric surfaces using Photometric Stereo

25

19. E. Prados and O. Faugeras. Shape from shading: a well-posed problem? In Proceedings of the IEEE Int. Conf. Comp. Vis. and Patt. Recog (vol. 2), pages 870–877, San Diego, California, USA, 2005. 20. I. Shimshoni, Y. Moses, and M. Lindenbaum. Shape reconstruction of 3d bilaterally symmetric surfaces. International Journal of Computer Vision, 39(2):97–110, 2000. 21. K. E. Torrance and E. M. Sparrow. Theory for off-specular reflecion from roughened surfaces. Journal of the Optical Society of America, 57(9):1105–1114, 1967. 22. P.S. Tsai and M. Shah. Shape from shading using linear approximation. Image and Vision Computing Journal (Elsevier, Oxford), 12(8):487–498, 1994. 23. O. Vogel, M. Breuss, and J. Weickert. Perspective shape from shading with non-lambertian reflectance. Lecture Notes in Computer Science. G. Rigoll (Ed.): Pattern Recognition, Springer, Berlin, 5096:517–526, 2008.

Abstract The reconstruction of a 3D surface through one gray scale digital image does not admit a unique solution in the orthographic Shape from Shading (SfS) framework. With the aim to make this type of problem well-posed it is possible to use the Photometric Stereo (PS) technique. It allows to add information about the surface introducing other images of the object taken from the same point of view but modifying, for each photo, the direction of the light source. The methods that use the PS technique with the orthographic model of SfS need of, at least, three images. However, even if three images are used, there is the possibility that the SfS-PS problem continues to be ill-posed. This is the case when the three images are taken using three coplanar light vectors. This work analyses this kind of illposedness in order to understand how it is possible to establish a connection among the images that do not guarantee uniqueness. A further result in this paper is given by a geometrical characterization of the surfaces for which it is possible to solve the classic SfS problem.

1 Introduction to the Shape from Shading - Photometric Stereo model and symmetric surfaces The Shape from Shading (SfS) problem consists on the 3D reconstruction of an object starting from only one gray scale photo. In the framework of the shape reconstruction the classical orthographic SfS problem fails its aims using only one digital

Roberto Mecca Dipartimento di Matematica “G. Castelnuovo”, Sapienza - University of Rome e-mail: [email protected] Silvia Tozza Dipartimento di Matematica “G. Castelnuovo”, Sapienza - University of Rome e-mail: [email protected]

1

2

Roberto Mecca and Silvia Tozza

image for the reconstruction of the corresponding surface. Our contribution is given by the SfS model where we are considering the following hypotheses: 1. the light sources are at infinity; 2. the surface is Lambertian; 3. the optical point is sufficiently far from the surface so that perspective deformations can be neglected. In this case the SfS model is given by the following irradiance equation: ρ(x, y)(n(x, y) · ω) = I(x, y)

(1)

where ρ(x, y) is the albedo (supposed constant and equal to one), n(x, y) = (n1 (x, y), n2 (x, y), n3 (x, y)) is the unit outgoing normal to the unknown surface e ω3 ) is the direction of the light source and z = u(x, y), ω = (ω1 , ω2 , ω3 ) = (ω, I : Ω → [0, 1] is the image function. Under the previous assumptions this inverse problem is ill-posed due to the concave/convex ambiguity. Many papers show the impossibility to have a unique solution for this SfS model (see [2, 6]). In order to avoid this ambiguity we can add information of the surface considering other data of the same type. The Shape from Shading problem with the Photometric-Stereo techniques (SfSPS) uses more than one image of the same surface taken from the same point of view, but using a different light source for each image, that is: n(x, y) · ω i = Ii (x, y).

(2)

Let us clarify by now the meaning of coplanar light sources: since we consider each light source as a unit vector in R3 , we will talk about coplanar light sources with the same meaning of coplanar unit vectors. That is, we will say that three or more light sources are coplanar if the representing vectors are coplanar. There are many possible approaches to model the SfS-PS problem (e.g. [11, 12, 15]). The main two ways that allow to reconstruct the surface are the normal vector field approach [4] and the differential approach [13] which we briefly describe later. The main difference between these approaches is that the first one is a local approach that works pixel by pixel in order to find the normal vector field of the unknown surface and then it reconstructs the surface z = u(x, y) all over the domain using the gradient field [4]. The second approach instead works globally: it is based on a differential problem that, using only two images, can recognize the shape of the object but, as we will see, this differential problem is well-posed only if the Dirichlet boundary condition (i.e. the height of the surface on the boundary) is known. In order to use it we need to approximate the boundary condition and another image must be added making the problem well-posed with three images. We emphasize one more time that both approaches need three image functions and our goal is to characterize this kind of functions with the aim to reduce the number of images needed considering the symmetry of the surface. There are several previous studies [7, 3] where this geometric constraint is used together with

Shape reconstruction of symmetric surfaces using Photometric Stereo

3

the photometric features. In particular there are works that focus on this field detecting the symmetries in the images [14, 8]. Other papers in which a single image of a symmetric object can be regarded as two images of another object taken from two different points of view and two different illumination directions [20] or other similar approaches as [9, 10]. In Section 2 we recall the local and the differential approach. We explain the impossibility to solve the SfS-PS problem with both approaches using three images obtained with three coplanar light source vectors. In this section we give a definition of linear dependent images used in all the subsequent sections. In Section 3 we introduce an algorithm for the “artificial” reconstruction of linear dependent images. Finally in Section 4 we study how it is possible to reduce the number of necessary images for the shape reconstruction of a symmetric surface. The main result in this section is the theorem that explains how to solve the SfS problem (with only one image) through the SfS-PS formulation when the surface under observation has at least four axes of symmetry. In the last section we present some numerical tests related to the algorithms presented.

2 Condition of linear independent images for the SfS-PS reconstruction We recall briefly the mathematical formulation of the approaches introduced in the previous section.

2.1 Normal vector approach As we have anticipated, the local approach solves the problem finding the outgoing normal vector to the surface locally, that is pixel by pixel [4, 5]. This means that, once we discretize the image domain Ω with an uniform lattice Ω d = Ωd ∪ ∂ Ωd , the following system is solved for every (xi , y j ) ∈ Ω d : 0 0 0 ω1 n1 (xi , y j ) + ω2 n2 (xi , y j ) + ω3 n3 (xi , y j ) = I1 (xi , y j ) 00 00 00 ω1 n1 (xi , y j ) + ω2 n2 (xi , y j ) + ω3 n3 (xi , y j ) = I2 (xi , y j ) (3) 000 000 000 ω1 n1 (xi , y j ) + ω2 n2 (xi , y j ) + ω3 n3 (xi , y j ) = I3 (xi , y j ). As explained before, next step for this method is to approximate the values of the height of the surface starting from the knowledge of its gradient field (i.e. knowing ∇u(xi , y j ) ∀(xi , y j ) ∈ Ω d ) [4]. In fact, using the definition of the outgoing normal vector: ∂u ∂u − (x, y), − (x, y), 1 ∂x ∂y (−∇u(x, y), 1) p n(x, y) = p = , (4) 1 + ||∇u(x, y)||2 1 + ||∇u(x, y)||2

4

Roberto Mecca and Silvia Tozza

it follows: ∇u(xi , y j ) =

! n1 (xi , y j ) n2 (xi , y j ) ,− − . n3 (xi , y j ) n3 (xi , y j )

(5)

In the local method it is very easy to see that the uniqueness of the normal field can be reached using the non-singularity of the matrix with respect to the linear system (3), namely: 0 ω1 ω20 ω30 A = ω100 ω200 ω300 , (6) ω1000 ω2000 ω3000 that is supposing the non-coplanarity of the three light sources. An important example that explains why the case of coplanar light sources is not trivial, is represented by the natural phenomenon of the sun illumination of the solar system planets (like Earth). In fact, since these planets are moving on a elliptic trajectory around the sun, we can see the sun as a light source moving on a plane (see Fig. 1).

Fig. 1 Schematic representation of the Earth lighted up in three different times by the Sun.

Shape reconstruction of symmetric surfaces using Photometric Stereo

5

2.2 PDE approach The differential approach is based on a global method to solve the problem given by the following non-linear PDEs system of the Hamilton-Jacobi type: e 0 + ω30 −∇u(x, y) · ω p = I1 (x, y), ∀(x, y) ∈ Ω 1 + ||∇u(x, y)||2 e 00 + ω300 y) · ω −∇u(x, p = I2 (x, y), ∀(x, y) ∈ Ω 1 + ||∇u(x, y)||2

(7)

with the Dirichlet boundary condition u(x, y) = g(x, y) known for all (x, y) ∈ ∂ Ω . This system can also be studied eliminating the non-linearity [13]. We can then consider the following hyperbolic problem: ( b(x, y) · ∇u(x, y) = f (x, y), a.e. (x, y) ∈ Ω (8) u(x, y) = g(x, y), ∀(x, y) ∈ ∂ Ω where and

b(x, y) = (I2 (x, y)ω10 − I1 (x, y)ω100 , I2 (x, y)ω20 − I1 (x, y)ω200 )

(9)

f (x, y) = I2 (x, y)ω30 − I1 (x, y)ω300 .

(10)

For the differential method some other passages have to be done in order to see the necessity to have non-coplanar light sources. It is clear that only taking pictures is not possible to get the height of the surface on the boundary points (i.e. the boundary condition g(x, y)). This supplementary information is necessary to make (8) well-posed and that is why we need to compute it. In fact the approximation of the Dirichlet boundary condition g(x, y) binds us to use another image (we call it I3 , as in the previous approach). This step allows us to start solving the problem with the same starting data set (i.e. ω 0 , ω 00 , ω 000 and the respective images I1 , I2 , I3 ) for both approaches. It is clear that we can approximate the boundary condition g(x, y) using the previous local method that is using (3) for all pixels (xi , y j ) ∈ ∂ Ωd . With the aim to emphasize that also using the differential approach the condition of non-coplanarity is required, we use the combinations of the three images in the hyperbolic model in order to approximate the boundary condition. We consider the following system of hyperbolic equations: (1,2) (x, y) · ∇u(x, y) = f (1,2) (x, y), b b(1,3) (x, y) · ∇u(x, y) = f (1,3) (x, y), (2,3) b (x, y) · ∇u(x, y) = f (2,3) (x, y),

a.e. (x, y) ∈ Ω a.e. (x, y) ∈ Ω a.e. (x, y) ∈ Ω

where b(h,k) (x, y) = (Ik (x, y)ω1h − Ih (x, y)ω1k , Ik (x, y)ω2h − Ik (x, y)ω2k ) and

(11)

6

Roberto Mecca and Silvia Tozza

f (x, y)(h,k) = Ik (x, y)ω3h − Ih (x, y)ω3k where (h, k) is the combination of two of the first three natural integer without repetition. In other words (h, k) contains the indices (and the respective light sources) of the images we are using. Then, using these equations, a way to approximate the boundary condition is to choose two of the previous three equations and to consider them in a generic boundary point (xi , y j ) ∈ ∂ Ωd obtaining the following linear system: ∂u ∂u (1,2) (1,2) b1 (xi , y j ) (xi , y j ) + b2 (xi , y j ) (xi , y j ) = f (1,2) (xi , y j ) ∂x ∂y b(1,3) (xi , y j ) ∂ u (xi , y j ) + b(1,3) (xi , y j ) ∂ u (xi , y j ) = f (1,3) (xi , y j ). 2 1 ∂x ∂y

(12)

It permits to compute the gradient field of u on the boundary points and then we can approximate the value of the height of u using some well know algorithm [4, 5]. Considering the coefficient matrix of the latter linear system it is possible to find a relation between the non-singularity of the matrix A (from (6)) and the matrix of the linear system (12) proving the following proposition: Proposition 1. Let (1,2)

L(xi , y j ) =

b1 (xi , y j ) (1,3) b1 (xi , y j )

! (1,2) b2 (xi , y j ) . (1,3) b2 (xi , y j )

(13)

Then det(L(xi , y j )) = 0 ∀(xi , y j ) ∈ Ω d if and only if ω 0 , ω 00 , ω 000 are coplanar. Proof. (⇒) Let us prove that: det(L(xi , y j )) = 0 implies ω 0 , ω 00 , ω 000 coplanar. We generalize the proof by simplifying the notation instead of considering the dependence on a specific point (xi , y j ). We explicit the functions that appear in the matrix L in terms of the images and of the light sources and we obtain: det(L) =ω10 ω20 I2 I3 − ω10 ω2000 I1 I2 − ω100 ω20 I1 I3 + ω100 ω2000 (I1 )2 −

− [ω10 ω20 I2 I3 − ω1000 ω20 I1 I2 − ω10 ω200 I1 I3 + ω1000 ω200 (I1 )2 ]

=I1 [−ω10 ω2000 I2 − ω100 ω20 I3 + ω100 ω2000 I1 + ω1000 ω20 I2 + ω10 ω200 I3 − ω1000 ω200 I1 ]

=I1 [(−ω10 ω2000 + ω1000 ω20 )I2 + (−ω100 ω20 + ω10 ω200 )I3 + (−ω1000 ω200 + ω100 ω2000 )I1 ]. (14) Let us consider the case where the three images I1 , I2 and I3 are defined by the surface u and the respective light sources as follows:

Shape reconstruction of symmetric surfaces using Photometric Stereo

I1 (x, y) = I2 (x, y) = I3 (x, y) =

7

−ω10 ∂∂ ux (x, y) − ω20 ∂∂ uy (x, y) + ω30 p , 1 + ||∇u(x, y)||2

−ω100 ∂∂ ux (x, y) − ω200 ∂∂ uy (x, y) + ω300 p , 1 + ||∇u(x, y)||2

(15)

−ω1000 ∂∂ ux (x, y) − ω2000 ∂∂ uy (x, y) + ω3000 p . 1 + ||∇u(x, y)||2

Replacing the previous equalities in (14) we have: det(L) = ∂ u 0 00 000 1 ∂u I1 p (ω1 ω1 ω2 − ω20 ω100 ω1000 ) + (ω10 ω200 ω2000 − ω20 ω200 ω1000 )− 2 ∂x ∂y 1 + ||∇u(x, y)|| ∂u ∂u ω10 ω300 ω2000 + ω20 ω300 ω1000 + (ω20 ω100 ω1000 − ω10 ω200 ω1000 ) + (ω20 ω100 ω2000 − ω10 ω200 ω2000 )− ∂x ∂y ∂ u ∂ u ω20 ω100 ω3000 + ω10 ω200 ω3000 + (ω10 ω200 ω1000 − ω10 ω100 ω2000 ) + (ω20 ω200 ω1000 − ω20 ω100 ω2000 )− ∂x ∂y ω30 ω200 ω1000 + ω30 ω100 ω2000 . (16) Since the common denominator (1 + ||∇u(x, y)||2 ) is always different from zero, we can continue the computation considering only the numerator of (16). Explicit I1 using (15) we can write as follows: ∂u − ω10 − ω10 ω300 ω2000 + ω20 ω300 ω1000 − ω20 ω100 ω3000 + ω10 ω200 ω3000 − ω30 ω200 ω1000 + ω30 ω100 ω2000 + ∂x ∂u − ω20 − ω10 ω300 ω2000 + ω20 ω300 ω1000 − ω20 ω100 ω3000 + ω10 ω200 ω3000 − ω30 ω200 ω1000 + ω30 ω100 ω2000 + ∂y ω30 − ω10 ω300 ω2000 + ω20 ω300 ω1000 − ω20 ω100 ω3000 + ω10 ω200 ω3000 − ω30 ω200 ω1000 + ω30 ω100 ω2000 . (17) With the aim to prove the coplanarity of ω 000 with respect to ω 0 and ω 00 , we define the direction parameters of the plane generated by these last two vectors: s = ω20 ω300 − ω30 ω200 , l = ω30 ω100 − ω10 ω300 ,

(18)

m = ω10 ω200 − ω20 ω100 . The direction parameters (s, l, m) of a plane represent a vector orthogonal to it, therefore, orthogonal to every vector v = (v1 , v2 , v3 ) belonging to the plane. This means that: sv1 + lv2 + mv3 = 0. (19)

8

Roberto Mecca and Silvia Tozza

Now, considering det(L)(x, y) (without denominator) as a polynomial in the variables ξ = ∂∂ ux (x, y) and η = ∂∂ uy (x, y), we have: det(L)(ξ , η) = 0 000 0 00 0 00 000 0 00 0 00 000 0 00 0 00 + ξ − ω1 ω1 ω2 ω3 − ω3 ω2 + ω2 ω3 ω1 − ω1 ω3 + ω3 ω1 ω2 − ω2 ω1 + η − ω20 ω1000 ω20 ω300 − ω30 ω200 + ω2000 ω30 ω100 − ω10 ω300 + ω3000 ω10 ω200 − ω20 ω100 ω30 ω1000 ω20 ω300 − ω30 ω200 + ω2000 ω30 ω100 − ω10 ω300 + ω3000 ω10 ω200 − ω20 ω100 . (20) Since det(L) = 0 ∀(ξ , η) then this means that all the coefficients are zero: ω1000 ω20 ω300 − ω30 ω200 + ω2000 ω30 ω100 − ω10 ω300 + ω3000 ω10 ω200 − ω20 ω100 = 0 (21) which implies the coplanarity. (⇐) Let us prove now that: ω 0 , ω 00 , ω 000 coplanar implies det(L) = 0. We start seeing how, considering two images I1 and I2 obtained by the light vectors ω 0 and ω 00 , it is possible to get all the images obtainable with a light vector that is coplanar to the first two (i.e. considering the vectorial equality ω 000 = αω 0 +β ω 00 ), that is I3 = αI1 + β I2 as α and β change. Considering the image definition of the orthographic SfS model we can write: I3 = =

− ∂∂ ux ω1000 − ∂∂ uy ω2000 + ω3000 p 1 + ||∇u||2

− ∂∂ ux (αω10 + β ω100 ) − ∂∂ uy (αω20 + β ω200 ) + (αω30 + β ω300 ) p 1 + ||∇u||2

=α

(22)

− ∂∂ ux ω10 − ∂∂ uy ω20 + ω30 − ∂∂ ux ω100 − ∂∂ uy ω200 + ω300 p p +β 1 + ||∇u||2 1 + ||∇u||2

= αI1 + β I2 .

We continue replacing in (14) the coplanarity equality ω 000 = αω 0 + β ω 00 and, in consequence, I3 with αI1 + β I2 . We obtain: det(L) = − ω10 (αω20 + β ω200 )I1 I2 − ω100 ω20 I1 (αI1 + β I2 ) + ω100 (αω20 + β ω200 )(I1 )2 +

(αω10 + β ω100 )ω20 I1 I2 + ω10 ω200 I1 (αI1 + β I2 ) − (αω10 + β ω100 )ω200 (I1 )2 = − αω10 ω20 I1 I2 − β ω10 ω200 I1 I2 − αω100 ω20 (I1 )2 − β ω100 ω20 I1 I2 +

αω100 ω20 (I1 )2 + β ω100 ω200 (I1 )2 + αω10 ω20 I1 I2 + β ω100 ω20 I1 I2 +

αω10 ω200 (I1 )2 + β ω10 ω200 I1 I2 − αω10 ω200 (I1 )2 − β ω100 ω200 (I1 )2 = 0.

Shape reconstruction of symmetric surfaces using Photometric Stereo

9

Now we want to remark why, from the point of view of the image data, this is a required condition in order to solve this PS formulation of the SfS problem. In particular we proved that, if we consider a third image I3 which is obtained by a light source ω 000 coplanar with ω 0 and ω 00 (that we use to obtain the first two images) we do not add any information. Rather, we will see that it is possible to compute such I3 “artificially”. Remark 1. The modelization of the orthographic SfS problem allows us to define the linear dependence of the images through their associated light vectors. In fact the way we consider the vector ω 000 coplanar with the other two allows us to define, in the same way, the concept of linearly dependent image functions. In order to give a better explanation of the linearly independent images and their importance for the resolution of the problem, we start giving the following definition: Definition 1. Three image functions I1 , I2 and I3 , obtained from the orthographic SfS model, are linearly independent if they are generated by three non coplanar light vectors. The idea behind this kind of definition uses the orthographic SfS model in order to compare the linear dependence/independence of vectors in R3 with respect to the image functions.

3 Linear dependent image reconstruction An interesting application regarding the orthographic SfS model is related to the possibility of generating the images of an object by varying the light source between two points ω 0 and ω 00 for which we already have the respective images. We suppose to have two images I1 and I2 obtained with two light sources ω 0 and ω 00 respectively, it is possible to deduce all the linearly dependent images. In the geometrical sense, if we consider these two light vectors like points of the upper hemisphere B((0, 0, 0), 1), then it is possible to calculate all the obtainable images with a light source in the geodesic identified by ω 0 and ω 00 . The main idea is related to the approximation of the points of a geodesic defined between the two light source vectors on the unit sphere. The method that we use to calculate the vector defined from the geodesic is based on the Remark 1. In fact, once we fix the extremal points, we discretize the geodesic through an affine transformation that allows us to see this curve as an arc of a circumference on the xy plane. Let us start by observing that the angle between the two unit light vectors θˆ is computable using the dot product between them. In fact: θˆ = arccos(ω 0 · ω 00 )

(23)

y

10

Roberto Mecca and Silvia Tozza

z

z ω!

ω !!

ω! θˆ

i=n

ω !!

ω !!

ω !! θˆ

(x1 , y1 )

x

j = 1, . . . , n

x

y

y

y

Fig. 2 Examples of geodesics (in blue) on the unit sphere between ω 0 and ω 00 . z

Now, we determine the affine transformation which allows the use of one angular parameter for the discretization of all the curve. For this purpose, we specify the linear mapping T like a matrix in R3x3 . In order to compute the coefficients of this matrix (i.e. the degree of freedom) we fix the following three conditions: T (1, 0, 0)T = ω 0 T (cos θˆ , sin θˆ , 0)T = ω 00

(24)

T

T (0, 0, 1) = q. y

z

z

ω!

ω

ω !!

!

ω !!

ω! θˆ θˆ

θˆ θˆ

x

x

x

θˆ

x

y

y

Fig. 3 Examples of affine transformations T between the area defined by the arc of the circumferz ence (in green) and the one z relates to the geodesic (in blue). !

We denote q asωthe unit vector (opportunely oriented) orthogonal to the plane where the vectors ω 0 and ω 00 lie (such that they are outgoing from the plane, i.e. q3 > 0). We consider the plane, generated by the vectors ω 0 and ω 00 , to have the following direction parameters:

Shape reconstruction of symmetric surfaces using Photometric Stereo

s = ω20 ω300 − ω30 ω200 l = ω30 ω100 − ω10 ω300

11

(25)

m = ω10 ω200 − ω20 ω100 . We consider the normalization of the vector (s, l, m) like a possible vector q depending on the sign of m. If m < 0 we take the opposite direction to the normalization. Let us make explicit the condition imposed by (24) writing the coefficients of the matrix that represents the transformation:

ω10

0 T = ω2 ω30

ω100 −ω10 cos θˆ sin θˆ ω200 −ω20 cos θˆ sin θˆ ω300 −ω30 cos θˆ sin θˆ

q1

q2 .

(26)

q3

If we consider an angle θ such that 0 < θ < θˆ , we obtain a third light source ω¯ using just the previous application: ¯ T (cos θ , sin θ , 0)T = ω.

(27)

Now, if we have built ω¯ coplanar to the previous two light vectors, then we have ω¯ = αω 0 + β ω 00 . The next step is to determine the coefficients α, β ∈ R. We write the previous equality and we obtain the following overdetermined linear system: 0 00 ω¯ 1 = αω1 + β ω1 ω¯ 2 = αω20 + β ω200 (28) ¯ 0 00 ω3 = αω3 + β ω3 With the aim to calculate the coefficients that determine the third light source, we solve the following linear system of two equations in two unknowns: α C =b (29) β where C ∈ R2×2 and clearly b ∈ R2 . Selecting two of the three equations of (28) in order to avoid singularity, the matrix C depends on the first two light vectors that we can express in spherical coordinates as follows: ω 0 = (cos θ1 sin ϕ1 , sin θ1 sin ϕ1 , cos ϕ1 ) and

ω 00 = (cos θ2 sin ϕ2 , sin θ2 sin ϕ2 , cos ϕ2 )

(with θ1 , θ2 ∈ [0, 2π] and ϕ1 , ϕ2 ∈ [0, π2 )). Our aim is to solve this overdetermined linear system considering only two of the three equations in order to avoid that one of the first two equations has all coefficients close or equal to zero (for example it

12

Roberto Mecca and Silvia Tozza

can happen if we have θ1 = θ2 = 0 with respect to the first equation or if we have θ1 = θ2 = π2 with respect to the second equation). Then, in order to avoid a bad conditioned linear system, if we consider the following set: π 3π 5π 7π ¯ A = 0, ∪ , ∪ , 2π (30) 4 4 4 4 we take

while

0 00 ω2 ω2 α ω¯ 2 = ω30 ω300 β ω¯ 3 0 00 ω1 ω1 α ω¯ 1 = ω30 ω300 β ω¯ 3

if

θ1 , θ2 ∈ A¯ otherwise.

(31)

(32)

Once we determine the coefficients α and β , we can compute the linear dependent image as follows: I¯ = αI1 + β I2 . (33) It is important to note that there is no approximation when we compute the linear ¯ In fact, the only computation that can introduce an error in the dependent image I. computation of the third image is related to the approximation of the linear system (29) for which we consider two different formulations ((31) and (32)) in order to make unimportant the error committed.

4 Reduction of the number of the images using symmetries Supposing that we are able to determine univocally the surface u(x, y) of the SfSPS problem using three images obtained with three non coplanar light sources, now we want to consider some classes of particular surfaces for which it is possible to resolve the SfS-PS problem using less than three images.

4.1 Symmetric surfaces We will prove a theorem that gives us the possibility to produce an image (of a particular surface) by using another image of the same surface, obtained from a light source appropriately chosen. Let us start to define the type of surface we can use in order to apply that. Definition 2. Let z = u(x, y) be a function defined in Ω ⊂ R2 . Let πs (symmetry plane) be a plane passing through the z axis and let rs (symmetry straight line) be its intersection with the xy plane. We say that u is symmetric with respect to the plane πs if, for every point (x0 , y0 ) ∈ Ω , it results:

Shape reconstruction of symmetric surfaces using Photometric Stereo

13

u(x0 , y0 ) = u(x1 , y1 )

(34)

where (x1 , y1 ) ∈ Ω is the symmetric point i =1, . . .of , n(x0 , y0 ) with respect to the straight line rs (Fig. 4). j=0

y rs

ω !!

(x0 , y0 ) α Fig. 4 Orthogonal view of the plane (x, y): example of a positioning of the points (x0 , y0 ) and (x1 , y1 ) with respect to the symmetry straight line rs .

(x1 , y1 )

x

With regard to the orthographic SfS problem, since the model has not a direct dependence on the value of the function u in a point, but only on its gradient, it is essential to determine a symmetry relation with respect to the partial derivatives of the surface u. We identify rs like a particular straight line passing through the origin of the reference system xyz, that is x sin α = y cos α, with α ∈ [0, π].

(35)

We consider this straight line like a subspace generated by the vector (cos α, sin α). It is possible to identify the orthogonal straight line to rs from its base vector. Then, let π π v = cos α + , sin α + = (− sin α, cos α) (36) 2 2 be this orthogonal straight line. The main relation about the symmetry that we use in the SfS model is the following: v · ∇u(x0 , y0 ) = −v · ∇u(x1 , y1 ).

(37)

4.2 Uniqueness theorem for the symmetric surfaces We aim at obtaining further information relative to the symmetric surfaces by deducing the grey scale values of an image corresponding to a particular light source

y

14

Roberto Mecca and Silvia Tozza

whose position depends on the symmetry straight line rs . These information, together with a third image, allow us to resolve the SfS-PS problem using three images linearly independent. Theorem 1. Let u(x, y) be a symmetric surface with respect to any straight line rs 00 be two light vectors such that: such that (0, 0) ∈ rs . Let ω 0 and ωded 1. 2.

e 0 is orthogonal with respect to the symmetry straight line rs ; ω 00 =−ω 00 = (−ω e e 0 (that is ωded e 0 , ω30 ) since the light vectors are constrained to ωded be taken in the upper hemisphere).

00 only by Then it is possible to deduce the image associated to the light source ωded I1 as follows: (38) I2ded (x0 , y0 ) = I1 (x1 , y1 ).

Proof. Considering the image definition of the SfS problem, we have: 00 · ∇u(x , y ) + ω 00 eded e 0 · ∇u(x0 , y0 ) + ω30 ∗ −ω ω 0 0 ded 3 ∗∗ p p = = 1 + ||∇u(x0 , y0 )||2 1 + ||∇u(x0 , y0 )||2 e 0 · ∇u(x1 , y1 ) + ω30 e 0 · ∇u(x1 , y1 ) + ω30 −ω −ω = p = p = I1 (x1 , y1 ) (39) 1 + ||∇u(x0 , y0 )||2 1 + ||∇u(x1 , y1 )||2

I2ded (x0 , y0 ) =

where in ∗∗ we use the assumption 2 of the theorem, whereas in ∗ , as the surface is symmetric, we exploit (37). In the last equality we consider furthermore that the normal vector in the point (x0 , y0 ) has the same norm as in (x1 , y1 ). t u Let us explain with an example the previous theorem considering the symmetric surface of Fig. 5 that we will use again in the numerical test section. n3

ω�

SR ω

Pˆ n2

Pˆ � χ

R χ

n2

n1

ω

n(x, y) (x, y, u(x, y))

(a)

(b)

y

Fig. 5 On the left the 3D shape of the symmetric surface taken into account. On the right the orthographic view of the surface is shown in order to make clear the symmetry (the red line). The color is proportional to surface height, starting from the cold color (for the lowest values) to the hotter color.

Shape reconstruction of symmetric surfaces using Photometric Stereo

15

It is clear that its symmetry straight line rs drown in Fig. 5 (b) defined by the angle α = 3π 4 , permits us to apply the previous theorem starting with a first light e 0 orthogonal to rs . For example we source ω 0 having its projection on the xy plane ω π can choose θ1 = 4 and, with the only aim to avoid shadows, we choose ϕ1 = 0.1 obtaining the image I1 shown in Fig. 6 (on the left). Theorem 1 says that we can deduce (and not approximate) a second image I2ded shown in Fig. 6 (on the right) 00 = (−ω 0 , ω 0 ). In this example the angles that associated to the light source ωded 3 ded = ϕ . 00 define the polar coordinates of ωded are θ2ded = θ1 + π = 5π 1 4 and ϕ2 As we will remark the deduced image will help to have the sufficient information (about the three linear independent images) in order to solve the SfS-PS problem using less than three images.

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

00 : ϕ ded = 0.1, θ ded = ωded 2 2

5π 4

→ (I2ded )

Fig. 6 On the left the starting image obtained by the light source ω 0 and on the right the deduced 00 . image computed using the light source ωded

16

Roberto Mecca and Silvia Tozza

Once we choose (x0 , y0 ) ∈ Ω , the problem consists in determining the coordinates of the symmetric point (x1 , y1 ) with respect to rs . The point (x1 , y1 ) belongs to the straight line rs⊥ , orthogonal to rs and passing through the point (x0 , y0 ), whose equation is (x − x0 ) cos α = −(y − y0 ) sin α. (40) Finally, observing that the symmetric straight line rs passes through the origin, it results necessarily that the distances with respect to (0, 0) of both points are the same.

i =1, . . . , n j=0

y

y rs

(x0 , y0 )

(x0 , y0 ) α

(x1 , y1 )

x Fig. 7 In this schematic example is shown the circumference containing the two symmetric points (x0 , y0 ) and (x1 , y1 ) with respect to the straight line rs .

That is, both points belong to the same circumference centered in the origin of the axis as shown in Fig. 7. This formulation of the problem allows that the following non-linear system ( (x − x0 ) cos α = −(y − y0 ) sin α (41) x2 + y2 = x02 + y20 admits only two solutions, that is (x0 , y0 ) and (x1 , y1 ). The procedure we use to solve the system (41) is related to an algorithm that takes into account the numerical instability. To this purpose, some particularities appear. For first in the computation of the solutions of (41) we have to take into account the angle α. In order to avoid numerical instability, we consider two possible sets to which this angle could belong: π 3π π 3π A = 0, ∪ , π and B = , . (42) 4 4 4 4 Let us suppose that α ∈ A; we use the first equation of (41), and we write the first coordinate of the point (since the division by cos α does not involve numerical

Shape reconstruction of symmetric surfaces using Photometric Stereo

17

instability). We have then: x = (y0 − y)

sin α + x0 cos α

which, replaced in the second equation of (41), gives a second order equation (in the general case written as ay2 + by + c = 0) where the coefficients are: a = 1 + (tan α)2 b = −2(y0 (tan α)2 + x0 tan α)

(43)

c = y20 (tan α)2 + 2y0 x0 tan α − y20 .

The two solutions of this equation are y0 and y1 . It is possible to determine the exact solution since y0 is known. In the same way, if α ∈ B it is possible to divide by sin α. We write then the second coordinate of the point: y = (x0 − x)

cos α + y0 sin α

and, by substitution, we solve the second order equation in the variable x with the following coefficients: a = 1 + (cot α)2 b = −2(x0 (cot α)2 + y0 cot α) c

(44)

= x02 (cot α)2 + 2x0 y0 cot α − x02 .

The choice of the right solution can be done exactly as explained before, that is excluding the value x0 already known. In the end, if we consider an image I1 obtained with a light source that respects the hypotheses of the Theorem 1, it is possible to deduce another image I2 in the way we explained before. Corollary 1. Let u(x, y) be a symmetric surface with respect to a straight line rs passing through the point (0, 0). e 0 , ω30 ) such that ω e 0 is 1. Let I1 be an image obtained with a light vector ω 0 = (ω orthogonal to the symmetry straight line of u(x, y), e 000 , ω3000 ) such that 2. let I3 be an image obtained with a light source ω 000 = (ω 0 000 0 000 0 000 e ·ω e 6= ±||ω e ||||ω e || (that is ω e and ω e have not the same direction). ω Then it is possible to univocally determine the surface u(x, y) using the SfS-PS model. The use of the image I2 defined in the Theorem 1 is necessary to the application of the reconstruction technique using three images obtained with non coplanar light vectors.

18

Roberto Mecca and Silvia Tozza

Other simplifications in the calculus of the point (x1 , y1 ) are possible when the surface u(x, y) has a particular kind of symmetry. As a first particular case we consider now a symmetry with respect to the orthogonal straight line: let rs0 and rs00 be two symmetry straight lines of u orthogonal between them. Thus we have the following relation: ∇u(x0 , y0 ) = −∇u(−x0 , −y0 ) (45) for which we note the facility of calculus in the point previously called (x1 , y1 ) = (−x0 , −y0 ).

4.3 Surfaces with four symmetry straight lines A very interesting particular case is related to the surfaces with four symmetry straight lines. In relation to these surfaces, in fact, it is possible to solve the classical SfS problem considering only one image generated from any light source. For this case the constraint concerning the light sources direction vanishes. As before, we follow the procedure of reconstructing the missing information (that is the three linear independent images). Now, because of the particular surface geometry, we are able to obtain the necessary three images only from one. That is we can reconstruct the two missing images (generated with non coplanar light vectors) starting from only one image of the surface lighted up by any light source. Let us consider the image I1 and let ω 0 be its relative light source. The other two images can be calculated in the following way: e 0 , ω30 ) I2 (x, y) = I1 (−x, −y) with ω 00 = (−ω

I3 (x, y) = I1 (−y, x) with ω 000 = (ω20 , −ω10 , ω30 ).

(46)

Relating to the third image, we use (v1 , v2 ) · ∇u(x, y) = −(v2 , −v1 ) · ∇u(−y, x), where (v1 , v2 ) is any vector of the plane (x, y). The choice of the light placing, together with the way of determining the images, is not unique. The aim of this construction is to obtain linearly independent reconstructed images for which the uniqueness of solution for the SfS problem stands. We can adapt to this particular situation based on symmetries the uniqueness theorem in [13] that proves the uniqueness of a Lipschitz solution (i.e. a solution in the Sobolev space W 1,∞ (Ω )) relatively to the differential approach: Theorem 2. Let I(x, y) be an image of the surface u(x, y) ∈ W 1,∞ (Ω ), with four straight lines of symmetry obtained by the SfS model such that 0 < I(x, y) ≤ 1 ∀(x, y) ∈ Ω with the light source ω. Then it is possible to determine this surface univocally solving the SfS classic problem. The procedure followed until now is a constructive proof of the previous theorem and allows to solve, for a particular class of surfaces, the classical SfS problem passing through the SfS-PS one. Furthermore, our analysis of the problem allows us to give also a numerical method for the approximation of this unique weak solution.

Shape reconstruction of symmetric surfaces using Photometric Stereo

19

Remark 2. Following this procedure it is possible to prove the uniqueness of a weak solution of the classic SfS problem even the surface is radial by symmetries. Remark 3. For both cases (one or four axes) the construction of the other images can be done only if the light source of the first image is not vertical. This means that in this case, even using this image deduction technique, the concave/convex ambiguity cannot be solved.

5 Numerical Tests This section is divided in two parts. The first one is related to the computation of the linear dependent images explained in Section 3. In the second part we present some shape recovery tests relative to the local and differential method for the SfSPS using symmetric surfaces in order to compute the error of both methods using less than three images.

5.1 Numerical computation of linear dependent images The tests on the linear independent images concept is carried out on synthetic and real surfaces. For the synthetic case we start using the surface (with the respective initial images) shown in Fig. 8. In Fig. 9 it is possible to see the images of the same surface for which the light sources are equidistant in the geodesic defined by ω 0 and ω 00 .

3D surface

ω 0 : ϕ1 = 0.26, θ1 = 0.0 → (I1 ) ω 00 : ϕ2 = 0.19, θ3 = π → (I2 )

Fig. 8 From left to right: the 3D shape of the surface and the starting images with the respective light sources.

For the real case we start from the images shown in Fig. 10, namely I1real and In Fig. 11 it is possible to see three linear dependent images obtained with the

I2real .

20

Roberto Mecca and Silvia Tozza

I1

I2

Fig. 9 Between I1 and I2 there are three reconstructed images with the intermediate light sources taken equidistant in the geodesic.

algorithm explained in the previous section. Also in this case we are considering 0 00 . three light sources equidistant in the geodesic defined by ωreal and ωreal

00 : ϕ = 0.281, θ = 3.502 → (I real ) 0 : ϕ1 = 0.263, θ1 = −0.305 → (I1real ) ωreal ωreal 2 2 2

Fig. 10 Real starting data.

I1real

I2real

Fig. 11 The results of the computation using our algorithm. The three images on the middle are lighted up from equidistant light sources in the geodesic defined by the starting images.

Shape reconstruction of symmetric surfaces using Photometric Stereo

21

5.2 Shape reconstruction for symmetric surfaces Let us start the numerical tests by studying several surfaces. Every single one is of a typology that allows us to apply each variant of the previous results explained. These surfaces have symmetries that allows us to apply all the previous techniques that solve the SfS-PS problem.

vsymm (x, y)

vlip (x, y)

vrad (x, y)

Fig. 12 Set of surfaces used for the numerical tests.

The different typologies of the three surfaces represented in Fig. 12 are such that their geometrical characteristics differentiate the number of necessary images for their tridimensional reconstruction. In fact: • vsymm (x, y), used also in the previous numerical test, is a surface with one symmetry axis (the straight line y + x = 0) for which we need two images; • vlip (x, y) is a Lipschitz surface with four axis of symmetry (that is y = ±x, y = 0, x = 0). For its reconstruction we need only one image; • vrad (x, y) is a particular case of the previous type of surfaces. In this radial case the curves where it is not differentiable are concentric circles. For these starting data, summarized in Fig. 13, we make different tests modifying the size of the grid ∆ , uniform for all the domain Ω . In order to solve the differential problem (8), once boundary condition g(xi , y j ) = u(xi , y j ) ∀(xi , y j ) ∈ ∂ Ωd is computed, we use the two following semi-lagrangian fix point schemes [13]: n un+1 i, j = u (xi − hγ1 (xi , y j ), y j − hγ2 (xi , y j )) +

fi, j h, ||bi, j ||

∀(xi , y j ) ∈ Ωd

(47)

n un+1 i, j = u (xi + hγ1 (xi , y j ), y j + hγ2 (xi , y j )) −

fi, j h, ||bi, j ||

∀(xi , y j ) ∈ Ωd

(48)

and

22

Roberto Mecca and Silvia Tozza

i =1, . . . , n j=1

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

j=0

ω 00 : ϕ2 = 0.1, θ2 = 34 π → (I2 )

(x0 , y0 )

ω 0 : ϕ1 = 0.1, θ1 =

π 4

→ (I1 )

Fig. 13 Set of images used with the respective light sources described by their spherical coordinates. For each surface it is possible to see the straight line of symmetry (in red) and the curve where the surface is not differentiable (in blue). In the first line we have the only two images used to compute the shape of vsym . On the second line we have the only image needed to compute the surfaces vlip and vrad (respectively from left to right).

b(x,y) , fi, j is the discrete version of (10) and where γ(x, y) = (γ1 (x, y), γ2 (x, y)) = ||b(x,y)|| h is a small parameter that discretizes the directional derivative. Using a bilinear interpolation with respect to the points of the lattice Ωd at point (xi ±hγ1 (xi , y j ), y j ± hγ2 (xi , y j )) ∈ / Ωd , we can reach the first order of consistency choosing h = ∆ . The previous two semi-lagrangian numerical schemes will be respectively called in the sequel as forward (47) and backward (48). For these iterative schemes we use the same initial guess, that is u0i, j = 0 in Ωd and u0i, j = gi, j in ∂ Ωd . The numerical approximation is carried out with both semi-lagrangian schemes and we calculate the errors in L∞ (Ωd ) norm by varying the uniform discretization step ∆ . In particular we consider a domain Ω d = [−1, 1]2 using synthetic images of size starting from 100 × 100 pixels (∆ = 0.02) to 800 × 800 pixels (∆ = 0.0025). In all the following tables we compute also the convergence order with the aim to verify the convergence of the methods. In the last column we show the error committed using the approximation by the local normal approach [4] always in the same L∞ (Ωd ) norm. In all the tables it is possible to see how both the numerical schemes work well for all the cases considered. In particular it is possible to see how the semi-lagrangian iterative schemes have an order of about one for all the tests. This means that the

Shape reconstruction of symmetric surfaces using Photometric Stereo

23

∆ SL forward error SL forward order SL backward error SL backward order 0.02 5.254 × 10−2 6.856 × 10−2 1.0055 1.0021 0.01 2.618 × 10−2 3.428 × 10−2 1.0033 1.0008 0.005 1.307 × 10−2 1.713 × 10−2 1.0011 0.9992 −3 −3 0.0025 6.530 × 10 8.570 × 10

Local Normal 1.63 × 10−2 1.74 × 10−2 6.75 × 10−2 1.289 × 10−1

∆ SL forward error SL forward order SL backward error SL backward order 0.02 2.731 × 10−2 2.731 × 10−2 0.8975 0.8975 −2 0.01 1.466 × 10 1.466 × 10−2 0.8630 0.8630 −3 −3 0.005 8.060 × 10 8.060 × 10 0.9165 0.9165 −3 −3 0.0025 4.270 × 10 4.270 × 10

Local Normal 1.01 × 10−2 5.1 × 10−3 3.0 × 10−3 4.32 × 10−2

SL forward error SL forward order SL backward error SL backward order 2.5491 × 10−1 3.7057 × 10−1 1.0313 1.0423 1.2472 × 10−1 1.7993 × 10−1 1 0.9833 −2 −2 6.2360 × 10 9.1010 × 10 0.9963 0.9878 3.1280 × 10−2 4.5890 × 10−2

Local Normal 1.136 × 10−1 4.53 × 10−2 2.0 × 10−2 8.3 × 10−3

Table 1 Error in L∞ (Ωd ) norm and order of convergence for vsym with angles: θ1 = ϕ1 = ϕ2 = 0.1.

Table 2 Error in L∞ (Ωd ) norm and order of convergence for vlip with angles: θ1 =

∆ 0.02 0.01 0.005 0.0025

Table 3 Error in L∞ (Ωd ) norm and order of convergence for vrad with angles: θ1 =

π 4

π 4

π 4

and θ2 = 34 π;

and ϕ1 = 0.1.

and ϕ1 = 0.1.

approximation of the symmetric surfaces can be done without loss of precision considering less than three images.

6 Conclusion and perspectives We have seen how it is possible to resolve the orthographic SfS problem with the PS technique uniquely, getting over the concave/convex ambiguity. We have shown how the linear independence of images is fundamental to obtain a unique reconstruction of the surface. Also an interesting application is explained: the possibility to generate linear dependent images of an object changing the light source between two points ω 0 and ω 00 for which we already have the relative images. The numerical tests give good results for both the two approaches: normal vector field approach and differential one. The main perspective is to apply these two results to a weaker model that take into account a less number of assumptions. For example, the well-know perspective SfS model [18, 19, 17], represented by the following equation e · ∇u(x)) − ω3 (x · ∇u(x) + u(x)) f (ω p = I(x) f 2 |∇u(x)|2 + (x · ∇u(x) + u(x))2

24

Roberto Mecca and Silvia Tozza

(where f is the focal distance), contains some degree of ill-posedness. A further step forward can be done solving the perspective problem with the same technique explained before instead of the orthographic projection assumption. With the aim to eliminate the Lambertian surface assumption the specular model can be considered (see [16, 21, 1, 22, 23] ).

References 1. J. F. Blinn. Models of light reflection for computer synthesized pictures. Computer Graphics, 11(2):192–198, 1977. 2. M. J. Brooks, W. Chojnacki, and R. Kozera. Impossible and ambiguous shading patterns. International Journal of Computer Vision, 7(2):119–126, 1992. 3. T.-J. Cham and R. Cipolla. Symmetry detection through local skewed symmetries. Image and Vision Computing, 13(5):439 – 450, 1995. 4. J.-D. Durou, J.-F. Aujol, and F. Courteille. Integrating the Normal Field of a Surface in the Presence of Discontinuities. In Proceeding of Energy Minimization Method in Computer Vision and Pattern Recognition 2009, volume 5681 of Lecture Notes in Computer Science, pages 261–273. 5. J.-D. Durou and F. Courteille. Integration of a normal field without boundary condition. In Proceedings of the First International Workshop on Photometric Analysis For Computer Vision, 2007. 6. J.-D. Durou, M. Falcone, and M. Sagona. Numerical methods for shape from shading: a new survey with benchmarks. Computer Vision and Image Understanding, 109(1):22–43, 2008. 7. S. A. Friedberg. Finding axes of skewed symmetry. Computer Vision, Graphics, and Image Processing, 34(2):138 – 155, 1986. 8. Y. Gofman and N. Kiryati. Detecting symmetries in grey level images: The global optimization approach. In Proceedings of the International Conference on Pattern Recognition, 1996. 9. G.G. Gordon. Shape from symmetry. In Proceedings of SPIE Intellent Robots and Computer Vision VIII, 1989. 10. G.G. Gordon. 3d reconstruction from a single view of an object ant its image in a plane mirror. In Proceedings of International Conference on Pattern Recognition, 1998. 11. R. Kozera. Existence and uniqueness in photometric stereo. Appl. Math. Comput., 44(1):103, 1991. 12. R. Kozera. Uniqueness in shape from shading revisited. Journal of Mathematical Imaging and Vision, 7(2):123–138, 1997. 13. R. Mecca. Uniqueness for shape from shading via photometric stereo technique. In Proceedings of the 18th IEEE International Conference on Image Processing, pages 2994–2997, 2011. 14. D.P. Mukherjee, A. Zisserman, and J.M. Brady. Shape from symmetry - detecting and exploiting symmetry in affine images. Technical Report, Report No. OUEL 1988/93. University of Oxford, Department of Engineering Science. 15. R. Onn and A. M. Bruckstein. Integrability disambiguates surface recovery in two-image photometric stereo. International Journal of Computer Vision, 5(1):105–113, 1990. 16. B. T. Phong. Illumination for computer generated pictures. Communications of the ACM, 18(6):311–317, 1975. 17. E. Prados, F. Camilli, and O. Faugeras. A viscosity solution method for shape-from-shading without image boundary data. Math. Model. Numer. Anal. 40, pages 393–412, 2006. 18. E. Prados and O. Faugeras. Perspective shape from shading and viscosity solutions. In Proceedings of the 9th IEEE International Conference on Computer Vision (vol. 2), pages 826– 831, 2003.

Shape reconstruction of symmetric surfaces using Photometric Stereo

25

19. E. Prados and O. Faugeras. Shape from shading: a well-posed problem? In Proceedings of the IEEE Int. Conf. Comp. Vis. and Patt. Recog (vol. 2), pages 870–877, San Diego, California, USA, 2005. 20. I. Shimshoni, Y. Moses, and M. Lindenbaum. Shape reconstruction of 3d bilaterally symmetric surfaces. International Journal of Computer Vision, 39(2):97–110, 2000. 21. K. E. Torrance and E. M. Sparrow. Theory for off-specular reflecion from roughened surfaces. Journal of the Optical Society of America, 57(9):1105–1114, 1967. 22. P.S. Tsai and M. Shah. Shape from shading using linear approximation. Image and Vision Computing Journal (Elsevier, Oxford), 12(8):487–498, 1994. 23. O. Vogel, M. Breuss, and J. Weickert. Perspective shape from shading with non-lambertian reflectance. Lecture Notes in Computer Science. G. Rigoll (Ed.): Pattern Recognition, Springer, Berlin, 5096:517–526, 2008.