CN1975323A - Method for making three-dimensional measurement of objects utilizing single digital camera to freely shoot - Google Patents

Abstract

This invention disclosed a three-dimension measurement method which consists of preparation, picture intussusception, coding spot identification, camera location confirmation, target curve pick-up, homonymic curve matching and target curve reconstruction. It has the following characteristics: mark on the target curve before measuring to make the picture identification much more easy; place a guage and a set of coding spots around the target object; screen randomly to obtain a set of image by the camera; calculate automatically the position and gesture of the camera according to the images; pick up the marked curve and match the homonymic curve in different images so that the information of three-dimensional spot line was automatically calculated.

Description

Method for carrying out three-dimensional measurement on object by using single digital camera to freely shoot

One, the technical field

Three-dimensional object measurement belongs to measurement, test technical field. The corresponding code in the international classification table is G01B.

Second, background Art

The rapid development of the three-dimensional measurement technology is greatly promoted by increasingly wide application requirements of the three-dimensional measurement technology in the fields of reverse engineering, industrial detection, quality control and the like, and various measurement methods based on optical, acoustic, electromagnetic and mechanical contact principles, such as a three-coordinate measuring machine, a laser scanner, a structured light measuring instrument and the like, appear. The three-coordinate measuring machine adopts mechanical contact type sensing, has high measuring precision, generally needs a special measuring chamber and a special measuring table, has high requirement on measuring environment, has limited measuring range and low measuring efficiency, and is not suitable for measuring soft objects; laser line scanning measurement and structured light irradiation measurement are the mainstream methods for measuring the three-dimensional geometric shape at present, dense point cloud data on the surface of a model can be quickly obtained by irradiating laser or structured light on the surface of an object, but the method is limited by a scanning range, high light reflection on the surface of the object and the like, and a laser scanner and a structured light measuring instrument are expensive. Especially for mechanical products, because of the obvious structural features, it is often necessary to mainly acquire data of angular points, edges, and some control lines on the surface of the model, which play a key role in reconstructing a digital model of a measured object, and the structured light irradiation measurement method, the laser line scanning measurement method, and the like acquire integral point cloud or grid data of the surface, so that on one hand, the output data amount is huge, and on the other hand, the required data of the edge features and some control lines on the key sections cannot be directly and explicitly acquired. Moreover, the measurement data output by these methods are generally more effective in smooth flat areas of the model, while the measurement results are less effective at the critical corners and edges.

In order to realize the measurement of three-dimensional geometric information in a simpler hardware condition and in a more flexible and convenient manner, it has become a research focus in recent years to accurately reconstruct the position and shape of an object from a plurality of images taken by a digital camera. Among them, TriTop from company Gom Germany^®The system can carry out high-precision three-dimensional coordinate positioning by using a single digital camera free shooting mode. The system is characterized in that a group of coding points and length scales are placed in a scene, mark points which are easy to identify are pasted on an interested part, then a user holds a digital camera to freely shoot a plurality of images, a certain overlap of the images is required, and after all the images are input into a corresponding software system, the system automatically calculates the position and the posture of the camera and the space coordinates of all the mark points at one time when shooting. The system is commercialized and sold in China, but the system can only carry out space coordinate positioning of a specific marked point target (consisting of a black circular ring surrounding a white circular point) at present, is generally used for matching with other measuring methods to carry out splicing of multi-view angle measuring data, and cannot carry out three-dimensional measurement on a curve targetAnd the method cannot be used for three-dimensional digital model reconstruction of products with complex geometric shapes.

Third, the invention

The invention aims to realize a practical measuring method which is oriented to industrial product measurement modeling, is convenient to implement, has higher precision and lower cost by using simple hardware conditions. Therefore, the characteristic lines of the object to be measured and certain key section control lines on the surface of the object required for digital model reconstruction are marked, so that the characteristic lines are obviously different from the color of the object to be measured in color brightness, and the image recognition is facilitated; a ruler and a group of specially designed coding points are arranged around a measured object; then, holding a digital camera by hand to obtain a group of images of the measured object in a free shooting mode; according to the group of images, the position and the posture of the camera during each shooting are automatically and accurately calculated, meanwhile, a convenient interaction means is provided for a user, the semi-automatic extraction of the identification curve and the optimized matching of homonymous curves in different images are realized, and further, the three-dimensional point list information of the identified curve structure is automatically calculated.

According to the scheme, a practical system which is flexible and convenient to use and suitable for measuring three-dimensional curve structures on objects with different sizes is developed, and the practical system can be well used for building a three-dimensional digital model of a mechanical product based on a real object. The object three-dimensional measurement method is characterized in that only one digital camera and one common personal computer are used for measurement, a group of coding points and a scale are used as auxiliary materials, a complex measurement hardware system is not needed, and complicated calibration of the measurement system is not needed; curve data required for digital model reconstruction of an object are directly and explicitly generated, data redundancy is avoided, and high-efficiency model reconstruction is facilitated; all the measured data are automatically located under a world coordinate system, so that the difficult problem that the data measured for many times in other measuring methods need to be spliced is avoided, and the accumulated error caused by splicing a plurality of data sets is also avoided. The method comprises the main steps of measurement preparation, image shooting, camera pose determination, target curve extraction, automatic matching optimization of homonymous curves, three-dimensional reconstruction of target curves and the like.

Measurement preparation and image capture

Measurement preparation mainly involves three aspects of work: 1) marking a target curve (generally a boundary line of a natural curved surface sheet, a key section control line and the like) according to the requirement of digital model reconstruction, so that the target curve is obviously different from a measured object in color and brightness, and the image identification is facilitated; 2) several code points are arranged in the measurement area. Each coding point has a unique identity code, and quick and reliable identification in different images is easy to realize. Both the marking points and the coding points can be generated simply by generating the corresponding patterns on a computer and then printing. The coding points are attached to the upper surfaces of the hard board, the wood chip and the like, and can be repeatedly used; 3) a scale with two code points is placed in the measuring scene, and the distance between the centers of the two code points is known. The purpose of placing the scale is to obtain the actual size of the object to be measured, otherwise, only three-dimensional structures differing by a proportionality coefficient can be obtained. After the preparation work is finished, the digital camera can be held by hand to shoot the measured object at multiple angles. It is required that there is some overlap between the images, i.e. that one image has at least more than 5 coding points and some object curves in common with the other image.

