Intel Patent | Method And System Of Point Cloud Registration For Image Processing

Patent: Method And System Of Point Cloud Registration For Image Processing

Publication Number: 20200388004

Publication Date: 20201210

Applicants: Intel

Abstract

A system, article, and method of point cloud registration using overlap regions for image processing.

BACKGROUND

[0001] Camera arrays, which may be provided on tablets or smartphones for example, may be provided to capture multiple images of the same scene except from different angles. This same effect is achieved with a single camera that is moved to capture images of the same scene from different angles. These images can then be used to generate a depth map, and in turn, 3D geometric or semantic models in order to accurately locate objects in a 3D space of the scene. This is performed so that objects in the scene, or distances between objects in the scene (or from the camera to the object) can be measured, for computer vision, artificial intelligence, object recognition, head mounted 3D displays (HMDs) or point of view displays, and otherwise whenever it is desirable to know the size, position, or identity of an object in a scene.

[0002] Many such conventional systems generate a point cloud for each image and then match the points from different images of the scene except from the different angles. Once such points are matched, triangulation or other algorithms can be used to measure depth from the camera to the points and other measurements. This is often performed using an iterative closest point (ICP) method. These methods, however, are typically inaccurate for systems that compare images along a sequence of frames generated from a moving camera. When the motion from one camera position to the next is too large, the system cannot accurately determine the correct corresponding point pairs, resulting in a high error rate. Such systems also are too slow since each iteration usually includes computations with all of the points in the point clouds being registered, which may result in pauses when a user views the images or delays when automatic actions are to be performed based on location and/or identification of 3D objects in the scene.

DESCRIPTION OF THE FIGURES

[0003] The material described herein is illustrated by way of example and not by way of limitation in the accompanying figures. For simplicity and clarity of illustration, elements illustrated in the figures are not necessarily drawn to scale. For example, the dimensions of some elements may be exaggerated relative to other elements for clarity. Further, where considered appropriate, reference labels have been repeated among the figures to indicate corresponding or analogous elements. In the figures:

[0004] FIGS. 1A and 1B is a pair of pictures of a scene with objects formed by point clouds and taken from a moving camera;

[0005] FIG. 2 is a flow diagram of a method of point cloud registration in accordance with at least one of the implementations disclosed herein;

[0006] FIG. 3 is a flow chart of a method of point cloud registration for image processing in accordance with at least one of the implementations disclosed herein;

[0007] FIG. 4 is a flow diagram of a method of principal axis extraction in accordance with at least one of the implementations disclosed herein;

[0008] FIGS. 5A-5G is a detailed flow chart of a method of point cloud registration for image processing in accordance with at least one of the implementations disclosed herein;

[0009] FIG. 5H is a schematic diagram to explain principal axes alternatives on a cuboid overlap region in accordance with at least one of the implementations disclosed herein;

[0010] FIG. 6 is an image of a resulting reconstructed scene using point cloud matching of a conventional method;

[0011] FIG. 7 is an image of a resulting reconstructed scene using point cloud matching of a conventional method;

[0012] FIG. 8 is an image of a resulting reconstructed scene using point cloud matching of a conventional method;

[0013] FIG. 9 is an image of a resulting reconstructed scene using point cloud matching of a conventional method;

[0014] FIG. 10** is an image of a resulting reconstructed scene using point cloud matching of a conventional method**

[0015] FIG. 11 is an image of a resulting reconstructed scene using the method of point cloud registration for image processing in accordance with at least one of the implementations disclosed herein;

[0016] FIG. 12 is an illustrative diagram of an example system;

[0017] FIG. 13 is an illustrative diagram of another example system;* and*

[0018] FIG. 14 illustrates another example device, all arranged in accordance with at least some implementations of the present disclosure.

DETAILED DESCRIPTION

[0019] One or more implementations are now described with reference to the enclosed figures. While specific configurations and arrangements are discussed, it should be understood that this is performed for illustrative purposes only. Persons skilled in the relevant art will recognize that other configurations and arrangements may be employed without departing from the spirit and scope of the description. It will be apparent to those skilled in the relevant art that techniques and/or arrangements described herein also may be employed in a variety of other systems and applications other than what is described herein.

[0020] While the following description sets forth various implementations that may be manifested in architectures such as system-on-a-chip (SoC) architectures for example, implementation of the techniques and/or arrangements described herein are not restricted to particular architectures and/or computing systems and may be implemented by any architecture and/or computing system for similar purposes. For instance, various architectures employing, for example, multiple integrated circuit (IC) chips and/or packages, and/or various computing devices and/or consumer electronic (CE) devices such as imaging devices, digital cameras, smart phones, webcams, video game panels or consoles, set top boxes, tablets with single or multiple cameras, wearables such as HMDs, robots, and so forth, may implement the techniques and/or arrangements described herein. Further, while the following description may set forth numerous specific details such as logic implementations, types and interrelationships of system components, logic partitioning/integration choices, and so forth, claimed subject matter may be practiced without such specific details. In other instances, some material such as, for example, control structures and full software instruction sequences, may not be shown in detail in order not to obscure the material disclosed herein. The material disclosed herein may be implemented in hardware, firmware, software, or any combination thereof.

[0021] The material disclosed herein may also be implemented as instructions stored on a machine-readable medium or memory, which may be read and executed by one or more processors. A machine-readable medium may include any medium and/or mechanism for storing or transmitting information in a form readable by a machine (for example, a computing device). For example, a machine-readable medium may include read-only memory (ROM); random access memory (RAM); magnetic disk storage media; optical storage media; flash memory devices; electrical, optical, acoustical or other forms of propagated signals (e.g., carrier waves, infrared signals, digital signals, and so forth), and others. In another form, a non-transitory article, such as a non-transitory computer readable medium, may be used with any of the examples mentioned above or other examples except that it does not include a transitory signal per se. It does include those elements other than a signal per se that may hold data temporarily in a “transitory” fashion such as RAM and so forth.

[0022] References in the specification to “one implementation”, “an implementation”, “an example implementation”, and so forth, indicate that the implementation described may include a particular feature, structure, or characteristic, but every implementation may not necessarily include the particular feature, structure, or characteristic. Moreover, such phrases are not necessarily referring to the same implementation. Further, when a particular feature, structure, or characteristic is described in connection with an implementation, it is submitted that it is within the knowledge of one skilled in the art to affect such feature, structure, or characteristic in connection with other implementations whether or not explicitly described herein.

[0023] Systems, articles, and methods to provide point cloud registration for image processing.

[0024] As mentioned above, point cloud registration may be used to match the same point on an object in multiple captured images either from multiple cameras in a camera array or a single camera being moved about the scene being captured. The point cloud registration then may be used for rectification for example. The rectification is performed to transform the images from multiple cameras (or multiple images from a single camera), and therefore multiple angles, into a common plane to simplify the feature matching process. The result is that the multiple images are now each viewed from the same perspective or angle so that the images can be placed in the same common plane, the points in the images are much more easily matched by parallel matching lines. Whether or not rectified depending on the application being used, the captured images then may be used for 3D reconstruction to form a 3D depth map, model, or other 3D space, or for other specific 3D reconstruction algorithms such as simultaneous localization and mapping (SLAM) to name one example. The captured images also may be used for object recognition to identify, recognize, and/or measure objects in the 3D space, and so forth for many different applications such as artificial intelligence (AI), robotics and computer vision such as with autonomous vehicles and residential service robots, and point of view (POV) systems such as with virtual, augmented, and mixed reality HMDs, to name a few examples.

