Three dimensional face recognition

FIELD AND BACKGROUND OF THE INVENTION

0001The present invention relates to a method and apparatus for three-dimensional face recognition and, more particularly, but not exclusively to such a method and apparatus that both obtains a 3D representation of a face and uses that representation for matching purposes.

0002Face recognition has recently become an important task of computer vision, and is needed in a wide range of biometric and security applications. However, most existing face recognition systems are reported to be sensitive to image acquisition conditions such as head position, illumination, etc. and can therefore be inaccurate and easily fooled. Reference is made to American Civil Liberties Union (ACLU) of Florida, <i>Press Release</i>, 14 May 2002, Available: http://www.aclufl.org/pbfaceitresults051402.html.

0003In general, modem face recognition approaches can be divided into two wide categories: 2D approaches, using only image information (which can be either grayscale or color), and 3D approaches, which incorporate three-dimensional information as well.

0004While simpler in data acquisition (which permits real-time surveillance applications, such as face recognition from a video-taped crowd in pubic places), the 2D approach suffers from sensitivity to illumination conditions and head rotation. Since the image represents the light reflected from the facial surface at a single observation angle, different illumination conditions can result in different images, which are likely to be recognized as different subjects (see <figref idref="f0003">Figure 3</figref>). One of the classical 2D face recognition algorithms is the Turk and Pentland eigenfaces algorithm. For a full discussion see <nplcit id="ncit0001" npl-type="s"><text>M. Turk and A. Pentland, Face recognition using eigenfaces, CVPR, May 1991, pp. 586-591</text></nplcit>, and M. Turk and A. Pentland, "Face recognition system, <patcit id="pcit0001" dnum="US5164992A"><text>US Patent 5,164,992, 17 Nov 1990</text></patcit>. The Eigenfaces algorithm works as follows: Given a set of faces arising from some statistical distribution, the principal components of this distribution form a set of features that characterize the variation between faces. "Eigenfaces" are the eigenvectors of the set of "all" faces. The eigenface approach treats face recognition as a 2D classification problem, without taking into consideration the significant difference in the images resulting from illumination conditions and head rotation. For this reason, eigenfaces usually produce mediocre results when faced with real life rather than laboratory conditions.

0005The 3D approach provides face geometry information, and face geometry information is independent of viewpoint and lighting conditions. Thus such information is complementary to the 2D image. 3D information, however, not only carries the actual facial geometry, but includes depth information which allows easy segmentation of the face from the background.

0006Gordon showed that combining frontal and profile views can improve the recognition accuracy and reference is made to <nplcit id="ncit0002" npl-type="s"><text>G. Gordon, "Face recognition from frontal and profile views", Proc. of the International Workshop on Face and Gesture Recognition, Zurich, Switzerland, pp. 47-52, Jun 1995</text></nplcit>.

0007Beumier and Acheroy show the adequacy of using geometric information in the rigid face profile for subject identification in a system using structured light for 3D acquisition, and reference is made to <nplcit id="ncit0003" npl-type="s"><text>C. Beumier and M. P. Acheroy, Automatic Face Identification, Applications of Digital Image Processing XVIII, SPIE, vol. 2564, pp. 311-323, July 1995</text></nplcit>.

0008The above described approach may be generalized to the whole surface, and reference is made to <nplcit id="ncit0004" npl-type="b"><text>C. Beumier and M. P. Acheroy, Automatic Face Authentication from 3D Surface, British Machine Vision Conference BMVC 98, University of Southampton UK, 14-17 Sep, 1998, pp 449-458.", 1998</text></nplcit>, , who describe such a generalization using global surface matching. However, surface matching is sensitive to facial expressions and cannot be considered a comprehensive solution.

0009There is thus a widely recognized need for, and it would be highly advantageous to have, a facial recognition system that uses three-dimensional information but is devoid of the above limitations such as being sensitive to facial expressions, lighting of the subject, or to angle.

SUMMARY OF THE INVENTION

0010The present invention obtains a canonical form representation of 3-D facial data and uses that canonical representation to compare with other faces.

0011According to one aspect of the present invention there is provided Apparatus for processing 3-dimensional data of a geometric body for matching, wherein said geometric body is a face, having soft geometric regions, being regions of said face susceptible to short term geometric changes, and hard geometric regions, being regions substantially insusceptible to said short term geometric changes, said apparatus comprising: <ul id="ul0001" list-style="none" compact="compact"><li>a geodesic converter, for receiving an input comprising 3-dimensional topographical data of said geometric body, and for converting said data into a series of geodesic distances between pairs of points of said data;</li><li>a triangulator connected prior to said geodesic converter to receive said 3-dimensional topographical data and to form therefrom a manifold by triangulation, said manifold being used as said input of said geodesic converter;</li><li>a multi-dimensional scaler, connected subsequently to said geodesic converter, for forming a low dimensional Euclidean representation of said series of geodesic distances, said low dimensional Euclidean representation providing a bending invariant representation of said geometric body suitable for matching with other geometric shapes;</li><li>a subsampler located after said triangulator and prior to said geodesic converter, configured to subsample said triangulated manifold, and to provide to said geodesic converter a subsampled triangulated manifold: <ul id="ul0002" list-style="none" compact="compact"><li>a preprocessor, located prior to said subsampler, for removing said soft geometric regions from said face, wherein said preprocessor is operable to carry out removal of said soft regions by application of a geodesic mask.</li></ul></li></ul>

0012The apparatus preferably comprises a subsampler located prior to said geodesic converter, configured to subsample said triangulated manifold, and to provide to said geodesic converter a subsampled triangulated manifold.

0013Preferably, said manifold comprises a plurality of vertices and said subsampler is operable to select a first vertex and to iteratively select a next vertex having a largest geodesic distance from vertices already selected, until a predetermined number of vertices has been selected.

0014Preferably, said subsampler is operable to use the fast marching method for triangulated domains to calculate geodesic distances between vertices for said iterative selection.

0015Preferably, said geometric body is a face, having soft geometric regions, being regions of said face susceptible to short term geometric changes, and hard geometric regions, being regions substantially insusceptible to said short term geometric changes, said apparatus comprising a preprocessor, located prior to said subsampler, for removing said soft geometric regions from said face.

0016Preferably, said preprocessor is operable to identify said soft regions by identification of an orientation point on said face.

0017Preferably, said orientation point is at least one of a nose tip, a pair of eyeball centers and a mouth center.

0018Preferably, said preprocessor is further operable to center said face.

0019Preferably, said preprocessor is further operable to crop said face.

0020Preferably, said preprocessor is operable to carry out removal of said soft regions by application of a geodesic mask.

0021Preferably, said geodesic converter is operable to use the fast marching method for triangulated domains to calculate said geodesic distances.

0022The apparatus preferably comprises a triangulator for forming said triangulated manifold from scan data of a geometric body.

0023The apparatus may further be operable to embed said triangulated manifold into a space of higher than two dimensions, thereby to include additional information with said topographical information.

0024Preferably, said additional information is any one of a group comprising texture information, albedo information, grayscale information, and color information. In the case of texture information, the texture can be extracted from the initial scan, and subsequently mapped onto the low dimensional Euclidian representation. Subsequently separate coefficients for the topographical and texture information can be used to contribute to a weighted distance to another surface in the matching process.

0025Preferably, said subsampler comprises an optimizer for allowing a user to select an optimum subsampling level by defining a trade-off between calculation complexity and representation accuracy.

0026Preferably, said multi-dimensional scalar is configured such that said Euclidean representation comprises a predetermined number of eigenvalues extractable to be used as co-ordinates in a feature space.

0027Preferably, said predetermined number of eigenvalues is at least three, and said feature space has a number of dimensions corresponding to said predetermined number.

0028According to an embodiment of the present invention there is provided apparatus for matching between geometric bodies based on 3-dimensional data comprising: <ul id="ul0003" list-style="none" compact="compact"><li>an input for receiving representations of geometric bodies as Euclidean representations of sets of geodesic distances between sampled points of a triangulated manifold, said Euclidean representations being substantially bending invariant representations,</li><li>a distance calculator for calculating distances between respective geometric bodies based on said Euclidean representation and</li><li>a thresholder for thresholding a calculated distance to determine the presence or absence of a match.</li></ul>

0029Preferably, said distance calculator comprises: <ul id="ul0004" list-style="none" compact="compact"><li>an eigenvalue extractor for extracting a predetermined number of eigenvalues from said Euclidean representations, and</li><li>a plotter for plotting said predetermined number of eigenvalues as a point on a feature space having a dimension for each of said predetermined number of Eigenvalues,</li><li>and wherein said thresholder is configured to be sensitive to clustering within said feature space, thereby to determine said presence or absence of said match.</li></ul>

0030Preferably, said predetermined number is three.

0031Preferably, said Euclidean representation is based upon geodesic distances between a subsampling of points of said triangulated manifold.

0032Preferably, said geometric body is a face, having soft geometric regions, being regions susceptible to short term geometric change and hard geometric regions, being regions substantially insusceptible to short term geometric changes, and wherein said Euclidian representation is substantially limited to said hard geometric regions.

0033Preferably, said distance calculator is configured to use the Hausdorff metric.

0034According to a third aspect of the present invention there is provided apparatus for obtaining 3-dimensional data of a geometric body for matching, and using said data to carry out matching between different bodies wherein said geometric body is a face, having soft geometric regions and hard geometric regions, said apparatus comprising: <ul id="ul0005" list-style="none" compact="compact"><li>a three dimensional scanner for obtaining three-dimensional topographical data of said body;</li><li>a triangulator for receiving said three-dimensional topographical data of said geometric body and forming said data into a triangulated manifold;</li><li>a geodesic converter, connected subsequently to said triangulator, for converting said triangulated manifold into a series of geodesic distances between pairs of points of said manifold;</li><li>a multi-dimensional scaler, connected subsequently to said geodesic converter, for forming a low dimensional Euclidean representation of said series of geodesic distances, said low dimensional Euclidean representation providing a bending invariant representation of said geometric body;</li><li>a distance calculator, connected subsequently to said multi-dimensional scaler, for calculating distances between geometric bodies based on said Euclidean representation;</li><li>a thresholder, connected subsequently to said distance calculator, for thresholding a calculated distance to determine the presence or absence of a match;</li><li>a subsampler located between said triangulator and said geodesic converter, configured to subsample said triangulated manifold, and to provide to said geodesic converter a subsampled triangulated manifold; and</li><li>a preprocessor, located between said triangulator and said subsampler, for removing said soft geometric regions from said face; wherein said preprocessor is operable to carry out removal of said soft regions by application of a geodesic mask.</li></ul>

0035Preferably, said distance calculator comprises: <ul id="ul0006" list-style="none" compact="compact"><li>an eigenvalue extractor for extracting a predetermined number of eigenvalues from said Euclidean representations, and</li><li>a plotter for plotting said predetermined number of Eigenvalues as a point on a feature space having a dimension for each of said predetermined number of eigenvalues,</li><li>and wherein said thresholder is configured to be sensitive to clustering within said feature space, thereby to determine said presence or absence of said match.</li></ul>

0036Preferably, said predetermined number is three.

0037Alternatively, said predetermined number is greater than three.

0038The apparatus preferably comprises a subsampler located between said triangulator and said geodesic converter, configured to subsample said triangulated manifold, and to provide to said geodesic converter a subsampled triangulated manifold.

0039Preferably, said subsampler is operable to use geodesic distances in selecting points from said triangulated manifold to include in said subsampled triangulated manifold.

0040Preferably, said subsampler is configured to take an initial point and then iteratively to select points by taking points furthest away in terms of a geodesic distance from already selected points.

0041Preferably, said geometric body is a face, having soft geometric regions and hard geometric regions, said apparatus comprising a preprocessor, located between said triangulator and said subsampler, for removing said soft geometric regions from said face.

0042Preferably, said geometric body is a face, having soft geometric regions and hard geometric regions, said apparatus comprising a preprocessor, located between said triangulator and said geodesic converter, for removing said soft geometric regions from said face.

0043Preferably, said preprocessor is operable to identify said soft regions by identification of an orientation point on said face.

0044Preferably, said orientation point is a nose tip.

0045Preferably, said preprocessor is further operable to center said face.

0046Preferably, said preprocessor is further operable to crop said face.

0047Preferably, said preprocessor is operable to carry out removal of said soft regions by application of a geodesic mask.

0048Preferably, said geodesic converter is operable to use the fast marching method for triangulated domains to calculate said geodesic distances.

0049Preferably, said subsampler comprises an optimizer for allowing a user to select an optimum subsampling level by defining a trade-off between calculation complexity and representation accuracy.

0050Preferably, said distance calculator is configured to use the Hausdorff metric.

0051According to a fourth aspect of the present invention there is provided a method of image preprocessing of three-dimensional topographical data for subsequent classification, wherein said geometric body is a face, having soft geometric regions and hard geometric regions, the method comprising: <ul id="ul0007" list-style="none" compact="compact"><li>providing said three-dimensional topographical data as a three-dimensional triangulated manifold,</li><li>subsampling said triangulated manifold,</li><li>generating a matrix of geodesic distances to selected vertices of said manifold,</li><li>removing said soft geometric regions from said face by application of a geodesic mask; and</li><li>using multi-dimensional scaling to reduce said matrix to a canonical representation in a low-dimensional Euclidean space, thereby to provide a representation suitable for subsequent classification.</li></ul>

0052The method may further comprise selecting said vertices for generating said matrix of geodesic distances by a procedure comprising iteratively selecting a next vertex having a largest geodesic distance from vertices already selected, until a predetermined number of vertices has been selected.

0053The method may further comprise determining geodesic distances for said procedure using the fast marching method for triangulated domains.

0054Preferably, said three-dimensional topographical data is a representation of a face, and further comprising cropping said face to exclude parts of said face being susceptible to short term changes, thereby to render said subsequent classification substantially invariant to said short term changes.

0055According to a fifth aspect of the present invention there is provided a method of classifying images of three-dimensional bodies, , wherein said three dimensional bodies comprise faces having soft geometric regions and hard geometric regions, the method comprising: <ul id="ul0008" list-style="none" compact="compact"><li>obtaining representations of said three dimensional bodies as canonical form representations derived from geodesic distances between selected sample points taken from surfaces of said three dimensional bodies;</li><li>removing said soft geometric regions from said face by application of a geodesic mask;</li><li>from each representation deriving co-ordinates on a feature space; and</li><li>classifying said bodies according to clustering on said feature space.</li></ul>

0056Preferably, n is at least three.

0057Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. The materials, methods, and examples provided herein are illustrative only and not intended to be limiting.