Identification and localization of code points

The pattern of code dots is designed to consist of a central white dot, a middle black ring, and an outer ring, wherein the outer ring is divided into 15 equal parts, each equal part is black or white, black represents a binary code "0", and white represents a binary code "1", which is called a "code strip". The code for each point in a set of code points is different. According to the coding, the identity of different images can be reliably identified, and the corresponding relation of the same-name coding points among the images is automatically established. It is based on the positions of these homonymous code points in the image and the correspondence between multiple images that the automatic calculation of the respective positions and attitudes of the camera at the time of free shooting can be realized.

Because the central circle of the coding point is elliptical after being imaged by the CCD, in order to identify and position the coding point in the image, the invention firstly adopts the Canny operator to segment the image, extracts contour information representing different areas in the image, and then gradually filters candidate coding point targets according to 5 constraint conditions of the size, the shape, the ellipse fitting residual error, the mean value of the area gray scale, the variance of the area gray scale and the like of each contour, thereby realizing the extraction of the coding point target.

After locking a code point target, it is decoded, i.e. it is determined which code point it is. The decoding is based on the gray scale of each sector on the encoded band. The invention adopts a method of fitting the middle ellipse of the coding band and performing median filtering on each pixel on the middle ellipse, comprehensively considers the gray values of most pixels in the coding band, and can eliminate the influence of isolated noise. The processing result of a large number of field shot images shows that the method is very effective for improving the robustness of the identification of the coding point.

And finally, determining the center coordinates of the coding points with sub-pixel positioning accuracy according to the gray values of all pixel points in the identified central circular area of the coding points.

Camera position and pose determination

The fundamental matrix between the two images is first calculated based on the pixel coordinates of at least 5 homonymous code point centers in the two images. The camera poses corresponding to the two images and the three-dimensional coordinates of the centers of the coding points which are visible together in the two images can be further recovered by the camera internal parameters and the basic matrix. And then, according to the homonymy corresponding relation between the obtained three-dimensional space point and the coding point on the third image, solving the camera attitude corresponding to the third image, further obtaining more central three-dimensional coordinates of the coding point, then solving the camera attitude corresponding to the next image, and increasing progressively in this way until all the camera attitudes and the central three-dimensional coordinates of the coding point are obtained. And finally, performing integral optimization on all camera parameters and the three-dimensional coordinates of the center of the coding point by adopting a light beam adjustment method so as to further improve the precision. The strategy of combining the incremental method and the global optimization method not only has higher efficiency, but also can ensure that the camera positioning reaches higher precision.

And reading a group of shot images at one time, and automatically calculating and recording the position and the posture of the camera during each shooting by the measuring system. The camera position and the camera pose during each shooting are determined, which means that the position and the pose of each image under a unified world coordinate system are determined, so that each target curve reconstructed by a subsequent algorithm according to different image pairs is directly positioned in the same coordinate system without data splicing.

Semi-automatic extraction of target curves

For a currently selected image pair displayed in the two image windows, points are respectively taken on the two images by a mouse near the same marking curve which is visible together, the connecting line of the points approximately reflects the corresponding image curve contour, and then the measuring system automatically and optimally fits the approximate image curve contour to the image curve according to the principle of energy optimization. In the half automatic extraction process, a user only needs to sequentially pick points near the image curve, the method is simple and easy to implement, the stability of target curve extraction and subsequent automatic matching optimization is greatly improved due to the good initial search position, and meanwhile, the optimal fitting method ensures the extraction precision of the target curve.

Automatic matching optimization of homonymous curves

After a pair of homonymous target curves is extracted from one image pair, a key problem is to establish the corresponding relation of each pixel point on the homonymous target curves between two images. Based on the correspondence of the image points, the spatial coordinates of these points can be reconstructed by stereo triangulation.

According to the basic theory of stereo vision, the epipolar constraint should be satisfied by the homonymous points between a pair of images obtained by the camera at different positions and angles. For two candidate homonymous matching points v in an image pair₁And v₂The invention is used on the second image v₂To v₁Of (2)Distance of lines, and v on the first image₁To v₂Is measured by the sum of the distances of the epipolar lines of₁And v₂The degree of matching. Meanwhile, the matching of the point columns on the two homonymous target curves meets the spatial correlation constraint, namely the point columns arranged in sequence on one curve are certain to correspond to the point columns arranged in sequence on the other curve. Based on such analysis, the present invention first employs dynamic programming to obtain initial matches of discrete pixel points on the homonymous curve. And after initial matching of the point pairs on the curves of the images with the same name is obtained, further performing matching optimization on the curves. Setting two homonymous image curves to use parameter equation c respectively₁(l)、c₂(l) Showing that the present invention optimizes the following objective function to achieve c₁、c₂Exact matching of points on a curve

<math> <mrow> <mi>min</mi> <munderover> <mo>&Integral;</mo> <mn>0</mn> <msub> <mi>L</mi> <mn>1</mn> </msub> </munderover> <mfrac> <mrow> <mo>|</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>F</mi> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <mi>eF</mi> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>F</mi> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>F</mi> <msup> <mi>e</mi> <mi>T</mi> </msup> </mrow> </mfrac> <mi>dl</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow> </math>

Wherein, $e=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right];$ f is a fundamental matrix between the pairs of images, which is known when the camera position and attitude at each shot have been determined; σ (l) is the mapping function to be solved, representing c₁The point on the curve with the parameter l is on the curve c₂The parameter values of (1).

Three-dimensional reconstruction of target curves

After matching of all pixel points on the homonymous curve is completed, the spatial coordinates of the points can be reconstructed by a mature triangulation method in binocular stereo vision because the camera positions and postures of the image pairs during respective shooting are automatically calculated, and thus the reconstruction of the whole curve is completed.