[0025] The task of point cloud registration is to calculate a transformation matrix by minimizing an alignment error between two point clouds. The iterative closest point (ICP) is the most popular technique on point cloud registration, see P. J. Best et al., “A method for registration of 3-D shapes,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 14, no. 2, pp. 239-25 (February 1992). This conventional method iteratively conducts a two-step process: a point-to-point correspondence across the point clouds, and a minimization of a least square distance metric describing the misalignment. It assumes that the Euclidean nearest point pair finds the correct correspondence in most cases. To satisfy this assumption, a common method is to limit the magnitude of sensor changes, for instance to move the camera slowly during visual mapping so there is only a small distance of motion from one image to the next image for the same point. However, this limitation is often violated in real scenarios. For example, when a user is holding the camera while video recording a scene, the user often simply cannot hold the camera sufficiently still for such image processing due to human physical limitations. This may be especially true when a user is attempting to record moving objects and is intentionally moving the camera to attempt to match the motion of the objects, often moving the camera too quickly for accurate 3D processing.

[0026] The known improvements of classic ICP can be classified into two categories as follows. First, some solutions select a subset of feature pairs to remove incorrect matching. One technique uses a normal vector to obtain more accurate point pairs. See J. Serafin et al., “NICP: Dense normal based point cloud registration”, IROS, pp. 742-749 (2015). Also, a robust point matching (RPM) algorithm was proposed to use soft assignment (softassign) instead of binary assignment for better point selection and matching. See Gold et al., “New algorithms for 2D and 3D point matching: Pose estimation and correspondence”, Advances in Neural Information Processing Systems, pp. 957-964 (1995). Other recent methods related to feature pairs combined two strategies to remove matching errors: pre-filtering and overlap filtering. See S. Nobili et al., “Overlap-based ICP Tuning for Robust Localization of a humanoid Robot”, ICRA (2017). In each filtering technique, a set of error matches are removed. Some techniques attached descriptors to the points to eliminate ambiguities, such as with signature of histograms of orientations (SHOT). See Serafin et al.; Y. Zhuo et al., “Automatic registration of partial overlap three-dimensional surfaces”, Mechanic Automation and Control Engineering, proc. Of the International Conf. on. Pp. 299-302 (2010); Nobili et al.; and Tombari et al., “Unique signatures of histograms for local surface description”, European conference on computer vision, Springer, Berlin, Heidelberg, pp: 356-369 (2010). The main drawback of all of these methods is that they still are easily affected by noise, large viewpoint changes, and occlusions. In addition, many of the complex descriptors mentioned above are time consuming.

[0027] Second, the other solution to improve ICP is to use various error metrics to improve robustness. Some techniques used multiple error metrics such as point-to-plane and plane-to-plane matching See respectively, Grant et al., “Point-to-plane registration of terrestrial laser scans J. ISPRS Journal of Photogrammetry and Remote Sensing”, 72: 16-26 (2012); and A. Segal et al., “Generalized-ICP”, in Proceedings of Robotics: Science and Systems (RSS) (2009). Another example, Iteratively Reweighted Least Squares (IRLS)” was proposed using various cost functions to provide robustness to noise. See P. Bergstrom et al., “Robust registration of point sets using iteratively reweighted least squares”, Computational Optimization and Applications, 58(3): 543-561 (2014). Another technique introduced a trimmed ICP algorithm, which extended the application of ICP to point clouds with partial overlap, See D. Chetverikov et al., “Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm”, Image and Vision Computing, 23(3): 299-309 (2005). These methods are robust for a small portion of matching errors to some extent. However, these techniques still cannot handle very small overlapped images of a scene or large rotation from image to image of a scene, e.g., when the overlap ratio (between two images) is lower than 50%, or the rotation between two images is larger than 60.degree., these techniques still result in inadequate point cloud matching, in turn resulting in many inaccuracies and errors.

[0028] The disadvantages of those previous solutions can be grouped into three main aspects. First, those techniques that use low level feature representations are often very noisy in real applications. For example, a conventional point detector used in these techniques are often very sensitive to various factors, including lighting, viewpoints, occlusions, and so on, resulting in many erroneous point matches. Second, the conventional low level feature matching typically is performed in a dense matching space, and is conducted for every iteration in ICP. For instance, in the classic ICP (see Best et al., cited above), the points in the source and target are set to correspond in every iteration, which is the main reason why ICP is slow. This may include brute matching where every point in one point cloud is compared to every, or numerous, points in the other point cloud. This is very time consuming, such that the point cloud registration may cause noticeable delays in other processes, and may consume an unwarranted percentage of a processor’s processing capacity. This alone is detrimental to the functioning of the computer, computing device, or camera performing the cloud point registration. Third, an error metric of the point cloud registration may be based on local descriptors which also are sensitive to large rotations and small overlaps. Thus, although some improvements have been proposed, these improvements do not work well with very large rotation or very small overlap from image to image, e.g., where the overlap ratio is lower than about 50%. The point registration becomes mostly unusable when overlap between point clouds is smaller than 40% or the rotation angle is larger than 60.degree.. For the inadequate improvements, see Grant et al., Bergstrom et al., and Chetverikov et al. cited above.

[0029] To resolve these issues, an efficient and accurate method and system of point cloud registration for image processing may include an iteratively modifying the position or dimensions or both of three-dimensional overlap regions that enclose some part of point clouds forming a scene in two images. The iterations may be used to modify the position and dimensions of the overlap regions until a convergence is reached. During the iterations, the overlap regions also are adjusted to align to principal axes detected in the point clouds of the two images. The best overlap regions respectively of the two images so far may be used as the initial overlap regions in the next iteration. The transformation including a rotation part (such as a rotation matrix) and a translation part (such as a translation vector) then may be established by using the overlap regions formed during the convergence iteration. This avoids the heavier necessity of storing and modifying data such as weights for many points from iteration to iteration, thereby reducing the memory capacity, processing time, and power capacity functions of the computing device.

[0030] Now in more detail, given two point clouds P and Q, the task of registration is to estimate the rotation matrix R and the translation vector T from P to Q. To accomplish this, an overlap region B may be used, and here shaped as a cuboid, with B=(x, y, z, w, h, l) where the variables (x, y, z) is the center point coordinate of the overlap region, and (w, h, l) are the width, height, and length dimensions of the overlap region. For simplicity, the dimensions may be set at w=h=l. Overlap regions B.sup.P and B.sup.Q are denoted as the overlap in P and Q respectively. The joint estimation of the transformation and the overlap regions is formulated to maximize their joint posterior probability:

R ^ , T ^ , B ^ P , B ^ Q = arg max R , T , B P , B Q p ( R , T , B P , B Q | Z ) , ( 1 ) ##EQU00001##