0058Implementation of the method and system of the present invention involves performing or completing selected tasks or steps manually, automatically, or as a combination thereof. Moreover, according to actual instrumentation and equipment of preferred embodiments of the method and system of the present invention, one or more steps could be implemented by hardware or by software on any operating system of any firmware or a combination thereof. For example, as hardware, selected steps of the invention could be implemented as a chip or a circuit As software, selected steps of the invention could be implemented as a plurality of software instructions being executed by a computer using any suitable operating system. In any case, selected steps of the method and system of the invention could be described as being performed by a data processor, such as a computing platform for executing a plurality of instructions.

BRIEF DESCRIPTION OF THE DRAWINGS

0059The invention is herein described, by way of example only, with reference to the accompanying drawings. With specific reference now to the drawings in detail, it is stressed that the particulars shown are by way of example and for purposes of illustrative discussion of the preferred embodiments of the present invention only, and are presented in the cause of providing what is believed to be the most useful and readily understood description of the principles and conceptual aspects of the invention. In this regard, no attempt is made to show structural details of the invention in more detail than is necessary for a fundamental understanding of the invention, the description taken with the drawings making apparent to those skilled in the art how the several forms of the invention may be embodied in practice.

0060In the drawings: <ul id="ul0009" list-style="none" compact="compact"><li><figref idref="f0001">FIG. 1</figref> is a simplified block diagram showing a preferred embodiment of a device for gathering 3D topographical data of a body and processing the data into a canonical form representation for efficient matching, according to a first preferred embodiment of the present invention;</li><li><figref idref="f0002">FIG. 2</figref> is a simplified diagram of apparatus for receiving 3D topographical data in canonical form representation and for carrying out matching or classification, according to a second preferred embodiment of the present invention;</li><li><figref idref="f0003">FIG. 3</figref> is a series of 2D representations of a single face, each taken under different lighting conditions. The faces are clearly the same to the human observer but are very difficult to match using conventional image analysis techniques;</li><li><figref idref="f0003">FIG. 4a</figref> is a simplified schematic diagram showing a first stage of 3D data gathering using depth code illumination;</li><li><figref idref="f0003">FIG. 4b</figref> is a simplified schematic diagram showing a second stage of 3D data gathering using depth code illumination;</li><li><figref idref="f0004">FIG. 5</figref> is a photograph showing a 3D scanner arrangement comprising a video camera and a moving laser projector;</li><li><figref idref="f0004">FIG. 6</figref> is a simplified schematic diagram illustrating a photometric stereo acquisition scheme;</li><li><figref idref="f0005">FIG. 7</figref> is a triangulated manifold representation using data points gathered by 3D scanning of a face;</li><li><figref idref="f0005">FIG. 8</figref> is a simplified diagram showing the manifold of <figref idref="f0005">Fig. 7</figref> after subsampling;</li><li><figref idref="f0006">FIG. 9</figref> is a simplified block diagram showing in greater detail the subsampler of <figref idref="f0001">Fig. 1</figref>;</li><li><figref idref="f0007">FIG. 10</figref> is a simplified block diagram showing in greater detail the geodesic converter of <figref idref="f0001">Fig. 1</figref>;</li><li><figref idref="f0008">FIG. 11</figref> is a simplified flow chart showing operation of a preferred embodiment of the present invention;</li><li><figref idref="f0009">FIG. 12</figref> is a simplified flow chart showing in greater detail the subsampling stage of <figref idref="f0008">FIG. 11</figref>;</li><li><figref idref="f0010">FIG. 13</figref> is a database of six faces, used as the subject of the first experiment;</li><li><figref idref="f0010">FIG. 14</figref> shows one of the faces of <figref idref="f0001">FIG. 1</figref>, a variation of the face to simulate change in expression, and a further variation to illustrate a change in a major feature;</li><li><figref idref="f0011">FIG. 15</figref> is a graph showing results of the first experiment plotted onto a 3-dimensional feature space. Faces differing only by a change in expression form a cluster on the feature space whereas faces differing by a change in a major feature are distant;</li><li><figref idref="f0012">FIG. 16</figref> illustrates nine different illuminations of a given face and their reconstruction, using photometric stereo, into a manifold by solution of the Poisson equation;</li><li><figref idref="f0013">FIG. 17</figref> is a set of ten faces each reconstructed using least squares from nine facial illuminations of a different subject for use in the second experiment;</li><li><figref idref="f0014">FIG. 18</figref> is a set of images of a single subject each with a different pose, for use in experiment II;</li><li><figref idref="f0014">FIG. 19</figref> is a graph showing results of the second experiment, plotted on a 3-dimensional feature space. Different poses of the same subject form clusters;</li><li><figref idref="f0015">FIG. 20</figref> is a schematic diagram illustrating an alignment process of two facial surfaces;</li><li><figref idref="f0015">FIG. 21</figref> is a graph showing the results of the third experiment without alignment, plotted on a three-dimensional feature space. No recognizable clusters appear;</li><li><figref idref="f0015">FIG. 22</figref> is a graph showing the results of the third experiment carried out with optimal alignment, plotted on a three-dimensional feature space. Clustering is present but indistinct;</li><li><figref idref="f0016">FIG. 23</figref> is a set of three faces from a further database, the faces used as subjects for the fourth experiment;</li><li><figref idref="f0016">FIG. 24</figref> is a set of three poses of a face of <figref idref="f0016">FIG. 23</figref>, also for use in the fourth experiment;</li><li><figref idref="f0016">FIG. 25</figref> shows a texture map of a face before application of a mask as part of preprocessing;</li><li><figref idref="f0016">FIG. 26</figref> shows the texture map of <figref idref="f0016">FIG. 25</figref> after application of a geodesic mask;</li><li><figref idref="f0016">FIG. 27</figref> shows a depth map of a face before application of a mask as a part of preprocessing;</li><li><figref idref="f0016">FIG. 28</figref> shows the depth map of <figref idref="f0016">FIG. 27</figref> after application of a geodesic mask;</li><li><figref idref="f0017">FIG. 29</figref> shows a triangulated manifold as obtained directly from 3D topographical image data;</li><li><figref idref="f0017">FIG. 30</figref> shows the manifold of <figref idref="f0017">FIG. 29</figref> following subsampling;</li><li><figref idref="f0017">FIG. 31</figref> shows the manifold of <figref idref="f0017">FIG. 29</figref> reduced to canonical form;</li><li><figref idref="f0018">FIG. 32</figref> is a simplified diagram illustrating viewing and illumination angles for 3D data acquisition;</li><li><figref idref="f0019">FIG. 33</figref> is a simplified diagram illustrating the relationship between surface gradient and local stripe direction;</li><li><figref idref="f0020">FIG. 34</figref> is a simplified diagram illustrating mapping of a Cartesian co-ordinate system onto a 3-dimensional manifold;</li><li><figref idref="f0021">FIG. 35</figref> is a simplified diagram illustrating the process of selection of additional vertices from which to measure geodesic distances;</li><li><figref idref="f0021">FIG. 36</figref> shows a series of the same face illuminated from different angles;</li><li><figref idref="f0022">FIG. 37A</figref> is a simplified block diagram illustrating apparatus for processing 3D representations of faces and performing matching of the processed results without forming the 3D surface as an intermediate step, according to a preferred embodiment of the present invention;</li><li><figref idref="f0023">FIG. 37B</figref> is a simplified flow chart illustrating a sequence of operations involved in use of the apparatus of <figref idref="f0022">Fig. 37A</figref>;</li><li><figref idref="f0024">FIG. 38</figref> is a simplified diagram illustrating three stages in the processing of gathered face information, initially gathered data, extraction of a 3D representation and preprocessing to concentrate on the non-changing regions of the face;</li><li><figref idref="f0025">FIG. 39</figref> shows results obtained using the embodiment of <figref idref="f0022">FIG. 37A</figref>;</li><li><figref idref="f0026">FIGs. 40</figref> and <figref idref="f0027">41</figref> are a single flow chart illustrating a procedure for incorporating texture information into the matching process of the embodiment of <figref idref="f0001">FIG. 1</figref>; and</li><li><figref idref="f0028">FIG. 42</figref> is a simplified diagram illustrating isolation and subsequent mapping onto the canonical form of the texture data according to the embodiment of <figref idref="f0026">FIGs. 40</figref> and <figref idref="f0027">41</figref>.</li></ul>

DESCRIPTION OF THE PREFERRED EMBODIMENTS

0061The present embodiments show a face recognition approach based on 3D geometric or topographical information. Given a 3D facial surface, a geometric bending-invariant canonical form can be obtained by sampling the surface, computing the geodesic distances between points on it (using the Fast Marching method) and applying Multidimensional scaling (MDS). MDS allows representation of the surface in a low-dimensional Euclidean space, the bending invariant canonical form, and the canonical form is a representation which can be used for classification in general and matching in particular.

0062In one embodiment a triangulated manifold of the facial surface is formed from which the geodesic distances are computed. In another embodiment the geodesic distances are obtained from surface gradient data without forming the manifold.

0063Facial recognition using 3D geometry may be used on its own or to augment conventional 2D imaging. In one embodiment, 2D or texture data is processed together with 3D range or geometry data and matching is then carried out based on a weighted Euclidean distance using both the texture and the range data. As the 3D geometry is independent of viewpoint and lighting conditions, accurate face classification can be achieved.

0064Referring now to the drawings, <figref idref="f0001">Figure 1</figref> is a simplified diagram showing apparatus for obtaining 3-Dimensional data of a geometric body for classification, including matching, according to a first preferred embodiment of the present invention. The preferred embodiments relate specifically to matching of faces but the skilled person will be aware that the principles of the present invention are applicable to any geometric body having a three-dimensional structure.

0065Apparatus 10 comprises a three dimensional scanner 12 for obtaining three-dimensional topographical data of the body. Several types of scanner are described hereinbelow together with brief discussions of how to process data therefrom in some of the cases.

0066Data from the three-dimensional scanner 12 is passed to a triangulator 14. The triangulator may perform triangulation on the data received from the scanner in order to generate a three-dimensional triangulated manifold to represent the topological features of the body. The exact operation of the triangulator to derive the manifold may vary depending on the way in which the 3D information is gathered. In certain cases the manifold may be formed directly from the gathered data without the need for any intervening processing stage. The manifold preferably represents all of the three-dimensional topology of the body and therefore is in theory sufficient for allowing matching. However, in practice, direct comparisons using the triangulated manifold have a number of disadvantages, as will be demonstrated in experiment 3 hereinbelow. They require a large amount of calculation. Matching does not distinguish reliably between different faces. Moreover matching generally fails when the same face has a different expression and matching is unreliable even when the same face is posed at a different angle.

0067Embodiments of the present invention therefore preferably include four additional processing stages, the first of which is a preprocessor 16. Preprocessor 16 takes a reference point in order to determine an orientation around the manifold. The reference point can be any point on the face provided it is relatively easy to find automatically from a manifold of a face. A suitable reference point as used in the preferred embodiments is the tip of the nose. Other possible reference points include centers of eyeballs and the center of the mouth. Once the preprocessor has found the tip of the nose it is able to orientate itself with respect to the rest of the face and then parts of the face whose geometry is particularly susceptible to expressions, hereinafter referred to as soft regions, can be ignored. Parts of the face that are invariant with change of expression and the like, hereinafter hard regions, can be retained or even emphasized. As will be explained in greater detail below, the definition of soft regions is not fixed. For some methods and in some circumstances soft regions to be excluded may include all of the lower region of the face around the mouth. In other cases less drastic exclusions may be considered. In one embodiment, described in greater detail below, soft regions are removed using a geodesic mask. The mask may be applied separately to a texture map of the face and a depth map of the face.

0068Following the preprocessor is a subsampler 18. The subsampler 18 takes the preprocessed manifold and removes points so as to produce a less well defined manifold, but one which still defines the essential geometry of the face it is desired to match. In preferred embodiments, the user is able to select the number of points to trade off between accurate matching - a large number of points - and faster processing - a smaller number of points. As will be discussed in greater detail below, a preferred embodiment of the sub-sampler uses the Voronoi subsampling technique which begins at an initial point or vertex on the manifold and then adds the point or vertex having the greatest distance therefrom. The procedure is repeated iteratively until the selected number of points are included. Preferably the technique uses geodesic distances, which may be obtained using the fast marching method for the triangulated domain (FMM-TD), as described below.

0069Following the subsampler is a geodesic converter 20. The geodesic converter 20 receives the list of points of the subsampled manifold and calculates a vector for each pair of points. The vectors are expressed as geodesic distances, and the fast marching algorithm for the triangulated domain is again used to obtain the geodesic distances in an efficient manner as possible.

0070Following the geodesic converter is a multi-dimensional scaler 22, which takes the matrix of the geodesic distances calculated by the geodesic converter 20, referred to below as the distance matrix, and forms a low dimensional Euclidean representation of the series of geodesic distances, using multi-dimensional scaling. Multi-dimensional scaling is discussed in detail below. The low dimensional Euclidean representation provides a bending invariant representation of the geometric body, as will be explained in the discussion on multi-dimensional scaling below. The use of such a bending invariant representation ensures that the matching process is not fooled by, for example, scans of the head at different angles.

0071The output 24 of the multi-dimensional scalar is a representation of the 3D face in terms of Euclidean distances between surface points, referred to hereinbelow as the canonical form representation.

0072Reference is now made to <figref idref="f0002">Fig. 2</figref>, which is a simplified diagram showing a matching apparatus for matching two faces using the canonical form output as described above. The matcher 30 may be a continuation of the apparatus 10 or may be supplied as a separate unit. The matcher 30 comprises a distance calculator 32, which takes as input two canonical form representations 34 and 36, and calculates a distance therebetween. The distance calculation may use any suitable method for comparison of the canonical forms for each of the faces to be matched. A straightforward approach is to measure a distance between two sets of points, using, for example, the Hausdorff metric. However, the Hausdorf metric based method is computationally extensive.

0073An alternative approach, used in the present embodiments, takes the first <i>m</i> eigenvalues obtained from the MDS procedure to provide coordinates in a low-dimensional feature space. Although the dominant eigenvalues do not describe the canonical form entirely, it is reasonable that similar faces have similar eigenvalues (and thus form clusters in the feature space). A distance is calculated between the geometric bodies, or, as will be described below, plotted on a graph of the feature space and a thresholder 38, which is connected subsequently to the distance calculator, thresholds the calculated distance to determine the presence or absence of a cluster in the feature space, the cluster indicating a match. In the embodiments described in detail herein, the first three Eigenvalues are taken and are plotted in a three dimensional feature space.

0074Reference is now made to <figref idref="f0003">Fig. 3</figref>, which shows a series of three two-dimensional images. It will be apparent to any human observer that the three images are of the same face, however conventional automatic matching techniques generally find very large distances between the three images and consequently fail to match them. Thus two-dimensional facial matching is prone to errors introduced by simple variables such as lighting direction. In fact, areas of the face that have high reflectivity, such as the eyes, can change substantially for very minor changes in lighting.