The invention has the outstanding advantages of simple measurement hardware (a digital camera, a common personal computer, a staff gauge and a group of coding points generated by printing), flexible measurement mode (free shooting), unlimited measurement range, automatic combination of measurement data of all angles, no redundancy of measurement output data, convenient use, low cost and the like, can be used for positioning space points, can measure three-dimensional information such as edges, characteristic lines, key section control lines and the like on a measured object according to the requirement of mechanical product measurement modeling, and has wide application prospect in the fields of reverse engineering, product quality detection and the like.

Drawings

FIG. 1: the invention provides a basic flow chart of a measuring method.

FIG. 2: and (5) a schematic diagram of an encoding point. Fig. 2(a) depicts a coding dot structure, which is composed of a central white dot, a middle black ring, and an outer ring, wherein the outer ring is divided into 15 equal parts, the identity of each sector is determined according to the color distinguished by each sector, the black color represents a binary code "0", and the white color represents a binary code "1"; fig. 2(b) is an example of three encoding points.

FIG. 3: the embodiment is a partition diagram of a software graphical interface of a measurement system. 1. A menu area; 2. an icon toolbar; 3. displaying a list of image files; 4. displaying one image in the current moving image pair; 5. displaying the other image in the current moving image pair; 6. and a three-dimensional graphic display area of the reconstructed target curve. The images currently displayed in the two windows 4 and 5 are specified by clicking the image file list in the area 3, and after a user interactively outlines the approximate shape of a pair of homonymous curves in the two image windows 4 and 5 according to the identification curves on the images, the system automatically calculates the three-dimensional point column on the curves and displays the three-dimensional point column in the three-dimensional graphic area.

Detailed Description

The embodiment of the method for measuring the three-dimensional curve structure of the object provided by the invention is as follows: the digital camera adopts a Nikon manual focusing digital camera with a built-in flash lamp and the resolution of 4256 multiplied by 2848, the computer adopts a Pentium IV microcomputer with the main frequency of 2.8GHz and the memory of 512MB, and a measuring software system is realized on a Visual C + + 6.0 platform.

The specific embodiment and principle of the invention are described with reference to fig. 1: before measurement, certain measurement preparation is carried out, including marking a target curve (generally a characteristic line of an object and certain key section control lines of the surface of the object required for digital model reconstruction) to be measured on the measured object, so that the target curve is obviously different from the color of the measured object on the aspect of color brightness, and the image identification is facilitated; a ruler and a group of specially designed coding points are arranged around a measured object, and the codes of all the coding points in the group of coding points are different, namely each coding point has identity uniqueness. After the preparation work of the measurement is finished, a digital camera is held by hand to obtain a group of images of the measured object in a free shooting mode, and certain overlap is required among the images, namely, more than 5 coding points and certain target curves which are visible together at least between one image and the other image are required. Based on the set of images, the measurement system automatically identifies and precisely locates the encoded points in each image, and then automatically and precisely calculates the camera position and pose at each shot. The user finishes semi-automatic extraction of one identification curve in the current moving image pair through simple interaction, and the measurement system automatically performs optimization matching of the image curve with the same name, so that three-dimensional point list information of the identified curve structure is automatically calculated. If there are no unprocessed target curves, the semi-automatic extraction, automatic matching and three-dimensional curve reconstruction processes are repeated for the next target curve (which may appear in a different image pair) until the three-dimensional reconstruction of all target curves is completed. A detailed description of several embodiments of the main steps in fig. 1 is provided below.

Coding point identity recognition

The three-dimensional measurement method of the invention is based on the analysis of a set of images taken, first of all the identification of the coding points shown in fig. 2. The center of each code point is a circular target point, the periphery of the code point is an annular code band, the code band is divided into 15 parts according to the angle average, one part is arranged at every 24 degrees, the code band is equivalent to a binary bit, white is taken as a foreground color, the corresponding binary code is '1', black is taken as a background color, and the corresponding binary code is '0'. For each code point, there are 15 possible binary codes, and the decimal number corresponding to the smallest number of the 15 binary numbers is taken as the ID of the code point.

The automatic detection algorithm of the coding points mainly comprises the following three main processes: firstly, extracting a coding point target, namely finding a target point in an image; determining the unique identity of the coding point according to the information on the coding point 'coding strip', namely decoding the coding point; and locating the sub-pixel in the center of the coding point.

(1) Encoded point target extraction

The central circle of the coding point is elliptical after being imaged by the CCD. Therefore, firstly, Canny operator is adopted to carry out image segmentation, contour information representing different areas is extracted from the image, and then a step-by-step filtering method is adopted to carry out target extraction of the coding points. Firstly, preliminarily filtering possible target points according to the size and the shape of a marked point target, and entering a further identification process for a closed contour meeting the following conditions:

P_min≤P≤P_max (1)

1≤P²/4πA≤1.5 (2)

where P and A are the perimeter and area, respectively, of the closed contour, P_min、P_maxRespectively, a minimum and a maximum threshold for the perimeter of the contour. Equation (1) is a definition of the size of the closed contour, and equation (2) measures its proximity to a circle.

For the closed contour satisfying the formulas (1) and (2), ellipse fitting is carried out by adopting a least square method, and the residual error epsilon_eliSatisfies a given tolerance e_τHaving a square as the center circle of the candidate code point, i.e.

ε_eli≤ε_τ (3)

After least squares template matching, all elliptical contours in the image have been found. However, some outlines which are not mark point targets but have elliptical shapes or are close to the elliptical shapes are often mistaken for mark point targets in a real scene, but because the gray level of the foreground of the mark point targets adopted by the method is white,the background gray scale is black, and the contrast between the two is strong. This is a significant feature of the marked point object from other non-marked point objects. Therefore, non-coding point objects are further excluded according to such a gradation characteristic of the marker point. Since the ellipse fitting criterion shown by equation (3) has been passed, the areas of the inner central ellipse and the black circle can be determined. Recording the average value of the gray scale of the internal area of the central white point as M_IMean gray level of black circle region is M_OThen M is_IAnd M_OIt should satisfy:

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <msub> <mi>M</mi> <mi>I</mi> </msub> <msub> <mrow> <mo>&GreaterEqual;</mo> <mi>M</mi> </mrow> <mi>t</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>M</mi> <mi>O</mi> </msub> <msub> <mrow> <mo>&le;</mo> <mi>M</mi> </mrow> <mi>t</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>M</mi> <mi>I</mi> </msub> <mo>-</mo> <msub> <mi>M</mi> <mi>O</mi> </msub> <mo>&GreaterEqual;</mo> <mi>&Delta;</mi> <msub> <mi>M</mi> <mi>t</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein M is_tA threshold for distinguishing foreground gray from background gray; Δ M_tIs the minimum value that the difference between the foreground gray and the background gray should satisfy.

