Linear Differential Constraints for Photo-polarimetric Height Estimation
Silvia Tozza William A. P. Smith Dizhong Zhu Ravi Ramamoorthi Edwin R. Hancock







Linear Differential Constraints for Photo-polarimetric Height Estimation
    Silvia Tozza

Sapienza - Università di Roma


tozza@mat.uniroma1.it





William A. P. Smith

University of York


william.smith@york.ac.uk
Dizhong Zhu
 
University of York

dz761@york.ac.uk
Ravi Ramamoorthi

UC San Diego 
 
ravir@cs.ucsd.edu
Edwin R. Hancock
 
University of York

edwin.hancock@york.ac.uk

    December 30, 2023
=======================================================================================================================================================================================================================================================================================================================

empty




In this paper we present a differential approach to photo-polarimetric shape estimation. We propose several alternative differential constraints based on polarisation and photometric shading information and show how to express them in a unified partial differential system. Our method uses the image ratios technique to combine shading and polarisation information in order to directly reconstruct surface height, without first computing surface normal vectors. Moreover, we are able to remove the non-linearities so that the problem reduces to solving a linear differential problem. We also introduce a new method for estimating a polarisation image from multichannel data and, finally, we show it is possible to estimate the illumination directions in a two source setup, extending the method into an uncalibrated scenario. From a numerical point of view, we use a least-squares formulation of the discrete version of the problem. To the best of our knowledge, this is the first work to consider a unified differential approach to solve photo-polarimetric shape estimation directly for height. Numerical results on synthetic and real-world data confirm the effectiveness of our proposed method. 





§ INTRODUCTION


A recent trend in photometric <cit.> and physics-based <cit.> shape recovery has been to develop methods that solve directly for surface height, rather than first estimating surface normals and then integrating them into a height map. Such methods are attractive since: 1. they only need solve for a single height value at each pixel (as opposed to the two components of surface orientation), 2. integrability is guaranteed, 3. errors do not accumulate through a two step pipeline of shape estimation and integration and 4. it enables combination with cues that provide depth information directly <cit.>. In both photometric stereo <cit.> and recently in shape-from-polarisation (SfP) <cit.>, such a direct solution was made possible by deriving equations that are linear in the unknown surface gradient.

In this paper, we explore the combination of SfP constraints with photometric constraints (i.e. photo-polarimetric shape estimation) provided by one or two light sources. Photometric stereo with three or more light sources is a very well studied problem with robust solutions available under a range of different assumptions. Two source photometric stereo is still considered a difficult problem <cit.> even when the illumination is calibrated and albedo is known. We show that various formulations of one and two source photo-polarimetric stereo lead to the same general problem
(in terms of surface height), 
that illumination can be estimated and that certain combinations of constraints lead to an albedo invariant formulation. Hence, with only modest additional data capture requirements (a polarisation image rather than an intensity image), we arrive at an approach for uncalibrated two source photometric stereo.
We make the following novel contributions:




  * We show how to estimate a polarisation image from multichannel data such as from colour images, multiple light source data or both (Sec. <ref>).



  * We show how polarisation and photometric constraints (Sec. <ref>) can be expressed in a unified formulation (of which previous work <cit.> is a special case) and that various combinations of these constraints provide different practical advantages (Sec. <ref>).



  * We show how to estimate the illumination directions in two source photo-polarimetric data leading to an uncalibrated solution (Sec. <ref>).




 §.§ Related Work


The polarisation state of light reflected by a surface provides a cue to the material properties of the surface and, via a relationship with surface orientation, the shape. Polarisation has been used for a number of applications, including early work on material segmentation <cit.> and diffuse/specular reflectance separation <cit.>. However, there has been a resurgent interest <cit.> in using polarisation information for shape estimation.

Shape-from-polarisation   The degree to which light is linearly polarised and the orientation associated with maximum reflection are related to the two degrees of freedom of surface orientation. In theory, this polarisation information alone restricts the surface normal at each pixel to two possible directions. Both Atkinson and Hancock <cit.> and Miyazaki <cit.> solve the problem of disambiguating these polarisation normals via propagation from the boundary under an assumption of global convexity. Huynh <cit.> also disambiguate polarisation normals with a global convexity assumption but estimate refractive index in addition. These works all 

used a diffuse polarisation model while Morel <cit.> use a specular polarisation model for metals. Recently, Taamazyan <cit.> introduced a mixed specular/diffuse polarisation model. All of these methods estimate surface normals that must be integrated into a height map. Moreover, since they rely entirely on the weak shape cue provided by polarisation and do not enforce integrability, the results are extremely sensitive to noise.

