248 
with f(z,y) = f(2 Lz + z,y) = f (2,2 L, -- y). These new boundary conditions lead to an exact reconstruction 
of f provided that, at the surface boundaries, the derivatives of f(z, y) vanish orthogonal to the boundaries. 
Although this assumption may still be invalid for a real surface it is much weaker since it concerns derivatives 
rather than direct values and since it leads to local constraints in contrast to the global constraints for a periodic 
surface. 
In order to get this idea to work, the original data must be enlarged accordingly. The x-derivatives p(z, y) 
must be converted to the enlarged data set 5(x, y) = —P(2 Le — x, y) = 9(x, 2 Ly — y) = —-9(2 D; — 2,2 L, — y) 
with twice the number of rows and columns of the original p. The minus sign is due to the fact that the 
derivative of a mirror-symmetric function is odd in z- and even in y-direction. Accordingly, the y-derivatives 
q(z, y) are converted to the data set §(z,y) = §(2 DL; — z,y) 2 —4(z,2 L, — y) 2 —4d(2 L; — 2,2 L, — y) with 
twice the number of rows and columns of the original q. 
Instead of calculating a direct FFT of these data sets, their particular symmetries can be exploited by 
calculating Discrete Cosine Transforms (DCT) for even functions and Discrete Sine Transforms (DST) for odd 
functions (Press, 1987). Thus, p is calculated with a 1D DST in x-direction and a 1D DCT in y-direction, q 
is calculated with a 1D DCT in x-direction and a 1D DST in y-direction, and the back-transform is calculated 
with two 1D DCTs. This implementation results in a speedup factor on the order of 8 since it exploits the 
symmetries in x and y as well as the fact that the data is real. The computational expense of all these algorithms 
is on the order of L; L,(log(L;) 4- log(L,)). Thus, they are inherently fast. 
The resulting algorithm for photometric stereo is presented in the following subsection. 
3.1 Surface reconstruction from noisy derivatives 
Let the noisy derivative images p(z, y) and q(z, y) have size L;, by L,. The reconstructed surface is obtained 
by applying the following steps: 
1. Construct the enlarged images p(z, y) and ÿ(x, y), with size 2 L, by 2 L,, from p(z, y) and from q(z, y) 
so that p is odd in z and d is odd in y. 
2. Calculate the 2D FFT of f(z, y) and q(z, y), using 1D DCT and 1D DST, and denote the results P(v) 
and Q(v). 
3. Apply Eq. 8 (regularization) or Eq. 9 (Wiener filtering) by using Eq. 10 for the derivative operators and 
by using a problem specific model of the SNR. Avoid divisions by zero by replacing zero values by a small 
€, such as € = 1071, 
4. Calculate the reconstructed surface f (2, y) by applying the inverse FFT in the lower left quadrant with 
size L, by L,. 
4 SIMULATION 
The proposed algorithms, as described in Sec. 3.1, were tested on synthetic and real data. The synthetic test 
object was a pyramid which contains both sharp edges and extended surfaces. The two derivative images p(z, y) 
and q(z, y) of a top view were calculated in an image of size 128 by 128. The surface reconstruction algorithms 
were tested on noiseless and noisy data. Finally, the algorithms were tested on real data taken of a dummy. 
4.1 Reconstruction from noiseless data 
The Wiener filtering reconstruction, using a low SNR (c — 1079, p = —2) is shown in Fig. 1(a) in a side 
view and shaded according to a Lambertian reflectance law. The choice of the exponent p is not important for 
noiseless data. The pyramid is reconstructed without artifacts. As discussed in Sec. 2, Wiener filtering with 
these parameters corresponds to regularization using A; = 0 and A; = 1078, 
4.2 Reconstruction from noisy data 
In a second experiment, the derivative images p(z, y) and q(z, y) were contaminated with Gaussian white noise 
of variance 1. The Wiener filtering reconstruction, using (c — 1, p — —2) is displayed in Fig. 1(b) and shows 
good reconstruction quality. The SNR exponent p was set to -2 since the desired depth map contains only edges 
in the derivatives and the noise is white. 
Ignoring noise affects the restoration as shown in Fig. 1(c) which is reconstructed using c — 1073and p = —2. 
Ignored noise cannot be compensated by a stabilizer of the form 2+ Fi. A reconstruction, using A; = 0.2, 
and A; = 0 is shown in Fig. 1(d). The result shows a distortion of the pyramid compared to little improvement 
IAPRS, Vol. 30, Part 5W1, ISPRS Intercommission Workshop "From Pixels to Sequences", Zurich, March 22-24 1995