# Ascertaining the Ideality of Photometric Stereo Datasets under Unknown Lighting
## 2 A Review of the Hayakawa Procedure for Determining the Light Position
Equations (3) can be represented in the form of the matrix equation  
\[DN^{T}L=M, \tag{4}\]  
after assembling data and unknowns in the matrices  
\[D =\text{diag}(\rho_{1},\ldots,\rho_{p})\in\mathbb{R}^{p\times p}, \qquad L =\big{[}\boldsymbol{\ell}_{1},\ldots,\boldsymbol{\ell}_{q}\big{]} \in\mathbb{R}^{3\times q},\] \[N =\big{[}\mathbf{n}_{1},\ldots,\mathbf{n}_{p}\big{]}\in\mathbb{R}^ {3\times p}, M =\big{[}\mathbf{m}_{1},\ldots,\mathbf{m}_{q}\big{]}\in\mathbb{R}^{p \times q}.\]  
If the matrix \(L\), whose columns are the vectors aiming at the lights, is known, one usually sets \(\tilde{N}=ND\) in (4), and computes \(\tilde{N}^{T}=ML^{\dagger}\), where \(L^{\dagger}\) denotes the Moore-Penrose pseudoinverse of \(L\)[27]. This produces a unique solution if \(q\geq 3\) and the matrix \(L\) is full-rank. Once \(\tilde{N}\) is available, by normalizing its columns, one easily obtains \(N\) and \(D\) from which, by integrating the normal vectors (see [28]), the approximation of the surface representation \(u(x,y)\) at the discretization grid can be computed. When the light positions are not known, the Hayakawa procedure [14] can be used to estimate this information from the data, but this requires the availability of at least six images with different lighting conditions. Here, we briefly summarize this technique (see [15] for more details). To obtain an initial Rank-3 factorization of the data matrix \(M\), we start computing the singular value decomposition (SVD) [29]  
\[M=U\Sigma V^{T}, \tag{5}\]  
where \(\Sigma=\text{diag}(\sigma_{1},\ldots,\sigma_{q})\) contains the singular values \(\sigma_{1}\geq\cdots\geq\sigma_{q}\geq 0\), and \(U\in\mathbb{R}^{p\times q}\), \(V\in\mathbb{R}^{q\times q}\) are matrices with orthonormal columns \(\mathbf{u}_{i}\), \(\mathbf{v}_{i}\). Theoretically, only the first three singular values should be positive, that is, \(\sigma_{4}=\cdots=\sigma_{q}=0\). Anyway, because of experimental errors propagation and lack of ideality in the data, one usually finds \(\sigma_{i}>0\), \(i=4,\ldots,q\). It is known that the best Rank-\(k\) approximation of a matrix with respect to both the 2-norm and the Frobenius norm is produced by truncating the SVD to \(k\) terms [27]. So, we define \(W=[\sigma_{1}\mathbf{u}_{1},\sigma_{2}\mathbf{u}_{2},\sigma_{3}\mathbf{u}_{3}] ^{T}\) and \(Z=[\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3}]^{T}\), to obtain the Rank-3 approximation \(W^{T}Z\) of \(M\). Following the proof of Theorem 1 from [15], we consider the factorization \(W^{T}Z=M\) as a tentative approximation of \(\bar{N}^{T}L=M\). Assuming \(\|\boldsymbol{\ell}_{t}\|=1\), we seek a matrix \(B\) such that  
\[\|B\mathbf{z}_{t}\|=1,\qquad t=1,\ldots,q, \tag{6}\]  
where \(\mathbf{z}_{t}\in\mathbb{R}^{3}\), \(t=1,\ldots,q\), denote the columns of \(Z\). The Hayakawa procedure consists of writing Problem (6) in the form  
\[\text{diag}(Z^{T}GZ)=\mathbf{1}, \tag{7}\]  
where \(\mathbf{1}=[1,\ldots,1]^{T}\in\mathbb{R}^{q}\) and \(G=B^{T}B\) is a symmetric positive definite \(3\times 3\) matrix, which depends upon six parameters, namely, the entries \(g_{ij}\) with \(i=1,2,3\) and \(j=i,\ldots,3\). Each equation of System (7) reads  
\[\mathbf{z}_{t}^{T}G\mathbf{z}_{t}=\sum_{i,j=1}^{3}z_{i}z_{j}g_{ij}=1.\]  
Such equations can be assembled in the linear system \(H\mathbf{g}=\mathbf{1}\), where \(H\) is a \(q\times 6\) matrix with rows  
\[[z_{1t}^{2}\quad z_{2t}^{2}\quad z_{3t}^{2}\quad 2z_{1t}z_{2t}\quad 2z_{1t}z_{3t} \quad 2z_{2t}z_{3t}],\qquad t=1,\ldots,q. \tag{8}\]  
The solution \(\mathbf{g}=[g_{11},g_{22},g_{33},g_{12},g_{13},g_{23}]^{T}\) is unique if the system is overdetermined, that is, if the matrix \(H\) is full-rank and \(q\geq 6\). This shows that a necessary condition for the position of the light sources to be uniquely determined is that the dataset contains at least six images. The factor \(B\) of \(G\) is determined up to a unitary transformation, so to simplify the problem, we represent \(B\) by its QR factorization and substitute \(B=R\).