where Z( ) is the observations and that Z=[Z(P),Z(Q)]. This may be rewritten as:

p(R,T,B.sup.P,B.sup.Q|Z)=p(B.sup.P,B.sup.Q|Z)p(R,T|B.sup.P,B.sup.Q), (2)

where the first probability on the right side of the equals represents the posterior probability of the overlap, while the second probability stands for the conditional probability of the transformation when their overlap is known. It is assumed that the prior probability distribution of the overlap p(B.sup.P, B.sup.Q), and the prior probability distribution of observation p(Z) are both uniform. According to the Bayes theory, the above equation (2) can be rewritten as:

p ( R , T , B P , B Q | Z ) = p ( Z | B P , B Q ) * p ( B P , B Q ) p ( Z ) .varies. p ( Z | B P , B Q ) p ( R , T | B P , B Q ) ( 3 ) ##EQU00002##

Therefore, the estimation of R, T, B.sup.P, B.sup.Q is divided into two sub problems: (1) the problem of overlap B.sup.P and B.sup.Q estimation by maximizing p(Z|B.sup.P, B.sup.Q), and (2) the problem of transformation R and T estimation given B.sup.P and B.sup.Q by maximizing p(R,T|B.sup.P, B.sup.Q). These two sub problems are interconnected and may be solved in an iterative process where the conditional probability of the observation may be represented by matching principal axes of the point clouds to the overlap regions, while the posterior probability of the transformation may be represented by the selection of the best overlap region (and its matched principal axes) for an iteration. Performing the disclosed method then effectively performs the probability equation (3) above. The details are explained below.

[0031] To accomplish these tasks, first, a Manhattan world assumption is used. The Manhattan world assumption is an assumption about the scene statistics of city and indoor scenes. See Flint et al., “A dynamic programming approach to reconstructing building interiors”, ECCV (2010). The assumption is that such city or indoor room scenes were built on a Cartesian grid in which surfaces in the world are aligned with one of three dominant directions where all planes are either parallel or perpendicular to each other. It is based on the observation that most indoor and city scenes are designed on a Manhattan three-dimensional grid. The disclosed method makes use of the geometric priors in Manhattan scenes in order to define a principal coordinate system (or axes) for the point clouds, and introduces a registration method based on the principal coordinate axes and registration of the three-dimensional overlap regions, here cuboids, that represent the two point clouds to be registered to each other.

[0032] Based on the representation, the overlap regions can be estimated with great precision, which accurately removes irrelevant points, leading to robustness of point cloud registration of multiple images with small overlaps of the same scene. To accomplish this, principal coordinate axes are aligned during a transformation estimation. This is performed very fast due to the ability to effectively use a small search space that does not rely on traditional low level feature matching, and which also leads to robustness towards large rotations from image to image.

[0033] The disclosed method and system achieves a significantly better result than the state-of-the-art in terms of both effectiveness and efficiency. When the overlap between point clouds is smaller than about 40% or the rotation angle is larger than about 60.degree., the disclosed method significantly outperforms the state-of-the-art methods. This is especially true when the overlap ratio is less than about 30%, or the rotation angle is larger than about 90.degree.. Such increase in accuracy and reduction in computational load, improves the functioning of the computer, computing device, or camera by reducing power and processor time consumed by performing these tasks leaving the computing device with more power and time to perform other tasks. Other details are provided below.

[0034] Referring now to FIG. 1A, each image 100 providing a point cloud 102 to be registered to each other may have the point cloud limited with a bounding box 106. The bounding box 106 may be unique for each point cloud and may be computed automatically by finding the minimum and maximum coordinate of the points on each default axis (x, y, z) for example. BY one form, the bounding box is the smallest possible bounding box to enclose the point set of a point cloud.

[0035] As mentioned, the main unit that is being manipulated by the methods here is an overlap region 104. The overlap region is placed within the point cloud 102 of the image 100, and in turn the bounding box, at a location likely to overlap two point clouds of two images being compared. The overlap region then may enclose at least some of the points within the point cloud. Thus, the bounding box is larger than the overlap region,

[0036] Referring to FIG. 1B, an example scene 110 shown in two images 112 and 114 is shown from two different perspectives whether captured from a single camera that is moved from image to image or captured from multiple cameras each providing one of the images 112 or 114. As shown, each image 112 and 114 has its own overlap region 118 and 116 respectively. By one form, the overlap region is cuboid, and more particularly cubic where the edge length is the same in all three directions (width, height, and length of the overlap region). It will be appreciated, however, that the methods herein could be applied to overlap regions with other shapes than cuboid or cubic. Also, and as described below, the methods herein modify the position and edge length of the overlap regions. Herein the position refers to the center point of the overlap region and the edge length is the edge forming a corner of the overlap region where two planes of the overlap region meet. It will be understood, however, that other alternatives could be used such as a corner or other desired location in or on the overlap region to be used as the position discussed herein, and the edge length (also referred to herein as edge size) may be referred to as edge dimension as well, but could be replaced with any overlap dimension that conveniently and accurately indicates a size of the overlap region for the purposes explained herein.

[0037] As mentioned, the point cloud registration may be achieved by a principal axes alignment within or on the overlap region. Such principal axes alignment is represented by a principal axes pair having two detected principal axes 120 and 122 that are detected by using the points in the point clouds. Herein, the principal axes refers to a Cartesian coordinate system with three axes (x, y, z), and the principal axes may be referred to herein as a principal coordinate system or principal axes system. Thus, the term “principal axes pair” or “candidate principal axes pair” discussed below refers to two matched principal axes each with three axes (x, y, z) and where the x-axis of one principal axes in the pair is matched to the x-axis of the other principal axes of the pair, and so on for each of the three axes. The principal axes 120 and 122 here are shown aligned to each other as well as respectively aligned to the overlap regions 116 and 118. The details of such alignment is provided below.

[0038] Referring to FIG. 2, a disclosed method 200 provides an overview for a method and system of point cloud registration for image processing. In the illustrated implementation, process 200 may include one or more operations, functions or actions 202 to 210 numbered evenly. By way of non-limiting example, process 200 may be described herein with reference to example image capture processing system 1200 of FIG. 12, and where relevant.

[0039] Process 200 may first include first obtaining a pair 202 of two captured images 212 and 214 each with a point cloud to be registered to each other. This may be an automatic process where point clouds are registered to each other every two frames of a video sequence or could be some other interval whether registered one or more times per video sequence as desired, and may be decided differently depending on the application. For instance, in a SLAM 3D model application, usually every two adjacent frames are registered, but for other 3D scene reconstruction applications, the images may be registered less frequently. By one form, the process is applied every other frame for 3D object reconstruction if the camera moves slowly. Also, the point clouds may be registered once every continuous scene or depending on some other content related indicator. Otherwise, the point clouds may be registered from two images once or other desired number of times upon a user selecting a certain mode on a computing device or camera, such as a 3D-related mode, or when setting a desired quality level or other computer function level (such as battery, processing speed, or memory capacity related settings). It also will be appreciated that the images may be single shot or still camera images as well.