In addition, the gray variance V of the internal area of the central white point is restricted_IAnd the gray variance of the black circular ring area satisfies:

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <msub> <mi>V</mi> <mi>I</mi> </msub> <msub> <mrow> <mo>&le;</mo> <mi>&delta;</mi> </mrow> <mi>I</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mi>O</mi> </msub> <msub> <mrow> <mo>&le;</mo> <mi>&delta;</mi> </mrow> <mi>O</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein, delta_I，δ_OIs the maximum gray variance allowed. The condition (5) constrains that the center of the encoded point must satisfy a certain gray uniformity. If the above expressions (1) to (5) are satisfied, the encoding point decoding process is performed.

(2) Encoded point decoding

The specific implementation steps of the coding point decoding algorithm provided by the invention are as follows:

step 1: fitting the outer contour ellipse (marked as ellipse A) of the central dot of the coding point, the outer contour ellipse (marked as ellipse B) of the middle black ring and the outer contour ellipse (marked as ellipse C) of the ring in which each white sector is positioned. An ellipse is then fitted in the middle of the ellipse B, C, with the center and rotation angle being the same as the center and rotation angle of ellipse B, C, and the major and minor axes are averaged over the major and minor axes of B, C, respectively. And obtaining the position coordinates of each pixel point on the ellipse D by adopting an ellipse drawing algorithm.

Step 2: and calculating the median of all pixel grays in the area enclosed by the A as the foreground grayscale, and calculating the median of all pixel grays in the area between the A and the B as the background grayscale. And the average value of the foreground gray and the background gray is used as a threshold value for subsequently determining the code value of each binary bit of the coding point.

Step 3: for any pixel point TD on the ellipse D, a ray is made through the center of the ellipse, and the intersection points of the ray and the ellipse B, C are recorded as TB and TC. And sequencing the gray values of all pixels on the line segment TBTC, and taking the gray value of the middle pixel as the new gray value of the TD.

Step 4: and (3) carrying out inverse affine transformation on each point on the ellipse D according to a formula (6), so that the ellipse D corresponds to a unit circle, and the gray level of each point on the unit circle corresponds to the new gray level value of the ellipse B.

<math> <mrow> <msup> <mi>X</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <msup> <mi>a</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mtd> <mtd> </mtd> </mtr> <mtr> <mtd> </mtd> <mtd> <msup> <mi>b</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mtd> </mtr> </mtable> </mfenced> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mi>cos</mi> <mi>&alpha;</mi> </mtd> <mtd> <mi>sin</mi> <mi>&alpha;</mi> </mtd> </mtr> <mtr> <mtd> <mo>-</mo> <mi>sin</mi> <mi>&alpha;</mi> </mtd> <mtd> <mi>cos</mi> <mi>&alpha;</mi> </mtd> </mtr> </mtable> </mfenced> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msub> <mi>X</mi> <mi>o</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>

Wherein X' is the coordinate of a point on the unit circle corresponding to TD, X is the coordinate of point TD_oIs the coordinate of the center O of the ellipse D, a, b are the lengths of the major and minor axes, respectively, of the ellipse D, and α is the rotation angle of the ellipse D.

Step 5: and (4) carrying out binarization on the pixels on the unit circle, and taking one edge point as a starting point.

Step 6: starting from the starting point, every 24 ° on the unit circle is a binary bit, and the average gray value of all pixel points in each bit is calculated. If the average value of the gray scale of a certain bit is larger than the threshold value, the bit takes a binary code as '1'; otherwise, take "0". A binary coding of the encoded points is thus obtained. And finding the minimum decimal number corresponding to the binary number, wherein the decimal number is the ID of the code point.

Step3 corresponds to median filtering the point TD within a linear window, where the pixels included in the filtered window are the pixels on the segment TBTC. The new gray value on the ellipse D after median filtering is adopted to determine the binary code of each bit of the coding point, the gray values of all pixels in a coding band are considered, and the influence of isolated noise can be eliminated. The processing result of a large number of field shot images shows that the method is very effective for improving the robustness of the identification of the coding point.

(3) Marker center calculation

Center sub-pixel positioning of encoded points using equation (7)

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mi>j</mi> </munder> <munder> <mi>&Sigma;</mi> <mi>i</mi> </munder> <mi>i</mi> <mo>&CenterDot;</mo> <msub> <mi>I</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>/</mo> <munder> <mi>&Sigma;</mi> <mi>j</mi> </munder> <munder> <mi>&Sigma;</mi> <mi>i</mi> </munder> <msub> <mi>I</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mi>j</mi> </munder> <munder> <mi>&Sigma;</mi> <mi>i</mi> </munder> <mi>j</mi> <mo>&CenterDot;</mo> <msub> <mi>I</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>/</mo> <munder> <mi>&Sigma;</mi> <mi>j</mi> </munder> <munder> <mi>&Sigma;</mi> <mi>i</mi> </munder> <msub> <mi>I</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow> </math>

In the formula (x)_c，y_c) As center coordinates of the code points, I_i，jThe gray value of the pixel point (i, j) in the central circular area.

Camera pose automatic determination

In homogeneous coordinate representation, the projection X of a three-dimensional space point X on the camera imaging plane can be represented as:

x＝K[R|t]X＝PX (8)

wherein K is a camera intrinsic parameter matrix; r and t are respectively a rotation transformation matrix and a translation transformation vector from a world coordinate system to a camera coordinate system; p is a 3 × 4 projective transformation matrix.

The method comprises the steps of shooting two images of the same scene by a camera at two different positions and orientations, wherein a rotation matrix R exists between the two cameras₁₂And a non-zero translation vector t₁₂Then from the epipolar geometry, the following constraints exist between the two images