Photo-polarimetric methods   There have been a number of attempts to combine photometric constraints with polarisation cues. Mahmoud <cit.> used a shape-from-shading cue with assumptions of known light source direction, known albedo and Lambertian reflectance to disambiguate the polarisation normals. Atkinson and Hancock <cit.> used calibrated, three source Lambertian photometric stereo for disambiguation but avoiding an assumption of known albedo. Smith <cit.> showed how to express polarisation and shading constraints directly in terms of surface height, leading to a robust and efficient linear least squares solution. They also show how to estimate the illumination, up to a binary ambiguity, making the method uncalibrated. However, they require known or uniform albedo. We explore variants of this method by introducing additional constraints that arise when a second light source is introduced, allowing us to relax the uniform albedo assumption. 
We also give an explanation for why the matrix they consider is full-rank except in a unique case. 
Recently, Ngo <cit.> derived constraints that allowed surface normals, light directions and refractive index to be estimated from polarisation images under varying lighting. However, this approach requires at least 4 lights. All of the above methods operate on single channel images and do not exploit the information available in colour images.

Polarisation with additional cues   Rahmann and Canterakis <cit.> combined a specular polarisation model with stereo cues. Similarly, Atkinson and Hancock <cit.> used polarisation normals to segment an object into patches, simplifying stereo matching. Stereo polarisation cues have also been used for transparent surface modelling <cit.>. Huynh <cit.> extended their earlier work to use multispectral measurements to estimate both shape and refractive index. Drbohlav and Sara <cit.> showed how the Bas-relief ambiguity <cit.> in uncalibrated photometric stereo could be resolved using polarisation. However, this approach requires a polarised light source. 
Recently, Kadambi <cit.> proposed an interesting approach in which a single polarisation image is combined with a depth map obtained by an RGBD camera. The depth map is used to disambiguate the normals and provide a base surface for integration. 





§ REPRESENTING POLARISATION INFORMATION

We place a camera at the origin of a three-dimensional coordinate system (Oxyz) in such a way that Oxy coincides with the image plane and Oz with the optical axis. In Sec. <ref> we propose a unified formulation for a variety of methods, all of which assume a) orthographic projection, b) known refractive index of the surface.
Other assumptions will be given later on, depending on the specific problem at hand. 
We  denote by 𝐯 the  viewer direction, by 𝐬 a general light source direction with 𝐯≠𝐬. We only require the third components of these unit vectors to be greater than zero (all the vectors belong to the upper hemisphere). We will denote by 𝐭 a second light source where required. 
We parametrise the unknown surface  height by the function z(𝐱), where 𝐱=(x,y) is an image location, and 
the unit normal to the 
surface at the point 𝐱 is given by:

    𝐧(𝐱) = 𝐧̂(𝐱)/|𝐧̂(𝐱)| = [-z_x, -z_y, 1]^T/√(1 + |∇ z(𝐱)|^2),
 
where 𝐧̂(𝐱) is the outgoing normal vector and z_x, z_y denotes the partial derivative of z(𝐱) w.r.t. x and y, respectively, so that ∇ z(𝐱) = (z_x, z_y). 
We now  introduce relevant polarization theory, describing how we can estimate a polarisation image from multichannel data. 




 §.§ Polarisation image


When unpolarised light is reflected by a surface it becomes partially polarised <cit.>.
A polarisation image can be estimated by capturing a sequence of images  in which a linear polarising filter in front of the camera lens  is rotated through a sequence of P≥ 3 different angles ϑ_j, j∈{1, …, P}. The measured intensity at a pixel varies sinusoidally with the polariser angle:

    i_ϑ_j(𝐱)=i_un(𝐱)  ( 1 + ρ(𝐱) cos(2ϑ_j-2ϕ(𝐱))).

The polarisation image is thus obtained by decomposing the sinusoid at every pixel location  into three quantities <cit.>: the phase angle, ϕ(𝐱), the degree of polarisation, ρ(𝐱), and the unpolarised intensity, i_un(𝐱).  
The parameters of the sinusoid can be estimated from the captured image sequence using non-linear least squares <cit.>, linear methods <cit.> or via a closed form solution <cit.> for the specific case of P=3, ϑ∈{ 0^∘, 45^∘, 90^∘}. 



 §.§ Multichannel polarisation image estimation