[0040] Thereafter, a principal axes extraction operation 204 includes the extraction of three dominant orthogonal directions, and name them as the principal axes. This is performed by detecting planes 216 using the Manhattan world assumption so that all resulting planes are parallel or perpendicular. The normals 218 of the planes are used to pair perpendicular planes to construct a principal axes with each pairing of perpendicular planes.

Once an initial overlap between point clouds is determined as represented by initial overlap regions (such as overlap regions 116 and 118 in FIG. 1B), an iterative process may be used to determine refined or best overlap regions that can be used to establish a final transformation. The details for selecting the initial overlap regions is provided below. The iterative process may include an estimation of transformation 210 that is achieved by using aligned principal axes 222 and 224 of the two overlap regions representing the two point clouds being registered, and shown here formed during a principal axes alignment operation 208. The estimation of the transformation 210 represents the posterior probability described above. Also, an estimation of the overlap may be used when the overlap still is considered to be unknown or too coarse, and may be referred to as an overlap tuning operation 206, which represents the conditional probability of the observation described above, and is used to determine a best overlap region 226. The estimation of the overlap and the transformation is a coupled question, and an integrated framework is proposed to jointly optimize them.

[0041] Referring to FIG. 3, a process 300 is provided for a method and system of point cloud registration for image processing. In the illustrated implementation, process 300 may include one or more operations, functions or actions 302 to 308 numbered evenly. By way of non-limiting example, process 300 may be described herein with reference to example image capture processing system 1200 of FIG. 12, and where relevant.

[0042] Process 300 may include “obtaining image data of multiple images of the same scene from one or more cameras and including a first image and a second image each respectively being associated with a first point cloud and a second point cloud both separately forming the scene” 302. As mentioned above, such pairing of images with point clouds to be registered may be performed once per video sequence or continuous scene, every two images, or any other desired time interval.

[0043] Then, the process 300 also may include “generating three-dimensional first and second overlap regions each respectively enclosing points from, and having a position and dimensions in, the respective first and second point cloud” 304. This may include generating initial first and second overlap regions each with a center point, or other position indicator, located at a certain point relative to its point cloud. This also may include setting an edge length(s) of the two point clouds. By one form, the center points of both overlap regions starts at the center point of the bounding box of the point clouds, and the edge length is some factor of the size of the bounding box.

[0044] Separately, the process may include forming a set of candidate principal axes pairs comprising determining at least one principal axes formed by using the first point cloud to be paired with a principal axes formed by using the second point cloud and to generate candidate principal axes pairs. As mentioned, this may be accomplished by first using a Manhattan world assumption to detect planes used to form the principal axes. The principal axes may be formed from two normals respectively of detected perpendicular planes and computing the third axis by using the two normals to form three axes for each principal axes. The pairing of two principal axes may be performed by comparing color or luminance image data histograms along each axis of the principal axes being considered for pairing.

[0045] Process 300 then may include “iteratively modifying the position or dimensions or both of the first and second overlap regions to generate best first and second overlap regions so far of each iteration and until a convergence criteria is met” 306. Particularly, this operation may include generating a set of overlap region samples with variations in position or dimensions or both of that of the initial first and second overlap regions that are initial relative to an individual iteration. The overlap region samples are varied from the initial overlap regions by use of Gaussian distribution-based random sampling of variations of the position, dimensions, or both of the initial first or second overlap region.

[0046] A best principal axes pair of the set of candidate principal axes pairs then may be iteratively matched to individual overlap region samples, and this represents the posterior probability mentioned above. Determining the best match may comprise differencing an individual overlap region sample varied from the first or second overlap region, and differenced from the other of the initial first or second overlap region maintained with an initial center point position but where an edge length is varied to match the edge length of the overlap region sample being compared to. Before the differencing is performed, the overlap region sample may be adjusted to align with the principal axes of one of the candidate principal axes pair, and the other initial overlap region (with its center point fixed and edge length equal to the edge length of the overlap region sample) may be adjusted to align with the other principal axes of the one candidate principal axes pair. By one form, the differencing comprises determining a hamming difference between binary histograms point distributions within the adjusted overlap region sample and the initial overlap region. The candidate principal axes pair causing the smallest hamming distance associated with an overlap region sample, among all of the candidate principal axes pairs, may be considered a match to the overlap region sample. The hamming distance herein refers to a difference in bin values of the same bin location. Thus, if the overlap regions have bins 1 to 5, where one overlap region has values [1 0 0 1 0] and another overlap region has bin values [0 1 0 10}, then the absolute value result is (1-0)+(0-1)+(0-0)+(1-1)+(0-0)=2.

[0047] Also at each iteration, the best overlap region sample and matched best principal axes pair among all of the overlap region samples may depend on a difference in point distribution of (1) the overlap region sample of one of the first and second point clouds and (2) an overlap region of the other of the first and second point clouds, and this represents the conditional probability mentioned above. By one form, the hamming distances are already computed for the previous operation and may be compared to each other to determine which overlap region sample is associated with the minimum hamming distance, and is therefore, the best overlap region sample of the iteration. An updating test may be performed to determine if the best overlap region samples are better than the ones resulting from the prior iterations. If not, the iterations are continued by using the best overlap regions of the current iteration as the initial overlap regions of the next iteration. If at least one of the two best overlap region samples is better than the prior overlap region samples, then a convergence operation is performed. Convergence occurs when a difference in position or dimensions or both of the best overlap regions from one iteration to another meets a convergence criteria, such as a threshold. Once convergence is reached, iterations are stopped and the last best overlap region samples are used to establish a transformation.

[0048] The process 300 may include “establishing a transformation at least partially based on the best overlap regions, and the transformation comprising a rotation part, a translation part, or both to convert point locations between the first and second point clouds” 308. Thereby establishing a transformation that can be used to convert point locations of one of the point clouds into point locations on the other point cloud. This point cloud registration then can be used to generate a 3D map, 3D model, and forth as mentioned above.

[0049] Referring to FIG. 4, a process 400 is provided for generating a set of candidate principal axes pairs to be used to explain operations of the detailed process 500 described below. In the illustrated implementation, process 400 may include one or more operations, functions or actions 402 to 414 numbered evenly. By way of non-limiting example, process 400 may be described herein with reference to example image capture processing system 1200 of FIG. 12, and where relevant.

[0050] Process 400 may include plane 416 and normal extraction 402 from two point clouds to be registered to each other. Perpendicular planes are used to form 404 axes pairs 418, and principal coordinates 420 are extracted 406. These three operations are performed for each point cloud. The principal axes 420 from each point cloud are then matched or sorted 412 to form candidate principal axes pairs. By one form, all possible pairs are maintained as candidates for the iterations. By another form, only some number of the principal axes pairs with the least difference between the paired principal axes (the top K candidate principal axes) are used for the iterations.

[0051] Referring to FIG. 5, process 500 is provided for a method and system of point cloud registration for image processing. In the illustrated implementation, process 500 may include one or more operations, functions or actions 502 to 550 numbered evenly. By way of non-limiting example, process 500 may be described herein with reference to example image capture processing system 900 of FIG. 9, and where relevant.

