246 
where the vectors p, q, and f contain the discrete data values in scan-line order so that each vector contains 
L; L, data elements. The matrices D, and D, are discrete representations of the x- and y-derivative operators. 
The vectors n, and n, represent the noise, i.e. the unpredicatable deviation between the true and the measured 
derivative values. It is assumed that the noise is additive. 
Eqs. 1 and 2 can be written as a single matrix equation 
Big © 
where p and q were concatenated to a single vector with 2 L, L, elements, the pair n, and n, was concatenated 
accordingly, and the matrices D, and D, were combined to a single matrix of size 2 L, L, by L4 L,. 
D, 
D, 
The matrix in Eq. 3 is not square and therefore has no inverse. Traditionally, this problem was addressed in 
the frame-work of regularization theory, due mainly to Tikhonov (Tikhonov, 1977) and established by (Poggio, 
1985) as a general framework for solving early-vision problems. The stabilizing functional, favoring ‘well- 
behaved’ solutions, usually is 
MU + £304 Mal FE, 42 12, + 52) (4) 
where A; and A; are the regularization parameters. In the discrete version of photometric stereo the vector f is 
then determined by minimizing 
Hal LS: 
where the stabilizing functional is represented by the corresponding discrete matrices. Minimizing Eq. 5 by 
2 
t A(IDz£ + [Dy£[) + A2([D2D2f|? +2 [D,Dyf|? + [Dy Dy £7) (5) 
  
t 
using Eq. 15 (Appendix A) and | De | = [DI DÍ ] leads to 
y 
f= (DID, DID, 4 A(DÍD, 4 DÍD,) 
+ X(D!D!D.D. +2D!D!D.D, +D!D!D,D,)) [Di Di] | £ | (6) 
where 'i' denotes the conjugate matrix transpose which equals the ordinary matrix transpose for real matrices. 
The regularization approach makes no explicit reference to signal or noise statistics. 
Alternatively, we can take the Wiener filtering approach where a matrix W is constructed that minimizes 
E{|W | a | — f|?) so that the estimated vector is f = W | > | Assuming that p and q are corrupted by 
an equal amount of noise and that the noise vectors n, and n, are uncorrelated, the solution is obtained by 
applying Eq. 24 (Appendix B) 
w=c;[Di p}1(| 5e |er[ mi Di ]+ [5 on (7) 
where C; and C, are the covariance matrices of signal f and noise. 
Let us first ignore boundary effects which will be discussed in Sec. 3. Assuming that signal and noise 
are wide-sense stationary, the matrices C, and C, are block-Toeplitz and therefore diagonalizable with the 
2D discrete Fourier transform. The same holds for the matrices D, and D, since they represent convolution 
operations. Thus, all matrices are diagonal in the Fourier basis. 
In the Fourier basis Eq. 6 takes the form 
D;(v)P(v) + Dy (v)Q(v) 
(1+ A1)(I1Dz|*(v) + | Dy[2(v)) + A2(| Da |*(v) - 21Dz (v) DP (v) + | Dy |*(v)) 
  
F(v) (8) 
where P(v), Q(v), and F(v) denote the 2D Fourier transforms of p and q and f. D,(v) and D, (v) denote the 
Fourier transforms of the derivative operators and v — (v5, Vy) is the spatial frequency. 
Correspondingly, Wiener filtering (Eq. 7) leads to the 2 by 2 matrix equations 
ID.(v)- &&D  Dz(v)D,(v) "p Piv) 
Fw zin mm {| 270 em ml) [^ | 
IAPRS, Vol. 30, Part 5W1, ISPRS Intercommission Workshop "From Pixels to Sequences", Zurich, March 22-24 1995