A polarisation image is usually computed by fitting the sinusoid in (<ref>) to observed data in a least squares sense. Hence, from P≥ 3 measurements we estimate i_un, ρ and ϕ. In practice, we may have access to multichannel measurements. For example, we may capture colour images (3 channels), polarisation images with two different light source directions (2 channels) or both (6 channels). Since ρ and ϕ depend only on surface geometry (assuming that, in the case of colour images, the refractive index does not vary with wavelength), then we expect these quantities to be constant over the channels. On the other hand, i_un will vary between channels either because of a shading change caused by the different lighting or because the albedo or light source intensity is different in the different colour channels. Hence, in a multichannel setting with C channels, we have C+2 unknowns and CP observations. If we use information across all channels simultaneously, the system is more constrained and the solution  will be more robust to noise. Moreover, we do not need to make an arbitrary choice about the channel from which we estimate the polarisation image. This idea shares something in common with that of Narasimhan <cit.>, though their material/shape separation was not in the context of polarisation.

Specifically, we can express the multichannel observations in channel c with polariser angle ϑ_j as

    i^c_ϑ_j(𝐱) =  
        i^c_un(𝐱)(1 + ρ(𝐱)cos(2ϑ_j-2ϕ(𝐱))).

The system of equations is linear in the unpolarised intensities and, by a change of variables, can be made linear in ρ and ϕ <cit.>. Hence, we wish to solve a bilinear system and do so in a least squares sense using interleaved alternating minimisation. Specifically, we a) fix ρ and ϕ and then solve linearly for the unpolarised intensity in each channel and b) then fix the unpolarised intensities and solve linearly for ρ and ϕ using  all channels simultaneously. Concretely, for a single  pixel, we obtain the unpolarised intensities across channels by solving:

    min_i^1_un(𝐱), …, i^C_un(𝐱) C_I [
      i^1_un(𝐱),
      …,
      i^C_un(𝐱)
    ]^T -  d_I
    ^2 ,

where C_I∈ℝ^CP× C is given by


    C_I =
    [  (1+ ρ(𝐱)cos(2ϑ_1-2ϕ(𝐱))) I_C;                             ⋮; ( 1+ ρ(𝐱)cos(2ϑ_P-2ϕ(𝐱))) I_C ],


with I_C  denoting the C× C identity matrix, 
and d_I∈ℝ^CP is given by


    d_I=
    [
      i^1_ϑ_1(𝐱),
      …,
      i^C_ϑ_1(𝐱),
      i^1_ϑ_2(𝐱),
      …,
      i^C_ϑ_P(𝐱)
    ]^T.


Then, with  the unpolarised intensities fixed, we solve for ρ and ϕ using the following linearisation:

    min_a,b C_ρϕ[ a; b ]-
     d_ρϕ^2,

where
[a    b]^T = [ρ(𝐱)cos(2ϕ(𝐱)), ρ(𝐱)sin(2ϕ(𝐱))]^T, and
C_ρϕ∈ℝ^CP× 2 is given by


    C_ρϕ=
        [ i^1_un(𝐱)cos(2ϑ_1) i^1_un(𝐱)sin(2ϑ_1);                  ⋮                  ⋮; i^1_un(𝐱)cos(2ϑ_P) i^1_un(𝐱)sin(2ϑ_P); i^2_un(𝐱)cos(2ϑ_1) i^2_un(𝐱)sin(2ϑ_1);                  ⋮                  ⋮; i^C_un(𝐱)cos(2ϑ_P) i^C_un(𝐱)sin(2ϑ_P);                    ],


and d_ρϕ∈ℝ^CP is given by:


    d_ρϕ=
        [ i^1_ϑ_1(𝐱)-i^1_un(𝐱);                    ⋮; i^1_ϑ_P(𝐱)-i^1_un(𝐱); i^2_ϑ_1(𝐱)-i^2_un(𝐱);                    ⋮; i^C_ϑ_P(𝐱)-i^C_un(𝐱) ].


We estimate  ρ and ϕ from the linear parameters using ϕ(𝐱)=1/2atan2(b,a) and ρ(𝐱)=√(a^2+b^2). 



We initialise by computing a polarisation image from one channel using linear least squares, as in <cit.>, and then use the estimated ρ and ϕ to begin alternating  interleaved optimisation by solving for the unpolarised intensities across channels.
We  interleave and alternate the two steps until convergence. In practice, we find that this approach not only dramatically reduces noise in the polarisation images but also removes the ad hoc step of choosing an arbitrary channel to process. We show an example of the results obtained in Figure <ref>. The multichannel result is visibly less noisy than the single channel performance.




§ PHOTO-POLARIMETRIC HEIGHT CONSTRAINTS
 

In this section we describe the different constraints provided by photo-polarimetric information and then  show how to combine them to arrive at linear equations in the unknown  surface height. 


 §.§ Degree of polarisation constraint