$x_2^{T}F{x}_1=0---\left(9\right)$

In the formula x₁And x₂Respectively is a projection point of a three-dimensional space point X on the first image and the second image; f is a 3 x 3 basis matrix that maps a point x on the right image₂To the corresponding epipolar line Fx of the left image₁The above.

According to the corresponding relation of the coding points between the two images 1, L, N, N is more than or equal to 5, and a basic matrix F between two images is calculated by using an MLESAC (maximum Likelihood Estimation SAmple consensus) method.

According to the internal parameters such as the focal length and the like marked by the camera, the initial value of the internal parameter matrix K can be constructed (the product marking parameters are only used as the initial values, and the measurement system can further optimize the product marking parameters in the subsequent process). Thus, from the basis matrix F, the essential matrix E between the two images can be further calculated

E＝K^TFK (10)

And according to the definition of the essence matrix, E ═ t]_×R (in the formula [ g ]]_×An anti-symmetric matrix representing a vector) can be easily derived using the orthogonality of the rotation matrices

${\hat{E}}^T\hat{E}=\left[\begin{array}{ccc}1-{\hat{t}}_x^2& -{\hat{t}}_x{\hat{t}}_y& -{\hat{t}}_x{\hat{t}}_z\\ -{\hat{t}}_y{\hat{t}}_x& 1-{\hat{t}}_y^2& -{\hat{t}}_y{\hat{t}}_z\\ -{\hat{t}}_z{\hat{t}}_x& -{\hat{t}}_z{\hat{t}}_y& 1-{\hat{t}}_z^2\end{array}\right]---\left(11\right)$

In the formula $\hat{E}=E/\sqrt{\mathrm{Tr}\left(E^{T}E\right)/2},$ Tr (g) represents the trace of the matrix, ${\hat{t}}_{12}=t_{12}/\left|\right|t_{12}\left|\right|$ is the normalized translation vector. Thus, the normalized translation vector can be easily obtained from the E matrix obtained by expression (10) and expression (11) ${\hat{t}}_{12}={({\hat{t}}_x,{\hat{t}}_y,{\hat{t}}_z)}^T.$ Due to the fact that ${(-\hat{E})}^T(-\hat{E})={\hat{E}}^T\hat{E},$ Thus, the resulting matrix is normalized

May be compared with realityDiffering by one symbol. In addition, due to the matrix

Is about a vector

The quadratic term of (a), the vector thus calculated according to equation (11)Also has ambiguity, i.e.All satisfy the formula (11). The solution will be given laterAnd

the method of (4).

For calculating a rotation matrix R between a first image and a second image₁₂Definition of

<math> <mrow> <msub> <mi>w</mi> <mi>i</mi> </msub> <mo>=</mo> <msub> <mover> <mi>E</mi> <mo>^</mo> </mover> <mi>i</mi> </msub> <mo>&times;</mo> <msub> <mover> <mi>t</mi> <mo>^</mo> </mover> <mn>12</mn> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>=</mo> <mn>1,2,3</mn> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow> </math>

In the formulaRepresentative matrixEach row vector of (2). Let r_iIs a rotation matrix R₁₂Each row vector of (1), then

r_i＝w_i+w_j×w_k (13)

Where (i, j, k) is a cyclic combination of (1, 2, 3). The camera poses of the first two views are thus determined.

Establishing a world coordinate system on a first camera, and deriving a Z-direction coordinate of a space point X in the first camera coordinate system according to the imaging geometric relation

${\mathrm{<figure-callout\; id="1"\; label="Z"\; filenames="A20061016127400182.png"\; state="\{\{state\}\}">Z</figure-callout>}}_1=f\frac{{({\mathrm{fr}}_1-x_2{r}_3)}^T{\hat{t}}_{12}}{{({\mathrm{fr}}_1-x_2{r}_3)}^T{x}_1}---\left(14\right)$

Further, two other coordinate components can be obtained

X₁＝x₁Z₁/f，Y₁＝y₁Z₁/f (15)

X in the second camera coordinate system is

X₂＝R₁₂(X₁-t₁₂) (16)

Due to the fact that

Andambiguity of (D) may be producedFour pairs of differentAccording to the actual situation when shooting, only when a certain shooting is adoptedAll points for reconstruction are in front of both cameras simultaneously, i.e. Z for all points (here referring to the respective coded point centers visible in common for both views)₁And Z₂All are positive, the reconstruction result is correct, and a corresponding groupThe correct solution is obtained.

Because the length of the base line between the two cameras is unknown in the reconstruction algorithm, only normalized translation vectors can be obtainedAs can be easily seen from equation (14), the reconstructed scene differs from the actual scene by a fixed scale factor. Therefore, a scale is placed in the scene, and the distance between two marking points on the scale is known, so that a scale factor can be determined, and the actual size of the measured object can be obtained.

On the basis of the two-view camera pose determination and the encoding point center three-dimensional coordinate reconstruction, the camera poses of other shooting times need to be further determined in sequence. When the jth image is processed, at least 6 coded points whose three-dimensional coordinates have been reconstructed in the previous step are required to exist in the image, i.e. the corresponding X of the known image point and the space point_i*x_iI is 1, L L, L is equal to 6. These constraints are introduced into projection equation (8) to obtain

x_i＝P_jX_i，i＝1，L L，L≥6 (17)

Since each group X_i*x_iGenerates two linear equations, so that the projection matrix P of the jth image can be solved by the least square method according to the formula (17)_j11 unknown elements.

A 3 × 4 projection matrix P_jIs shown as

P_j＝K[R_j|t_j]＝[KR_j|Kt_j]＝[M|p₄] (18)

Where M is a matrix P_jFirst 3 × 3 sub-matrix, p₄Representation matrix P_jThe fourth column of (2). The translation vector t is easily determined from equation (18)_j＝K^-1p₄。