[0052] Process 500 may include “determine if Manhattan space is to be processed” 501. By one form, the user may have the process initiated only when the scene to be captured and used for point cloud registration has a structure suitable for a Manhattan world assumption such as with an indoor room or city scene with many parallel and/or perpendicular planes. An indoor room with four walls, a floor, a ceiling, and many tables or counter tops is one example. This also may be performed automatically by a computing device such as detecting the computing device is indoors due to GPS, lighting, or other sensed characteristic of the environment to automatically initiate the process.

[0053] Process 500 may include “obtain image data of multiple images of the same scene” 502, and this may include obtaining images as frames of a video sequence or a series of still photographs. The scene being captured may be a stationary scene or may have moving objects, and the camera may be held still or may be moving as well. With the present process 500, better quality and better functionality of the computing device can be achieved despite a relatively low overlap and/or large difference in viewing angle from image to image as described below.

[0054] Process 500 may include “perform pre-processing” 504, and pre-processing at least sufficient to perform the point cloud registration when raw image data is received from one or more camera sensors for example. Such pre-processing may include noise reduction, pixel linearization, shading compensation, resolution reduction, Bayer demosaic, vignette elimination, and/or 3A related operations including automatic white balance (AWB), automatic focus (AF), and/or automatic exposure (AE), and so forth.

[0055] Process 500 may include “obtain image pair both with points of the same scene to be matched to each other” 506. The process 500 may be performed on two point clouds. Thus, the process 500 may be applied on just two RGBD images for example, but could also be applied on every two images in a video sequence as already described above. This can be an option provided to a user or switched on or off automatically also as described above.

Process 500 may include “generate a point cloud for each image of the pair within a bounding box” 508. Particularly, many different algorithms for point cloud generation exist, and are often based on some feature detection where edges or boundaries of objects are detected by algorithms that use gradients of color and/or luminance in the image data. Point cloud generation techniques used here for RGBD images may simply convert coordinates to 3D points for the point cloud and may follow (X, Y)=((x-cx)*Z/f, (y-cy)*Z/f), where (X, Y, Z) is the 3D coordinate in the camera coordinate system, (x, y) is the image coordinate, (cx, cy, f) are the camera intrinsic parameters. See for example, M Langer, “World vs. Camera Coordinates”, University of McGill, Center of Intelligent Machines, COMP 558, Lecture 4, www.cim.mcgill.ca/.about.langer/558/4-cameramodel.pdf (2010); R. Collins, “Lecture 12: Camera Projection Reading: T&V Section 2.4”, CSE 486 Computer Vision I, Introduction to Computer Vision, Penn. State, www.cse.psu.edu/.about.rtc12/CSE486/lecture12.pdf (2007). For an RGB image, depth is estimated first and then the point cloud. For instance, estimate of depth using stereo images is disclosed by K Zhang, et al., “Cross-Scale Cost Aggregation for Stereo Matching[C]//Computer Vision and Pattern Recognition”, IEEE, 2014:1590-1597; estimate of depth using video sequence is disclosed by K Karsch et al., “Depth Extraction from Video Using Non-parametric Sampling[C]//”, European Conference on Computer Vision, Springer, Berlin, Heidelberg, 2012:775-788, and estimate of depth using a single image is disclosed by D. Eigen et al., “Depth Map Prediction from a Single Image using a Multi-Scale Deep Network[J]”, Dept. of Computer Science, Courant Institute, New York University, 2014:2366-2374. The point clouds for registration are designated as P and Q for the remainder of the process 500.

[0056] While referring to process 400 (FIG. 4) as well, process 500 may include “extract candidate rotationally-aligned principal axes pairs” 510. This may first involve an operation to “detect planes in both point clouds” 512 (also operation 402 (FIG. 4)). For this operation, first all possible planes may be detected no matter the orientation of the plane. This may be performed using techniques such as Hough transform or random sample consensus (RANSAC) to name a few examples. Then, to “use Manhattan assumption planes” 514, only planes perpendicular to at least one other plane is maintained.

[0057] Process 500 may include “form perpendicular axis pairs on crossing planes” 516. Where it is detected that planes are perpendicular, the normals of those two planes are maintained together as the first two of three axes to form a principal axes. The normals here are unit normals of each plane (where length is designated as 1). This is repeated for each perpendicular pair of planes that are detected, as represented by operation 404 (FIG. 4), where axis v.sub.1 is paired with perpendicular axis v.sub.2 for each pair.

[0058] Process 500 may include “form third axis for each axis pair” 518, and specifically to form a principal axes (or principal coordinate system) S. The principal axes may comprise three orthogonal axes selected from the normal vectors of the dominant detected planes. The three principal axis may be denoted as v.sub.1, v.sub.2, v.sub.3 where v.sub.1.perp.v.sub.2, v.sub.1.perp.v.sub.3, v.sub.2.perp.v.sub.3 which comply with the right hand rule. Thus, the principal axes may be defined as S=[v.sub.1v.sub.2v.sub.3] where v.sub.1, v.sub.2, v.sub.3 are each 3.times.1 column vectors so that S is a 3.times.3 matrix. By one form, the three values that form each axis (or vector v) are values that indicate the direction of the normal that forms that axis, or in other words, each axis is a normal vector [a, b, c] from the equation of a plane (ax+by+cz+d=0), although other forms could be used.

[0059] For each pair of normal vectors (or axes), if they are perpendicular to each other, they are selected as candidates for v.sub.1 and v.sub.2 as mentioned above, while the third principal axis is obtained by:

v.sub.3=v.sub.1.times.v.sub.2 (4)

where “x” is the cross product of the two vectors. This operation is also shown on process 400 (see operation 406, FIG. 4). The generation of the principal axes for the two point clouds P and Q also may be referred to as S.sup.P={S.sub.i.sup.P, i=1 to C.sup.P} and S.sup.Q={S.sub.i.sup.Q, i=1 to C.sup.Q} respectively shown at operations 408 and 410 of process (FIG. 4).

[0060] Process 500 may include “select top K axes candidate pairs” 520. A point cloud may have many different principal axes. Thus, the process 500 may generate all or many individual candidate principal axes for P and Q. In order to determine the best overlap regions during the iterations, a principal axes form one point cloud is matched to a principal axes from the other point cloud to form a candidate principal axes pair. This is performed to eventually generate an aligned principal axes (or aligned coordinate system).

[0061] If (S.sup.P, S.sup.Q) is the aligned principal axes for P and Q, the corresponding points in P and Q will only have translation relations in the created coordinate system. In other words, each principal axes establishes its own coordinate system, such as (x, y, z) coordinates in 3D space. Therefore, when the coordinates of points in P and Q are converted to new coordinate systems as indicated by aligned principal axes pair (S.sup.P, S.sup.Q) and where P’=S.sup.P.times.P and Q’=S.sup.Q.times.Q (where “x” here is matrix multiplication), the corresponding points in P’ and Q’ only have a difference in translation between corresponding points but no difference in rotation, thereby simplifying the determination of the transformation.