0075Returning now to <figref idref="f0001">Fig. 1</figref>, and as mentioned above, apparatus 10 preferably comprises a three-dimensional scanner 12. The face recognition described in the present embodiments treats faces as three-dimensional surfaces. It is therefore first of all necessary to obtain the facial surface of the subject that it is desired to recognize. Below is a short overview of currently available range finding techniques that are able to scan a face and generate three-dimensional data.

Leaser range camera (Zcam)

0076Reference is now made to <figref idref="f0003">Figs. 4a and 4b</figref>, which are simplified diagrams showing successive stages of operation of a laser range camera using depth encoding illumination. Currently, the fastest and most accurate, but at the same time most expensive, range cameras are those that are based on depth-encoding illumination.

0077Depth-encoded illumination is based upon the generation of a wall of light 40 which moves along the field of view. As the light hits the objects 42 in the scene, it is reflected back towards the camera 44. Due to the finite speed of the light, the reflection carries an imprint of the object depth which may be decoded using the appropriate electronics. The illumination is typically generated using IR laser diodes.

0078A 3D acquisition equipment of this type is known from <patcit id="pcit0002" dnum="WO0118563A"><text>WO Patent 01/18563, March 15, 2001</text></patcit>, the contents of which are hereby incorporated by reference. Such a technology allows real-time data acquisition at about 30fps frame rate. Depth resolution can be greatly improved by averaging frames in time. Typical prices range between US$2K - 50K.

3D scanner

0079A slower and cheaper version of a 3D scanner is based upon <nplcit id="ncit0005" npl-type="s"><text>J.-Y. Bouguet and P. Perona, "3D photography on your desk", in Proc. of the Int. Conf. on Computer Vision, Bombay, India, Jan. 1998</text></nplcit>. The scanner is based on a lamp and a pencil casting a line of shadow on a desk, and an implementation by Zigelman and Kimmel uses a narrow laser beam instead of shadow, see <nplcit id="ncit0006" npl-type="s"><text>G. Zigelman and R. Kimmel, "Fast 3D laser scanner", Dept. of Computer Science, Technion - Israel Institute of Technology, Tech. Rep. CIS-2000-07, 2000</text></nplcit>. A typical device is shown in <figref idref="f0004">Fig. 5</figref> and comprises a video camera 46 and moving laser projector 48. Using the device of <figref idref="f0004">Fig. 5</figref>, depth reconstruction is performed by analyzing the laser beam deformation as it illuminates the object. A related approach is discussed below under the heading "structured light".

0080Such a scanner can be constructed from cheap materials having a cost of between US$ 50-200. Typical scan speeds for faces achievable with such devices are within the range 15-30 sec, and the low scan speed limits the approach to laboratory and like controlled conditions.

Structured light

0081The idea of structured light relates to the projecting of a known pattern (e.g. parallel stripes) on to an object, and subsequently capturing depth information from the pattern deformation. Further details are available from <nplcit id="ncit0007" npl-type="b"><text>C. Rocchini, P. Cignoni, C. Montani, P. Pingi and R. Scopigno, A low cost 3D scanner based on structured light, EUROGRAPHICS 2001, A. Chalmers and T.-M. Rhyne (Guest Editors), Volume 20 (2001), Number 3</text></nplcit>, the contents of which are hereby incorporated by reference.

0082Typical data acquisition setup includes a CCD camera and a projector and is relatively inexpensive with current costs being in the region of US$1K - 5K. Such a data acquisition device was used by Beumier and Acheroy for 3D face recognition, and reference is made to <nplcit id="ncit0008" npl-type="s"><text>C. Beumier and M. P. Acheroy, Automatic Face Identification, Applications of Digital Image Processing XVIII, SPIE, vol. 2564, pp. 311-323, July 1995</text></nplcit>, and to <nplcit id="ncit0009" npl-type="b"><text>C. Beumier and M. P. Acheroy, Automatic Face Authentication from 3D Surface, British Machine Vision Conference BMVC 98, University of Southampton UK, 14-17 Sep, 1998, pp 449-458, 1998</text></nplcit>, the contents of which are hereby incorporated by reference.

0083The disadvantage of the structured light method is the need for controlled illumination conditions, again rendering the method inapplicable to many real life applications.

Photometric stereo

0084An alternative way of facial surface acquisition, which does not require any dedicated hardware, is surface reconstruction from photometric stereo. Photometric stereo requires the acquiring of several images of the same subject in different illumination conditions and extracting the 3D geometry from the images by assuming a Lambertian reflection model. A surface is said to exhibit Lambertian reflection if it behaves as a dull or matte surface. That is to say, for incident light from any direction, the reflection is entirely diffuse, meaning that light is reflected in equal intensity in all directions. For any given surface, the brightness depends only on the angle between the illumination source and the surface normal.

0085Reference is now made to <figref idref="f0004">Fig. 6</figref>, which is a simplified diagram showing a facial surface 50, represented as a function, viewed from a given position 52 on the <i>z</i>-axis. The object is illuminated by a source of parallel rays 54 directed along /<i>'</i>.

0086In the following, Lambertian reflection is assumed and the observed picture is thus given by <maths id="math0001" num="(1)"><math display="block"><msup><mi mathvariant="italic">I</mi><mi>i</mi></msup><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>=</mo><mi>ρ</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>〈</mo><mi>n</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>,</mo><msup><mi mathvariant="italic">l</mi><mi>i</mi></msup><mo>〉</mo><mo>,</mo></math><img file="EP1550082B1_D0001.tif" /></maths> where <i>p</i>(<i>x,y</i>) is the object albedo or reflective power, and <i>n</i>(<i>x,y</i>) is the normal to the object surface, expressed as <maths id="math0002" num="(2)"><math display="block"><mi>n</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>=</mo><mfrac><mfenced open="[" close="]"><mo>-</mo><msub><mi>z</mi><mi>x</mi></msub><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>,</mo><mo>-</mo><msub><mi>z</mi><mi>y</mi></msub><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>,</mo><mn>1</mn></mfenced><msqrt><mn>1</mn><mo>+</mo><msubsup><mrow><mo>‖</mo><mo>∇</mo><mi>z</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>‖</mo></mrow><mn>2</mn><mn>2</mn></msubsup></msqrt></mfrac><mn>.</mn></math><img file="EP1550082B1_D0002.tif" /></maths>

0087Using matrix-vector notation, equation (1) can be rewritten as <maths id="math0003" num="(3)"><math display="block"><mi>I</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>=</mo><mi mathvariant="normal">L</mi><mo>⁢</mo><mi>v</mi></math><img file="EP1550082B1_D0003.tif" /></maths> where <maths id="math0004" num="(4)"><math display="block"><mi mathvariant="normal">L</mi><mo>=</mo><mo>-</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msubsup><mi>l</mi><mn>1</mn><mn>1</mn></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>2</mn><mn>1</mn></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>3</mn><mn>1</mn></msubsup></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msubsup><mi>l</mi><mn>1</mn><mi>N</mi></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>2</mn><mi>N</mi></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>3</mn><mi>N</mi></msubsup></mtd></mtr></mtable></mfenced><mo>;</mo><mmultiscripts><mrow><mi>I</mi><mfenced><mi>x</mi><mi>y</mi></mfenced></mrow><mprescripts /><mspace width="1em" /><none /></mmultiscripts><mo>=</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msup><mi>I</mi><mn>1</mn></msup><mfenced><mi>x</mi><mi>y</mi></mfenced></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msup><mi>I</mi><mi>N</mi></msup><mfenced><mi>x</mi><mi>y</mi></mfenced></mtd></mtr></mtable></mfenced></math><img file="EP1550082B1_D0004.tif" /></maths> and <maths id="math0005" num="(5)"><math display="block"><msub><mi>v</mi><mn>1</mn></msub><mo>=</mo><msub><mi>z</mi><mi>x</mi></msub><mo>⁢</mo><msub><mi>v</mi><mn>3</mn></msub><mo>;</mo><msub><mi>v</mi><mn>2</mn></msub><mo>=</mo><msub><mi>z</mi><mi>y</mi></msub><mo>⁢</mo><msub><mi>v</mi><mn>3</mn></msub><mo>;</mo><msub><mi>v</mi><mn>3</mn></msub><mo>=</mo><mfrac><mrow><mi>ρ</mi><mfenced><mi>x</mi><mi>y</mi></mfenced></mrow><msqrt><mn>1</mn><mo>+</mo><msubsup><mrow><mo>‖</mo><mo>∇</mo><mi>z</mi><mo>‖</mo></mrow><mn>2</mn><mn>2</mn></msubsup></msqrt></mfrac><mn>.</mn></math><img file="EP1550082B1_D0005.tif" /></maths>

0088Given at least 3 linearly independent illuminations <maths id="math0006"><math display="inline"><msubsup><mfenced open="{" close="}"><msup><mi>l</mi><mi>i</mi></msup></mfenced><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mo>,</mo></math><img file="EP1550082B1_D0006.tif" /></maths> and the corresponding observations <maths id="math0007"><math display="inline"><msubsup><mfenced open="{" close="}"><msup><mi>I</mi><mi>i</mi></msup></mfenced><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mo>,</mo></math><img file="EP1550082B1_D0007.tif" /></maths> one can reconstruct the values of ∇<i>z</i> by pointwise least-squares solution <maths id="math0008" num="(6)"><math display="block"><mi>v</mi><mo>=</mo><msup><mi>L</mi><mi>†</mi></msup><mo>⁢</mo><mi>I</mi><mfenced><mi>x</mi><mi>y</mi></mfenced></math><img file="EP1550082B1_D0008.tif" /></maths> where L⁺ = (L^TL)^-1 L^T denotes the Moore-Penrose pseudoinverse. It is noted that the Moore-Penrose pseudoinverse is a substitute for the inverse matrix, and is useful in circumstances in which a standard inverse matrix does not exist. A function for obtaining the Moore-Penrose pseudoinverse of a matrix is provided for example as a standard feature in Matlab.

0089Having the gradient ∇<i>z</i>, the surface 50 can be reconstructed by minimization of the following function: <maths id="math0009" num="(7)"><math display="block"><mover><mi>z</mi><mo>‾</mo></mover><mo>=</mo><munder><mi>arg min</mi><mover><mi>z</mi><mo>‾</mo></mover></munder><mspace width="1em" /><mo mathvariant="normal">∬</mo><mfenced open="[" close="]"><msup><mfenced><msub><mover><mi>z</mi><mo>‾</mo></mover><mi>x</mi></msub><mo>-</mo><msub><mi>z</mi><mi>x</mi></msub></mfenced><mn>2</mn></msup><mo>+</mo><msup><mfenced><msub><mover><mi>z</mi><mo>‾</mo></mover><mi>y</mi></msub><mo>-</mo><msub><mi>z</mi><mi>y</mi></msub></mfenced><mn>2</mn></msup></mfenced><mo>⁢</mo><mi mathvariant="italic">dxdy</mi></math><img file="EP1550082B1_D0009.tif" /></maths>

0090The Euler-Lagrange conditions of equation (7), in turn, allow rewriting of the Poisson equation thus, <maths id="math0010" num="(8)"><math display="block"><msub><mover><mi mathvariant="italic">z</mi><mo>‾</mo></mover><mi mathvariant="italic">xx</mi></msub><mo mathvariant="italic">+</mo><msub><mover><mi mathvariant="italic">z</mi><mo>‾</mo></mover><mi mathvariant="italic">yy</mi></msub><mo mathvariant="italic">=</mo><msub><mi mathvariant="italic">z</mi><mi mathvariant="italic">xx</mi></msub><mo mathvariant="italic">+</mo><msub><mi mathvariant="italic">z</mi><mi mathvariant="italic">yy</mi></msub></math><img file="EP1550082B1_D0010.tif" /></maths><ul id="ul0010" list-style="none" compact="compact"><li>[1] the solution of which yields a reconstructed surface <i>z̃</i>. See R. Kimmel, <i>Numerical geometry of images</i>, Lecture notes.</li></ul>

0091One of the obvious problems in the surface reconstruction from photometric stereo approach arises from deviations from the Lambertian model in real faces. In face recognition applications, some face features (such as eyes, hair, beard etc.) have a strongly non-Lambertian nature.

0092To reduce the irregularities in the reconstructed surface, one can add a penalty term on surface non-smoothness to the least-squares solution, for example to give a total variation: <maths id="math0011" num="(8a)"><math display="block"><mover><mi>z</mi><mo>‾</mo></mover><mo>=</mo><munder><mi>arg min</mi><mover><mi>z</mi><mo>‾</mo></mover></munder><mo>∬</mo><mfenced open="[" close="]"><msup><mfenced><msub><mover><mi>z</mi><mo>‾</mo></mover><mi>x</mi></msub><mo>-</mo><msub><mi>z</mi><mi>x</mi></msub></mfenced><mn>2</mn></msup><mo>+</mo><msup><mfenced><msub><mover><mi>z</mi><mo>‾</mo></mover><mi>y</mi></msub><mo>-</mo><msub><mi>z</mi><mi>y</mi></msub></mfenced><mn>2</mn></msup></mfenced><mo>⁢</mo><mi mathvariant="italic">dxdy</mi><mo>+</mo><mi>λ</mi><mo>∬</mo><msub><mrow><mo>‖</mo><mo>∇</mo><mi>z</mi><mo>‖</mo></mrow><mn>2</mn></msub><mo>⁢</mo><mi mathvariant="italic">dxdy</mi></math><img file="EP1550082B1_D0011.tif" /></maths>

Geometric Face Recognition Scheme

0093The face recognition scheme of the present embodiments is based on the geometric features of the three-dimensional facial surface. The geometric features may be obtained either directly from a 3D scan and incorporated as points or vertices into a triangulated manifold of the kind shown in <figref idref="f0005">Fig. 7</figref>, or the vertices may be derived from photometric stereo and then arranged to form the triangulated manifold.

0094Preliminary processing, such as centering and cropping, is preferably performed on the manifold prior to recognition. Centering and cropping actions can be carried out by simple pattern matching, which may use the eyes, or the tip of the nose 60, as the most recognizable feature of the human face. The facial contour may also be extracted in order to limit processing to the surface belonging to the face itself, thus to exclude foreground items such as spectacles and the like. Such tasks are preferably performed by the preprocessor 16 referred to above.

0095Preprocessing preferably emphasizes those sections of the face less susceptible to alteration, the so-called hard regions, including the upper outlines of the eye sockets, the areas surrounding one's cheekbones, and the sides of the mouth. On the other hand, sections, which can be easily changed (e.g. hair), the soft regions, are preferably excluded from the recognition process. As will be discussed below, regions that change significantly according to expression can also be included as soft regions and excluded, and in certain embodiments processing may be limited just to the upper part of the face.

