tains 
ators. 
sured 
(3) 
1ated 
ed in 
ggio, 
‘well- 
(4) 
r f is 
rices. 
nizes 
d by 
d by 
noise 
1 the 
ıtion 
(8) 
e the 
7995 
247 
} Ca¥) pr (yy) Ca(V)p- P(v) 
ID: Pv v) + |ID4J?(v)) sm | us EC 4 | | Q(v) | 
Dz (v)P(v) + D;(v)Q(v) 
SNR; + 1Ps[?(v) + |Dy[*(v) 
  
(9) 
where C;(v) and C,(v) denote the expected power spectra of signal and noise. SNR is the signal-to-noise 
power ratio at each spatial frequency. These power spectra represent the diagonal elements of the covariance 
matrices in the Fourier basis since the covariance matrices describe autocorrelations and the Fourier transform 
of an autocorrelation is the power spectrum of the corresponding function (Pratt, 1991). 
Eq. 8 and Eq. 9 have a similar structure and, in fact for a particular choice of the SNR power ratio 
C+(v)/Cn(v) the two approaches become equivalent. In order to demonstrate this, let us now model the 
derivative operators and the SNR power ratio. 
The Fourier transforms of the partial derivative operators 0/0x and 0/0y, is (Pratt, 1991) 
D,(v)- —i27v, and Dy(v) = —i27y, (10) 
where 7 is the imaginary unit. 
Traditional implementations use discrete derivative operators in the spatial domain which may be represented 
as 6(1/2) —6(—1/2) so that first order derivatives are calculated between grid points and second order derivatives 
are calculated at grid points. The connection to discrete derivative operators is made by representing the Fourier 
transforms of the partial derivatives by 
D,(v) 2 —i2sin(vv,) and D,(v) = —i2sin(xvy). (11) 
Numerical experiments showed that the inverse filtering results obtained by using D, and D, are almost 
indistinguishable from the results obtained by using D, and Dy: 
The signal-to-noise power ratio is problem specific. Many real world image processing applications are 
approximatively invariant to orientation so that the power spectra are also approximatively rotation invariant. 
The decay-law of the high frequency components depends on the worst type of singularities that are expected 
to be found in the surface f. If the surface is continuous but contains step edges in the first derivatives, its 
power spectrum decays as 1/v*. If the surface contains step edges, its power spectrum decays as 1/v?. Typically, 
therefore, we can assume a power law 
1 1 > 
SNR(v) = Ju = = (+0) a) 
where p € [-2,0] is the exponent and c depends on the noise level. For a white noise power spectrum and a 
1/v* signal spectrum (surface without step edges but with sharp bends) we get p — —2. Thus, in this case, the 
Wiener filtering approach is equivalent to regularization using À1 — 0 and A5 = c. An arbitrary combination of 
À, and A» corresponds to a SNR model 1/SNR - 2, : |v|? + Xz |v]: 
In more general cases, Wiener filtering will lead to better results as soon as a good model of the SNR is 
available. Regularization, on the other hand, only makes implicit assumptions about the SNR. This aspect 
was ignored, so far. Usually, in fact, not even the selected regularization parameters A; and A; are reported in 
photometric stereo publications. 
3 IMPLEMENTATION ON A DISCRETE GRID 
Regularization and Wiener filtering, as discussed in Sec. 2, lead to especially simple and fast solutions in the 
Fourier domain. In order to construct a numerical solution, the desired functions have to be reconstructed 
on a discrete grid. Fast algorithms are obtained if the 2D Fourier transform can be realized as a discrete 2D 
Fourier transform with its fast implementation as a Fast Fourier Transform (FFT), assuming that the image 
dimensions are powers of 2. Using the FFT, however, implies periodic boundary conditions at the borders of 
the reconstructed surface, 4.e., f(z, y) = f(z,y + Ly) = f(z + Ly, y) where L, and L, are the x- and y-size of 
the height map. In practice, such an assumption is unrealistic and can lead to severe reconstruction errors. 
Much better results are obtained if the desired function f(z, y) is reconstructed within a larger height map 
F(z, y) of size 2 Lg by 2 L, with the condition f(z, v) - /(-z,y) = /(z,-y) = f2 Le — 2,9) = f(z, 2 L, — y) 
so that the f(z,y) is mirror symmetric at the original boundaries and the solution f(z, y) corresponds to the 
lower left quadrant of f (z,y). As a consequence of the symmetry conditions, f (2, y) is automatically periodic 
IAPRS, Vol. 30, Part 5W1, ISPRS Intercommission Workshop "From Pixels to Sequences", Zurich, March 22-24 1995