[0062] Thus, process 500 may include “perform axes matching” 522, and the principal axes of one point cloud may be matched to another principal axes of the other point cloud to form the set of candidate principal axes pairs to be used during the iteration and to determine the aligned principal axes pair and best overlap regions. The matching of principal axes may be performed by using a set of heuristic rules. By one form, process 500 optionally may include “use color histograms” 524. By this option, the image data of the axes are projected into color histograms on each axis, and the principal axes being compared and that have the least difference are designated as candidate principal axes pairs. Other methods rather than the color histograms may be used for matching the principal axes as well. The principal axes matching is shown as operation 412 of process 400 (FIG. 4).

[0063] In this example, there may be a threshold cutoff number of candidate principal axes pairs that will be used for the analysis. The selected candidate set may be denoted as S.sup.PQ={(S.sub.1.sup.P,S.sub.1.sup.Q),(S.sub.2.sup.P,S.sub.2.sup.Q), … (S.sub.K.sup.P, S.sub.K.sup.Q)} for k=1 to K candidate principal axes pairs, and this is also shown as operation 414 of process 400 (FIG. 4). By one example as shown on FIG. 5H, since there are six faces (such as 502-H and 504-H) in a cuboid overlap region 500-H, there are at most 24 candidates, and therefore 24 possible different orientations, for the principal axes in a cuboid including one principal axes 506-H at each corner 508-H, where one axis (such as the x-axis) may extend in any of three different directions (thus, 8 corners.times.3 directions for the axes=24). Therefore, the maximum value of K here is 24. Note it does not matter where the candidate axes is located on or in the point cloud or overlap region; the placement of the principal axes at a corner of the cuboid overlap region still may represent the same principal axes extending in the same directions.

[0064] By another approach, the process 500 may use all or many individual candidate principal axes pairs rather than limiting the set to only the top K, although such a process may unnecessarily increase the computational load of the point cloud registration.

[0065] The process 500 then turns to the iterations. Each iteration has two stages. Generally, the first stage is to estimate (or re-estimate) the overlap region B.sup.P. In the first stage, overlap regions B.sup.Q and B.sup.P are initialized (and the initial overlap regions of each iteration are denoted below as {circumflex over (B)}.sup.P and {circumflex over (B)}.sup.Q). B.sup.Q is fixed in the point cloud Q such that its center point position does not change during the first stage. N overlap region samples are generated to represent the candidate overlap regions B.sup.P in point cloud P are then generated as variations in position or edge length or both of the initial overlap region B.sup.P. The edge length of B.sup.Q is changed to match the edge length of each overlap region sample it is to be compared to throughout the first stage and still in a single iteration. The best overlap region B.sup.P is then determined during the first stage by finding the maximum observation probability as described below.

[0066] In the second stage, the best overlap region B.sup.P from point cloud P and determined in the first stage now has its center point fixed in point cloud P, and then N overlap region samples are drawn to represent the candidate overlap region B.sup.Q in point cloud Q. The overlap region B.sup.P and The iteration proceeds to select the best overlap region B.sup.Q by finding the maximum observation probability as with the first stage. These two stages include both estimation of transformation and overlap tuning processes that are interactively conducted until the overlap region does not change (convergence is reached). When no convergence is reached, the best overlap regions B.sup.P and B.sup.Q in both stages of the iteration are used as the initial overlap regions B.sup.P and B.sup.Q in the next iteration.

[0067] Now in more detail, process 500 may include “obtain initial selected first and second overlap regions” 526. This preliminary operation refers to obtaining any parameters set manually by a user or automatically by the point cloud registration or 3D unit, or any other application. By one example, this may include setting the desired edge length, as some fraction of the length of the bounding box of a point cloud for example. If the user, or a program automatically, wants to use the method in a video sequence, different decimal numbers may be input such as “0.2”, “0.3”, and “0.5” as fractions of the bounding box length to be set as the edge length of the overlap regions. The best quality or efficient result then may be selected. Otherwise, the input may be a setting such as the decimal number “0.5”, which does not ensure the best result, but still provides acceptable results. Many other parameters may be set this way as well.

[0068] Process 500 may include “initialize first and second overlapping regions B.sup.P, B.sup.Q” 528. This may include setting both the position (here, the center point) and the edge length (and in turn, the edge location) of the overlap regions. Thus, process 500 may include “determine center points of first and second 3D overlapping regions” 530. By one form, the overlap regions {circumflex over (B)}.sup.P and {circumflex over (B)}.sup.Q are initialized to have their center points located at the center of point clouds P and Q respectively, which refers to the center coordinates of the three-dimensional bounding box of P and Q.

[0069] Process 500 also may include “determine initial length and location of edges of the first and second overlapping regions” 532. When the overlap region is cubic as in the example provided here, only one number is needed as the input for all three dimensions (width, height and length of the overlap region), and it will be understood that the center point of the three-dimensional cubic overlap region has coordinates one-half edge length from the edge of the overlap region and to each edge in any elevational view of the overlap region. As mentioned, the initial length of the edge of the overlap region may be set by a user or automatically set by a 3D unit or some other application. By one example, the bounding box for the point clouds P and Q are obtained. Then, the longest edge (denoted as a) of the two bounding boxes is determined. The initial edge length of the overlap region then may be set at a/2, so an example input may be a decimal number “0.5”. The user or application may set the initial edge length according to experimentation or prior knowledge of the point clouds. If no prior information is known, the initial decimal number may be set to some default such as “0.5”. A good initial edge length can speed up the process and also helps to improve the registration accuracy, but the number does not significantly affect the final result due to the updating in each iteration that occurs no matter the input edge length.

[0070] Process 500 may include an iteration counter that sets the time or iteration as “set t=0” 534 to start on a first iteration. Then, process 500 may include “re-estimate first overlapping region {circumflex over (B)}.sup.P” 536 to begin the first stage of the iteration. The first stage is considered to continue until operation 589 in the present example. In the first stage, overlap region samples are generated with modified positions, dimensions or both to be used to determine a best overlap region by performing a transformation estimation and overlap tuning as introduced above and explained in greater detail below. The generation of the overlap regions samples may be explained as follows.