0096Following treatment by the preprocessor, the manifold typically still comprises too many points for efficient processing and more points than are really necessary to convey the underlying facial geometry. Hence subsampling is carried out using subsampler 18 in order to form the subsampled manifold of <figref idref="f0005">Fig. 8</figref>.

0097Reference is now made to <figref idref="f0006">Fig. 9</figref>, which is a simplified block diagram showing the subsampler 18 in greater detail. The subsampler comprises an initializer 70, and a Voronoi sampler 72, and takes as inputs the preprocessed full manifold and a desired number of points, or any other user-friendly way of defining a trade-off between accuracy and computational efficiency.

0098In the subsampling procedure, a subset of <i>n</i> uniformly distributed vertices is selected from the triangulated surface within the facial contour. The sub-sampling is performed using the iterative Voronoi sampling procedure, where on each iteration a vertex with the largest geodesic distance from the already selected ones is selected. The procedure is initialized by a constant vertex selected by the initializer and the geodesic distances needed are preferably computed using FMM-TD, the mathematics of which are discussed in greater detail below. For further information on Voronoi sampling, reference is made to <nplcit id="ncit0010" npl-type="b"><text>CISM Journal ACSGC Vol. 45 No. 1. Spring 1991 pp 65-80</text></nplcit>, Problems with handling spatial data, the Voronoi approach, Christopher M. Gold. FMM-TD is discussed in greater detail below.

0099Reference is now made to <figref idref="f0007">Fig. 10</figref>, which is a simplified block diagram showing in greater detail the geodesic converter 20 of <figref idref="f0001">Fig. 1</figref>. In the geodesic converter 20, sets of points 80 are received from the subsampler 18. An <i>n</i> × <i>n</i> distance matrix of geodesic distances is then created by applying FMM-TD from each of the <i>n</i> selected vertices to each other vertex. One of the principles ensuring the low computational complexity of FMM-TD is the fact that already calculated distances from the Voronoi sampling can be reused in subsequent computations such as computing the geodesic distance matrix.

0100Using the n x n distance matrix, multi-dimensional scaling (MDS) is applied to produce a dissimilarity matrix, a map of similarities or dissimilarities that provides a canonical representation of the face in a low-dimensional Euclidean space. Since isometric transformations do not affect the geodesic distances, it may be expected that a face undergoing such a transformation yields the same canonical form as that of the original face. Studies carried out by the inventors show that indeed slight non-isometric transformations produce small but insignificant deviations from the original canonical form. MDS is discussed in greater detail below.

0101The last stage of the facial matching operation consists of matching or classification, i.e. comparison of the canonical forms. A straightforward approach is to measure a distance between two sets of points, using, for example, the Hausdorff metric. However, the Hausdorf metric based method is computationally extensive and should be therefore avoided.

0102An alternative approach, used in the present embodiments, takes the first <i>m</i> eigenvalues obtained from the MDS procedure to provide coordinates in a low-dimensional feature space. Although the dominant eigenvalues do not describe the canonical form entirely, it is reasonable that similar faces have similar eigenvalues (and thus form clusters in the feature space). Experiments, cited below, show that the described comparison is sufficiently accurate and has low computational complexity. In the present embodiments the first three eigenvalues are taken and plotted onto a three-dimensional feature space.

Fast Marching on Triangulated Manifolds

0103Face recognition according to the present embodiments uses geodesic distances between points on the facial surface, both for subsampling and for creating a distance matrix from the subsampled points. Initially, one computes distances between pairs of points on a triangulated manifold representing the facial surface, and then an efficient numerical method is required to obtain geodesic distances therefrom.

0104A method known as The Fast Marching Method (FMM), is disclosed by <nplcit id="ncit0011" npl-type="s"><text>J. A. Sethian, A fast marching level set method for monotonically advancing fronts, Proc. Nat. Acad. Sci., 93, 4, 1996</text></nplcit>. The fast marching method was subsequently extended to triangulated domains (FMM-TD) as disclosed in R. Kimmel and J. A. Sethian, Computing geodesic paths on manifolds. FMM-TD is an efficient numerical method to compute a first-order approximation of the geodesic distances.

0105Given a set of source points {<i>s_i</i>} on a manifold such as that of <figref idref="f0005">Fig. 7</figref>, the distance map <i>T</i>(<i>x</i>,<i>y</i>) from these points to other points on the manifold is obtained as the solution of the Eikonal equation <maths id="math0012" num="(9)"><math display="block"><mo>‖</mo><mo>∇</mo><mi>T</mi><mo>‖</mo><mo>=</mo><mn>1</mn><mo>;</mo><mi>T</mi><mfenced><msub><mi>s</mi><mi>i</mi></msub></mfenced><mo>=</mo><mn>0.</mn></math><img file="EP1550082B1_D0012.tif" /></maths>

0106FMM-TD allows computing of the distance map with <i>O</i>(<i>N</i>log<i>N</i>) complexity, where <i>N</i> is the number of points on the manifold.

0107When the face geometry is obtained from photometric stereo additional efficiency can be achieved as follows. There is no actual need to compute the surface itself from equation (8) above. It is sufficient to find the gradient ∇<i>z</i>, and use its values to construct the metric used by the FMM.

Multidimensional Scaling

0108Multidimensional scaling (MDS) is a method that maps measured similarities or dissimilarities among a set of objects into a representation of the pattern of proximities in a low-dimensional Euclidean space, and in this context, reference is made to <nplcit id="ncit0012" npl-type="s"><text>G. Zigelman, R. Kimmel, and N. Kiryati, Texture mapping using surface flattening via multi-dimensional scaling, Accepted to IEEE Trans. on Visualization and Computer Graphics, 2001</text></nplcit>, and <nplcit id="ncit0013" npl-type="s"><text>R. Grossmann, N. Kiryati, and R. Kimmel. Computational surface flattening: A voxel-based approach. Accepted to IEEE Trans. on PAMI, 2001</text></nplcit>.

0109Given a set of <i>n</i> objects, their mutual similarities <maths id="math0013"><math display="inline"><msubsup><mfenced open="{" close="}"><msub><mi>d</mi><mi mathvariant="italic">ij</mi></msub></mfenced><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup></math><img file="EP1550082B1_D0013.tif" /></maths> and the desired dimensionality <i>m</i>, MDS finds a set of vectors in <i>m</i>-dimensional space (each vector corresponding to an object) such that the matrix of Euclidean distances among them corresponds as closely as possible to a function of the input matrix D according to a certain criterion function.

0110In the present embodiments, proximity values are obtained by measuring the geodesic distances between points on the facial surface using the FMM-TD method. <nplcit id="ncit0014" npl-type="s"><text>A. Elad and R. Kimmel, Bending invariant representations for surfaces, Proc. of CVPR'01 Hawaii, Dec. 2001</text></nplcit>, showed that applying MDS to geodesic distances on a manifold produces a bending-invariant canonical form.

0111There exists a variety of different algorithms for solving the MDS problem; in each of them a tradeoff between computational complexity and algorithm accuracy is made. The present embodiments make use of the classical scaling algorithm introduced by Young et al. however the skilled person will appreciate the applicability of other methods.

0112Classical scaling finds the coordinates <i>x_i</i> of <i>n</i> points in a <i>k</i>-dimensional Euclidean space, given their mutual distances <maths id="math0014"><math display="inline"><msubsup><mfenced open="{" close="}"><msub><mi>d</mi><mi mathvariant="italic">ij</mi></msub></mfenced><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup><mn>.</mn></math><img file="EP1550082B1_D0014.tif" /></maths> The Euclidean distance between the points <i>i</i> and <i>j</i> is expressed by <maths id="math0015" num="(10)"><math display="block"><msubsup><mi>d</mi><mi mathvariant="italic">ij</mi><mn>2</mn></msubsup><mo>=</mo><msubsup><mrow><mo>‖</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mi>x</mi><mi>j</mi></msub><mo>‖</mo></mrow><mn>2</mn><mn>2</mn></msubsup><mo>=</mo><msup><mfenced><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mi>x</mi><mi>j</mi></msub></mfenced><mi mathvariant="normal">T</mi></msup><mo>⁢</mo><mfenced><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mi>x</mi><mi>j</mi></msub></mfenced><mn>.</mn></math><img file="EP1550082B1_D0015.tif" /></maths>

0113Given the squared-distance matrix D with elements as in (10) one can construct the inner-product matrix <maths id="math0016" num="(11)"><math display="block"><mi mathvariant="normal">B</mi><mo>=</mo><mo>-</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mo>⁢</mo><mfenced><mi mathvariant="normal">I</mi><mo>-</mo><mfrac><mn>1</mn><mi>n</mi></mfrac><mo>⁢</mo><msup><mn>11</mn><mi mathvariant="normal">T</mi></msup></mfenced><mo>⁢</mo><mi mathvariant="normal">D</mi><mo>⁢</mo><mfenced><mi mathvariant="normal">I</mi><mo>-</mo><mfrac><mn>1</mn><mi>n</mi></mfrac><mo>⁢</mo><msup><mn>11</mn><mi mathvariant="normal">T</mi></msup></mfenced><mo>=</mo><msup><mi>XX</mi><mi mathvariant="normal">T</mi></msup></math><img file="EP1550082B1_D0016.tif" /></maths> where 1=[1,...,1]^T and X = [<i>x</i>₁,...,<i>x_n</i>]^T. B is a positive semi-definite matrix with at most <i>k</i> positive eigenvalues, which can be expressed as <maths id="math0017" num="(12)"><math display="block"><mi mathvariant="normal">B</mi><mo>=</mo><msup><mi>UΛU</mi><mi mathvariant="normal">T</mi></msup><mn>.</mn></math><img file="EP1550082B1_D0017.tif" /></maths>

0114The coordinate matrix X is therefore obtained by computing the square root matrix of B <maths id="math0018" num="(13)"><math display="block"><mi mathvariant="normal">X</mi><mo>=</mo><msup><mi>UΛ</mi><mfrac><mn>1</mn><mn>2</mn></mfrac></msup><mn>.</mn></math><img file="EP1550082B1_D0018.tif" /></maths>

0115It is noted that, from a statistical point of view, the above-described approach is equivalent to principal component analysis (PCA), which finds an orthogonal basis that maximizes the variance of the given <i>n</i> vectors projected to the basis vectors.

0116Empirical observations show that three dimensions usually suffice for representation of most of the geometric structure contained in the geodesic distances measured on a smooth surface.

0117Reference is now made to <figref idref="f0008">Fig. 11</figref>, which is a simplified flow chart showing the process flow of embodiments of the present invention, and showing how the various procedures and functions described above may work together. In a first stage S1, 3D data is collated from the object it is desired to match. Collation may use any of the scanning methods described above or any other suitable method of gathering 3D topological information of a body. In an optional stage S2, albedo or reflectivity power information may be gathered. One of the possibilities of incorporating albedo information is by embedding the two-dimensional face manifold into a 4D or 6D space and measuring distances on the manifold using a combined metric, thus the 4^th coordinate in the 4D-embedding and the 4^th-6^th coordinates in the 6D-embedding represent the gray level or the RGB channels of the albedo information, respectively.

0118The data gathered is then presented in a stage S3 as a series of points on a triangulated manifold such as that of <figref idref="f0004">Fig. 6</figref>. The manifold may be constructed directly from the 3D data or may require auxiliary calculations, depending on the data gathering method used.

0119Stage S4 is a preprocessing stage. The preprocessing stage removes soft - that is often changing regions of the face geometry and retains hard regions, that is those regions that remain constant. Preprocessing may also include cropping the image and like operations. Preprocessing preferably involves locating a reference point on the face and then using the general facial form to determine the locations of the soft regions that are to be ignored.

0120Stage S5 is a subsampling stage. The preprocessed manifold is reduced by selecting only those points which are most crucial for retaining geometric information. As discussed above, efficient techniques such as Voronoi sampling may be used. Stage S5 is shown in greater detail in <figref idref="f0009">Fig. 12</figref> in which a stage S5.1 comprises input of the preprocessed manifold. A selection stage S5.2 allows a user to set a desired number of points for the subsampling in a tradeoff between accuracy and computational complexity. Finally a stage 5.3 involves Voronoi sampling, which looks for the minimum number of points having the maximum topological information. As discussed above, an initial point is taken and then a farthest point therefrom is added until the desired number of points is reached. Geodesic distances are used to determine furthest points and FMM - TD can be used to find the geodesic distances.

0121In a stage S6, a distance matrix is computed of geodesic distances between each subsampled vertex and each other point. As explained above, FMM-TD is preferably used as an efficient method of computing geodesic distances.

0122In a stage S7, MDS is used to measure dissimilarity in the distance matrix and thereby express the geodesic distances in a low dimensional Euclidean representation. Finally in a stage S8, the face as represented in low dimensional Euclidean representation can be compared with other faces by various distance calculation techniques. As discussed above, a preferred method uses the first <i>m</i> eigenvalues obtained from the MDS procedure to provide coordinates in a low-dimensional feature space. In the present embodiments the first three eigenvalues are taken and plotted on a three-dimensional feature space. Although the dominant eigenvalues do not describe the canonical form entirely, similar faces nevertheless form recognizable clusters in the feature space and thus allow matching.

Experimental Results

0123Four experiments were performed in order to evaluate the approach and validate the present embodiments. Three databases were used: <ul id="ul0011" list-style="none" compact="compact"><li>I. A database ¹ of 3D facial surfaces obtained by a simple 3D scanner of the kind described above and shown in <figref idref="f0010">Fig. 13</figref>. The database consisted of 6 subjects a..f including one artificial face (subject b) and 5 human faces. ¹The database is available from http://www.cs.technion.ac.il/∼zgil, courtesy of Gil Zigelman.</li><li>II. The Yale Face Database B ². The database consists of high-resolution grayscale images of 10 subjects of both Caucasian and Asian type, taken in controlled illumination conditions. Each subject was depicted in several poses (slight head rotations). 3D surface reconstruction was carried out using photometric stereo. ²The database is available from http://cvc.yale.edu/projects/yalefacesB/ yalefacesB.html</li><li>III. A Database of high-resolution facial surfaces and textures of human subjects with different facial expressions, obtained using the 3D scanner.</li></ul>

Scanned surfaces with artificial facial expressions

0124The first experiment was performed on database I (<figref idref="f0010">Fig 13</figref>) and was intended to test algorithm sensitivity to artificial facial expressions.

0125Reference is now made to <figref idref="f0010">Fig. 14</figref> which shows how one of the faces of <figref idref="f0010">Fig. 13</figref> was deformed to simulate facial expressions. Face a is undeformed face d of <figref idref="f0010">Fig. 13</figref>, and faces b and c are different deformations applied thereto. In b, the forehead, cheek and chins of subject d were deformed, simulating different facial expressions of the same subject. In c, a feature, the nose, was substituted with that of a different face. The latter was applied to test algorithm sensitivity to major feature alteration It is emphasized that these deformations were not isometric.