A polarisation image provides a constraint on the surface normal direction at each pixel. 
The exact nature of the constraint depends on the polarisation model used.
In this paper we will consider diffuse polarisation, due to subsurface scattering (see <cit.> for more details). 
The degree of diffuse polarisation ρ_d(𝐱) at each point 𝐱 can be expressed in terms of the refractive index η and the surface zenith angle θ∈ [0, π/2] as follows (Cf. <cit.>):

    ρ_d(𝐱) =  
          (η-1/η)^2 sin^2(θ)/2 + 2η^2 - (η+ 1/η)^2 sin^2(θ)  + 4 cos(θ) √(η^2 - sin^2 (θ)).

Recall that the zenith angle is the angle between the unit surface normal vector 𝐧(𝐱) and the viewing direction 𝐯. If we 
 know the degree of polarisation ρ_d(𝐱) and the refractive index η (or have good  estimates of  them at hand), 
equation (<ref>) can be rewritten  
with respect to the cosine of the zenith angle, and expressed in terms of the function, f(ρ_d(𝐱),η), that depends on the measured 
degree of polarisation and the refractive index:


    cos(θ) = 𝐧(𝐱)·𝐯 = f(ρ_d(𝐱),η) =  
       √(η^4 (1-ρ_d^2) + 2η^2 (2ρ_d^2 +ρ_d -1) +ρ_d^2 + 2ρ_d -4 η^3 ρ_d √(1-ρ_d^2)+1/(ρ_d+1)^2   (η^4 + 1) + 2η^2(3ρ_d^2 + 2ρ_d -1))


where we drop the dependency of ρ_d on (𝐱) for brevity.




 §.§ Shading constraint

The unpolarised intensity provides an additional constraint on the surface normal direction via an appropriate reflectance model. We assume that pixels have been labelled as diffuse or specular dominant and restrict consideration to diffuse shading. In practice, we deal with specular pixels in the same way as <cit.> and simply assume that they point in the direction of the halfway vector between 𝐬 and 𝐯. For the diffuse pixels, we therefore assume that light is reflected according to the Lambert's law. 
Hence, the unpolarised intensity is related to the surface normal by:

    i_un(𝐱) = γ(𝐱) cos(θ_i) = γ(𝐱) 𝐧(𝐱) ·𝐬,

 where 
 γ(𝐱) is the albedo.  
Writing 𝐧(𝐱) in terms of the gradient of z as reported in (<ref>), (<ref>) can be rewritten as follows:

    i_un(𝐱) = γ(𝐱) -∇ z(𝐱)·𝐬̃ + s_3/√(1 + |∇ z(𝐱)|^2),

with 𝐬̃ = (s_1, s_2). 
This is a non-linear equation, but we will see in Sec. <ref> and <ref> how it is possible to remove the non-linearity by using the ratios technique.



 §.§ Phase angle constraint

An additional  constraint comes from the phase angle, which determines the azimuth angle of the surface normal α(𝐱) ∈ [0, 2π] up to a 180^∘ ambiguity. This constraint can be rewritten as a collinearity condition <cit.>, that is satisfied by either of the two possible azimuth angles implied by the phase angle measurement. 
Specifically, for diffuse pixels we require the projection of the surface normal into the x-y plane, [n_x n_y], and a vector in the image plane pointing in the phase angle direction, [sin(ϕ) cos(ϕ)], to be collinear. This corresponds to requiring 

    𝐧(𝐱)· [cos(ϕ(𝐱)) -sin(ϕ(𝐱)) 0]^T = 0.

In terms of the surface gradient, using (<ref>), it is equivalent to 

    (-cosϕ,sinϕ) ·∇ z = 0.

A similar expression can be obtained for specular pixels, substituting in the π/2-shifted phase angles.
The advantage of doing this will become  clear in Sec. <ref>.



 §.§ Degree of polarisation ratio constraint

Combining the two constraints illustrated in Sec. <ref> and <ref>, we can arrive at a linear equation, that we refer to as the DOP ratio constraint. 
Recall that cos(θ) = 𝐧(𝐱) ·𝐯 and that we can express 𝐧 in terms of the gradient of z by using (<ref>), then isolating the non-linear term in (<ref>) we obtain

    √(1 + |∇ z(𝐱)|^2) = -∇ z(𝐱)·𝐯̃ + v_3/f(ρ_d(𝐱),η),

where  𝐯̃ = (v_1, v_2). 
On the other hand, considering the shading information contained in (<ref>), and again isolating the non-linearity  
we arrive at the following

    √(1 + |∇ z(𝐱)|^2) = γ(𝐱) -∇ z(𝐱)·𝐬̃ + s_3/i_un(𝐱).