Because the internal parameter matrix is triangular and the rotation matrix is orthogonal, QR decomposition is carried out on the matrix to obtain the rotation matrix R_j. Estimating external attitude parameter R of camera relative to world coordinate system when j image is shot_jAnd t_jOn the basis, the three-dimensional space point coordinates of the encoding points which newly appear in the j image and can be matched with the same name in the previous j-1 image are reconstructed by an optical triangulation method. Thus, the camera pose determination and the encoding point three-dimensional coordinate calculation corresponding to the jth image are completed. And then continuing to process the next image until all images are processed.

(3) Camera pose optimization

Three-dimensional space point X due to image noise and other factors_iProjected matrix P_jTransformed image point, with the actually recognized X_iCoordinates x of image points in the jth image_ijAnd do not coincide. In order to further improve the system precision, the invention establishes an objective function with minimum reprojection error based on a beam adjustment algorithm

<math> <mrow> <munder> <mi>&Sigma;</mi> <mi>ij</mi> </munder> <mi>d</mi> <msup> <mrow> <mo>(</mo> <msub> <mi>P</mi> <mi>j</mi> </msub> <msub> <mi>X</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>x</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>&RightArrow;</mo> <mi>min</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow> </math>

And carrying out global optimization on the camera parameters and the three-dimensional space point coordinates obtained in the previous step. The specific solution adopts LM (Levenberg-Marquardt) algorithm. Since it is X which has been found previously and is relatively close to the true value_iAnd P_jAs an initial value, therefore, the global optimization can converge faster.

Semi-automatic extraction of target curves

The graphical interface schematic diagram of the three-dimensional measurement software system of the object curve structure developed by the invention is shown in fig. 3, the left side of the diagram displays a list of all the shot image files, the lower right side of the diagram is two image display windows of the current active image pair, the currently displayed images in the two windows are specified by clicking the image file list on the left side, and the upper right side of the diagram is a three-dimensional graph display area of the reconstructed target curve.

The invention adopts the basic idea of energy optimization to automatically 'fit' the approximate contour of a target curve mutually outlined on two moving images to the corresponding target curve through a background algorithm.

In the concrete implementation, firstly, the points on the target curve inputted by user are connected to form a broken line (polygon under closed condition), then the invention adopts DDA algorithm of linear segment raster scan conversion in computer graphics to quickly obtain all the pixel points passed by broken line, and makes them pass through according to fixed interval (in the embodiment, every two pixels are taken one)

Sampling, denoted v_iAnd i is (0, 1, L, n). Here, since it is necessary to further automatically detect image edge information by using a Canny edge detection operator, the detected edge point set is denoted as P. Note additionally v_ijJ is (1, L, 8) is a point v_iEight neighborhood pixel points, and for descriptive convenience, note v_i0＝v_i。

The invention establishes each point v_iAnd the following energy functions of the eight neighborhood points:

E(v_ij)＝[αE_tension(v_ij)+βE_bend(v_ij)+γE_img(v_ij)+δE_attr(v_ij)] (20)

wherein E is_tension(v_ij)，E_bend(v_ij)，E_img(v_ij)，E_attr(v_ij) Are each v_ijThe stretching energy, the bending energy, the image energy and the energy generated by the attraction of the edge points at the points, and alpha, beta, gamma and delta are respectively the weight of each energy item and are used for adjusting the specific gravity of each energy item. To balance the effects of the terms, the energy terms are all normalized to the interval [0, 1 ]]：

<math> <mrow> <msub> <mi>E</mi> <mi>tension</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mover> <mi>d</mi> <mo>&OverBar;</mo> </mover> <mo>-</mo> <mo>|</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mrow> <munder> <mi>max</mi> <mrow> <mn>0</mn> <mo>&le;</mo> <mi>j</mi> <mo>&le;</mo> <mn>8</mn> </mrow> </munder> <mo>{</mo> <mo>|</mo> <mover> <mi>d</mi> <mo>&OverBar;</mo> </mover> <mo>-</mo> <mo>|</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <mo>|</mo> <mo>}</mo> </mrow> </mfrac> <mo>,</mo> </mrow> </math> <math> <mrow> <mover> <mi>d</mi> <mo>&OverBar;</mo> </mover> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <mo>|</mo> <msub> <mi>v</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> </mrow> </math>

<math> <mrow> <msub> <mi>E</mi> <mi>bend</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <msup> <mrow> <mo>|</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <mn>2</mn> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>+</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mrow> <munder> <mi>max</mi> <mrow> <mn>0</mn> <mo>&le;</mo> <mi>j</mi> <mo>&le;</mo> <mn>8</mn> </mrow> </munder> <mo>{</mo> <msup> <mrow> <mo>|</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mrow> <mn>2</mn> <mi>v</mi> </mrow> <mi>ij</mi> </msub> <mo>+</mo> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> </mrow> <mn>2</mn> </msup> <mo>}</mo> </mrow> </mfrac> </mrow> </math> <math> <mrow> <msub> <mi>E</mi> <mi>img</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mi>min</mi> <mrow> <mn>0</mn> <mo>&le;</mo> <mi>j</mi> <mo>&le;</mo> <mn>8</mn> </mrow> </munder> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mi>img</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>E</mi> <mi>img</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <munder> <mi>max</mi> <mrow> <mn>0</mn> <mo>&le;</mo> <mi>j</mi> <mo>&le;</mo> <mn>8</mn> </mrow> </munder> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mi>img</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mo>-</mo> <munder> <mi>min</mi> <mrow> <mn>0</mn> <mo>&le;</mo> <mi>j</mi> <mo>&le;</mo> <mn>8</mn> </mrow> </munder> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mi>img</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </math>

<math> <mrow> <msub> <mi>E</mi> <mi>attr</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>-</mo> <msub> <mi>p</mi> <mi>ij</mi> </msub> <mo>|</mo> </mrow> <mrow> <munder> <mi>max</mi> <mrow> <mn>0</mn> <mo>&le;</mo> <mi>j</mi> <mo>&le;</mo> <mn>8</mn> </mrow> </munder> <mo>{</mo> <mo>|</mo> <msub> <mi>v</mi> <mi>ij</mi> </msub> <mo>-</mo> <msub> <mi>p</mi> <mi>ij</mi> </msub> <mo>|</mo> <mo>}</mo> </mrow> </mfrac> </mrow> </math>