0126The surfaces were scaled to the size 60×95 and downsampled to 400 points using Voronoi sampling. Afterwards, the geodesic distance from each point was computed using FMM-TD and multidimensional scaling was applied to the distance matrix.

0127Reference is now made to <figref idref="f0011">Fig. 15</figref>, which is a graph of a 3-dimensional feature space for matching. More specifically, the graph of <figref idref="f0011">Fig. 15</figref> depicts the feature-space representation on which the first 3 eigenvalues obtained by MDS may be plotted. Each point in the feature space corresponds to one of the faces in <figref idref="f0010">Figs 13 and 14</figref>, viz A (dot), B (circle), C (cross), D (pentagram), E (triangle) and F (hexagram). It can be seen that face deformations, though not purely isometric, produce closely located points that form clusters. On the other hand, faces obtained by feature substitution are far from the original face of subject D.

Photometric stereo

0128The second experiment was performed on database II. The aims of the experiment were to determine the feasibility of using photometric stereo for 3D surface acquisition and the sensitivity of the algorithm to head rotations.

0129It appears that faithful surface reconstruction from photometric stereo demands a large number of different illumination directions. The reason is primarily due to deviations in practice from the Lambertian model (especially notable in forehead, hair and eyes) and slight subject displacement during image acquisition (e.g. eye blinking). Reference is now made to <figref idref="f0012">Fig. 16</figref>, which shows nine different illuminations of the same face and a least squares surface reconstruction made by surface reconstruction from all nine of the images. Nine is simply a number found empirically to be an efficient number. Smaller numbers of images resulted in less accurate surfaces and therefore less accurate recognition.

0130Reference is now made to <figref idref="f0013">Fig. 17</figref>, which depicts a series of ten facial surfaces a-j each reconstructed from photometric stereo from nine images of respective faces. It is noted that relatively high noise appears to be inevitable when reconstructing 3D geometry from photometric stereo.

0131Reference is now made to <figref idref="f0014">Fig. 18</figref>, which shows subject G of <figref idref="f0013">fig. 17</figref> taken at several poses. For some subjects, several poses were taken to test the recognition sensitivity to small head rotations. Since the algorithm is bending invariant, the influence of rotation (an isometric transformation) should be negligible, and all the possible deviations found seem mainly to be due to inaccurate surface reconstruction.

0132In the experiment, faces were manually centered according to the eye and lip locations and processing was performed in the region bounded by the facial contour. The surfaces were scaled to the size 56×80 and sampled at 400 points using Voronoi sampling. Afterwards, the distance from each point was computed using fast marching, and multidimensional scaling was applied to the distance matrix.

0133Reference is now made to <figref idref="f0014">Fig. 19</figref>, which is a simplified graph depicting the feature-space representation using the first 3 eigenvalues obtained by singular value decomposition of the faces of <figref idref="f0013">Fig. 17</figref>, some of them in several poses as exemplified by <figref idref="f0014">Fig. 18</figref>. A (dot), B (circle), C (cross), D (plus), E (star), F (pentagram), G (hexagram), H (diamond), I (square) and J (triangle). Different poses of subjects A, B, D, E and G form clusters. It is clear from the results of the second experiment that different poses result in relatively narrow clusters.

Surface matching

0134The third experiment was performed on database I, the aim of the experiment being to compare the geometric approach to a straightforward surface comparison.

0135The experiment uses a simple L₂-norm-based measure for surface similarity. Reference is now made to <figref idref="f0015">Fig. 20</figref>, which shows two pairs of surfaces 90 and 92. Each surface has a plane 94, 94' defined by the two eyeball centers and the center of the mouth. Pair of surfaces 90 are not aligned and pair of surfaces 92 are the same as those in 90 but after alignment. In the experiment, the surfaces were aligned in two different ways: <ol id="ol0001" compact="compact"><li>(a) as shown in <figref idref="f0015">Figure 20</figref>, so that the planes 94, 94' formed by the two eyeball centers and the center of the mouth in the two surfaces coincided, and</li><li>(b) by finding the Euclidean transformation <maths id="math0019" num="(14)"><math display="block"><mi mathvariant="italic">Rx</mi><mo>+</mo><mi>b</mi></math><img file="EP1550082B1_D0019.tif" /></maths></li></ol> in which <i>R</i> is a rotation matrix, and <i>b</i> is a translation vector, so as to minimize the sum of Euclidean distances between the points on the two surfaces, that is the rotation is selected that minimizes the surface dissimilarity (in this case, the approximated volume in the gap between two aligned surfaces). Optimal alignment was found by the gradient descent algorithm with numeric gradient computation.

0136Once an optimal alignment is found, the proximity of surfaces is computed by integrating the volume of the gap between the two surfaces. Such an approach follows the global surface matching idea disclosed in <nplcit id="ncit0015" npl-type="b"><text>C. Beumier and M. P. Acheroy, Automatic Face Authentication from 3D Surface, British Machine Vision Conference BMVC 98, University of Southampton UK, 14-17 Sep, 1998, pp 449-458</text></nplcit>.", 1998.

0137For visualization purposes and for displaying the subjects as points in 3D space, a matrix D is formed of the distances between the surfaces according to <maths id="math0020" num="(15)"><math display="block"><msub><mi>d</mi><mi mathvariant="italic">ij</mi></msub><mo>=</mo><munderover><mi mathvariant="normal">Σ</mi><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><msub><mrow><mo>‖</mo><msubsup><mi>x</mi><mi>k</mi><mi>i</mi></msubsup><mo>-</mo><msubsup><mi>x</mi><mi>k</mi><mi>j</mi></msubsup><mo>‖</mo></mrow><mn>2</mn></msub></math><img file="EP1550082B1_D0020.tif" /></maths> where <maths id="math0021"><math display="inline"><msubsup><mfenced open="{" close="}"><msubsup><mi>x</mi><mi>k</mi><mi>i</mi></msubsup></mfenced><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup></math><img file="EP1550082B1_D0021.tif" /></maths> are the N points forming the <i>i</i>-th surface, and applying MDS to D. The first three dominant eigenvalues were taken in order to obtain the representation of the proximity pattern in a three-dimensional Euclidean space.

0138Reference is now made to <figref idref="f0015">Figs 21 and 22</figref>, which are two graphs summing up the results of this third experiment. It is apparent that the straightforward comparison of manifolds after alignment is sensitive to facial expressions and fails to distinguish between different expressions of the same subject on the one hand and different subjects on the other hand.

0139More particularly, <figref idref="f0015">Figs. 21 and 22</figref> show that unlike the geodesic distances, L₂ distances between the surfaces do not distinguish between faces, even using an optimal surface alignment. Different expressions and head poses of the same subject do not form recognizable clusters. The reason is the fact that the straightforward comparison does not take into consideration the geometry of the facial surfaces.

0140<figref idref="f0015">Fig. 21</figref> shows the results of an L₂-based algorithm for alignment according to 3 points for the set of faces of <figref idref="f0013">Fig. 17</figref> and some of the faces showing different expressions: A (dot), B (circle), C (cross), D (plus), E (star), F (pentagram), G (hexagram), H (diamond), I (square) and J (triangle). No clusters are recognizable.

0141<figref idref="f0015">Fig. 22</figref> uses an L₂-based algorithm with optimal alignment for the same set of faces: A (dot), B (circle), C (cross), D (plus), E (star), F (pentagram), G (hexagram), H (diamond), I (square) and J (triangle). Again, no distinct clusters appear.

Human subjects with real expressions

0142A fourth experiment was performed, this time using database III. The main goal was testing the algorithm in field conditions, on subjects with real facial expressions. Three subjects from database III are shown in <figref idref="f0016">Fig. 23, and Fig. 24</figref> shows a single subject with three different expressions.

0143Reference is now made to <figref idref="f0016">Figs 25 - 28</figref> which show preprocessing of the faces to extract only features that are relevant. The faces were processed in order to extract the relevant features only. A geodesic mask was applied both to the texture and the depth maps. The mask was obtained by computing a ball of constant radius -in the geodesic sense-around the tip of the nose. In the experiment the nose location was found manually. <figref idref="f0016">Figs. 25 and 26</figref> are texture maps, and <figref idref="f0016">Figs 27 and 28</figref> are depth maps. <figref idref="f0016">Figs. 25 and 27</figref> show texture and depth maps respectively before application of the mask and <figref idref="f0016">Figs. 26 and 28</figref> show texture and depth maps respectively after application of the mask.

0144Reference is now made to <figref idref="f0017">Figs. 29, 30 and 31</figref>, which respectively show mesh processing stages of the present embodiments. <figref idref="f0017">Fig. 29</figref> shows the high-resolution triangulated manifold, <figref idref="f0017">Fig. 30</figref> shows the downsampled version of the manifold and <figref idref="f0017">Fig. 31</figref> shows the surface in canonical form.

0145The minimal mean squared distance measure yielded much better results applied to the canonical surface rather than to the original manifold. Currently, the minimal mean squared distance measure appears to be computationally expensive. Canonical form comparison, based on central moments which are much simpler to compute, yielded a bit less accurate results than the minimal squared distance, but still provided a good measure for face similarity.

0146The classical scaling algorithm was used for creating the canonical forms. The present inventors have seen no significant improvement over the classical scaling algorithm in using the least-square MDS algorithm. 250-500 subsamples have been found to give sufficient accuracy while keeping the computational complexity relatively low.

0147It is noted that the present embodiments yield poor results when considering strong non-isometric deformations of the facial surface such as inflated cheeks or open mouth. Limiting the region of interest to the upper part of the face, being a hard region and containing substantially rigid features, on the preprocessing stage yields significantly better results under the above conditions.

Geodesic Distance Calculations Using a Metric Tensor.

0148The above-described embodiment describes a 3D face recognition approach based on geometric invariant signatures. The key idea of the embodiment is a representation of the facial surface, invariant to isometric deformations, such as those resulting from facial expressions. One of the stages in the construction of the geometric invariants is the measurement of geodesic distances on triangulated surfaces, which, as described above, is carried out by fast marching on triangulated domains (FMTD). The embodiment to be described below uses only the metric tensor of the surface for geodesic distance computation. When combined with photometric stereo for facial surface acquisition, the embodiment permits construction of a bending-invariant representation of the face without reconstructing the 3D the surface.

0149In the present embodiment, a variation is provided of the above-described fast marching on triangulated domains (FMTD). As mentioned in the previous paragraph, the variation is capable of computing geodesic distances given only the metric tensor of the surface, and does not require an actual 3D reconstruction of the facial surface, thus making the method simpler and less prone to errors. The embodiment is preferably combined with the photometric stereo method for facial surface acquisition. Photometric stereo is a cheap and simple approach, producing the metric tensor without reconstructing the surface. As a result, a simple and fast face recognition method is obtained.

Surface Acquisition

0150The face recognition algorithm of the present embodiment treats faces as three-dimensional surfaces. It is therefore necessary to obtain first the facial surface of the subject that we are trying to recognize. The present embodiment focuses on methods that reproduce the surface gradient. As will be shown hereinbelow, the actual surface reconstruction is not needed, saving computational effort and reducing numerical errors.

A. Photometric stereo

0151Reference is now made to <figref idref="f0018">Fig. 32</figref>, which shows a face 100, a viewing position 102 and a lighting source 104. The photometric stereo technique consists of obtaining several pictures of the same subject in different illumination conditions and extracting the 3D geometry by assuming the Lambertian reflection model. We assume that the facial surface, represented as a function, is viewed from given viewing position 102 along the z-axis. The object is illuminated by a source of parallel rays directed along <i>_il</i> from 104.

0152In the following, the equation numbering is started again from 1.

0153We assume a Lambertian reflection model, i.e. the observed image is given by <maths id="math0022" num="(1)"><math display="block"><msup><mi>I</mi><mi>i</mi></msup><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>=</mo><mi>ρ</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>⁢</mo><mi>n</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>⋅</mo><msup><mi>l</mi><mi>i</mi></msup></math><img file="EP1550082B1_D0022.tif" /></maths> where <i>ρ(x,y)</i> is the object albedo, and <i>n</i>(<i>x,y</i>) is the normal to the object surface, expressed as <maths id="math0023" num="(2)"><math display="block"><mi>n</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>=</mo><mfrac><mfenced><mo>-</mo><msub><mi>z</mi><mi>x</mi></msub><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>,</mo><mo>-</mo><msub><mi>z</mi><mi>y</mi></msub><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>,</mo><mn>1</mn></mfenced><msqrt><mn>1</mn><mo>+</mo><msubsup><mrow><mo>‖</mo><mo>∇</mo><mi>z</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>‖</mo></mrow><mn>2</mn><mn>2</mn></msubsup></msqrt></mfrac></math><img file="EP1550082B1_D0023.tif" /></maths>

0154Using matrix-vector notation, (1) can be rewritten as <maths id="math0024"><math display="block"><mi mathvariant="italic">I</mi><mfenced><mi mathvariant="italic">x</mi><mi mathvariant="italic">y</mi></mfenced><mo mathvariant="italic">=</mo><mi mathvariant="normal">L</mi><mo mathvariant="normal">⁢</mo><mi mathvariant="italic">v</mi></math><img file="EP1550082B1_D0024.tif" /></maths> where <maths id="math0025" num="(4)"><math display="block"><mi mathvariant="normal">L</mi><mo>=</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msubsup><mi>l</mi><mn>1</mn><mn>1</mn></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>2</mn><mn>1</mn></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>3</mn><mn>1</mn></msubsup></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msubsup><mi>l</mi><mn>1</mn><mi>N</mi></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>2</mn><mi>N</mi></msubsup></mtd><mtd><msubsup><mi>l</mi><mn>3</mn><mi>N</mi></msubsup></mtd></mtr></mtable></mfenced><mo>;</mo><mi>I</mi><mfenced><mi>x</mi><mi>y</mi></mfenced><mo>=</mo><mfenced open="[" close="]"><mtable><mtr><mtd><msup><mi>I</mi><mn>1</mn></msup><mfenced><mi>x</mi><mi>y</mi></mfenced></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msup><mi>I</mi><mi>N</mi></msup><mfenced><mi>x</mi><mi>y</mi></mfenced></mtd></mtr></mtable></mfenced></math><img file="EP1550082B1_D0025.tif" /></maths> and <maths id="math0026" num="(5)"><math display="block"><msub><mi>v</mi><mn>1</mn></msub><mo>=</mo><mo>-</mo><msub><mi>z</mi><mi>x</mi></msub><mo>⁢</mo><msub><mi>v</mi><mn>3</mn></msub><mo>;</mo><msub><mi>v</mi><mn>2</mn></msub><mo>=</mo><mo>-</mo><msub><mi>z</mi><mi>y</mi></msub><mo>⁢</mo><msub><mi>v</mi><mn>3</mn></msub><mo>;</mo><msub><mi>v</mi><mn>3</mn></msub><mo>=</mo><mfrac><mrow><mi>ρ</mi><mfenced><mi>x</mi><mi>y</mi></mfenced></mrow><msqrt><mn>1</mn><mo>+</mo><msubsup><mrow><mo>‖</mo><mo>∇</mo><mi>z</mi><mo>‖</mo></mrow><mn>2</mn><mn>2</mn></msubsup></msqrt></mfrac><mn>.</mn></math><img file="EP1550082B1_D0026.tif" /></maths>