Note that we are supposing 𝐬𝐯, 
and i_un(𝐱)≠ 0, f(ρ_d(𝐱),η) ≠ 0. 
Inspecting  Eqs. (<ref>) and (<ref>) we obtain

    -∇ z(𝐱)·𝐯̃ + v_3/f(ρ_d(𝐱),η) = γ(𝐱) -∇ z(𝐱)·𝐬̃ + s_3/i_un(𝐱).

We thus  arrive at  the following partial differential equation (PDE):

    𝐛(𝐱)·∇ z(𝐱)=h(𝐱),

where

    𝐛(𝐱) := 
    𝐛^(f,i_un) = 
    i_un(𝐱) 𝐯̃ - γ(𝐱)
    f(ρ_d(𝐱),η)  𝐬̃,

and

    h(𝐱) := 
    h^(f,i_un)=
    i_un(𝐱)v_3 - γ(𝐱)
    f(ρ_d(𝐱),η)   s_3.
 



 §.§ Intensity ratio constraint


Finally, we construct an intensity ratio constraint by considering two unpolarised images, i_un,1,i_un,2, taken from two different light source directions, 𝐬,𝐭. 
We construct our constraint equation by applying (<ref>) twice, once for each light source. We can remove the non-linearity as before and take a ratio, arriving at the following equation:

    i_un,2(-∇ z(𝐱)·𝐬̃ +s_3) = i_un,1(-∇ z(𝐱)·𝐭̃ + t_3).

The above  equation is  independent of albedo, light source intensity and non-linear normalisation term.
Again as before, we can rewrite (<ref>) as a PDE in the form  of (<ref>) with

    𝐛(𝐱) :=
    𝐛^(i_un,1,i_un,2) = 
    i_un,2(𝐱) 𝐬̃ - i_un,1(𝐱)  𝐭̃,

where 𝐭̃ = (t_1, t_2),  
and

    h(𝐱) := 
    h^(i_un,1,i_un,2) = 
    i_un,2(𝐱)s_3 -
    i_un,1(𝐱)   t_3.
 



§ A UNIFIED PDE FORMULATION


Commencing  from the constraints introduced in Sec. <ref>, in this section we show how to solve  several different  problems in photo-polarimetric shape estimation. The common feature is that these are all linear in the unknown height, and are expressed in a unified formulation in terms of 
a system of PDEs in the same general form: 

    𝐁(𝐱) ∇ z(𝐱)=𝐡(𝐱),

where 𝐁: Ω̅→ℝ^J× 2, 𝐡: Ω̅→ℝ^J× 1, denoting by Ω the reconstruction domain and being J=2,3 or 4 depending on the cases. 
(<ref>) 
is a compact and general equation, suitable for  describing several cases in a unified differential formulation that  solves directly for surface  height. 

Different combinations of the three constraints described in Sec. <ref> that are linear in the surface gradient can be combined in the formulation of (<ref>). Each corresponds to different assumptions and have different pros and cons. We explore three variants and show that <cit.> is a special case of our formulation. We summarise the alternative formulations in Tab. <ref>. 





 §.§ Single light and polarisation formulation


This case has been studied in <cit.>. It uses a single polarisation image, requires known illumination (though <cit.> show how this can be estimated if unknown) and assumes that albedo is known or uniform. This last assumption is quite restrictive, since it can only be applied to objects with  homogeneous surfaces. With just a single illumination condition, only the phase angle and DOP ratio constraints are available.
This  thus becomes a special case of our general unified formulation (<ref>), 
where 𝐁 and 𝐡 are defined as

    𝐁 =  
    [ b^(f,i_un)_1 b^(f,i_un)_2;        -cosϕ         sinϕ;              ],
      𝐡 = [h^(f,i_un), 0]^T,

with 𝐛^(f,i_un) and h^(f,i_un) defined by (<ref>) and (<ref>), with uniform γ(𝐱) and 𝐯 = [0,0,1]^T. 



 §.§ Proposed 1: Albedo invariant formulation


Our first proposed method uses
the phase angle constraint (<ref>) and two unpolarised images, taken from two different light source directions, 
obtained through (<ref>) and combined as in (<ref>). 
In this case the problem studied  is described by the system of PDEs (<ref>) with


    𝐁(𝐱) = 
    [ b^(i_un,1,i_un,2)_1 b^(i_un,1,i_un,2)_2;               -cosϕ                sinϕ;                     ], 𝐡(𝐱)= [ h^(i_un,1,i_un,2);                 0 ],