Wherein E is_img＝-|▽I(x，y)|²，E_attr(v_ij) P in (1)_ijE P is and v_iThe closest edge point. The invention adds the edge point gravitational energy E outside the stretching energy, the bending energy and the image energy_attr(v_i) The purpose of (1) is to further promote convergence of the point array towards the target curve. p is a radical of_ijIs limited to v_iWithin a window at the center. If there are no edge points in the window, E_attr(v_ij)＝0，j＝0，L，8。

By iterative method making the dot sequence v_iAnd i is moved to a position where the energy of the equation (20) is extremely small (0, 1, L, n), and finally locked near the image feature, thereby forming a smooth target curve point row. The points are then fitted with a uniform B-spline curve for subsequent curve matching. For the open curve, the two end points of the open curve are constrained to be kept unchanged all the time, thereby preventing the open curve from being openedThe contraction of the row of points of the curve to a point produces degradation.

Automatic matching optimization of homonymous curves

(1) Measure of degree of match

One basic constraint that a camera should satisfy for a same-name point between a pair of images obtained at different positions and angles is the epipolar constraint described by equation (9). In the case where the intra-camera parameters and relative positions and poses are known, the fundamental matrix F between the pairs of images is known. For two candidate matching points v in an image pair₁And v₂V can be used on the second image₂To v₁Distance of polar line of

$D_2({\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="A20061016127400182.png"\; state="\{\{state\}\}">v</figure-callout>}}_1,v_2)=\frac{\left|v_2^{\mathrm{<figure-callout\; id="1"\; label="T\; F\; v"\; filenames="A20061016127400182.png"\; state="\{\{state\}\}">T</figure-callout>}}F{v}_{}{}_1\right|}{\left|\right|\mathrm{<figure-callout\; id="1"\; label="eF\; v"\; filenames="A20061016127400182.png"\; state="\{\{state\}\}">eF</figure-callout>}{v}_{}{}_1\left|\right|}---\left(21\right)$

And v on the first image₁To v₂Distance of polar line of

$D_1({\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="A20061016127400182.png"\; state="\{\{state\}\}">v</figure-callout>}}_1,v_2)=\frac{\left|v_2^{\mathrm{<figure-callout\; id="1"\; label="T\; F\; v"\; filenames="A20061016127400182.png"\; state="\{\{state\}\}">T</figure-callout>}}F{v}_{}{}_1\right|}{\left|\right|v_2^{T}F{e}^T\left|\right|}---\left(22\right)$

Measure v₁And v₂The degree of matching of (a), wherein, $e=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right].$ the invention is based on the epipolar constraint conditions of the formula (21) and the formula (22) to establish the optimization objective function of point pair matching between homonymous curves.

(2) Image curve resampling

To obtain a fitted uniform B-spline curve p (u), u ∈ [0, 1 ]]Taking the discrete parameter interval delta u as 1/L and L as the accumulated chord length of each type value point to obtain the discrete pixel points v on the uniform B spline curve_iI is 0, 1,.. times.n, constructing a piecewise linear interpolation curve satisfies:

<math> <mrow> <mfenced open='{' close=''> <mtable> <mtr> <mtd> <mi>c</mi> <mrow> <mo>(</mo> <msub> <mi>l</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>v</mi> <mi>i</mi> </msub> <mo>,</mo> </mtd> <mtd> <mi>i</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <mi>N</mi> </mtd> </mtr> <mtr> <mtd> <mi>c</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>l</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <mi>l</mi> </mrow> <mrow> <msub> <mi>l</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>l</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>v</mi> <mi>i</mi> </msub> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <msub> <mi>l</mi> <mi>i</mi> </msub> </mrow> <mrow> <msub> <mi>l</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>l</mi> <mi>i</mi> </msub> </mrow> </mfrac> <msub> <mi>v</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> </mtd> <mtd> <msub> <mi>l</mi> <mi>i</mi> </msub> <mo>&le;</mo> <mi>l</mi> <mo>&lt;</mo> <msub> <mi>l</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow> </math>

wherein l₀＝0.0， <math> <mrow> <msub> <mi>l</mi> <mi>i</mi> </msub> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>i</mi> </munderover> <mo>|</mo> <msub> <mi>v</mi> <mi>j</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mrow> <mi>j</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <mo>.</mo> </mrow> </math> Note the book <math> <mrow> <mi>L</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mo>|</mo> <msub> <mi>v</mi> <mi>j</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mrow> <mi>j</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <mo>.</mo> </mrow> </math> Without loss of generality, one of the two homonymous image curves with more pixel points is marked as c₁(l) The less one is denoted as c₂(l)。c₁(l) Has a discrete pixel point of v_k ⁽¹⁾，0≤k≤N₁；c₂(l) Has a discrete pixel point of v_j ⁽²⁾，0≤j≤N₂. Seeking c₁(l)，l∈[0，L₁]Each pixel point on c₂(l)，l∈[0，L₂]Can realize that it is at c₂(l) The sub-pixels are matched, so that high-precision three-dimensional curve reconstruction is realized.

(3) Matching optimization

Dynamic programming is first employed to initially match discrete pixel points on the corresponding curve. The cumulative cost function of a corresponding point pair of the homonymous curve is defined as

<math> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mi>k</mi> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>v</mi> <mi>j</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mi>D</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mi>k</mi> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>v</mi> <mi>j</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <munder> <mi>min</mi> <mrow> <mi>m</mi> <mo>&Element;</mo> <msub> <mi>G</mi> <mi>kj</mi> </msub> </mrow> </munder> <mi>C</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>v</mi> <mi>m</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>24</mn> <mo>)</mo> </mrow> </mrow> </math>

Wherein $D(v_k^{\left(1\right)},v_j^{\left(2\right)})=D_1(v_k^{\left(1\right)},v_j^{\left(2\right)})+D_2(v_k^{\left(1\right)},v_j^{\left(2\right)}),$ G_kjIs shown at v_k ⁽¹⁾And v_j ⁽²⁾All possible values of m in the case of already matching. Because the number of pixel points on the two homonymous image curves is unequal, the result of dynamic programming matching can have many-to-one condition, namely longer c₁A plurality of pixel points on the curve correspond to shorter c₂The same pixel point on the curve.