0155Given at least 3 linearly independent illuminations <maths id="math0027"><math display="inline"><msubsup><mfenced open="{" close="}"><msup><mi>l</mi><mi>i</mi></msup></mfenced><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mo>,</mo></math><img file="EP1550082B1_D0027.tif" /></maths> and the corresponding observations <maths id="math0028"><math display="inline"><msubsup><mfenced open="{" close="}"><msup><mi>I</mi><mi>i</mi></msup></mfenced><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mo>,</mo></math><img file="EP1550082B1_D0028.tif" /></maths> and <i>v</i>=L^†<i>I</i>(<i>x,y</i>) one can reconstruct the values of ∇<i>z</i> by pointwise least-squares solution where L^† = (L^TL)^-1L^T denotes the Moore-Penrose pseudoinverse. When needed, the surface can be reconstructed by solving the Poisson equation <maths id="math0029"><math display="block"><msub><mover><mi mathvariant="italic">z</mi><mo>‾</mo></mover><mi mathvariant="italic">xx</mi></msub><mo mathvariant="italic">+</mo><msub><mover><mi mathvariant="italic">z</mi><mo>‾</mo></mover><mi mathvariant="italic">yy</mi></msub><mo mathvariant="italic">=</mo><msub><mi mathvariant="italic">z</mi><mi mathvariant="italic">xx</mi></msub><mo mathvariant="italic">+</mo><msub><mi mathvariant="italic">z</mi><mi mathvariant="italic">yy</mi></msub></math><img file="EP1550082B1_D0029.tif" /></maths> with respect to ∇<i>z</i>. In the present embodiment we adopt the photometric stereo approach due to its simplicity, however it will be appreciated that the invention is in no way so limited.

B. Structured light

0156<nplcit id="ncit0016" npl-type="s"><text>M. Proesmans, L. Van Gool and A. Oosterlinck, One-Shot Active Shape Acquisition, Proc. Internat. Conf. Pattern Recognition, Vienna, Vol. C, pp. 336- 340, 1996</text></nplcit>; and <nplcit id="ncit0017" npl-type="b"><text>S. Winkelbach and F. M. Wahl, Shape from 2D Edge Gradients, Pattern Recognition, Lecture Notes in Computer Sciences 2191, pp. 377-384, Springer, 2001</text></nplcit>, proposed a shape from 2D edge gradients reconstruction technique, which permits a reconstruction of the surface normals (gradients) from two stripe patterns projected onto the object. The reconstruction technique is based on the fact that directions of the projected stripes in the captured 2D images depend on the local orientation of the surface in 3D. Classical edge-detecting operators, such as Sobel, Canny, etc can be used to find the direction of the stripe edges.

0157Reference is now made to <figref idref="f0019">Fig. 33</figref> which is a simplified diagram illustrating the relationship between surface gradient and local stripe direction. A pixel in the image plane defines a viewing vector <i>s</i>. The stripe direction determines the stripe direction vector <i>v'</i>, lying in both the image plane and in the viewing plane. The real tangential vector of projected stripe <i>v</i>₁ is perpendicular to the normal <i>v</i> = <i>s</i>' × <i>c</i> of the viewing plane and to the normal <i>p</i> of the stripe projection plane. Assuming parallel projection, we obtain <maths id="math0030" num="(8)"><math display="block"><msub><mi>v</mi><mn>1</mn></msub><mo>=</mo><mi>c</mi><mo>×</mo><mi>p</mi></math><img file="EP1550082B1_D0030.tif" /></maths>

0158Acquiring a second image of the scene with a second stripe illumination rotated relative to the first image, allows a second tangential vector <i>v</i>₂ to be calculated. Next, the surface normal is computed according to <maths id="math0031" num="(9)"><math display="block"><mi>n</mi><mo>=</mo><msub><mi>v</mi><mn>1</mn></msub><mo>×</mo><msub><mi>v</mi><mn>2</mn></msub><mn>.</mn></math><img file="EP1550082B1_D0031.tif" /></maths>

0159In <nplcit id="ncit0018" npl-type="b"><text>S. Winkelbach and F. M. Wahl, Shape from Single Stripe Pattern Illumination, L. Van Gool (Ed): Pattern Recognition (DAGM 2002), Lecture Notes in Computer Science 2449, pp. 240-247, Springer, 2002</text></nplcit>, the authors use a single lighting pattern to estimate the surface normal from the local directions and <i>widths</i> of the projected stripes.

Bending-Invariant Representation

0160Classical surface matching methods, based on finding an Euclidean transformation of two surfaces which maximizes some shape similarity criterion, are suitable mainly for rigid objects. Use of such methods is detailed in <nplcit id="ncit0019" npl-type="s"><text>O. D. Faugeras and M. Hebert, A 3D recognition and positioning algorithm using geometrical matching between primitive surfaces, Proc. 7th Int'l Joint Conf. on Artificial Intelligence, pp. 996-1002, 1983</text></nplcit>; further in <nplcit id="ncit0020" npl-type="b"><text>P. J. Besl, The free form matching problem. Machine vision for threedimensional scene, Freeman, H. (ed.) New York Academic, 1990</text></nplcit>, and also in <nplcit id="ncit0021" npl-type="s"><text>G. Barequet and M. Sharir, Recovering the position and orientation of free-form objects from image contours using 3D distance map, IEEE Trans. PAMI, 19(9), pp. 929-948, 1997</text></nplcit>.

0161As clear from the present disclosure, the human face can by no means be considered as a rigid object since it undergoes deformations resulting from facial expressions. On the other hand, the class of transformations that a facial surface can undergo is not arbitrary, and empirical observations show that facial expressions can be modeled as <i>isometric</i> (or length preserving) transformations. Such transformations do not stretch and nor do they tear the surface. More rigorously, it is possible to say that the transformations that the human face undergoes preserve the surface metric. The surfaces resulting from such transformations are called <i>isometric surfaces</i>. The requirement of a deformable surface matching algorithm is to find a representation, which is the same for all isometric surfaces.

0162<nplcit id="ncit0022" npl-type="s"><text>E. L. Schwartz, A. Shaw and E. Wolfson, A numerical solution to the generalized mapmaker's problem: flattening non-convex polyhedral surfaces, IEEE Trans. PAMI, 11, pp. 1005-1008, 1989</text></nplcit>, were the first to introduce the use of multidimensional scaling (MDS) as a tool for studying curved surfaces by planar models. <nplcit id="ncit0023" npl-type="s"><text>G. Zigelman, R. Kimmel and N. Kiryati, Texture mapping using surface flattening via multi-dimensional scaling, IEEE Trans. Visualization and Comp. Graphics, 8, pp. 198-207, 2002</text></nplcit>; and <nplcit id="ncit0024" npl-type="s"><text>R. Grossman, N. Kiryati and R. Kimmel, Computational surface flattening: a voxel-based approach, IEEE Trans. PAMI, 24, pp. 433-441, 2002</text></nplcit> extended some of these ideas to the problem of texture mapping and voxel-based cortex flattening.

0163A generalization of this approach was introduced in <nplcit id="ncit0025" npl-type="s"><text>A. Elad, R. Kimmel, Bending invariant representations for surfaces, Proc. CVPR, pp. 168-174, 2001</text></nplcit>, as a framework for object recognition. They showed an efficient algorithm for constructing a representation of surfaces, invariant under isometric transformations. The method, described therein as bending-invariant canonical forms, is utilized as a convenient representation in the present embodiment.

0164Given a polyhedral approximation of the facial surface, S. one can think of such an approximation as having been obtained by sampling the underlying continuous surface on a finite set of points <i>p_i</i>(<i>i</i> = 1,...,<i>n</i>), and discretizing the metric δ associated with the surface <maths id="math0032" num="(10)"><math display="block"><mi>δ</mi><mfenced><msub><mi>p</mi><mi>i</mi></msub><msub><mi>p</mi><mi>j</mi></msub></mfenced><mo>=</mo><msub><mi>δ</mi><mi mathvariant="italic">ij</mi></msub><mn>.</mn></math><img file="EP1550082B1_D0032.tif" /></maths>

0165Writing the values of δ<i>_ij</i> in matrix form, we obtain a matrix of mutual distances between the surface points. For convenience, we define the <i>squared</i> mutual distances, <maths id="math0033" num="(11)"><math display="block"><msub><mfenced><mi mathvariant="normal">Δ</mi></mfenced><mi mathvariant="italic">ij</mi></msub><mo>=</mo><msubsup><mi>δ</mi><mi mathvariant="italic">ij</mi><mn>2</mn></msubsup><mn>.</mn></math><img file="EP1550082B1_D0033.tif" /></maths>

0166The matrix Δ is invariant under isometric surface deformations, but is not a unique representation of isometric surfaces, since it depends on arbitrary ordering and the selection of the surface points. We would like to obtain a geometric invariant, which would be unique for isometric surfaces on one hand, and which, on the other hand allows use of simple rigid surface matching algorithms to compare such invariants. Treating the squared mutual distances as a particular case of <i>dissimilarities</i>, one can apply the dimensionality reduction technique, mentioned in respect of the previous embodiments and called <i>multidimensional scaling</i> (MDS), in order to embed the surface into a low-dimensional Euclidean space <b>R</b><i>_m</i>. This is equivalent to finding a mapping between two metric spaces, <maths id="math0034" num="(12)"><math display="block"><mi>ϕ</mi><mo>:</mo><mfenced><mi mathvariant="normal">S</mi><mi>δ</mi></mfenced><mo>→</mo><mfenced><msup><mi mathvariant="normal">R</mi><mi>m</mi></msup><mi>d</mi></mfenced><mo>;</mo><mi>ϕ</mi><mfenced><msub><mi>p</mi><mi>i</mi></msub></mfenced><mo>=</mo><msub><mi>x</mi><mi>i</mi></msub><mo>,</mo></math><img file="EP1550082B1_D0034.tif" /></maths> that minimizes the embedding error, <maths id="math0035" num="(13)"><math display="block"><mi>ε</mi><mo>=</mo><mi>f</mi><mo>⁢</mo><mfenced><mfenced open="|" close="|"><msub><mi>δ</mi><mi mathvariant="italic">ij</mi></msub><mo>-</mo><msub><mi>d</mi><mi mathvariant="italic">ij</mi></msub></mfenced></mfenced><mo>;</mo><msub><mi>d</mi><mi mathvariant="italic">ij</mi></msub><mo>=</mo><msub><mrow><mo>‖</mo><msub><mi>x</mi><mi>i</mi></msub><mo>-</mo><msub><mi>x</mi><mi>j</mi></msub><mo>‖</mo></mrow><mn>2</mn></msub><mo>,</mo></math><img file="EP1550082B1_D0035.tif" /></maths> for some monotonous function <i>f</i>.

0167The obtained <i>m</i>-dimensional representation is a set of points <i>x_i</i> ∈ <b>R</b><i>_m</i> (<i>i</i> = 1,...,<i>n</i>), corresponding to the surface points <i>p_i</i>.

0168Different MDS methods can be derived using different embedding error criteria, thus <nplcit id="ncit0026" npl-type="b"><text>I. Borg and P. Groenen, Modem multidimensional scaling - theory and applications, Springer, 1997</text></nplcit>. A particular case is <i>classical scaling</i>, introduced by <nplcit id="ncit0027" npl-type="s"><text>G. Young and G. S. Householder, Discussion of a set of points in terms of their mutual distances, Psychometrika 3, pp. 19-22, 1938</text></nplcit>. The embedding in <b>R</b><i>_m</i> is performed by double-centering the matrix <maths id="math0036" num="(14)"><math display="block"><mi mathvariant="normal">B</mi><mo>=</mo><mo>-</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mo>⁢</mo><mi>JΔJ</mi></math><img file="EP1550082B1_D0036.tif" /></maths> where J = I-_{1 <i>n</i>}U; I is a <i>n</i>×<i>n</i> identity matrix, and U is a matrix consisting entirely of ones. The first <i>m</i> eigenvectors <i>e_i</i>, corresponding to the <i>m</i> largest eigenvalues of B, are used as the embedding coordinates <maths id="math0037" num="(15)"><math display="block"><msubsup><mi>x</mi><mi>i</mi><mi>j</mi></msubsup><mo>=</mo><msubsup><mi>e</mi><mi>i</mi><mi>j</mi></msubsup><mo>⁢</mo><msub><mo>;</mo><mspace width="1em" /></msub><mo>⁢</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>…</mo><mo>,</mo><mmultiscripts><mi>n</mi><mprescripts /><mspace width="1em" /><none /></mmultiscripts><mo>;</mo><mi>j</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>…</mo><mo>,</mo><mi>m</mi><mo>,</mo></math><img file="EP1550082B1_D0037.tif" /></maths><i>where</i><maths id="math0038"><math display="inline"><msubsup><mi>x</mi><mi>i</mi><mi>j</mi></msubsup></math><img file="EP1550082B1_D0038.tif" /></maths> denotes the <i>j</i>-th coordinate of the vector <i>x_i</i>. We refer to the set of points <i>x_i</i> obtained by the MDS as the <i>bending invariant canonical form</i> of the surface. When <i>m</i>=3, the set of points can be plotted as a surface. Standard rigid surface matching methods can be used in order to compare between two deformable surfaces, using their bending-invariant representations instead of the surfaces themselves. Since the canonical form is computed up to a translation, rotation, and reflection transformation, to allow comparison between canonical forms, they must be <i>aligned</i>. This can be done, for instance, by setting the first-order moments, typically centers of mass, and the mixed second-order moments to zero, as described in <nplcit id="ncit0028" npl-type="s"><text>A. Tal, M. Elad and S. Ar, Content based retrieval of VRML objects - an iterative and interactive approach", EG Multimedia, 97, pp. 97-108, 2001</text></nplcit>.

Measuring Geodesic Distances

0169One of the crucial steps in the construction of the canonical form of a given surface, is an efficient algorithm for the computation of geodesic distances on surfaces, that is, δ<i>_ij</i>. A numerically consistent algorithm for distance computation on triangulated domains, henceforth referred to as <i>fast marching on triangulated domains</i> (FMTD), was used by Elad and Kimmel and has been discussed hereinabove in respect of the previous embodiments. FMTD was proposed by Kimmel and Sethian as a generalization of the <i>fast marching method</i>.