where 𝐛^(i_un,1,i_un,2) and h^(i_un,1,i_un,2) defined as in (<ref>) and (<ref>). 
The phase angle does not depend on albedo and the intensity ratio constraint is invariant to albedo. As a result, this formulation is particularly powerful because it allows albedo invariant height estimation. Moreover, the light source directions in the two images can be estimated (again, in an albedo invariant manner) using the method in Sec. <ref>. 

Once surface height has been estimated, we can compute the surface normal at each pixel and it is then straightforward to estimate an albedo map using (<ref>). Where we have two diffuse observations, we can compute albedo from two equations of the form of (<ref>) in a least squares sense. In real data, where we have specular pixel labels, we use only the diffuse observations at each pixel. To avoid artifacts at the boundary of specular regions, we introduce a gradient consistency term into the albedo estimation. We encourage the gradient of the albedo map to match the gradients of the intensity image for diffuse pixels.




 §.§ Proposed 2: Phase invariant formulation


Our second proposed method uses only the  DOP ratio and the intensity ratio constraints. This means that phase angle estimates are not used. The advantage of this is that phase angles are subject to a shift of π/2 at specular reflections when compared to diffuse reflections. So, the phase angle constraint relies upon having accurate per-pixel specularity labels, which  classify reflections as either dominantly  specular  or diffuse  (or  alternatively use a mixed polarisation model <cit.> with a four way ambiguity). 
In this case we need a) two unpolarised intensity images, taken with two different 
light source directions, 𝐬 and 𝐭, obtained through (<ref>), b)   polarisation information  from the function f(ρ,η) and c) knowledge of the albedo map. We need 𝐬, 𝐭, 𝐯 non-coplanar in order to have the matrix field 𝐁 not singular. 
Note that the function f, obtained from polarization information (as in (<ref>)), is the same for the two required   images. The reason for this is that it  does not depend on the light source directions but only on the viewer direction  𝐯  which does not change. 
This formulation can be deduced starting from (<ref>) and (<ref>),
arriving at a PDE system as in (<ref>) with

    𝐁 =[𝐛^(f,i_un,1) , 𝐛^(f,i_un,2) , 𝐛^(i_un,1,i_un,2)]^T,

and
𝐡 =[h^(f,i_un,1) , h^(f,i_un,2) , h^(i_un,1,i_un,2)]^T,
using  (<ref>), (<ref>), 
 (<ref>), (<ref>) to define the vector fields 𝐛 and the scalar fields h that appear in 𝐁 and 𝐡.



 §.§ Proposed 3: Most constrained formulation


Our final proposed method combines all of the previous constraints, leading to a problem of  the form 
(<ref>)    
with

    𝐁=[      b^(f,i_un,1)_1      b^(f,i_un,1)_2;      b^(f,i_un,2)_1      b^(f,i_un,2)_2; b^(i_un,1,i_un,2)_1 b^(i_un,1,i_un,2)_2;               -cosϕ                sinϕ;                     ],
      𝐡 = 
    [      h^(f,i_un,1);      h^(f,i_un,2); h^(i_un,1,i_un,2);                 0;                   ].

This formulation uses the most information and so is potentially the most robust method. However, it requires known albedo in order to use the DOP ratio constraint. Nevertheless, it is possible to first apply proposed method 1, estimate the albedo and then re-estimate surface height using the maximally constrained formulation and the estimated albedo map. In fact, the best performance is obtained by iterating these two steps, alternately using the surface height estimate to compute albedo and then using the updated albedo to re-compute surface height.



 §.§ Extension to colour images


We now consider how to extend the above systems of equations when colour information is available. If a surface is lit by a coloured point source, then each pixel provides three equations of the form in (<ref>). In principle, this provides no more information than a grayscale observation since the surface normal and light source direction are fixed across colour channels. However, in the presence of noise using all three observations improves robustness. In particular, if the albedo value at a pixel is lower in one colour channel, the signal to noise ratio will be worse in that channel than the others. For a multicoloured object, it is impossible to choose a single colour channel that provides the best signal to noise ratio across the whole object. For this reason, we propose to use information from all colour channels where available. 

We already exploit colour information in the estimation of the polarisation image in Sec. <ref>. Hence, the phase angle estimates have already benefited from the improved robustness. Both the DOP ratio and intensity ratio constraints can also exploit colour information by repeating each constraint three times, once for each colour channel. In the case of the intensity ratio, the colour albedo once again cancels if ratios are taken between the same colour channels under different light source directions.







§ HEIGHT ESTIMATION VIA LINEAR LEAST SQUARES 

We have seen that each of  the  variants illustrated in the previous section, each with different advantages, can be written as a PDE system (<ref>). 
Denoting by M the number of pixels, we discretise the gradient in (<ref>) via finite differences, arriving at the following linear system in 𝐳

    𝐀𝐳 = 𝐡̅,

