A Mumford-Shah Diffusion Process for Shape-from-Shading Antonio Robles-Kelly Edwin R. Hancock Department of Computer Science, University of York,York YO1 5DD, UK. arobkell,erh @minster.cs.york.ac.uk Abstract
In this paper we show how the Mumford-Shah functional can be used to derive diffusion kernels that can be employed in the recovery of surface height in shape-from-shading. We commence from an initial field of surface normals which are constrained to fall on the Lambertian reflectance cone and to point in the direction of the local Canny edge gradient. We aim to find a path through this field of surface normals which can be used for surface height reconstruction. We demonstrate that the Mumford-Shah functional leads to a diffusion process which is a Markov chain on the field of surface normals. Moreover, the diffusion kernel for the Markov process depends on the rate of change of surface normal direction. As a result, we can find the steady state of the Markov chain using the leading eigenvector for the transition probability matrix computed from the diffusion kernels. We show how the steady state path can be used for height recovery and also for smoothing the initial field of surface normals.
1 Introduction Shape-from-shading is a process which allows surface shape to be recovered from shading variations using the physics of light. Unfortunately the problem is an ill-posed one [1] since at each image location it requires the two degrees of freedom of the surface normal, i.e. the local slant and tilt angles of a surface, to be recovered from a single measured brightness value. To overcome this problem the process is frequently posed in a variational setting, in which the goal is to find the field of surface normals so as to minimise a regularised energy functional [2]. The main problem with this approach is that the smoothness constraints tend to dominate the data-closeness constraints, resulting in a recovered field of surface normals that is devoid of fine surface detail. To overcome this problem, Worthington and Hancock [3] have recently developed a new geometric framework for shape-from-shading. Datacloseness is guaranteed by constraining the surface normals to fall on the local irradiance cone defined by Lamberts law. That is to say they must reside on a cone whose axis is the light source direction and whose opening angle is the inverse cosine of the measured image intensity. Although effective, this new framework leaves considerable scope for further development. In particular, it is concerned with the recovery of fields of surface normals and does not address the problem of recovering surface height. In this paper our aim is to develop a height recovery method in which the surface normals are constrained to fall on
Supported by CONACYT, under grant No. 146475/151752. 708
the Lambertian reflectance cone. We pose the height recovery process as that of finding an integration path through a field of surface normals which satisfy Lambert’s equation. Our starting point is the Mumford-Shah functional [4]. This introduces penalties associated with curve length and total squared curvature of the path. The required curvature may be computed from the change in surface normal direction between locations on the image plane. We use the Mumford-Shah energy to define a probability distribution for the integration path. We demonstrate that this is a Markov chain on the sites of the field of surface normals. The diffusion kernels for the Markov chain hence depend on the sectional curvature associated with the trace of the integration path on the surface. We make use of the Markov property and develop a graph-spectral method for recovering an integration path for surface height recovery. To do this we borrow ideas from spectral graph theory concerning random walks on graphs [5, 6, 7]. The steady state random walk on the graph can be located by finding the leading eigenvector of the associated transition probability matrix. Moreover, the eigenstructure of the transition matrix also allows us to locate surface patches. We use this property to find a curvature minimizing path across the field of available surface normals associated with each surface patch. We adjust the initial field of surface normals by traversing these paths and adjusting the surface normal direction. Each surface normal is rotated on the cone so that it points in the direction of minimum curvature along the path. This approach offers a number of advantages. First, we avoid smoothing across patch boundaries. Second, smoothing only proceeds in the direction of minimum curvature, and is hence regulated. Third, the propagation process is less computationally demanding since it is uni-dimensional.
2 Lambertian Reflectance In the case of Lambertian reflectance from a matte surface of constant albedo illuminated with a single point light-source, the observed intensity is independent of the viewing direction. The observed intensity depends only on the quantity of absorbed light, and this in turn is proportional to the cosine of the incidence angle. Suppose that is the is the unit-vector in the surface unit-vector in the direction of the light source and that normal direction at the pixel . According to Lambert’s law, the observed image intensity . at the pixel indexed is Lambert’s equation provides insufficient information to uniquely determine the surface normal direction. However, as recently observed by Worthington and Hancock [3], the equation does have a simple geometric interpretation which can be used to constrain the direction of the surface normal. The equation specifies that the surface normal must fall on the surface of a right-cone whose axis is aligned in the light-source direction and . whose apex angle is Worthington and Hancock [3] exploit this property to develop a two-step iterative process for SFS. The process commences from a configuration in which the surface normals are placed on the position on the irradiance cone where their projections onto the image plane are aligned in the direction of the local (Canny) image gradient.
3 Diffusion Kernels Stated formally, our goal is the recovery of height information from the field of surface normals. ¿From a computational standpoint the aim is to find a path along which simple trigonometry may be applied to increment the estimated height function. To be more
709
formal suppose that the surface under study is and that the field of surface normals is sampled on the plane . Our aim here is to find a curve across the plane that can be used as an integration path to reconstruct the height-function of the surface . The projection of the curve onto the surface is denoted by . Further, suppose that is the sectional curvature of the curve at the point with parametric co-ordinate . that minimises the Mumford-Shah [4] functional We seek the path
)!
'!
(
"!
* % !+ -,/.10 2"35476 #89$& ;:=/$
# %$& $
(1)
.0 3 6 where * and are constants. The probability of the path can be written as ? @ACB8DFE % !+ HG . The field of unit surface normals for the surface on the plane is denoted by . Accordingly, and following do Carmo [8], we let I8JK9 L represent the tangent plane to the surface at the point ( which belongs to the curve ! . To compute the sectional cur vature # %$& we require the differential of the surface or Hessian matrix >M JON1IPJK% L RQ I J % S . The maximum and minimum eigenvalues TVU and T : of >W J are the principal curvatures at the point ( . The corresponding eigenvectors X UZY[I J % L and X Y[I J 9 L : normal . the tangent plane I8JK9 L . At the point ( the unit form an orthogonal basis on \ vector to the curve is and the unit tangent vector is ]^J_Y`IPJK% L . The sectional curvature of at ( is given by 4 E X # %$& ]; Jba U : \ %T . U T : T : (2) a J . 4 E X \ where ]; Jba U : %T U T T is the normal curvature and c de af J is the angle : : between the curve normal and the surface normal. In practice, we will be dealing with points which are positioned at discrete positions on the pixel lattice. Suppose that and are the pixel indices of neighbouring points sampled on the pixel lattice along the path . With this discrete notation, the path energy and its probability are given by
g '!
* % ! h . 0 2 3o4p6 #=k : j j k and ? .10 q . 0 @ACBZrsE 2 3o4p6 #fk : j tj ku (3) < $ < $ i j kmlHn i j kmlHn j k is an estimate of the curvature based on the surface normal directions at the where # j k is the path distance between these points. Hence, we pixel locations and g , and $ can view the integration path as Markov chain across the field of surface normals. The diffusion kernel for the process, i.e. probability of migrating between sites and g is vw @ACB8DFE * j k G @ACB r E 2)3o4p6 # : j k j k u (4) sÒÓ%Ì m Ì : aFa¶mÌ Fa a¶aÔ Ê
i Fl
£ Ê Õ Ç Õ ±Ç 4 h T Õ Õ ± (9) ¡ : , Õ Ç is the leading eigenvector and the where the leading eigenvalue is unity i.e. T U Õ . To identify the patches, eigenvectors are normalised to be of unit length, i.e. x x Ê i l Ê we use the following iterative li U Õ iU U l Õ i U l ± . iÕ U l procedure. We initialise Êthei U lalgorithm by letting Ç Ç Further suppose that Ç is the Ê leading eigenvector of . The matrix Ì represents the first block of . The nodes with non-zero entries belong to the patch. These nodes may be identified and removed i U l consideration. To do this we Ê i l from Ê i U l l further E Ì in which the elements of the compute the residual transition matrix : Ê i: l Õl Çi : l of the residual transition matrix first patch are nulled. The leading eigenvector i l Õ Çi : Õ Çi : ± is used to compute the second block Ì : i È l is repeated iteratively Ê . The Õ process \ Ç , is the leading to identify all of the principal blocks of . Ati iteration l iÈ l Õ iÈ l± i È l eigenvector Ê Õ ¼ Ö È \ Ç Ç . , and the block is Ì of the residual transition probability matrix \ is the set of nodes for which the components of The index set of the patch Õ Çi È l indexed S ¼ Ö F l are non-zero. Hence, the index-set for the patch is the leading eigenvector i ¨&× x Õ Ç × 5 ¬ ® . It is important to stress that the patches are non-overlapping, i ¤ l Õ i l i.e. the Õ Ç , where inner product of the block eigenvectors for different patches is zero a Ç × ¬ ¸ . 6 Normal Correction and Height Recovery Our surface height recovery algorithm proceeds along the sequence of pixel sites defined by the order of the co-efficients of the leading eigenvector associated with the separate patch, the path is where the order of the patches. For the co-efficients of the leading eigenvector for this patch is such that . As we move from pixel-site to pixel-site defined by this path we perform two computations. First, we adjust the position of the surface normals on the irradiance cone so that they are smoothed in the local direction of minimum curvature. Second, we recover the surface from the Gauss map by incrementing the height function using the local slope parameters. At step of the algorithm, we make a transition from the pixel to the pixel with path-index . To perform smoothing and height with path-index and at the two pixel-sites. recovery, we make use of the surface normals Turning our attention first to the smoothing process, our aim is to adjust the surface normal directions so that they are consistent with the direction of the curvature minimis-
× ¼ Ö
°P¤Ç Fg ¤ ´ ¿ a¶a¶a
¤ Fg ¤ U 9 g ¤ : Hg ¤ ´ a¶a¶aFa¶ ° ¤Ç Fg ¤ U ¿ P° ¤Ç Fg ¤ : ¿
g È PU
\
kHØ g È kHØdÙ1Ú
713
300
300
300
300
250
250
250
250
200
200
200
200
150
150
150
150
100
100
100
100
50
50
50
0 40
0 40
40
30
0
40
20 10
10 0
0
40
40
35
35
35
30
30
30
25
25
25
25
20
20
20
20
15
15
15
15
10
10
10
10
5
5
5
0
0
0
30
0
5
10
15
20
25
30
35
0
40
5
10
15
20
25
30
35
0
5
10
15
20
25
30
35
0
40
300
300
300
300
250
250
250
250
200
200
200
200
150
150
150
150
100
100
100
50
50
50
0 40
0 40
0 40
40
30
30
20
10
10
10
300
300
300
250
250
250
200
200
200
150
150
150
150
100
100
100
50
50
50
0 40
0 40
0 40
30
20 10
40
40
40
35
35
30
30
30
25
25
25
25
20
20
20
20
15
15
15
15
10
10
10
10
5
5
5
0
0
0
0
5
10
15
20
25
30
35
0
40
5
10
15
20
25
30
35
0
5
10
15
20
25
30
35
0
40
300
300
300
250
250
250
250
200
200
200
200
150
150
150
150
100
100
100
50
50
50
0 40
40
30
30
20
40
10
10
10
300
300
300
250
250
250
200
200
200
150
150
150
100
100
100
100
50
50
50
40
30
40
30
20
20 10
10 0
30
20
40
30
30
20
20 10
10 0
0
0 40
40
30
20 10
50
40
0
50
0 40
30
10 0
0
40
10 0
150
20
35
20
10
200
30
30
30
20
0
250
0 40
25
30
10 0
0
300
0 40
20
40
20
10 0
0
15
50
30
20
20
10 0
10
0 50
40
30
30
20
20
5
100
0 40
30
10
5
0
40
300
0 40
40
30
20 10
0
35
30
0
20
10 0
0
40
35
0
30
20 10
10 0
0
0
50
30
20
20 10
10 0
40
30
10 0
0 40
40
30
30
20
40
100
40
30
35
20 10
200
20
30
20
0
250
40
25
30
10 0
0
300
30
20
20
10 0
0
15
50
30
20
20
10 0
10
0 40
40
30
30
20
20
5
100
40
30
0
5
0
40
10 0
0
40
35
40
30
20
20 10
10
0
30
30
20
20 10
10 0
40
30
30
20
20 10
0 40
40
30
30
20
50
0 40
0
20 10
10 0
0
Figure 2: Ground-truth, needlemaps, error plots and surface-height recovery results for four basic shapes under noise-controlled conditions.
ing integration path and also remain on their respective irradiance cone. We therefore parameterise the surface normal directions using two angles. The first of these is the polar between the surface normal and the light source direction. This is simply the angle opening angle of the irradiance cone, and this angle must be kept fixed in order to ensure that the recovered surface normal satisfies Lambert’s law. The second is the azimuthal , which measures the position of the surface normal on the irradiance cone. The angle and can be defined in terms of the surface normal components. Using angles
kHØ
Û Hk Ø Hk Ø Û
kHØ
714
some simple trigonometry we can write
} Ø Ø iFá l Zk Ø E Þ ß y} if âC k Ø de1SZ k Ø %ä1 m Û k Ø ÝÜ d à de8 } Ø Ø iF á l 1 y} otherwise ã kØ Zk Ø onto the image plane given by where ~ projection of the surface normal k~ HØ mÃFå disethe H k Ø M ä/ ^ .
These angles are adjusted in the smoothing step. This smoothing step uses the following equations to update the tilt and azimuthal angles of the surface normals
ç kHØ
Û)æ kHØ Û k Ø Eèç k Ø 4 Û k سÙ1Ú ç k Ø
ã
æk Ø k Ø
Û æ k Ø æ kHØ
where is a weight which measures the consistency of the surface normal directions are the updated anwith the direction of the curvature minimising path and gles. For the smoothing step, we make use of the thresholded change in surface normal direction if (10) otherwise
é
ksØ 2 x ksØ E kHسÙ1Ú x x ksØ E ksسÙ1Ú xêìë where the threshold ë is a constant that has been foundØ empirically to have aÔí kç
as its using the hyperbolic tangent optimum value. We compute the smoothing weight function as suggested by robust statistics and as used by Worthington and Hancock in their work on needle-map smoothing. The weight is given by
é ç kHØ ÂFî£mïVmx³Û kHØ E Û kHØdÙ1Ú x Hk Ø m :
(11)
and takes the product of the change in azimuthal angle and surface normal direction as its arguement. Once the updated angles have been computed, then the surface normals . With may be updated using the equation the smoothed surface normals to hand, we can compute the height increment. Using the trapezium rule the increment may be computed using the components of the surface normals and is given by
kHæ Ø %~ sk Ø ðÛ æ Hk Ø ~ sk Ø ^ öåKÛ æ sk Ø m æ kHØ ±
ñ
>È\
kHØdÙ1Ú kØ È > È 2 æ kHØ t ò 4 ksسÙ1Ú tò < æ âC âð
(12)
ñ
ä kHó
where is the Euclidean distance on the x-y plane between the two pixel sites associated , then the with the th transition. If the height-function is initialised by setting is . We merge abutting centre-height for the pixel with path-index patches by assuring they have the same mean boundary height.
gÈ
ä k Øô/Ú ä k Ø 4 È
7 Experiments The experimental evaluation of the new surface reconstruction method is divided into two parts. We commence with a sensitivity study aimed at evaluating the method on synthetic data. In the second part of the study, we focus on real-world data.
715
7.1 Synthetic Data In this section we provide some experiments on synthetic data. The aim here is to determine the accuracy of the surface reconstruction method. To this end we have generated synthetic surfaces. From the surfaces, we have computed the field of surface normal directions. We have then applied the graph-spectral method to the 2D field of surface normals to recover an estimate of the surface height. In Figure 2 we show the results obtained for a series of different surfaces under conditions of controlled noise. To do this we have added random measurement errors to the surface height. The measurement errors have been sampled from a Gaussian distribution with zero mean and known variance. In the top row we show the original noise-free synthetic surface, i.e. ground-truth. In the second row of the figure we show the field of surface normals for the Figure 3: Error versus variance noise-free surface. The third row shows the surface reconstructed from the field of surface normals when no Gaussian noise has been added. The fourth row shows the absolute error between the ground-truth and reconstructed surface height. From left-to-right the surfaces studied are a dome, a sharp ridge, a torus and a volcano. In all four cases the surface reconstructions are qualitatively good. For the dome the height errors are greater at the edges of the surface where the slope is largest. In the case of the ridge, there are errors at the crest. For the volcano, there are some problems with the recovery of the correct depth of the “caldera”, i.e. the depression in the centre. For the reconstructed surfaces, the relative mean-squared errors are 5.6% for the dome, 10.8% for the ridge, 7.8% for the torus and 4.7% for the volcano. Hence, the method seems to have greater difficulty for surfaces containing sharp creases. In the fifth row of Figure 2 we show the field of surface normals for the noise-corrupted surface. In the sixth row, we show the result of reconstructing the surfaces shown in the top row, when random height errors have been added. The seventh , i.e. bottom, row shows the difference between the height of the surface reconstructed from the noisy surface normals and the ground-truth height function. In the case of all four surfaces, the gross structure is maintained. However, the recovered height is clearly noisy. The height difference plots are relatively unstructured. These are important Figure 4: Results of the surface integration process. ERROR AS A FUNCTION OF THE VARIANCE
35
DOME RIDGE TORUS VOLCANO
30
25
20
15
10
5
0
0
0.2
0.4
0.6
0.8
1
70
80
90
60
80
70
70
50
60
60
50
40
50
40
30
40
30
30
20
20
20
10
10
10
0
0
0
10
20
30
716
40
50
60
70
80
0
0
10
20
30
40
50
60
70
0
10
20
30
40
50
60
70
80
90
observations. They mean that our graph-spectral method simply transfers errors in surface normal direction into errors in height, without producing structural noise artifacts. To investigate the effect of noise further, Figure 4 we plot the mean-squared error for the reconstructed surface height as a function of the standard deviation of the added Gaussian noise. From the plots for the different surfaces shown in Figure 3, it is clear that the mean-squared error grows slowly with increasing noise standard deviation.
7.2 Real World Experiments We have performed experiments using real world imagery. The images used in our study have been selected because they have proved problematic when conventional shape-fromshading algorithms are used. In the top panel of Figure 4 we show the three gray-scale images studied. These are a wooden block, an urn handle and a marble toe. The wooden block image illustrates the usefulness of the new method when dealing with ambiguous initial normal estimates (i.e. the planar faces of the block) while the other two images show the ability of the new method for smoothing on the minimum curvature direction keeping patch boundaries fixed. In the second row, we show the final needle-maps. The smoothing process improves the quality of the field of surface normals, but does not oversmooth across the patch boundaries. In the third row, we show the surfaces reconstructed from the smoothed fields of surface normals. It is evident that the smoothing process leads to surfaces which are in good subjective agreement with the image contents.
8 Conclusions In this paper, we have described a graph-spectral algorithm shape-from-shading. The method commences from an initial surface normal estimate and performs smoothing along a curvature minimising path defined by the leading eigenvector of a curvature-based transition probability matrix. By traversing this path and applying some simple trigonometry we reconstruct the surface height function. Results on real world images reveal that the method is able to deliver accurate height reconstructions for complex surfaces.
References [1] B. K. P. Horn and M. J. Brooks. Height and gradient from shading. International Journal of Computer Vision, 5(1):37–75, 1986. [2] B. K. P. Horn and M. J. Brooks. The variational approach to shape from shading. CVGIP, 33(2):174–208, 1986. [3] P. L. Worthington and E. R. Hancock. Needle map recovery using robust regularizers. Image and Vision Computing, 17:545–557, 1999. [4] D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated variational problems. Comm. in Pure and Appl. Math., 42(5):577–685, 1989. [5] L. Lov´asz. Random walks on graphs: a survey. Bolyai Society Mathematical Studies, 2(2):1– 46, 1993. [6] Fan R. K. Chung. Spectral Graph Theory. American Mathematical Society, 1997. [7] M. Doob D. Cvetkovi´c and H. Sachs. Spectra of Graphs:Theory and Application. Academic Press, 1980. [8] M. P. Do Carmo. Differential Geometry of Curves and Surfaces. Prentice Hall, 1976. [9] R. S. Varga. Matrix Iterative Analysis. Springer, second edition, 2000.
717