[0071] The process 500 may include “generate N overlap region samples for n=0 to N samples as variations of {circumflex over (B)}.sup.P” 538. By one form, N=50 for n=0 to 49 samples, but N may be set to other numbers as desired. This operation of process 500 may include “generate individual samples B.sub.t.sup.Pn with variations of edge dimensions, center point position, or both relative to {circumflex over (B)}.sup.P” 540. By the present example, each overlap region sample has both a change in position and a change in edge length. It will be understood, however, that other examples are possible such as changing the edge length for 25 samples (or some other number) and changing the center position for 25 samples (or some other number. Another option may include having samples generated with edge length changes on some iterations, such as odd iterations t=1, 3, 5, … , for example, while samples are generated with center position changes at even iteration t=2, 4, 6, … , for example. Many other variations are contemplated.

[0072] Process 500 may include “generate the samples B.sub.t.sup.Pn using a Gaussian distributed-based random sampling algorithm” 542. Thus, the overlap region samples may be generated completely randomly or at some interval as a ratio of the edge length and center position of the bounding box as described above. Instead, a simplified form of sequential importance resampling (SIR) particle filtering is applied, and this involves using a more efficient technique by using a Gaussian function as a likelihood for sampling. The other details of the particle filtering are provided below with explanation of the overlap tuning. Relevant here, a Gaussian function can be applied to obtain the overlap region samples n.

[0073] For instance, the position and the size (edge length) are sampled independently. Here, let edge length (or edge size) of an overlap region sample be “s” and the edge length of the initial overlap region {circumflex over (B)}.sup.P is as an example, so that to generate one sample is to select an “s” from all possible “s”. The Gaussian function or distribution for edge length can then be stated as:

N ( s , , .sigma. s 2 ) = 1 2 .pi. .sigma. s exp ( - ( s - ) 2 2 .sigma. s 2 ) ( 5 ) ##EQU00003##

where .sigma. is the deviation (and in turn .sigma..sup.2 is the variance) that may be predetermined by experimentation and provided by the 3D unit or other application. In the present implementation, the variance for position may be set to 0.25*edge length, and the variance for edge length may be set to 0.2*edge length although other values may be used.

[0074] The Gaussian function indicates which “s” is more likely to be selected by providing a higher value of “s” when “s” is near to and a lower value of “s” when “s” is farther from , thereby establishing a Gaussian function or distribution whose mean is the edge length of the best overlap region selected by the previous iteration. The mean of the best overlap region {circumflex over (B)}.sup.P is favored because we assume that the and obtained at the t-th iteration has higher or equal conditional observation probability than that in the (t-1)-th iteration. Many different techniques may be used to implement the Gaussian sampling. See for example, www.alanzucconi.com/2015/09/16/how-to-sample-from-a-gaussian-distribution- .

[0075] For the change in position (or center point), the position may be represented as a vector c=(x, y, z) which is the center point of the cuboid,* its Gaussian function is*

N ( c , , ) = 1 ( 2 .pi. ) 3 2 1 2 exp ( - 1 2 ( c - ) T - 1 ( c - ) ) ( 6 ) ##EQU00004##

where is the center point of {circumflex over (B)}.sup.P, and E is the covariance (identity) matrix:

= [ 1 0 0 0 1 0 0 0 1 ] ( 7 ) ##EQU00005##

[0076] Still in the first stage, process 500 may include “generate individual samples {circumflex over (B)}.sup.Qn with variations of edge dimensions with the same edge dimensions as a corresponding B.sub.t.sup.Pn, and a center point position fixed as the same as {circumflex over (B)}.sup.Q” 544. Fixed here refers to maintaining the center point of each overlap region {circumflex over (B)}.sup.Qn to be the same as the initial (or latest) {circumflex over (B)}.sup.Q. Thus, in the first stage, when overlap region samples B.sub.t.sup.Pn are generated, each sample changes in both its center position and edge dimensions. Then, when the n overlap region samples {circumflex over (B)}.sup.Qn are generated in the first stage, the {circumflex over (B)}.sup.Qn are generated each with a same edge size as B.sub.t.sup.Pn but with a center position that is the same as {circumflex over (B)}.sup.Q so that the two overlap regions used as input in equation (X) below in order to align (or match) the principal axes with the overlap regions have the same size. In other words, comparing overlap regions with different sizes, and in turn different sizes in point distributions, is more likely to be less accurate.

[0077] In the second stage, operation 540 includes (in the 2.sup.nd stage, only modify center points to generate individual samples B.sub.t.sup.Qn). Specifically, similar operations may be performed except in the second stage where overlap region samples for B.sub.t.sup.Qn are generated using the Gaussian distribution. Thus, in the second stage, the center points are modified using equation (6) above except for switching each P with Q. However, the edge length of all of the overlap region samples B.sub.t.sup.Qn are fixed to be the same as the edge length of the latest best {circumflex over (B)}.sup.P, the center point of B.sub.t.sup.Qn is otherwise independent of the center point of {circumflex over (B)}.sup.P. The overlap region samples B.sub.t.sup.Qn will be compared to the latest best {circumflex over (B)}.sup.P with its edge length and center point unchanged (544).

[0078] Returning to the first stage description, the transformation estimation may be performed by using principal axes or principal coordinate alignment. The overlap regions may be set to have cuboid faces that are each parallel to one of the principal axes of the point cloud. Once the center points and edge lengths of the overlap regions B.sup.P and B.sup.Q are set in the point clouds P and Q respectively, the problem of transformation estimation is to select the best transformation rotation part R and translation part T that maximize the probability p(R,T|B.sup.P,B.sup.Q), and which represents the conditional probability of the transformation given their overlap regions as explained above.

[0079] Given points in an overlap region B, its observation Z(B) may be described by a binary point distribution histogram Z(B)=H.sub.p(B). The overlap region B is divided into M cubic bins where each bin size may be set to about 0.05.sup.3 m.sup.3 for example in this disclosure, the function H.sub.p( ) is defined as a binary histogram where the value in each bin H.sub.p[m] (m=1, 2, … M) is 1 if the number of points in this bin is larger than a threshold and 0 otherwise. The cost of matching two overlap regions B.sup.P and B.sup.Q then may be defined as:

cost(Z(B.sup.P),Z(B.sup.Q))=d(H.sub.p(B.sup.P),H.sub.p(B.sup.Q)) (8)

where function d( ) is the Hamming distance. In order to select the aligned coordinate system pair for P and Q, we define a likelihood that measures the coordinate system alignment probability when given the overlap regions:

p ( S i P , S i Q | B P , B Q ) = exp - cost ( Z ( S k P .times. B P ) , Z ( S k Q .times. B Q ) ) = exp - d ( H p ( S k P .times. B P ) , H Q ( S k Q .times. B Q ) ) ( 9 ) ##EQU00006##

and the aligned coordinate system pair may be selected by:

( S ^ P , S ^ Q ) = argmax S k P , S k Q .di-elect cons. S PQ p ( S k P , S k Q | B P , B Q ) = exp - d ( H p ( S k P .times. B P ) , H Q ( S k Q .times. B Q ) ) ( 10 ) ##EQU00007##

where a principal axes S.sub.k.sup.P,S.sub.k.sup.Q belongs to one of the k=1 to K candidate principal axes pairs S.sup.PQ. As mentioned, K may be limited to the 24 possible orientations on a cuboid. Due to such a sparse space for principal axes alignment rather than testing every possible principal axes pair in the point clouds, the method can be very fast.

[0080] To implement this principal axis alignment equation (or matching of principal axes to overlap region sample), process 500 may include starting a sample counter to “set n=0” 545 to initiate the process. Thereafter, process 500 may include “set B.sup.P=B.sub.t.sup.Pn” 546, “set B.sup.Q={circumflex over (B)}.sup.Qn” 548, and “set k=0” 550.

[0081] Process 500 then may include “determine adjusted overlap regions for B.sub.t.sup.Pn and {circumflex over (B)}.sup.Qn aligned to principal axes k as (S.sub.k.sup.P.times.B.sub.t.sup.Pn) and (S.sub.k.sup.Q.times.{circumflex over (B)}.sup.Qn)” 552. This refers to matrix multiplication, and adjusts the coordinates of the points in the overlap region samples B.sub.t.sup.Pn and overlap region {circumflex over (B)}.sup.Qn respectively to the coordinates of the principal axes S.sub.k.sup.P and S.sub.k.sup.Q to thereby align the overlap region samples and the overlap region to the respective principal axes.

[0082] Process 500 may include “determine difference between adjusted overlap regions” 554, and by one form this includes an operation to “determine hamming distance d between binary point distribution histograms of each adjusted overlap region” 556, and according to equation (7) above. It will be understood that other equations may be used for determining hamming distance or a different representative difference between the overlap regions aligned to candidate principal axes.

[0083] PQ Process 500 may include “save S.sub.k.sup.PQ as best for n if hamming distance exp.sup.l/d is a maximum so far” 558, where the candidate principal axes pair is saved as the potential aligned principal axes pair that has the maximum value for alignment equation (7). The first time the alignment equation (7) is used for an iteration, the first value is automatically saved.

[0084] Process 500 proceeds with an inquiry “k=K?” 560 to determine if all candidate principal axes have been tested yet, and if not the process 500 may include “set k=k+1” 562, and the process loops back to operation 552 to test the next principal axes pair S.sub.k.sup.PQ for alignment and matching to an overlap region sample B.sub.t.sup.Pn Once all of the candidate principal axes pairs K are tested, process 500 may include “set last saved S.sub.k.sup.PQ as best aligned principal axes pair (S.sub.t.sup.Pn,S.sub.t.sup.Qn) for n” 564 such that a principal axes pair is now aligned (or matched to) an overlap region sample B.sub.t.sup.Pn.

[0085] Process 500 then may include an inquiry “n=N?” 566 to determine if all overlap region samples N have been tested already. If not, process 500 may include “set n=n+1” 567 to obtain the next overlap region sample B.sub.t.sup.Pn and its corresponding overlap region {circumflex over (B)}.sup.Qn. This next overlap region sample is then tested with all of the candidate principal axes S.sub.k.sup.PQ until all (or some desired amount less than all) of the overlap region samples B.sub.t.sup.Pn are matched or aligned to a principal axes pair.

[0086] Once the principal axes alignment is complete for each overlap region sample, process 500 proceeds with an overlap tuning operation to determine the single best overlap region for the iteration, and which represents the posterior probability as explained above. The basis for the overlap tuning is as follows.

[0087] A recursion of Bayesian filtering may be employed to estimate the overlapped region. In the Bayesian predict stage, the overlap regions in P and Q are predicted, and in the update stage, the overlap regions are updated based on the observation:

Predict : p ( B t + 1 P , B t + 1 Q | Z t ) = .intg. .intg. p ( B t + 1 P , B t + 1 Q | B t P , B t Q ) p ( B t P , B t Q | Z t ) dB t P dB t Q ( 11 ) Update : p ( B t + 1 P , B t + 1 Q | Z t + 1 ) = p ( Z t + 1 | B t + 1 P , B t + 1 Q ) p ( B t + 1 P , B t + 1 Q | Z t ) p ( Z t + 1 | Z t ) .varies. p ( Z t + 1 | B t + 1 P , B t + 1 Q ) p ( B t + 1 P , B t + 1 Q | Z t ) ( 12 ) ##EQU00008##

where t is the iteration time. The above predict and update steps are iteratively conducted until convergence. This iterative process is theoretically appealing, but it is intractable due to the large search space of the combination of B.sup.P and B.sup.Q.

[0088] Alternatively, it has been found that an accurate approximation of the above Bayes filtering may be achieved with particle filtering. Particularly, the disclosed method may be a simplified version of a sequential importance resampling (SIR) particle filter. The particle filter represents the likelihood P(X) by random sampling. It computes the weight of each sample and selects those samples with the largest weights. See for example, P. Abbeel, “Particle Filters”, UC Berkeley EEC S (November 2012) (at people.eecs.berkeley.edu/pabbeel/cs287-fa12/slides/ParticleFilters.pdf). The particle filter may be applied to the present point cloud registration method as follows.

[0089] The first part of the particle filtering is prediction, and here the disclosed method predicts the position and size of the overlap regions based on a Gaussian distribution:

p(B.sub.t+1.sup.Pn|{circumflex over (B)}.sup.P).about.N({circumflex over (B)}.sup.P,.sigma..sup.2) (13)

which is applied to the present cloud registration method when determining the overlap region samples N as already described above, and that is described above with operation 542.

[0090] Next, weighting is performed to compute the importance and/or weights of each sample. This is accomplished by the hamming distance equation (24) and (25) formed below and to select the best overlap region samples of the iteration.

[0091] Lastly, the SIR particle filtering includes resampling based on the computed weights to update the samples at a next iteration. The difference between the present method and the basic SIR particle filter is at this stage, the samples are resampled based on the distribution of weights in the original SIR method. In contrast, the present point registration method merely reuses the best overlap region sample with the maximum weight. Thus, the present method’s resampling generates only one sample, or it can be stated that the method generates N replicate samples by copying the best sample from the previous iteration.

[0092] To these ends, the probability density function p(B.sub.t.sup.P, B.sub.t.sup.Q|Z.sub.t) may be represented by a set of particles {.omega..sub.t.sup.k, B.sub.t.sup.Pk, B.sub.t.sup.Qk}.sub.k=1.sup.N. However, the joint sampling space of B.sub.t.sup.P and B.sub.t.sup.Q is still too large, and a common solution to reduce the sampling space is to sample B.sub.t.sup.P and B.sub.t.sup.Q one by one. The formulation can be approximated by:

[0093] Predict B.sub.t+1.sup.P:

p(B.sub.t+1.sup.P|Z.sub.t)=.intg.p(B.sub.t+1.sup.P|B.sub.t.sup.P)p(B.sub- .t.sup.P|Z.sub.t)dB.sub.t.sup.P (14)

[0094] Update B.sub.t+1.sup.P:

p(B.sub.t+1.sup.P|Z.sub.t+1).varies.p(Z.sub.t+1|B.sub.t+1.sup.P,B.sub.t.- sup.Qbest)p(B.sub.t+1.sup.P|Z.sub.t) (15)

[0095] Predict B.sub.t+1.sup.Q:

p(B.sub.t+1.sup.P|Z.sub.t)=.intg.p(B.sub.t+1.sup.Q|B.sub.t.sup.Q)p(B.sub- .t.sup.Q|Z.sub.t)dB.sub.t.sup.Q (16)

[0096] Update B.sub.t+1.sup.Q:

p(B.sub.t+1.sup.Q|Z.sub.t+1).varies.p(Z.sub.t+1|B.sub.t+1.sup.Pbest,B.su- b.t+1.sup.Q)p(B.sub.t+1.sup.Q|Z.sub.t) (17)

where p(B.sub.t+1.sup.P|B.sub.t.sup.P) and p(B.sub.t+1.sup.Q|B.sub.t.sup.Q) are the proposal prediction function from a current state to the next state which are modeled by Gaussian distributions centered at B.sub.t.sup.P and B.sub.t.sup.Q and already applied by the overlap region sampling N, and P(B.sub.t.sup.P|Z.sub.t) and P(B.sub.t.sup.Q|Z.sub.t) are the probability density function which will be modeled by the particle filter.

……
……
……