0170Using FMTD, the geodesic distances between a surface vertex and the rest of the <i>n</i> surface vertices can be computed in O(<i>n</i>) operations. Measuring distances on manifolds has been carried out before for graphs of functions and implicit manifolds, as discussed in <nplcit id="ncit0029" npl-type="s"><text>J. Sethian and A. Vladimirsky, Ordered upwind methods for static Hamilton-Jacobi equations: theory and applications. Technical Report PAM 792 (UCB), Center for Pure and Applied Mathematics, 2001 Submitted for publication to SIAM Journal on Numerical Analysis, July 2001</text></nplcit>; and also <nplcit id="ncit0030" npl-type="s"><text>F. M'moli and G. Sapiro, Fast computation of weighted distance functions and geodesics on implicit hyper-surfaces. Journal of Computational Physics, 173(2), pp. 730-764, 2001</text></nplcit>.

0171By contrast with the above citations however, the focus of the present embodiment is to avoid making a surface reconstruction. Thus, the present embodiment makes use of a modified version of FMTD, which computes the geodesic distances on a surface, using the values of the surface gradient Vz only. The surface gradient values can be obtained, for example, from photometric stereo or structured light.

0172Reference is now made to <figref idref="f0020">Fig. 34</figref>, which is a simplified diagram illustrating an orthogonal grid U 120 projected via transform X onto a manifold 122 to form a non-orthogonal grid X(U) 124.

0173The facial surface can be thought of as a parametric manifold, represented by a mapping ₂₃X: → <b>R R</b> from the parameterization plane U=(<i>u</i>¹,<i>u</i>²)=(<i>x,</i>y) to the parametric manifold <maths id="math0039" num="16"><math display="block"><mi mathvariant="normal">X</mi><mfenced><mi mathvariant="normal">U</mi></mfenced><mo>=</mo><mfenced><msup><mi>x</mi><mn>1</mn></msup><mfenced><msup><mi>u</mi><mn>1</mn></msup><msup><mi>u</mi><mn>2</mn></msup></mfenced><mo>,</mo><msup><mi>x</mi><mn>2</mn></msup><mfenced><msup><mi>u</mi><mn>1</mn></msup><msup><mi>u</mi><mn>2</mn></msup></mfenced><mo>,</mo><msup><mi>x</mi><mn>3</mn></msup><mfenced><msup><mi>u</mi><mn>1</mn></msup><msup><mi>u</mi><mn>2</mn></msup></mfenced></mfenced><mo>;</mo></math><img file="EP1550082B1_D0039.tif" /></maths> which, in turn, can be written as <maths id="math0040" num="(17)"><math display="block"><mi mathvariant="normal">X</mi><mfenced><mi mathvariant="normal">U</mi></mfenced><mo>=</mo><mfenced><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mfenced><mi>x</mi><mi>y</mi></mfenced></mfenced><mn>.</mn></math><img file="EP1550082B1_D0040.tif" /></maths>

0174The derivatives of X with respect to <i>u_i</i> are defined as <maths id="math0041"><math display="block"><msub><mi>X</mi><mi>i</mi></msub><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msup><mi>u</mi><mi>i</mi></msup></mrow></mfrac><mo>⁢</mo><mi mathvariant="normal">X</mi></math><img file="EP1550082B1_D0041.tif" /></maths> and they constitute non-orthogonal coordinate system 124 on the manifold as illustrated. In the particular case of equation (17) above, <maths id="math0042" num="(18)"><math display="block"><msub><mi>X</mi><mn>1</mn></msub><mo>=</mo><mfenced><mn>1</mn><mn>0</mn><msub><mi>z</mi><mi>x</mi></msub></mfenced><mo>,</mo><msub><mi>X</mi><mn>2</mn></msub><mo>=</mo><mfenced><mn>0</mn><mn>1</mn><msub><mi>z</mi><mi>y</mi></msub></mfenced><mn>.</mn></math><img file="EP1550082B1_D0042.tif" /></maths>

0175The distance element on the manifold is <maths id="math0043" num="(19)"><math display="block"><mi mathvariant="italic">ds</mi><mo>=</mo><msqrt><msub><mi>g</mi><mi mathvariant="italic">ij</mi></msub><mo>⁢</mo><msup><mi>u</mi><mi>i</mi></msup><mo>⁢</mo><msup><mi>u</mi><mi>j</mi></msup></msqrt><mo>,</mo></math><img file="EP1550082B1_D0043.tif" /></maths> where we use Einstein's summation convention, and the metric tensor <i>g_ij</i> of the manifold is given by <maths id="math0044" num="(20)"><math display="block"><mfenced><msub><mi>g</mi><mi mathvariant="italic">ij</mi></msub></mfenced><mo>=</mo><mfenced><mtable><mtr><mtd><msub><mi>g</mi><mn>11</mn></msub></mtd><mtd><msub><mi>g</mi><mn>12</mn></msub></mtd></mtr><mtr><mtd><msub><mi>g</mi><mn>21</mn></msub></mtd><mtd><msub><mi>g</mi><mn>22</mn></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced><mtable><mtr><mtd><msub><mi>X</mi><mn>1</mn></msub><mo>⋅</mo><msub><mi>X</mi><mn>1</mn></msub></mtd><mtd><msub><mi>X</mi><mn>1</mn></msub><mo>⋅</mo><msub><mi>X</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><msub><mi>X</mi><mn>2</mn></msub><mo>⋅</mo><msub><mi>X</mi><mn>1</mn></msub></mtd><mtd><msub><mi>X</mi><mn>2</mn></msub><mo>⋅</mo><msub><mi>X</mi><mn>2</mn></msub></mtd></mtr></mtable></mfenced><mn>.</mn></math><img file="EP1550082B1_D0044.tif" /></maths>

0176The classical fast marching method, discussed in <nplcit id="ncit0031" npl-type="s"><text>J. A. Sethian, A review of the theory, algorithms, and applications of level set method for propagating surfaces. Acta numerica, pp. 309-395, 1996</text></nplcit>, calculates distances in an orthogonal coordinate system. The numerical stencil for the update of a grid point consists of the vertices of a right angled triangle. In the present case, g12 <img file="EP1550082B1_D0045.tif" /> and the resulting triangles are not necessarily right angled. If a grid point is updated by a stencil which is an obtuse triangle, a problem may thus arise. The values of one of the points of the stencil might not be set in time and cannot be used. There is a similar problem with fast marching on triangulated domains which include obtuse triangles. The problem, and a possible solution, is discussed further in <nplcit id="ncit0032" npl-type="s"><text>R. Kimmel, and J. A. Sethian, Computing geodesic on manifolds. Proc. US National Academy of Science 95, pp. 8431-8435, 1998</text></nplcit>.

0177Our solution is similar to that of the above reference. In the present embodiment, we perform a preprocessing stage for the grid, in which we split every obtuse triangle into two acute ones, as illustrated in <figref idref="f0021">Fig. 35</figref>. <figref idref="f0004">Fig. 5</figref> illustrates two original triangles labeled 1 and 2. Triangle 1 gives a proper numerical support, yet triangle 2 is obtuse. Thus triangle 2 is replaced by triangle 3 and triangle 4.The split is performed by adding an additional edge, connecting the updated grid point with a non-neighboring grid point. The distant grid point becomes part of the numerical stencil. The need for splitting is determined according to the angle between the non-orthogonal axes at the grid point. It is calculated by <maths id="math0045" num="21"><math display="block"><mi>cos</mi><mspace width="1em" /><mi>α</mi><mo>=</mo><mfenced><mfrac><mtable><mtr><mtd><msub><mi>X</mi><mn>1</mn></msub><mo>⁢</mo><msub><mi>X</mi><mn>2</mn></msub></mtd></mtr></mtable><mtable><mtr><mtd><mo>‖</mo><msub><mi>X</mi><mn>1</mn></msub><mo>‖</mo><mo>‖</mo><msub><mi>X</mi><mn>2</mn></msub><mo>‖</mo></mtd></mtr></mtable></mfrac></mfenced><mo>=</mo><mfrac><msub><mi>g</mi><mn>12</mn></msub><msqrt><msub><mi>g</mi><mn>11</mn></msub><mo>⁢</mo><msub><mi>g</mi><mn>22</mn></msub></msqrt></mfrac></math><img file="EP1550082B1_D0046.tif" /></maths>

0178If cos α= 0, the axes are perpendicular, and no splitting is required. If cos α&lt;0, the angle α is obtuse and should be split. The denominator of (21) is always positive, so we need only check the sign of the numerator g12.

0179In order to split an angle, we connect the updated grid point with another point, located <i>m</i> grid points from the point in the <i>X</i>₁ direction, and <i>n</i> grid points in the <i>X</i>₂ direction. It is noted that <i>m</i> and <i>n</i> can be negative. The newly connected point is a proper supporting point if the obtuse angle is split into two acute ones. For cos α<0 this is the case if <maths id="math0046" num="(22)"><math display="block"><mtable columnalign="left"><mtr><mtd><mi>cos</mi><mspace width="1em" /><msub><mi>β</mi><mn>1</mn></msub></mtd><mtd><mo>=</mo><mfenced><mfrac><mrow><msub><mi>X</mi><mn>1</mn></msub><mo>⁢</mo><mfenced><mi>m</mi><mo>⁢</mo><msub><mi>X</mi><mn>1</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>X</mi><mn>2</mn></msub></mfenced></mrow><mrow><mo>‖</mo><msub><mi>X</mi><mn>1</mn></msub><mo>‖</mo><mo>‖</mo><mi>m</mi><mo>⁢</mo><msub><mi>X</mi><mn>1</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>X</mi><mn>2</mn></msub><mo>‖</mo></mrow></mfrac></mfenced><mo>=</mo></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mo>=</mo><mfrac><mrow><mi>m</mi><mo>⁢</mo><msub><mi>g</mi><mn>11</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>g</mi><mn>12</mn></msub></mrow><msqrt><msub><mi>g</mi><mn>11</mn></msub><mo>⁢</mo><mfenced><msup><mi>m</mi><mn>2</mn></msup><mo>⁢</mo><msub><mi>g</mi><mn>11</mn></msub><mo>+</mo><mn>2</mn><mo>⁢</mo><mi mathvariant="italic">mn</mi><mo>⁢</mo><msub><mi>g</mi><mn>12</mn></msub><mo>+</mo><msup><mi>n</mi><mn>2</mn></msup><mo>⁢</mo><msub><mi>g</mi><mn>22</mn></msub></mfenced></msqrt></mfrac><mo>></mo><mn>0</mn></mtd></mtr></mtable></math><img file="EP1550082B1_D0047.tif" /></maths> and <maths id="math0047" num="(23)"><math display="block"><mtable columnalign="left"><mtr><mtd><mi>cos</mi><mspace width="1em" /><msub><mi>β</mi><mn>2</mn></msub></mtd><mtd><mo>=</mo><mfenced><mfrac><mrow><msub><mi>X</mi><mn>2</mn></msub><mo>⁢</mo><mfenced><mi>m</mi><mo>⁢</mo><msub><mi>X</mi><mn>1</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>X</mi><mn>2</mn></msub></mfenced></mrow><mrow><mo>‖</mo><msub><mi>X</mi><mn>2</mn></msub><mo>‖</mo><mo>‖</mo><mi>m</mi><mo>⁢</mo><msub><mi>X</mi><mn>1</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>X</mi><mn>2</mn></msub><mo>‖</mo></mrow></mfrac></mfenced><mo>=</mo></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mo>=</mo><mfrac><mrow><mi>m</mi><mo>⁢</mo><msub><mi>g</mi><mn>12</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>g</mi><mn>22</mn></msub></mrow><msqrt><msub><mi>g</mi><mn>22</mn></msub><mo>⁢</mo><mfenced><msup><mi>m</mi><mn>2</mn></msup><mo>⁢</mo><msub><mi>g</mi><mn>11</mn></msub><mo>+</mo><mn>2</mn><mo>⁢</mo><mi mathvariant="italic">mn</mi><mo>⁢</mo><msub><mi>g</mi><mn>12</mn></msub><mo>+</mo><msup><mi>n</mi><mn>2</mn></msup><mo>⁢</mo><msub><mi>g</mi><mn>22</mn></msub></mfenced></msqrt></mfrac><mo>></mo><mn>0.</mn></mtd></mtr></mtable></math><img file="EP1550082B1_D0048.tif" /></maths>

0180Also, it is enough to check the sign of the numerators. For cos α>0, cos β₂ changes its sign and the constraints are <maths id="math0048" num="(24)"><math display="block"><mi>m</mi><mo>⁢</mo><msub><mi>g</mi><mn>11</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>g</mi><mn>12</mn></msub><mo>></mo><mn>0</mn><mo>,</mo><mi>m</mi><mo>⁢</mo><msub><mi>g</mi><mn>12</mn></msub><mo>+</mo><mi>n</mi><mo>⁢</mo><msub><mi>g</mi><mn>22</mn></msub><mo><</mo><mn>0.</mn></math><img file="EP1550082B1_D0049.tif" /></maths>

0181This process is done for all grid points. Once the preprocessing stage is done, we have a suitable numerical stencil for each grid point and we can calculate the distances.

0182The numerical scheme used is similar to that of Kimmel and Sethian as referenced above, with the exception that there is no need to perform the unfolding step. The supporting grid points that split the obtuse angles can be found more efficiently. The required triangle edge lengths and angles are calculated according to the surface metric <i>g_ij</i> at the grid point, which, in turn, is computed using the surface gradients <i>z_x</i>, <i>z_y</i>. A more detailed description appears in A. Spira and R. Kimmel, An efficient solution to the eikonal equation on parametric manifolds, submitted for publication, 2003, the contents of which are hereby incorporated by reference.

3d Face Recognition

0183Reference is now made to <figref idref="f0023">Fig. 37B</figref>, which is a simplified block diagram illustrating apparatus for carrying out the present embodiment. A 3-D scanning system 126 carries out 3D data acquisition in order to obtain surface gradient data. A preprocessor 128 carries out preprocessing of the raw data. Preliminary processing, such as centering and cropping can be carried out by simple pattern matching, which can use the eyes as the most recognizable feature of the human face, as discussed with the previous embodiments.

0184The facial contour is preferably also be extracted in order to limit processing to the surface belonging to the face itself. That is to say, the initial scanning may include regions not belonging to the face itself.

0185As discussed in respect of the previous embodiments, preprocessing preferably emphasizes those sections of the face less susceptible to alteration and excludes the parts that can be easily changed, such as the hair. In experimental work using the present embodiment, the preprocessing stage was limited to cropping the triangulated manifold by removing the parts lying outside an ellipse (in the geodesic sense) centered at the nose tip.