where 𝐀 = 𝐁̅𝐆, with 𝐆∈ℝ^2M× M the matrix of finite difference gradients.  𝐁̅∈ℝ^JM× 2M is the discrete per-pixel version of the matrix 𝐁(𝐱), hence 𝐀∈ℝ^JM× M, where J depends on the various proposed cases reported in Sec. <ref> (J=2 for (<ref>) and (<ref>), J=3 for (<ref>) and J=4 for (<ref>)). 𝐡̅ is the discrete per-pixel version of the function 𝐡(𝐱), 𝐡̅∈ℝ^JM× 1, and 𝐳∈ℝ^M× 1 the vector of the unknown height values. 
The resulting discrete system is large, since we have JM linear equations in M unknowns, but sparse, since 𝐀 has few non-zero values for each row, and has as unknowns the height values.  
The per-pixel matrix 𝐀 is a full-rank matrix, for each choice of 𝐁̅ that comes from the proposed formulations in Sec. <ref>, under the different assumptions specified for each case. 
The per-pixel matrix 𝐀 related to <cit.> is full-rank except in one case: when the first two components of the light vector 𝐬 are non-zero and s_1 = -s_2 and it happens that the phase angle is ϕ = π/4 at least in one pixel. In that case, the matrix has a rank-deficiency (though in practice ϕ assuming a value of exactly π/4, up to numerical tolerance, is unlikely).

We want to find a solution of  (<ref>) in the least-squares sense, find a vector 𝐳∈ℝ^M such that

    ||𝐀𝐳 - 𝐡̅||^2_2 ≤ ||𝐀𝐲 - 𝐡̅||^2_2,   ∀𝐲∈ℝ^M.

Considering the associated system of  normal equations

    𝐀^T (𝐀𝐳 - 𝐡̅) = 0,

it is well-known that if there exists 𝐳∈ℝ^M that satisfies (<ref>), then 𝐳 is also solution of the least-squares problem, 𝐳 satisfies (<ref>). 
Since 𝐀 is  a full-rank  matrix, then the matrix 𝐀^T 𝐀 is not singular, hence there exists a unique solution 𝐳 of (<ref>) for each data term 𝐡̅. 
Since neither 𝐁 nor 𝐡 depend on z in (<ref>), the solution can be computed only up to an additive constant (which is  consistent  with the orthographic projection assumption). 
To resolve the unknown constant, 
knowledge of z at just one pixel is sufficient.  
In our implementation, we remove the height of one pixel from the variables and substitute its zero value elsewhere.



§ TWO SOURCE LIGHTING ESTIMATION


Our three proposed shape estimation methods require knowledge of the two light source directions. Previously, Smith <cit.> showed that a single polarisation image can be used to estimate illumination conditions up to a binary ambiguity. 
However, to do so, they assumed that the albedo was known or uniform,  and they also worked only with a single colour channel. In a two source setting, we show that it is possible to estimate both light source directions simultaneously, and do so in an albedo invariant manner. Moreover, we can exploit information across different colour channels to improve robustness to noise. Hence, our three methods can be used in an uncalibrated setting.

The intensity ratio (<ref>) provides one equation per pixel relating unpolarised intensities, surface gradient and light source directions. Given two polarisation images with different light directions, we have one such equation per pixel and six unknowns in total.
We
assume that ambiguous surface gradient estimates are known from ρ and ϕ, and then  use (<ref>) to estimate the light source directions.



The intensity ratio (<ref>) is homogeneous in s and t and so has a trivial solution s= t=[0 0 0]^T. If we assume that the intensity of the light source remains constant in each colour channel across the two images, then this intensity divides out when taking an intensity ratio and so the length of the light source vectors is arbitrary. We therefore constrain them to unit length (avoiding the trivial solution), and represent them by spherical coordinates (θ_s,α_s) and (θ_t,α_t), such that [s_1,s_2,s_3]=[cosα_ssinθ_s,sinα_ssinθ_s,cosθ_s] and [t_1,t_2,t_3]=[cosα_tsinθ_t,sinα_tsinθ_t,cosθ_t].
This reduces the number of unknowns to four. We can now write the residual at each pixel given an estimate of the light source directions. There are two possible residuals, depending on which of the two ambiguous polarisation normals we use. From the phase angle and the zenith angle estimated from the degree of polarisation using (<ref>), we have two possible surface normal directions at each pixel and therefore two possible gradients: 
z_x(𝐱) ≈±cosϕ(𝐱)tanθ(𝐱),   z_y(𝐱) ≈±sinϕ(𝐱)tanθ(𝐱). 
Hence, the residuals at pixel 𝐱_j in channel c are given by either:


    r_j,c(θ_s,α_s,θ_t,α_t) = 
           i_un,1^c(𝐱_j)(-z_x(𝐱_j)t_1-z_y(𝐱_j)t_2+t_3) -
           i_un,2^c(𝐱_j)(-z_x(𝐱_j)s_1-z_y(𝐱_j)s_2+s_3),
       or
    
    q_j,c(θ_s,α_s,θ_t,α_t) = 
           i_un,1^c(𝐱_j)(z_x(𝐱_j)t_1+z_y(𝐱_j)t_2+t_3) -
           i_un,2^c(𝐱_j)(z_x(𝐱_j)s_1+z_y(𝐱_j)s_2+s_3).