And after initial matching of the point pairs on the image curve is obtained, matching optimization of the curve is carried out. Based on equations (21) and (22), the present invention optimizes the following objective function to achieve c₁、c₂Exact matching of points on a curve

<math> <mrow> <mi>min</mi> <munderover> <mo>&Integral;</mo> <mn>0</mn> <msub> <mi>L</mi> <mn>1</mn> </msub> </munderover> <mfrac> <mrow> <mo>|</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>F</mi> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <mi>eF</mi> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>F</mi> <msub> <mi>c</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>c</mi> <mn>2</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>F</mi> <msup> <mi>e</mi> <mi>T</mi> </msup> </mrow> </mfrac> <mi>dl</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>25</mn> <mo>)</mo> </mrow> </mrow> </math>

Where σ (l) is the mapping function to be solved, representing c₁The point on the curve with the parameter l is on the curve c₂The parameter values of (1). The form of integration of equation (25) is rewritten into the form of summation

(26) In the formula I₀，l₁，L，l_N1Are each v_k ⁽¹⁾，0≤k≤N₁Corresponds to c₁(l) The parameter values of (1). The coarse matching generated by the dynamic programming method is used as an iteration initial condition, and the minimization problem of the formula (26) can be solved by adopting a conjugate gradient method.

Three-dimensional reconstruction of curves

After matching of all pixel points on the homonymous curve is completed, because the parameters in the camera are known, and the positions and postures of the camera at the time of respective shooting of the image pair are automatically calculated, namely the rotation matrix and the translation vector relative to the world coordinate system are also known, the space coordinates of the points can be reconstructed by a mature triangulation method in binocular stereo vision, and a specific calculation formula can be expressed as

$\left.\begin{array}{c}{x}^{\left(j\right)}=f\frac{R_{11}^{\left(j\right)}X+R_{12}^{\left(j\right)}Y+R_{13}^{\left(j\right)}Z+T_x^{\left(j\right)}}{R_{31}^{\left(j\right)}X+R_{32}^{\left(j\right)}Y+R_{33}^{\left(j\right)}Z+T_z^{\left(j\right)}}\\ y^{\left(j\right)}=f\frac{R_{21}^{\left(j\right)}X+R_{22}^{\left(j\right)}Y+R_{23}^{\left(j\right)}Z+T_y^{\left(j\right)}}{R_{31}^{\left(j\right)}X+R_{32}^{\left(j\right)}Y+R_{33}^{\left(j\right)}Z+T_z^{\left(j\right)}}\end{array}\right\},j=\mathrm{1,2}---\left(27\right)$

Where f is the focal length of the camera, x^(j)And y^(j)Two components, R, of the pixel coordinates of the homonymous matching point on the jth image, respectively.^(j)Is the component, T, of the rotation matrix of the j-th camera with respect to the world coordinate system.^(j)The translation components of the j-th shooting phase relative to the world coordinate system are obtained, and 3 unknown components (X, Y, Z) of the space coordinate of the point can be solved by a least square method according to 4 equations in the formula (27). And (4) executing the solving process for each matching point pair on the homonymous curve, and completing the three-dimensional measurement of the whole target curve point column. After the three-dimensional measurement of all target curves is completed, a parameter equation of the curve and a curved surface equation of the model can be further constructed according to the point lists.

Claims (5)

1. A method for carrying out three-dimensional measurement on an object by utilizing free shooting of a single digital camera is characterized by comprising seven steps of measurement preparation, image shooting, identification and positioning of encoding points, camera pose determination, target curve extraction, automatic matching optimization of homonymous curves and three-dimensional reconstruction of target curves, and the specific method comprises the following steps: firstly, marking a characteristic line of a measured object and a key section control line on the surface of the object required for digital model reconstruction, so that the characteristic line is obviously different from the color of the measured object in color brightness, and the identification of images is facilitated; a ruler and a group of coding points are arranged around a measured object; then, holding a digital camera by hand to obtain a group of images of the measured object in a free shooting mode; according to the group of images, the position and the posture of the camera during each shooting are automatically and accurately calculated, meanwhile, a convenient interaction means is provided for a user, the semi-automatic extraction of the identification curve and the optimized matching of homonymous curves in different images are realized, and further, the three-dimensional point list information of the identified curve structure is automatically calculated.

2. The method of claim 1, wherein the identification of the code points is that the codes of each point in a set of designed code points are different, and the candidate code point targets are gradually filtered by adopting 5 constraint conditions of size, shape, ellipse fitting residual, area gray mean, and area gray variance; the gray values of most pixels in a coding band are comprehensively considered in the coding point decoding process, so that the effects of eliminating noise influence and improving the robustness of the identification of the coding points are achieved; and (4) positioning the sub-pixels in the center of the coding point by adopting a gray scale weighting method in the target area.

3. The method of claim 1, wherein the camera pose determination is based on identity uniqueness of the code points, and automatically establishes correspondence between the code points of the same name in each image; and recovering the camera postures corresponding to the two images and the three-dimensional coordinates of the centers of the coding points which are commonly visible in the two images according to the pixel coordinates of the centers of at least 5 coding points with the same name in the two images, and further adopting a strategy of combining an incremental method and a global optimization method to realize the automatic solution of the camera positions and postures corresponding to all the images.

4. The method as claimed in claim 1, wherein the target curve extraction is performed by a user by using a mouse to pick points on the two images near the same identification curve, and making the connecting line of the points approximately reflect the corresponding image curve contour, and then the measurement software system iteratively minimizes the energy of the image curve contour, thereby automatically and optimally fitting the interactively sketched approximate image curve contour to the image curve.

5. The method as claimed in claim 1, wherein the automatic matching optimization of the homonymous curve is performed by using a dynamic programming method to obtain an initial matching of pixel points on the corresponding curve, and then using a nonlinear optimization method to minimize the sum of distances from all matching points on the corresponding curve to respective polar lines, so as to achieve an optimal matching of the homonymous target curve point rows on the image, and further calculate the three-dimensional coordinates of the point rows.

Images