0186FMTD unit 130 applies FMTD to distances from each of n selected vertices to points on the face to form an nxn dissimilarity distance matrix. MDS unit 132 then applies multi-dimensional scaling to the matrix to form a low dimensional canonical form of the surface, again as discussed with the previous embodiments. The canonical form is then used as an input to comparator 134 to compare with other faces represented as canonical forms in database 136.

0187Reference is now made to <figref idref="f0023">Fig. 37B</figref> which is a simplified flow diagram illustrating the process of 3D Face recognition without facial reconstruction. As a first step S140, the 3D face recognition system according to the present embodiment acquires the surface gradient Vz, as discussed hereinabove. At stage S142, the raw data are preprocessed, as discussed in respect of <figref idref="f0022">Fig. 37A</figref>. Briefly, centering and cropping can be carried out by simple pattern matching, using the eyes as the most recognizable feature of the human face.

0188The facial contour may also be extracted in order to limit the processing to the surface belonging only to the face itself.

0189As discussed, preprocessing preferably emphasizes those sections of the face less susceptible to alteration and excludes the parts that can be easily changed, such as the hair. In experimental work using the present embodiment, the preprocessing stage was limited to cropping the triangulated manifold by removing the parts lying outside an ellipse (in the geodesic sense) centered at the nose tip.

0190Next, in stage S144 an <i>n</i>×<i>n</i> dissimilarity matrix of geodesic distances, hereinafter the distance matrix, is created by applying FMTD from each of the n selected vertices. Then, in stage S 146, MDS - multi-dimensional scaling - is applied to the distance matrix, producing a canonical form of the face in a low-dimensional Euclidean space. In the above referred to experimental work a three-dimensional space was used. The canonical form obtained in the above manner is then compared with a database of templates corresponding to other subjects in a one-to many type matching process, or to a single stored template in a one-to-one type match. If the correspondence of the compared canonical form falls within a certain statistical range of values, the match is considered to be valid.

0191Finding correspondence between the canonical forms brings us back to the surface-matching problem discussed in respect of the previous embodiments. An advantage of the present embodiments in general is the fact that matching is performed between the rigid canonical forms carrying the intrinsic object geometry, rather than the deformable facial surfaces themselves. Thus even the simplest rigid surface matching algorithms produce plausible results.

0192In the present embodiment, we adopted the moments system referred to above in reference to A Tal et al, due to its simplicity. The canonical form's (<i>p, q, r</i>)-th moment is given by <maths id="math0049" num="(25)"><math display="block"><msub><mi>M</mi><mi mathvariant="italic">pqr</mi></msub><mo>=</mo><munder><mi mathvariant="normal">Σ</mi><mi>n</mi></munder><mo>⁢</mo><msup><mfenced><msubsup><mi>x</mi><mi>n</mi><mn>1</mn></msubsup></mfenced><mi>p</mi></msup><mo>⁢</mo><msup><mfenced><msubsup><mi>x</mi><mi>n</mi><mn>2</mn></msubsup></mfenced><mi>q</mi></msup><mo>⁢</mo><msup><mfenced><msubsup><mi>x</mi><mi>n</mi><mn>3</mn></msubsup></mfenced><mi>r</mi></msup></math><img file="EP1550082B1_D0050.tif" /></maths> where <maths id="math0050"><math display="inline"><msubsup><mi>x</mi><mi>n</mi><mi>i</mi></msubsup></math><img file="EP1550082B1_D0051.tif" /></maths> denotes the <i>i</i>_th coordinate of the <i>n</i>_th point in the canonical surface samples. In order to compare between two canonical forms, the vector (<i>M_{p<sub2>1</sub2>q<sub2>1</sub2>r<sub2>1</sub2>}</i>,...,<i>M_{p<sub2>M</sub2>q<sub2>M</sub2>r<sub2>M</sub2>}</i>), termed the <i>moments signature</i>, is computed for each surface. The Euclidean distance between two moments signatures measures the dissimilarity between the two surfaces.

Experimental Results

0193As already referred to, experimental work was performed in accordance with the above-described embodiment, to demonstrate that comparison of canonical forms obtained without actual facial surface reconstruction is better than reconstruction and direct comparison of the surfaces. As in experimental work on the previous embodiments, the Yale Face Database B was used. The database consisted of high-resolution grayscale images of different instances of 10 subjects of both Caucasian and Asian type, taken in controlled illumination conditions. Reference is now made to <figref idref="f0021">Fig. 36</figref> which shows a single subject illuminated from five different angles. For the five cases, 36A..36E, illumination angles respectively of (0°,0°) (0°,-20°) (0°,+20°) (-25°,0°) (+25°,0°) are used. Numbers in brackets indicate the azimuth and the elevation angles, respectively, thereby determining the illumination direction. The experiments used instances of the kind shown in <figref idref="f0021">Fig. 36</figref> taken from seven different subjects taken from the database.

0194Direct surface matching consisted of the retrieval of the surface gradient according to equation (6) above using five different illumination directions, reconstruction of the surface according to equation (7), and alignment and computation of the surface moments signature according to equation (25). Canonical forms were computed from the surface gradient, aligned and converted into a moment signature according to equation (25).

0195Reference is now made to <figref idref="f0024">Fig. 38</figref>, which shows three successive stages in processing of the surface data. <figref idref="f0024">Fig. 38A</figref> shows a surface gradient field, <figref idref="f0024">Fig. 38B</figref> shows a reconstructed surface, and <figref idref="f0024">Fig. 38C</figref> shows the bending-invariant canonical form which is represented as a surface. It is stressed that <figref idref="f0024">Fig. 38</figref> is shown for explanatory purposes only and, as explained, the present embodiment is able to use the canonical form represented by <figref idref="f0024">Fig. 38C</figref> without actually reconstructing the surface at any of the processing stages.

0196Reference is now made to <figref idref="f0025">Fig. 39</figref>, which shows a set of faces in the center surrounded by 3D plots of the faces and distances between different faces according to the present embodiments. Visualization of the face recognition results is shown in terms of three-dimensional proximity patterns. Subjects from the face database are represented as points obtained by applying MDS to the relative distances between subjects. Shown in the graph are both straightforward surface matching (A) and canonical forms (B). That is to say, in order to obtain an impression of the accuracy of the embodiment, we converted the relative distances between the subjects produced by each algorithm into 3D proximity patterns. The proximity patterns, representing each subject as a point in R³, were obtained by applying MDS to the relative distances, producing a distortion less than 1%.

0197The entire cloud of dots was partitioned into clusters formed by instances of the seven subjects, presently referred to as C₁-C₇. Visually, the more C<i>_i</i> are compact and distant from other clusters, the more accurate is the algorithm. Quantitatively, we measured <ul id="ul0012" list-style="none" compact="compact"><li>(i) the variance σ<i>_i</i> of C<i>_i</i> and</li><li>(ii) the distance <i>d_i</i> between the centroid of C<i>_i</i> and the centroid of the nearest cluster.</li></ul>

0198Table I shows a quantitative comparison of the algorithms. Inter-cluster distances <i>di</i> are given in the units of of variance σ<i>_i</i>. The table is limited to the results of C₁ to C₄. <tables id="tabl0001" num="0001"><table frame="bottom"><title>TABLE I</title><tgroup cols="5" colsep="0"><colspec colnum="1" colname="col1" colwidth="20mm" /><colspec colnum="2" colname="col2" colwidth="17mm" /><colspec colnum="3" colname="col3" colwidth="16mm" /><colspec colnum="4" colname="col4" colwidth="22mm" /><colspec colnum="5" colname="col5" colwidth="22mm" /><thead><row><entry namest="col1" nameend="col5" align="center" valign="top">FACE RECOGNITION ACCURACY</entry></row><row><entry align="center" valign="top">CLUSTER</entry><entry align="center" valign="top">σ_DIRECT</entry><entry align="center" valign="top"><i>d</i>_DIRECT</entry><entry align="center" valign="top">σ_CANONICAL</entry><entry align="center" valign="top"><i>d</i>_CANONICAL</entry></row></thead><tbody><row rowsep="0"><entry align="center">C₁</entry><entry align="char" char="." charoff="11">0.1749</entry><entry align="char" char="." charoff="12">0.1704</entry><entry align="char" char="." charoff="9">0.0140</entry><entry align="char" char="." charoff="9">4.3714</entry></row><row rowsep="0"><entry align="center">C₂</entry><entry align="char" char="." charoff="11">0.2828</entry><entry align="char" char="." charoff="12">0.3745</entry><entry align="char" char="." charoff="9">0.0120</entry><entry align="char" char="." charoff="9">5.1000</entry></row><row rowsep="0"><entry align="center">C₃</entry><entry align="char" char="." charoff="11">0.0695</entry><entry align="char" char="." charoff="12">0.8676</entry><entry align="char" char="." charoff="9">0.0269</entry><entry align="char" char="." charoff="9">2.3569</entry></row><row><entry align="center">C₄</entry><entry align="char" char="." charoff="11">0.0764</entry><entry align="char" char="." charoff="12">0.7814</entry><entry align="char" char="." charoff="9">0.0139</entry><entry align="char" char="." charoff="9">4.5611</entry></row></tbody></tgroup></table></tables>

0199The use of canonical forms improved the cluster variance and the inter-cluster distance by about one order of magnitude, compared to direct facial surface matching.

0200The present embodiment with associated experiment demonstrates how to perform 3D face recognition without reconstructing the facial surface. 3D face recognition based on bending-invariant representations, unlike previously proposed solutions, makes face recognition robust to facial expressions, head orientations and illumination conditions. The present embodiment improves on the previous embodiments discussed herein by allowing an efficient use of simple 3D acquisition techniques (e.g. photometric stereo) for fast and accurate face recognition.

3D face recognition using eigenforms

0201Reference is now made to <figref idref="f0026">Fig. 40</figref>, which is a simplified flow diagram illustrating a further preferred embodiment of the present invention. In the embodiment of <figref idref="f0026">Fig. 40</figref>, as a first step, S 150, a range camera is used to acquire a 3D image. The acquired image in fact comprises two separate images, the range image or three-dimensional geometry of the face, and the 2D image or texture of the face. The range image is converted into a triangulated surface in step 152 and smoothed using splining in stage 154. Regions outside the facial contour are cropped in stage 156, and the surface is decimated to the size of approximately 2000-2500 vertices in stage 158. Next, in stage 160, the bending-invariant canonical form of the face is computed and aligned using the procedure described in respect of the previous embodiments.

0202It will be appreciated that there is full correspondence between the texture image pixels and the canonical surface vertices <maths id="math0051"><math display="inline"><mfenced><msubsup><mi>x</mi><mi>n</mi><mn>1</mn></msubsup><msubsup><mi>x</mi><mi>n</mi><mn>2</mn></msubsup><msubsup><mi>x</mi><mi>n</mi><mn>3</mn></msubsup></mfenced><mn>.</mn></math><img file="EP1550082B1_D0052.tif" /></maths> Thus the face texture image may be mapped onto the aligned canonical surface in the canonical form space, as shown in S162 - <figref idref="f0027">Fig. 41</figref>. In stage S164 an interpolation is carried out of the texture data and the range data <maths id="math0052"><math display="inline"><msubsup><mi>x</mi><mi>n</mi><mn>3</mn></msubsup></math><img file="EP1550082B1_D0053.tif" /></maths> onto a Cartesian grid in the X1X2 plane. The result is to obtain flattened texture <i>ã</i> and the canonical image <i>x̃</i> respectively. Both <i>ã</i> and <i>x̃</i> conserve the invariance of the canonical form to isometric transformations, and can be represented as images. Reference is made in this respect to <figref idref="f0028">Fig. 42</figref>, which illustrates texture flattening in association with the canonical image. More particularly, the figure shows texture flattening by interpolation onto the X1X2 plane as follows: texture mapping on the facial surface (A) and on the canonical form (B); the resulting flattened texture (C) and the canonical image (D).

0203Application of eigendecomposition is straightforward in this representation. As with eigenfaces, we have a training set, which is now a set of duplets of the form <maths id="math0053"><math display="block"><msubsup><mfenced open="{" close="}"><msub><mover><mi>x</mi><mo>˜</mo></mover><mi>n</mi></msub><msub><mover><mi>a</mi><mo>˜</mo></mover><mi>n</mi></msub></mfenced><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mn>.</mn></math><img file="EP1550082B1_D0054.tif" /></maths>

0204Returning to <figref idref="f0027">Fig. 41</figref>, and in stage 166, we apply eigendecomposition separately on the set of ã and <i>x̃</i>. The result is to produce two sets of eigenspaces corresponding to the flattened textures and the canonical images. We term the respective sets of eigenvectors _a e<i>_n</i> and _x e<i>_n eigenforms</i>.

0205For a new subject represented by (<i>x̃</i>',<i>ã</i>'), the decomposition coefficients are computed according to <maths id="math0054"><math display="block"><mi>α</mi><mo>=</mo><mfenced open="[" close="]"><msup><msub><mi>e</mi><mn>1</mn></msub><mi>a</mi></msup><mo>…</mo><msubsup><mi>e</mi><mi>N</mi><mi>a</mi></msubsup></mfenced><mo>⁢</mo><mfenced><mover><mi>a</mi><mo>˜</mo></mover><mo>⁢</mo><mi>ʹ</mi><mo>-</mo><mover><mi>a</mi><mo>‾</mo></mover></mfenced><mn>.</mn></math><img file="EP1550082B1_D0055.tif" /></maths><maths id="math0055"><math display="block"><mi>β</mi><mo>=</mo><mfenced open="[" close="]"><msubsup><msub><mi>e</mi><mn>1</mn></msub><mspace width="1em" /><mi>x</mi></msubsup><mo>…</mo><msubsup><mi>e</mi><mi>N</mi><mi>x</mi></msubsup></mfenced><mo>⁢</mo><mfenced><mover><mi>x</mi><mo>˜</mo></mover><mo>⁢</mo><mi>ʹ</mi><mo>-</mo><mover><mi>x</mi><mo>‾</mo></mover></mfenced><mn>.</mn></math><img file="EP1550082B1_D0056.tif" /></maths> where a and x denote the average of <i>ã_n</i> and <i>x̃_n</i> in the training set, respectively. The distance between two subjects represented by (<i>x̃</i>₁, <i>ã</i>₁) and (<i>x̃</i>₂',<i>ã</i>₂') may then be computed, stage 168, as the weighted Euclidean distance between the corresponding decomposition coefficients, (α₁, β₁) and (α₂, β₂).

0206It is appreciated that certain features of the invention, which are, for clarity, described in the context of separate embodiments, may also be provided in combination in a single embodiment. Conversely, various features of the invention, which are, for brevity, described in the context of a single embodiment, may also be provided separately or in any suitable subcombination.