We can now write a minimisation problem for light source direction estimation by summing the minimum of the two residuals over all pixels and colour channels:

    min_θ_s,α_s,θ_t,α_t∑_j,cmin [ r_j,c^2(θ_s,α_s,θ_t,α_t),q_j,c^2(θ_s,α_s,θ_t,α_t) ].

The minimum of two convex functions is not itself  convex and so this optimisation is non-convex. However, we find that, even with a random initialisation, it almost always converges to the global minimum. As in <cit.>, the solution is still subject to a binary ambiguity, in that if ( s, t) is a solution then ( Ts, Tt) is also a solution (with T=diag([-1,-1,1])), corresponding to the convex/concave ambiguity. We resolve this simply by choosing the maximal solution when surface height is later recovered.



§ EXPERIMENTS



#190#1


We begin by using synthetic data generated from the Mozart height map (Fig. <ref>). We differentiate to obtain surface normals and compute unpolarised intensities by rendering the surface using light sources s=[1,0,5]^T and t=[-1,-2,7]^T according to (<ref>). We experiment with both uniform albedo and varying albedo for which we use a checkerboard pattern. We simulate the effect of polarisation according to (<ref>), varying the polariser angle between 0^∘ and 180^∘ in 10^∘ increments. Next, we corrupt this data by adding Gaussian noise with zero mean and standard deviation σ, saturate and quantise to 8 bits. This noisy data provides the input to our reconstruction. First, we estimate a polarisation image using the method in Sec. <ref>, then apply each of the proposed methods or the state-of-the-art comparison method <cit.> to recover the height map. 

In Tab. <ref> we report Root-Mean-Square (RMS) error in the surface height (in pixels) and mean angular error (in degrees) in the surface normals obtained by differentiating the estimated surface height. In Fig. <ref> we show a sample of qualitative results from this experiment. In all cases, more than one 
of our proposed methods outperform <cit.>. When albedo is uniform, our phase invariant (Prop. 2) or maximally constrained solution (Prop. 3) provides the best results. When albedo is non-uniform, the albedo invariant method (Prop. 1) provides much better performance. Although the combination of the albedo invariant method followed by the maximally constrained method (Prop. 1+3) does not give quantitatively the best performance, we find that on real world data containing more complex noise and specular reflections, this approach is most robust.



In Fig. <ref> we show qualitative results on two real objects with spatially varying albedo. From left to right we show: an image from the input sequence; the surface normals of the estimated height map (inset sphere shows how orientation is visualised as colour); the estimated albedo map; a re-rendering of the estimated surface and albedo map under novel lighting with Blinn-Phong reflectance <cit.>; a rotated view of the estimated surface; and, for comparison, reconstructions of the same surfaces using <cit.>. The results of <cit.> are highly distorted in the presence of varying albedo. Our approach avoids transfer of albedo details into the recovered shape, leading to convincing relighting results.
 




§ CONCLUSIONS


In this paper we have introduced a unifying formulation for recovering height from photo-polarimetric data and proposed a variety of methods that use different combinations of linear constraints. We proposed a more robust way to estimate a polarisation image from multichannel data and showed how to estimate lighting from two source photo-polarimetric images. Together, our methods provide uncalibrated, albedo invariant shape estimation with only two light sources.  
Since our unified differential formulation does not depend on a specific camera setup or a chosen reflectance model, the most obvious target for future work is to move to a perspective projection, considering more complex reflectance models, exploiting better the information available in specular reflection and polarisation. 
In addition, since our methods directly estimate surface height, it would be straightforward to incorporate positional constraints, for example provided by binocular stereo. 



  §.§.§ Acknowledgements

This work was supported mainly by the “GNCS - INdAM”, in
part by ONR grant N000141512013 and the UC San Diego Center for Visual Computing. 
W. Smith was supported by 
EPSRC grant EP/N028481/1.





ieee



