Method for real-time measurement of length-width size distribution of crystal population in crystallization reactor using binocular telecentric cameras

Abstract

The present invention relates to the technical field of industrial process control and detection. Disclosed is a method for real-time measurement of length-width size distribution of a crystal population in a crystallization reactor using binocular telecentric cameras. In-situ environment calibration is performed for a binocular telecentric stereo vision system, a calibration rod capable of probing into an in-situ environment (a reactor / glass tube) is designed, a binocular telecentric stereo vision imaging model is established, and a simple calibration method of rotating the calibration rod suitable for an in-situ limited space and a calibration plate design scheme are proposed. In order to improve binocular matching efficiency, a simplified telecentric stereo epipolar rectification method is further provided, so as to ensure accurate matching of on-site snapshot image pairs. On the basis of the matched image pairs, a three-dimensional reconstruction method using analytical solution-based ray intersection is provided to measure the three-dimensional pose of particles within a crystallizer. Finally, the three-dimensional length and width of a crystal are quantitatively evaluated by means of statistical data of Euclidean distances of pairs of length and width feature points of the crystal. The present invention has high operability, and can achieve the effect of automatically measuring the three-dimensional size of crystals.

Description

A method for real-time detection of long and wide size distribution of crystal seed population in a crystallization kettle by using binocular telecentric camera TECHNICAL FIELD

[0001] The present application belongs to the field of industrial process control and detection technology, and relates to a method for real-time detection of long and wide size distribution of crystal seed population in a crystallization kettle by using binocular telecentric camera. BACKGROUND

[0002] The process and control technology of crystallization process plays a very important role and significance for the high-quality development of China's advanced manufacturing field. In the past few decades, imaging systems have been increasingly used to measure particle size in various crystallization processes in laboratory and industrial scale crystallizers. Few literatures and patents at home and abroad have studied the influence of refraction generated in the cross-medium optical system, such as a glass crystallizer filled with solution and a glass jacket for heat exchange.

[0003] In order to ensure the accuracy of on-site particle size measurement based on microscopic image analysis, it is crucial to establish an accurate camera calibration model. However, compared with the existing pinhole imaging system calibration method, the calibration of telecentric system has not been fully studied in the past few years. For the problem of deformation measurement using telecentric lens, a simple calibration method based on linear fitting is proposed, which actually does not consider lens distortion and can only get the internal parameter magnification. Li proposed a planar pattern calibration method for imaging using single telecentric lens in the international journal "Optics and Lasers in Engineering", but did not solve the ambiguity problem of external parameters. Further, Chen et al. developed a calibration method by providing Z-axis translation displacement for the calibration pattern by a micro-positioning stage, which can solve the sign ambiguity problem in the planar pattern-based calibration process. Liu et al. proposed a two-step calibration method based on a three-dimensional calibration board in the international journal "Optics Express", but due to the very limited accuracy of the three-dimensional calibration mode, a planar calibration mode is still needed. In order to adapt to the limited in-situ space of calibration target motion, Cheng et al. proposed a simplified telecentric camera calibration method without auxiliary equipment in the top journal "IEEE Transactions on Instrumentation and Measurement", which uses a Scheimpflug telecentric model similar to the specific pinhole camera model to solve the related sign problem. For telecentric systems without Scheimpflug properties, this strategy may be ineffective, resulting in larger reprojection errors. In addition, the above calibration methods require auxiliary equipment for on-site calibration. In the limited in-situ calibration motion space, how to achieve accurate calibration and reconstruction simply and quickly is still a challenging problem.

[0004] For the crystal size measurement aspect of the crystallization process, there are few reports on the practical application of in-situ detection image analysis methods. Gao and Cardona et al. proposed an ellipse fitting method for measuring different crystal shapes in the international top chemical journals "Crystal Growth & Design" and "Chemical Engineering Science". These methods are limited to measuring the crystal size distribution (CSD) based on two-dimensional image projection. Zhang and Huo et al. developed a simplified three-dimensional geometric model of crystal morphology to calculate the crystal size in the volume space and published it in an international chemical journal. However, the accuracy of these methods depends on the key feature points identified to construct the geometric model. Fan et al. proposed a size measurement technology based on the reconstruction of key corner points of the crystal image in the top journal "IEEE Transactions on Instrumentation and Measurement". However, how to accurately obtain the geometric information of the crystal in three-dimensional space to accurately measure the CSD and its surface area remains a problem to be solved. SUMMARY

[0005] The technical problem to be solved by the present application is to calibrate the refractive effect of the reaction kettle medium on the optical system in-situ and measure the three-dimensional length and width of the crystal in the crystallization process. To solve the above problem, an in-situ calibration rod design scheme and a binocular telecentric camera imaging calibration modeling method are systematically proposed, as well as a telecentric stereo vision three-dimensional reconstruction technology method based on analytical solution, to realize the in-situ online detection of the three-dimensional length and width of the crystal in the solution.

[0006] The technical scheme adopted by the present application is as follows:

[0007] A method for real-time detection of the length and width size distribution of the crystal population in the crystallization kettle by binocular telecentric cameras, which is realized based on a non-contact high-resolution backlit calibration rod, ensures that two telecentric cameras synchronously collect images of the crystal solution in the reaction kettle, and thus measures the three-dimensional size of the crystal. The backlit calibration rod is composed of a light-emitting diode (LED), a ceramic calibration plate printed with a checkerboard pattern, a reflecting prism, an optical fiber, and a stainless steel sheath. The light emitted by the LED at the top of the calibration rod is guided to the reflecting prism at the bottom of the calibration rod through the optical fiber, and the reflecting prism reflects the light vertically to the back of the ceramic calibration plate engraved with the checkerboard pattern. The reflected light provides illumination for the shooting of the two telecentric cameras. This device can provide sufficient imaging light intensity in a short exposure time (such as 36-100 μs), thereby reducing the motion blur caused by the turbulent flow of the solution in the crystallizer and other environmental disturbances. The calibration rod is inserted into the glass crystallizer and rotated clockwise or counterclockwise along its pivot axis to perform on-site calibration. The backlit calibration rod can be assembled into different sizes to adapt to the on-site calibration of crystallizers of different volumes.

[0008] The present application synchronously acquires images at two different angles outside the reactor, which are called binocular images, including left view and right view. A two-step in-situ stereo imaging calibration model is established for the pose of binocular camera, including calculating the closed-loop solution of camera parameters to obtain the initial value of model parameters and subsequently using optimization algorithm to seek optimal camera model parameters. For the acquired binocular images, a crystal image matching analysis method is proposed, including image preprocessing, interest point detection, feature matching and mismatch removal. After detecting the long and wide key angle points of the crystal by using the crystal contour feature, the three-dimensional reconstruction of the long and wide angle points of each matched crystal is carried out through the calibrated stereo imaging model, so as to quantitatively evaluate the three-dimensional size. Specifically, the following steps are included:

[0009] First step, stereo imaging calibration

[0010] The imaging model of telecentric camera, i.e. the projection of point P(xw, yw, zw) in the world coordinate system to pixel coordinate system P(u, v), is expressed as

[0011] where m is the effective magnification of telecentric lens, and R and T are the rotation and translation matrices between the two coordinate systems, respectively.

[0012] (1) Monocular camera parameter calibration:

[0013] The world coordinates of the planar calibration board can be related to the captured image coordinates through the homography matrix H, and the relationship is as follows:

[0014] For ease of analysis, Euler angles are introduced to describe the rotation matrix R(α, β, θ), and α, β and θ are the rotation angles around the X, Y and Z axes of the coordinate system, respectively. The elements of the corresponding rotation matrix are described as:

[0015] By equating the corresponding elements of the 3x3 matrix in (2), the following equations can be easily established:

[0016] wherein α, β and θ are rotation angles in the range of (-π / 2, π / 2). In addition, the effective magnification factor is defined as a positive value. The solution of the above equation is:

[0017] where c=(h 11 h 22 -h 12 h 21 ) 2 . The translation matrix T(t x ,ty ) are:

[0018] Therefore, the magnification factor m and the Euler angle θ can be uniquely determined. However, the Euler angles a and β have two possible signs, as shown in (5). Therefore, there can be two ambiguous poses related to the world coordinate system in the single camera coordinate system, and four ambiguous poses in the binocular camera system. Recovering the true Euler angles a and β is the key to the calibration method of the system.

[0019] Due to the orthographic projection property of the above-mentioned telecentric camera, the camera pose recovered from the two-dimensional pattern is ambiguous. The existing method relies on a three-dimensional calibration pattern with limited accuracy, or makes full use of external displacement, but cannot be used for on-site calibration of the stereo vision system because it cannot establish a physical connection between the imaging target and the stereo vision system. In order to overcome the above problems, a practical rotating template calibration method is proposed.

[0020] During the calibration process, the calibration pattern is rotated around the Y w axis at a specific angle sequence (clockwise or counterclockwise) so that the camera can obtain an ordered calibration shot image. Therefore, it is observed that the placement angle of the calibration pattern changes to be sequentially increased (or decreased). At the same time, the angle of the rotation Euler angle β of the camera coordinate system relative to the world coordinate system is also monotonically increasing (or decreasing). For ease of analysis, it is assumed that the calibration pattern is rotated clockwise around the Y w axis.

[0021] Considering that the actual rotation angle β is a monotonically increasing sequence, the absolute value of β may present a local minimum point or remain a monotonous trend. Calculating the absolute value of each function can check whether there is a local minimum point. If there are multiple local minimum points or there are local maximum points, the calibration condition is not met. However, if the specified condition is met, it can be inferred that the extreme value point will appear near the minimum value of the discrete sequence, resulting in

[0022] Through a number of discrete, the extreme value point can be determined to be near the minimum value, such as the left or right side of the minimum value. Since the rotation angle β is monotonically increasing, in addition to the sign of the minimum value point, the following equation is used to recover the sign:

[0023] The rotation angle a is recovered by the following formula:

[0024] (2) Stereo vision system parameter calibration:

[0025] The essence of binocular vision calibration is to obtain the relative position relationship between the two telecentric cameras. The left and right camera coordinate position relationship is represented by a rotation matrix R CTranslation vector T C Description, i.e.

[0026] where R L and T L are the rotation matrix and translation vector associated with the left and right camera coordinate system, respectively. R and T R are the rotation matrix and translation vector associated with the left and right camera coordinate system, respectively.

[0027] For each pair of calibrated images, the extrinsic matrix of the left and right cameras changes, while the relative rotation matrix R C between the left and right cameras remains unchanged. The Euler angle representation of R C obtained from the kth pair of calibrated images is given by

[0028] where α C , β C , θ C are the rotation angles around the X, Y, and Z axes of the left camera coordinate system. The rotation Euler angles and represent the extrinsic parameters R L and R R in the left and right camera coordinate system, respectively. The superscript k indicates that the data come from the kth set of calibrated images.

[0029] The Euler angles α C , β C , θ C obtained from N pairs of different calibrated images are optimized

[0030] where N k is the number of variables in the set of calibrated images.

[0031] The calculated rotation matrix Euler angle representation is converted to matrix representation R C , and T C is obtained according to (11), completing the stereo calibration.

[0032] Second step, simplified stereo rectification

[0033] Epipolar rectification plays an important role in accurate stereo matching, which simplifies the two-dimensional search of matching points to one-dimensional search. In telecentric stereo vision, existing methods use a 3x4 matrix for epipolar rectification. In fact, from the geometric and physical point of view, the rectification parameters can be further simplified to a 3x3 homography matrix, which is more convenient for computer-aided image analysis. The specific method is as follows:

[0034] The projection of point P from the world coordinate system to the left and right camera pixel coordinate systems p L and p R is represented as

[0035] The projection of point P from the world coordinate system to the corrected left and right camera pixel coordinate system is determined as follows:

[0036] According to formulas (14) and (15), the conversion from the original image coordinates to the corrected coordinates is established:

[0037] Thirdly, reconstruction and measurement

[0038] (1) Three-dimensional reconstruction:

[0039] An improved ray intersection three-dimensional reconstruction method based on analytical solution is proposed for telecentric stereo vision system.

[0040] Firstly, the position equation of each ray passing through the matching point pair on the left and right physical image planes is constructed respectively.

[0041] where (x l ,y l ,z l ) represents the coordinates of any point on the left image plane ray, represents the coordinates of the matching point on the left physical image plane, represents the direction unit vector of the left camera optical axis, τ L is the ray coefficient of the left camera, (x r ,y r ,z r ) represents the coordinates of any point on the right image plane ray, represents the coordinates of the matching point on the right physical image plane, represents the direction unit vector of the right camera optical axis, τ R is the ray coefficient of the right camera. If the calculated point on the physical image plane and the position equation of the above ray are accurate, there is only one intersection point between the left and right rays, i.e. the reconstructed point on the object surface. However, due to the existence of noise, the left and right rays may not intersect. In order to find the best intersection point, the two closest points to the left and right rays are calculated respectively, and the coefficient of the left ray when the closest point position is obtained by calculating the analytical solution:

[0042] where

[0043] Correspondingly, the coordinates of the reconstructed point are obtained by the following formula:

[0044] (2) Crystal image matching and measurement:

[0045] Firstly, the image is corrected by the correction method in step two, and the crystal image is quickly matched by one-dimensional search. Then, the length feature point pairs are obtained by the distance from the contour to the contour center point, and the crystal main axis is obtained by connecting the length feature point pairs. The projection of the contour on the main axis is calculated to obtain multiple width feature point pairs. Then, the matched feature point pairs are matched by polar constraint, and the three-dimensional coordinates of the feature points are calculated. Finally, the distance between the length feature points is taken as the length of the crystal, and the distance between the multiple width feature points is voted to select the width feature point distance with the highest occurrence frequency as the crystal width, thereby completing the measurement of the length and width of the crystal.

[0046] The matching and measurement method in the third step can also be realized by the following method:

[0047] (1) Crystal image matching:

[0048] Firstly, the image is corrected by the correction method in step two, and then the crystal segmentation (Segment) and oriented bounding boxes object detection (Oriented Bounding Boxes Object Detection) are performed by using a deep learning algorithm (YOLOv8). The crystal image is quickly matched by one-dimensional search. Then, the part with straight line feature on the crystal contour is fitted by using Bézier curve segmentation, and then the matching is performed on the y coordinates and the angle direction of the line segment between the left and right images. The matched contour line segment is simplified as a feature point with weight, and the three-dimensional reconstruction is realized by using the above reconstruction algorithm.

[0049] (2) Virtual imaging plane establishment:

[0050] According to the three-dimensional length of the feature contour represented by the feature point, the weight of the feature point is defined as

[0051] wherein and represent the weights of the starting point, the middle point and the end point of the i-th feature contour, respectively. and represent the starting point and the end point of the i-th three-dimensional feature contour.

[0052] According to the three-dimensional feature points after weight allocation, the virtual three-dimensional plane where the maximum projection plane of the crystal is located is fitted as the imaging plane of the virtual camera. The virtual plane can be determined by minimizing the following equation to obtain the spatial plane parameters a p , b p , c p and d p .

[0053] wherein (x i , y i , zi ) is the coordinate of the i-th three-dimensional planar point, w i is the weight of the three-dimensional planar point, and n is the number of the points. The normal vector of the fitted virtual plane can be represented as

[0054] (3) Virtual camera projection imaging:

[0055] The projection of a two-dimensional image from the camera imaging plane to the virtual imaging plane is a rotation and translation transformation process of the imaging plane in space. The conversion between the camera and the virtual imaging plane is represented as

[0056] where p v and p c are the coordinates of a three-dimensional point in the camera plane and the virtual plane in space, respectively. Rv is the rotation matrix between the normal vector of the virtual plane and the normal vector of the binocular camera imaging plane , and t v is the translation matrix between the two planes. In order to eliminate the perspective distortion and maintain the geometric shape and position of each snapshot crystal, the present application adopts an orthogonal projection transformation. Considering that p c is a point on the imaging plane, Z c = 0, the transformation from the camera imaging plane to the virtual imaging plane is constructed by eliminating the Z coordinate on both sides of equation (2-4),

[0057] where H v is the homography matrix, and λ is the scaling factor, which can be set to 1 since the projection scale is invariant; and are the four elements of R v ; and are the elements of the translation vector t v .

[0058] Let the intrinsic parameters of the virtual camera be the same as those of the real camera, so the projection transformation from the camera imaging plane to the virtual imaging plane can be represented by the homography matrix H v ,

[0059] wherein and x c = [u c , v c , 1] T represent the homogeneous pixel coordinates of points on the virtual plane and the camera plane, respectively.

[0060] (4) Reconstruction quality evaluation:

[0061] First, the homography transformation is used to project the left and right camera images to the same virtual camera imaging plane using equation (2-6). The real projection from the spatial object to the camera imaging plane is a three-dimensional orthogonal projection, and the virtual plane coincides with the maximum projection plane of the crystal. Therefore, the virtual images obtained by the homography transformation of the left and right camera images to the same virtual plane should be consistent.

[0062] The homography-transformed left and right camera images are overlapped and compared. In this way, the accuracy of the three-dimensional reconstruction is verified by the degree of overlap of the images. If the left and right images have a high degree of overlap on the virtual plane, it can be considered that the reconstruction result is more accurate.

[0063] The overlap ratio of the images projected onto the virtual plane by the left and right camera images can be calculated as follows

[0064] wherein, are the virtual images of the crystal recovered from the left and right camera images, respectively.

[0065] (5) Statistical measurement based on virtual images:

[0066] The crystal virtual image segmentation result is rotated so that the crystal principal axis is perpendicular to the horizontal direction of the image. The rotated segmentation image is decomposed and projected in the horizontal and vertical directions for analysis,

[0067] wherein P U is the projection data along the u-axis one-dimensional length, N V is the total number of rows of the rotated crystal image; P V is the 1-D width projection data along the v-axis, and N U is the total number of columns of the rotated crystal image.

[0068] The average, median or maximum value and area of the length and width of each crystal image are respectively statistically analyzed as follows L mean = mean(P V ) / m (2-10) L median = median(P V ) / m (2-11) L max = max(P V ) / m (2-12) W mean = mean(P U ) / m (2-13) W median = median(P U ) / m (2-14) W max= max (P U ) / m (2-15)

[0069] In the formula, m is the camera magnification obtained by camera calibration.

[0070] According to different crystal shape features, suitable statistical parameters of the size of the crystal are selected, and the three-dimensional size of the crystal is measured.

[0071] The application can realize in-situ calibration and reconstruction of the telecentric stereovision system, and realize three-dimensional length and width measurement through the reconstructed crystal feature points, so that the measurement error caused by refraction of light in the cross-medium optical system can be effectively avoided, and the growth state of the crystal in the crystallization process can be more accurately analyzed. The method has strong operability, low experience and technical requirements, can achieve the effect of fast system in-situ calibration and automatic measurement of three-dimensional size of the crystal, and is convenient for actual industrial application. BRIEF DESCRIPTION OF DRAWINGS

[0072] Fig. 1 is a schematic diagram of the application scene and the calibration rod, Fig. 1(a) is a schematic diagram of the whole, Fig. 1(b) is a schematic diagram of the structure of the calibration rod, and Fig. 1(c) is a sectional view of the calibration rod;

[0073] Fig. 2 is a schematic diagram of a telecentric imaging model of the application;

[0074] Fig. 3 is a schematic diagram of a binocular telecentric stereoscopic imaging model of the application;

[0075] Fig. 4 is a schematic diagram of a epipolar correction model of the application;

[0076] Fig. 5 is a schematic diagram of epipolar correction, Fig. 5(a) is an original collected crystal image, and Fig. 5(b) is a corrected crystal image;

[0077] Fig. 6 is a schematic diagram of width feature point acquisition of crystal main shaft projection of the application, wherein Fig. 6(a) is a schematic diagram of crystal profile projection on the main shaft and width feature point selection, and Fig. 6(b) is a schematic diagram of extracted crystal length feature matching point and width feature matching point;

[0078] Fig. 7 is a schematic diagram of a light ray intersection reconstruction model of the application;

[0079] Fig. 8 is a schematic diagram of crystal length and width size calculation of the application, Fig. 8(a) is a three-dimensional reconstructed crystal length and width feature point, and Fig. 8(b) is a measurement result of the width of the crystal;

[0080] Fig. 9 is a schematic diagram of epipolar correction, Fig. 9(a) is an original collected crystal image, and Fig. 9(b) is a corrected crystal image and a crystal detection result;

[0081] Figure 10 is a schematic diagram of the feature profile segment of the crystal matching and the corresponding feature point pairs of the present application;

[0082] Figure 11 is a schematic diagram of the matching point three-dimensional reconstruction and virtual imaging plane fitting of the present application;

[0083] Figure 12 is a schematic diagram of the virtual imaging process from the binocular camera to the virtual camera of the present application, Figure 12(a) is the virtual imaging result corresponding to the image pair captured by the left camera; Figure 12(b) is the virtual imaging result corresponding to the image pair captured by the right camera;

[0084] Figure 13 is a schematic diagram of the crystal reconstruction result consistency evaluation of the present application, the left figure is the superimposed synthetic view of the left camera and the right camera virtual image; the right figure is the synthetic view of the binary segmentation result of the left camera and the right camera virtual image;

[0085] Figure 14 is a diagram of the three-dimensional size measurement method of the rod-shaped crystal of the present application, Figure 14(a) is a virtual imaging image of the crystal; Figure 14(b) is a virtual imaging image of the crystal after rotation; Figure 14(c) is a 1-D length projection along the v-axis; Figure 14(d) is a virtual binary image after rotation; Figure 14(e) is a one-dimensional width projection along the u-axis;

[0086] Figure 15 is the distribution result of the crystal virtual imaging plane recovered by the measurement method of the present application in space. DETAILED DESCRIPTION

[0087] In order to better understand the technical solutions of the present application, the embodiments of the present application will be described in detail below with reference to the accompanying drawings.

[0088] The method first synchronously collects images in real time by placing two telecentric cameras at different angles outside the reaction kettle, inserts a calibration rod into the reaction kettle, and rotates the calibration rod clockwise in order to calibrate the binocular system in situ. Second, a stereovision epipolar rectification method expressed by a 3*3 matrix is established to simplify the two-dimensional search of matching points to one-dimensional search. Finally, the stereoscopic imaging model and the stereovision epipolar rectification method are used to match each pair of crystal images, and the three-dimensional length-width feature points are reconstructed by the ray intersection method based on analytical solution, and the three-dimensional size is quantitatively evaluated by multiple sets of data.

[0089] The method is based on a non-contact industrial binocular telecentric imaging system, as shown in Figure 1(a). An example adopts a 1L crystallization glass reactor, the reactor is provided with a 4-blade stirring paddle connected with a computer, and 500 mL of crystal water solution is injected into the reactor. A non-contact binocular imaging device is arranged outside the reactor, which includes two high-speed high-resolution industrial monochrome cameras connected with the computer, and two LED point light sources are provided to provide backlight illumination, and the cameras are located on the two sides of the reactor. In the calibration process, the calibration rod is inserted into the imaging target area through the device hole at the top of the reactor, and the system is calibrated by rotating the calibration rod. The design of the calibration rod is shown in Figures 1(b) and (c). The light emitted by the LED at the top of the calibration rod is guided to the reflection prism at the bottom of the calibration rod through the optical fiber, and the reflection prism reflects the light vertically to the back of the ceramic calibration plate engraved with a chessboard pattern. The reflected light provides illumination for the shooting of the two telecentric cameras, and the shooting of the calibration image set is completed. In the measurement process, the binocular images (left view and right view) of the crystallization process are synchronously collected by the two telecentric cameras, and after the epipolar constraint correction and matching, the crystals in any posture in the solution are reconstructed to realize the three-dimensional size measurement of the length and width of the crystals based on binocular vision.

[0090] The method specifically comprises the following steps:

[0091] First step, calibration:

[0092] In order to eliminate the influence of refraction caused by the propagation of light in the cross-medium optical system to some extent, the present application proposes a new rotation calibration method. The calibration of the stereoscopic telecentric system can be simply completed by manually rotating the calibration plate around its axis to image in situ. A backlit calibration rod composed of a light-emitting diode (LED), a ceramic calibration plate covered with a chessboard pattern, a reflection prism and an optical fiber is designed for in-situ calibration into target spaces such as reactors.

[0093] (1) Monocular camera parameter calibration:

[0094] According to the telecentric camera imaging model shown in Figure 2, the world coordinates of the plane calibration plate can be related to the captured image coordinates through the homography matrix H, and the relationship is as follows:

[0095] In order to facilitate analysis, Euler angles are introduced to describe the rotation matrix R(α,β,θ), and α,β,θ are the rotation angles around the X, Y and Z axes of the coordinate system, respectively. The elements of the corresponding rotation matrix are described as:

[0096] By equating the corresponding elements of the 3x3 matrix in (2), the following equations can be easily established:

[0097] where, where a, b and q are rotation angles in the range of (-p / 2, p / 2). In addition, the effective amplification factor is defined as a positive value. The solution of the above equation is:

[0098] where c=(h 11 h 22 -h 12 h 21 ) 2 The translation matrix T(t x ,t y ) is obtained as:

[0099] Therefore, the amplification factor m and the Euler angle q can be uniquely determined. However, the Euler angles a and b have two possible signs. Therefore, there may be two ambiguous poses related to the world coordinate system in the single camera coordinate system, and there may be four ambiguous poses in the binocular camera system. Recovering the true Euler angles a and b is the key to the calibration method of the system.

[0100] Due to the orthographic projection characteristics of the above-mentioned telecentric camera, the camera pose recovered from the two-dimensional pattern is ambiguous. The existing method relies on a three-dimensional calibration pattern with limited accuracy, or makes full use of external displacement, but due to the existence of the reactor wall, it cannot establish a physical connection between the imaged target and the stereo vision system, and cannot be used for on-site calibration of the stereo vision system. In order to overcome the above problems, the present application proposes a practical rotating template calibration method.

[0101] In the calibration process, the calibration pattern is rotated around the Y w axis at a specific angle sequence (clockwise or counterclockwise) so that the camera can obtain an ordered calibration shot image. Therefore, the change of the placement angle of the calibration pattern is observed to be sequentially increased (or decreased). At the same time, the angle of the rotation Euler angle b of the camera coordinate system relative to the world coordinate system is also monotonously increased (or decreased). For ease of analysis, it is assumed that the calibration pattern is rotated clockwise around the Y w axis.

[0102] For each image captured during the calibration process, the extrinsic matrix of the camera is changing. By taking the positive sign of b in (25), a set of values of the rotation matrix R k (a k , b k , q k ), translation vector and amplification factor m k of the kth calibration image can be obtained. The calculation formula of the amplification factor m is

[0103] Considering that the actual rotation angle β is a monotonically increasing sequence, the absolute value of β can present a local minimum point or remain monotonous. Calculating the absolute value of each function can check whether there is a local minimum point. If there are multiple local minimum points or there are local maximum points, the calibration condition is not met. However, if the specified condition is met, it can be inferred that the extreme value point will appear near the minimum value of the discrete sequence, thereby causing

[0104] Through several discrete, the extreme value point can be determined to be near the minimum value, such as the left or right side of the minimum value, as shown in FIG. 5(a). Since the rotation angle β is monotonically increasing, in addition to the sign of the minimum value point, the sign is recovered by the following equation:

[0105] The rotation angle α is recovered by the following equation:

[0106] (2) Calibration of stereo vision system parameters:

[0107] The essence of binocular vision calibration is to obtain the relative position relationship between the two telecentric cameras, as shown in FIG. 3. The left and right camera coordinate position relationship is described by a rotation matrix R C and a translation vector T C , that is:

[0108] where R L and T L , R R and T R are the rotation matrix and translation vector related to the world coordinate system and the left and right camera coordinate system, respectively.

[0109] For each pair of calibrated images, the external matrix of the left and right cameras changes, while the relative rotation matrix R C between the left and right cameras remains unchanged. The Euler angle representation of R C obtained from the kth pair of calibration images is:

[0110] where α C , β C , θ C are the rotation angles around the X, Y, and Z axes of the left camera coordinate system. The rotation Euler angles and represent the external parameters R L and R R in the left and right camera coordinate systems. The superscript k indicates that the data comes from the kth set of calibration images.

[0111] The Euler angles α C , β C , θ C obtained from N pairs of different calibration images are optimized

[0112] where N k is the number of variables in the set of calibration images.

[0113] The calculated rotation matrix in Euler angle representation is converted to matrix representation R C and T C is obtained according to (11), completing the stereo calibration.

[0114] Second step, correction:

[0115] In order to reduce the computational complexity of the matching process, the present application proposes a simplified telecentric stereo vision epipolar correction method, which completes the re-projection correction process through a 3*3 homography matrix between the images before and after correction. As shown in Figure 4, the coordinate systems of the corrected left and right camera pixels are aligned along the common y-axis.

[0116] The projections of point P from the world coordinate system to the left and right camera pixel coordinate systems are respectively represented as

[0117] Considering that the magnification factors of the left and right cameras after correction should be equal to ensure the alignment of the y-axes of the two cameras in the pixel coordinate system, the magnification factors of the two cameras are determined by finding the maximum magnification factor from them:

[0118] Therefore, the projections of point P from the world coordinate system to the corrected left and right camera pixel coordinate systems are respectively determined as

[0119] As shown in Figure 4 represent the X-, Y-, and Z-axis direction vectors in the left and right camera coordinate systems before and after correction, i.e. RL, RR, , respectively, which are the row vectors of the matrix. It should be noted that after correction, the y-axis of the left camera coordinate system is parallel to the y-axis of the right camera, and the z-axis representing the light direction remains unchanged. Therefore, the common axis is perpendicular to the optical axes of the left and right cameras. The direction vectors of the ideal corrected coordinate systems of the left and right cameras can be represented as:

[0120] In order to align the y-axes of the corrected left and right camera coordinate systems, the new translation vector is calculated as

[0121] Combining equation (34) and equation (36), the conversion from the original image coordinates to the corrected coordinates is established:

[0122] Third step, matching and reconstruction measurement:

[0123] (1) Matching

[0124] The original image collected is shown in Fig. 5(a). First, the image is corrected by the correction method in step two, and the crystal image is matched by one-dimensional search. The corrected image and the matched crystal are shown in Fig. 5(b). Then, the length feature point pairs are obtained by the distance from the contour to the center point of the contour, and the crystal main axis is obtained by connecting the length feature point pairs, as shown by the dashed line connected by the star-shaped markers in Fig. 6(b).

[0125] It should be noted that during the crystallization process, small crystals will adhere to the edges of large crystals, as shown in Fig. 6(b). In order to avoid mis-extraction when measuring the crystal width, the present application proposes a method for selecting effective width contour points for measurement as follows:

[0126] The measurement of the width of each particle refers to the calculation of the Euclidean distance between two contour points, and the connecting line between the contour points is perpendicular to the main axis of the measured particle. In order to ensure measurement accuracy, all contour point pairs perpendicular to the connecting line of the main axis are selected as candidate marker points for measuring the width of the particle. For this purpose, the particle contour is divided into upper and lower halves according to the main axis. The distance from each half-contour point to the main axis is calculated as

[0127] Considering that the crystal width distance should reflect the dominant distance of the above-mentioned contour point pairs, rather than their average distance, it is suggested to judge whether the distance deviation is less than half of the average absolute deviation defined respectively:

[0128] where N p is the number of points on each half-contour. If yes, the corresponding point is identified as a candidate point, and vice versa, as shown by the candidate points on the left and right contours in Fig. 6(a).

[0129] In order to match these candidate points, the contour coordinates in the left and right images are projected onto their main axes respectively by the following formula, and the normalized contour distance to the main axis is shown in Fig. 6(a).

[0130] where x u and x d are the x-coordinate values of p u and p d respectively.

[0131] The candidate points on the upper and lower contours of the left image with the same projection coordinates are regarded as a pair of particle width marker points, as shown by the circular point markers on the horizontal coordinate axis in Fig. 6(a). Similarly, the corresponding matching points in the right image also have the same projection coordinates. That is, the four candidate points with the same projection coordinates on the upper and lower contours of the left and right images constitute a set of particle width marker points.

[0132] (2) Reconstruction

[0133] For the telecentric stereo vision system, the application proposes a new telecentric stereo vision three-dimensional reconstruction method based on analytical solution of ray intersection as shown in Figure 7, which is as follows:

[0134] First, the position equation of each ray passing through the matching point pair on the left and right physical image planes is constructed respectively.

[0135] Where (x l ,y l ,z l ) represents the coordinates of any point on the left image plane ray, represents the coordinates of the matching point on the left physical image plane, represents the direction unit vector of the left camera optical axis, τ L is the ray coefficient of the left camera, (x r ,y r ,z r ) represents the coordinates of any point on the right image plane ray, represents the coordinates of the matching point on the right physical image plane, represents the direction unit vector of the right camera optical axis, τ R is the ray coefficient of the right camera. If the position equation of the calculated point on the physical image plane and the above ray is accurate, there is only one intersection point between the left and right rays, i.e. the reconstruction point on the object surface. However, due to the existence of noise, the left and right rays may not intersect. In order to find the best intersection point, the two points closest to the left and right rays are calculated respectively:

[0136] The distance between any two points on the left and right rays is represented as:

[0137] Substitute (45) and (46) into (47) to obtain a quadratic equation about τ L and τ R :

[0138] The minimum value of f(τ L , τ R ) is obtained by letting , and the values of τ L and τ R are determined as:

[0139] Where

[0140] Correspondingly, the coordinates of the reconstructed points are given by:

[0141] (3) Measurement

[0142] First, the three-dimensional coordinates of the matched feature points are calculated using the above three-dimensional reconstruction method, and the reconstructed feature points are shown in Fig. 8(a). Then, the Euclidean distance between the length feature points is taken as the length of the crystal. The Euclidean distance between multiple width feature points is voted, and the width feature point distance with the highest occurrence frequency is selected as the width of the crystal, i.e., the width median, as shown in Fig. 8(b), completing the measurement of the length and width of the crystal.

[0143] The three-dimensional coordinates of the length marker points P u , P d and the width marker point set P f , P l are calculated using the above three-dimensional reconstruction method. Therefore, the length and width of each crystal can be calculated as: L = |P u -P d | (52) W med = median(|P f -P l |) (53) W mean = mean(|P f -P l |) (54)

[0144] Step three can also be achieved by the following matching and reconstruction measurement method:

[0145] (1) Matching:

[0146] The original image collected is shown in Fig. 9(a). First, the image is corrected by the correction method in step two, and then the crystal segmentation (Segment) and oriented bounding boxes object detection (Oriented Bounding Boxes Object Detection) are performed using the deep learning algorithm (YOLOv8). The crystal image is quickly matched through one-dimensional search, and the corrected image and the matched crystal are shown in Fig. 9(b). Then, the part with straight line feature on the crystal contour is fitted by using the Bézier curve segmentation, and then its y coordinate and the line segment angle direction are matched, as shown by the directed line segment in Fig. 10, and simplified as a feature point with weight, as shown by the intersection point in Fig. 10.

[0147] (2) Reconstruction:

[0148] Referring to the above analytical solution-based telecentric stereo vision three-dimensional reconstruction method of light intersection.

[0149] (3) Establish a virtual imaging plane:

[0150] To convert the three-dimensional spatial pose of the crystal image into a virtual two-dimensional image for analysis, the present application proposes a virtual stereoscopic imaging method based on homographic projection transformation. Inspired by the offline microscope imaging scene, the virtual camera can capture the vertical projection of the maximum crystal surface, and by reconstructing the maximum projection surface, the length, width and area of the captured crystal can be accurately measured. In addition, a self-supervised verification strategy based on the virtual imaging consistency between left and right camera images is proposed to ensure the accuracy of the three-dimensional reconstruction of the crystal image. Specifically as follows:

[0151] The crystal feature contour points matched in (1) are reconstructed using the method in (2), and the three-dimensional coordinates of the start point, end point and midpoint of the matched crystal feature contour can be obtained, as shown by the black cross marks in FIG. 11.

[0152] Based on the above representative 3-D points, a virtual imaging plane is constructed for fitting. According to the three-dimensional length of the feature contour represented by the feature points, the feature point weight is defined as

[0153] where and represent the weights assigned to the start point, midpoint and end point of the i-th feature contour, respectively. and represent the start point and end point of the i-th feature contour.

[0154] The virtual three-dimensional imaging plane equation is a P x+b P y+c P z+d P =0 (2-18)

[0155] where (x, y, z) is a point on the fitted plane; a p ,b p ,c p and d p are equation coefficients.

[0156] After assigning the weights, the plane fitting is converted into a least squares fitting problem, which can be determined by minimizing the following equation for coefficients a p ,b p ,c p and d p .

[0157] where (x i ,y i ,zi ) is the coordinate of the i-th three-dimensional planar point, w i is the weight of the three-dimensional planar point, and n is the number of the points.

[0158] Take the partial derivatives of F with respect to a p ,b p ,c p and d p , and set them equal to zero, we can get

[0159] The pseudo-inverse of the above coefficient matrix is obtained by singular value decomposition (SVD), so as to determine the coefficients of a p ,b p ,c p and d p of the best fitting plane.

[0160] The above four coefficients are normalized by using the following formula to calculate the normal vector of the fitting plane as shown in Figure 11.

[0161] (4) Virtual camera projection imaging:

[0162] The projection of a two-dimensional image from the camera imaging plane to the virtual imaging plane is a rotation and translation transformation process of the imaging plane in space. The conversion between the camera and the virtual imaging plane is represented as

[0163] where p v and p c are the coordinates of a three-dimensional point in the camera plane and the virtual plane in space, respectively.

[0164] The normal vector of the virtual plane is determined by equation (2-21) The normal vector of the binocular camera imaging plane can be obtained by system calibration in the first step The rotation matrix between the two planes is obtained by the following formula.

[0165] The position of the first matching point in the camera plane image is designated as the center of the virtual imaging plane, denoted as p1, and the translation vector t v is

[0166] In order to eliminate perspective distortion and maintain the geometric shape and position of each snapshot crystal, the present application adopts an orthogonal projection transformation. Considering that p c is a point on the imaging plane, Zc=0, by eliminating the Z coordinates on both sides of (2-22) to construct the transformation from the camera imaging plane to the virtual imaging plane,

[0167] where H v is a homography matrix, and λ is a scaling factor, which can be set to 1 due to the scale invariance of the projection; and are the four elements of R v ; and are the elements of translation vector t v . Since the intrinsic parameters of the virtual camera are the same as the real camera, the projection transformation from the camera plane to the virtual plane can be represented by a homography matrix H v ,

[0168] where, and x c = [u c , v c , 1] T represent the homogeneous pixel coordinates of points on the virtual plane and the camera plane, respectively.

[0169] Fig. 12 shows the virtual imaging from the left and right cameras to the virtual camera. In Fig. 12(a), the right camera represents the position and orientation of the real left camera in the world coordinate system. The image taken by the left camera is represented by the image in the upper right corner. The image in three-dimensional space represents the projection of the left camera image on the virtual plane. The left camera represents the position and orientation of the virtual camera perpendicular to the crystal surface in the world coordinate system. The image taken by the virtual camera is represented by the image in the lower right corner.

[0170] (5) Reconstruction result consistency evaluation:

[0171] The homography transformation in the above (4) provides the mapping from the camera imaging plane to the virtual imaging plane. While the real projection from the spatial object to the camera imaging plane belongs to the three-dimensional orthogonal projection. In other words, the projection from the camera plane to the virtual plane can only accurately recover the part of the space that overlaps with the virtual imaging plane. The virtual images obtained by the above homography transformation of the left and right camera images to the same virtual plane should be consistent. In order to evaluate the accuracy of the three-dimensional reconstruction result, the present application proposes a method of overlapping comparison of the orthographic images obtained by the left and right cameras.

[0172] According to equation (2-27), the projection of the crystal image recovered from the left camera image on the virtual plane is

[0173] Similarly, the projection calculation of the right camera is

[0174] Therefore, the corrected projection of the crystal image on the virtual plane can be represented by

[0175] where, and are the images projected onto the virtual plane from the left and right cameras respectively.

[0176] The overlap ratio of the images of the left and right camera images projected onto the virtual plane can be calculated as follows

[0177] If the overlap ratio O V is less than a threshold value of 0.9, it is considered that the virtual projection of the three-dimensional reconstruction is incorrect, and therefore the projection result is invalid.

[0178] Figure 12(b) shows the virtual imaging from the right camera to the virtual camera. The camera on the right represents the right camera, and the image in the top right corner is captured by the right camera. The camera on the left represents the virtual camera, and the image in the bottom right corner represents the projection of the virtual camera. Note that in Figures 12(a) and (b), the position of the virtual camera is unique and the same.

[0179] Figure 13 is the combined result of the virtual camera virtual imaging constructed by the left and right cameras. It can be seen that there is a slight non-overlap on the left edge of the crystal image, which is caused by the error in the thickness of the crystal shape. Overall, the consistency is good, with an overlap ratio of 96.65%, and further three-dimensional size measurement of the crystal grains can be performed.

[0180] (6) Statistical measurement based on virtual images:

[0181] Since the outline of the rod-shaped crystal image obtained in the actual crystallization process is not strictly rectangular, the present application proposes a method for statistical measurement of the length and width of a single crystal, i.e., by in-situ binocular image segmentation of the crystallization process, and reconstructing the virtual imaging image of each crystal to measure the average, median or maximum value of its length and width, and the maximum projection area.

[0182] The rotation angle around the z-axis between the virtual left (or right) plane and the left (or right) camera coordinate system is defined as follows

[0183] Therefore, the angle between the principal axis of the rod-shaped crystal image projected onto the virtual plane and the vertical direction of the virtual image is represented as

[0184] where are the principal axis directions of the crystals detected by the YOLOv8-Obb network respectively.

[0185] The two-dimensional coordinate transformation of the virtual imaging image rotated by θ Z can be represented as

[0186] The binary image I Vrotating the crystal main direction to align with the v-axis of the binary image I V and then performing orthogonal decomposition projection on the rotated binary image, the sum of the pixel values along the rows and columns provides the basic information of the size distribution and crystal shape in the crystal binary image.

[0187] The orthogonal decomposition projections of the rotated image along the u-axis and the v-axis are

[0188] where P U is the one-dimensional length projection data along the u-axis, N V is the total number of rows of the rotated crystal image; P V is the one-dimensional width projection data along the v-axis, N U is the total number of columns of the rotated crystal image.

[0189] Therefore, the average, median or maximum value of the length and width of each crystal image and their areas are respectively calculated as follows L mean = mean(P V ) / m (2-37) L median = median(P V ) / m (2-38) L max = max(P V ) / m (2-39) W mean = mean(P U ) / m (2-40) W median = median(P U ) / m (2-41) W max = max(P U ) / m (2-42)

[0190] where m is the magnification of the camera obtained by camera calibration.

[0191] The present application proposes to take the maximum value in the length projection data P U as the crystal length, which represents the longest straight line distance between any two contour points parallel to the crystal main axis, i.e. the calculation result of formula (2-39). At the same time, the crystal width is characterized by the median value in the width projection data P V , which represents the average distance of multiple repeated straight lines connecting any two contour points perpendicular to the crystal main axis, i.e. the calculation result of formula (2-41).

[0192] Figure 14 illustrates the measurement of rod-like crystals based on orthogonal decomposition projection. The angle θ between the main axis of the crystal and the v-axis of the virtual image Z determined from (2-32) and (2-33) as shown in Figure 14(a). Rotating the image (2-34) will produce the rotated virtual image in Figure 14(b). The rotated binary image is shown in Figure 14(d), and the one-dimensional data P representing the length and width of the crystal calculated from (2-35), (2-36) U and P V are shown in Figures 14(c) and (e), respectively. In Figure 14(c), the dashed line represents the maximum length.

[0193] Figure 15 shows the virtual imaging results of each crystal in a set of detection images (as shown in Figure 9) coming from the left camera and the spatial position of the crystal, illustrating the basic principle and rationality of the reconstruction method based on virtual plane proposed by the present application.

Claims

1. A method for real-time detection of the size distribution of crystal population in a crystallizer by using binocular telecentric cameras, characterized in that: The method is realized by using a non-contact high-resolution backlit calibration rod, which ensures that two telecentric cameras synchronously collect images of the crystal solution in the reactor, thereby measuring the three-dimensional size of the crystal; the backlit calibration rod is composed of a light-emitting diode (LED), a ceramic calibration plate printed with a chessboard pattern, a reflecting prism, an optical fiber, and a stainless steel sheath; the light emitted by the LED at the top of the calibration rod is guided to the reflecting prism at the bottom of the calibration rod through the optical fiber, and the reflecting prism reflects the light vertically to the back of the ceramic calibration plate with a chessboard pattern; the reflected light provides illumination for the shooting of the two telecentric cameras; the calibration rod can be inserted into the glass crystallizer and rotated clockwise or counterclockwise along the pivot axis to perform on-site calibration; the backlit calibration rod can be assembled into different sizes to adapt to the on-site calibration of crystallizers of different volumes; By placing two telecentric cameras at different angles outside the reactor to synchronously collect images, the images are called binocular images, which include left and right views; a two-step in-situ stereo imaging calibration model is established for the pose of the binocular camera, including calculating the camera parameter closed loop to obtain the initial value of the model parameter and subsequently using an optimization algorithm to seek the optimal camera model parameter; for the collected binocular images, a crystal image matching analysis method is proposed, including image preprocessing, interest point detection, feature matching, and mismatch removal; after detecting the length and width key corner points of the crystal using the crystal contour feature, the three-dimensional reconstruction of the length and width corner points of each matched crystal is performed through the calibrated stereo imaging model, thereby quantitatively evaluating the three-dimensional size. 2.The method according to claim 1, characterized in that: Step 1, stereo imaging calibration The telecentric camera imaging model, i.e. the projection of a point P(xw,yw,zw) in the world coordinate system to the pixel coordinate system P(u,v) is represented as where m is the effective magnification of the telecentric lens, and R and T are the rotation and translation matrices between the two coordinate systems, respectively; (1) Monocular camera parameter calibration: The world coordinates of the planar calibration plate can be related to the captured image coordinates by a homography matrix H, the relationship being as follows: For the convenience of analysis, Euler angles are introduced to describe the rotation matrix R(a, β, θ), a, β, θ are the rotation angles around the X-axis, Y-axis and Z-axis of the coordinate system, respectively; the elements of the corresponding rotation matrix are described as: By equating the corresponding elements of the 3x3 matrix in (2), the following equations can be easily established: wherein, wherein a, β and θ are rotation angles in the range (-π / 2, π / 2); further, the effective magnification factor is defined as a positive value; the solution of the above equation is: wherein c = (h 11 h 22 -h 12 h 21 ) 2 ; the translation matrix T(t x ,t y ) is: Therefore, the magnification factor m and the Euler angle θ can be uniquely determined; however, the Euler angles α and β have two possible signs, as shown in (5); therefore, there may be two ambiguous poses related to the world coordinate system in the single camera coordinate system, and there may be four ambiguous poses in the binocular camera system; recovering the true Euler angles α and β is the key to the calibration method of the system; Due to the orthographic projection characteristics of the above-mentioned telecentric camera, the camera pose recovered from the two-dimensional pattern has ambiguity; existing methods rely on three-dimensional calibration patterns with limited accuracy, or make full use of external displacement, but due to the inability to establish a physical connection between the imaging target and the stereo vision system, it cannot be used for on-site calibration of the stereo vision system; in order to overcome the above problems, a practical rotating template calibration method is proposed; During the calibration process, the calibration pattern rotates around the Y w axis at a specific angle sequence, so that the camera can obtain an ordered calibration shot image; therefore, the change of the placement angle of the calibration pattern is observed to be sequentially increased or decreased; at the same time, the angle of the rotation Euler angle β of the camera coordinate system relative to the world coordinate system is also monotonously increased or decreased; for the convenience of analysis, it is assumed that the calibration mode is clockwise rotation around the Y w axis; Considering that the actual rotation angle β is a monotonically increasing sequence, the absolute value of β can present a local minimum point or remain monotonous; calculating the absolute value of each function can check whether there is a local minimum point; if there are multiple local minimum points or there is a local maximum point, the calibration condition is not met; however, if the prescribed condition is met, it can be inferred that the extreme value point will appear near the minimum value of the discrete sequence, thereby causing By several discrete, extreme points can be determined near the minimum, such as the left or right of the minimum; since the rotation angle β is monotonically increasing, in addition to the sign of the minimum point, the following equation is used to restore the sign: The rotation angle a is recovered by the following equation: (2) Stereo vision system parameter calibration: The essence of binocular vision calibration is to obtain the relative position relationship between two telecentric cameras; the coordinate position relationship of left and right cameras is described by a rotation matrix R C and a translation vector T C , that is: where R L and T L are the rotation matrix and translation vector associated with the world coordinate system and the left and right camera coordinate system, respectively. R and T R are the rotation matrix and translation vector associated with the world coordinate system and the left and right camera coordinate system, respectively. For each pair of calibrated images, the extrinsic matrix of the left and right cameras changes, while the relative rotation matrix R C between the left and right cameras remains unchanged; the Euler angle representation of R C from the kth pair of calibrated images is denoted as where α C , β C , θ C are the rotation angles around the X, Y and Z axes of the left camera coordinate system; the rotation Euler angles and represent the extrinsic parameters R L and R R in the left and right camera coordinate systems; the superscript k indicates that the data comes from the kth set of calibration images; Euler angles a C , b C , g C are optimized for N pairs of different calibration images where N k is the number of variables in the set of calibration images; Convert the calculated rotation matrix Euler angle representation into matrix representation R C and obtain T according to (11) C , complete the stereo calibration work; Step 2, simplified stereo rectification Epipolar rectification plays an important role in accurate stereo matching, which simplifies the two-dimensional search of matching points to one-dimensional search; in the context of telecentric stereo vision, the existing method uses a 3x4 matrix for epipolar rectification; in fact, from the geometric and physical point of view, the rectification parameters can be further simplified to a 3x3 homography matrix, which is more convenient for computer-aided image analysis; the specific method is as follows: The projection of point P from the world coordinate system to the left camera pixel coordinate system p L and p R are represented as The projection of point P from the world coordinate system to the corrected left and right camera pixel coordinate systems is determined as: According to the formulas (14) and (15), the conversion from the original image coordinates to the corrected coordinates is established: Third step, reconstruction and measurement (1) Three-dimensional reconstruction: An improved ray intersection three-dimensional reconstruction method based on analytical solution is proposed for telecentric stereo vision system; First, the left and right physical image planes are constructed by matching pairs of points on the left and right physical image planes the position equation of each ray of the light beam; where (x l ,y l ,z l ) represents the coordinates of an arbitrary point on the left image plane ray, representing the coordinates of the matching point on the left physical image plane, denotes the direction unit vector of the left camera optical axis, τ L is the optical ray coefficient of the left camera, (x r ,y r ,z r ) denotes the coordinates of an arbitrary point on the right image plane optical ray, representing the coordinates of the matching point on the right physical image plane, denotes the direction unit vector of the right camera optical axis, τ R is the optical ray coefficient of the right camera; if the position equation of the calculated point on the physical image plane and the above-mentioned ray is accurate, there is only one intersection point between the left and right rays, i.e. the reconstructed point on the object surface; however, due to the existence of noise, the left and right rays may not intersect; in order to find the best intersection point, the two points closest to the left and right rays are calculated respectively, and the coefficient of the left ray when the nearest point position is obtained by calculating the analytical solution: wherein Correspondingly, the coordinates of the reconstruction points are found by the following equation: (2) Crystal image matching and measurement: First, the image is corrected by the correction method in step two, and the crystal image is matched quickly through one-dimensional search; then the length feature point pairs are obtained by the distance from the contour to the contour center point, and the crystal main axis is obtained by connecting the length feature point pairs; then, the feature point pairs are matched by the epipolar constraint, and the three-dimensional coordinates of the feature points are calculated; finally, the distance between the length feature points is taken as the length of the crystal, and the distance between the width feature points is voted to select the width feature point distance with the most occurrences as the crystal width, and the length and width of the crystal are measured.

3. The method according to claim 1, wherein the matching and measurement method of the third step is also realized by the following method: (1) Crystal image matching: First, the image is corrected by the correction method in step two, and then the crystal segmentation and rotation axis frame detection are performed by using the deep learning algorithm YOLOv8, and the crystal image is matched quickly through one-dimensional search; then the contour on the crystal is fitted by using the Bezier curve segmentation, and then the y coordinates and the angle direction of the line segment between the left and right images are matched, the matched contour segment is simplified as a feature point with weight, and the three-dimensional reconstruction is realized by using the above reconstruction algorithm; (2) Virtual imaging plane establishment: The starting point and the ending point of the i-th three-dimensional feature contour are represented by (xi, yi) and (xi+1, yi+1), respectively; According to the three-dimensional length of the feature contour represented by the feature point, the feature point weight is defined as wherein and respectively denote the weights assigned to the start, middle and end points of the i-th feature profile; and (3) Virtual camera projection imaging: According to the three-dimensional feature points after the distribution weight, a virtual three-dimensional plane where the maximum projection plane of the crystal is located is fitted as the imaging plane of the virtual camera, and the spatial plane parameters a p , b p , c p and d p of the virtual plane are determined by minimizing the following equation where (x i ,y i ,z i ) are the coordinates of the i-th three-dimensional planar point, w i is the weight of the three-dimensional planar point, and n is the number of points; the normal vector of the fitted virtual plane is represented as (4) Reconstruction quality evaluation: The projection of a two-dimensional image from the camera imaging plane to the virtual imaging plane is a process of rotation and translation transformation of the imaging plane in space; the conversion between the camera and the virtual imaging plane is represented as where p v and p c are the coordinates of the three-dimensional point in the space in the camera plane and in the virtual plane, respectively; Rv is the normal vector of the virtual plane Normal vector to binocular camera imaging plane the rotation matrix between t v is the translation matrix between the two planes; in order to remove the perspective distortion and preserve the geometry and position of each snapshot crystal, an orthographic projection transformation is adopted; considering p c is a point on the imaging plane, Z c = 0, by eliminating the Z coordinate on both sides of equation (2-4) to construct the transformation from the camera imaging plane to the virtual imaging plane, wherein H v is a homography matrix, and λ is a scale factor, which can be set to 1 due to the scale invariance of the projection; and is R v four elements of R and is the translation vector t v an element of ; Let the internal parameters of the virtual camera be identical to the real camera, so the projection transformation from the camera imaging plane to the virtual imaging plane uses a homography matrix H v denotes, In the formulae, and x c = [u c , v c , 1] T respectively denote the homogeneous pixel coordinates of a point on the virtual plane and the camera plane; First, the homography transformation is used to project the left and right camera images to the same virtual camera imaging plane according to formula (2-6); the real projection from the space object to the camera imaging plane belongs to three-dimensional orthogonal projection, and the virtual plane coincides with the maximum projection plane of the crystal, so the virtual images obtained by the homography transformation of the left and right camera images to the same virtual plane should be consistent; The left and right camera images after homography transformation are overlapped and compared; in this way, the accuracy of three-dimensional reconstruction is verified by the degree of image overlap; if the left and right images have high coincidence on the virtual plane, it can be considered that the reconstruction result is more accurate; The crystal virtual images recovered from the left and right camera images are respectively The overlap ratio of the images projected onto the image of the virtual plane by the left and right camera images is calculated as follows wherein, (5) Statistical measurement based on virtual images: The average, median or maximum value of the length and width of each crystal image and its area are respectively counted as follows rotating the crystal virtual image segmentation result to make the crystal main axis perpendicular to the horizontal direction of the image; and performing decomposition projection analysis on the segmented image in the horizontal and vertical directions, where P U is the projection data along the u-axis one-dimensional length, N V is the total number of rows of the crystal image after rotation; P V is the projection data along the v-axis 1-D width, N U is the total number of columns of the crystal rotation image; In the formula, m is the magnification of the camera obtained by camera calibration; L mean = mean(P V ) / m (2-10) L median = median(P V ) / m (2-11) L max = max(P V ) / m (2-12) W mean = mean(P U ) / m (2-13) W median = median(P U ) / m (2-14) W max = max(P U ) / m (2-15) ​ According to different crystal shape features, suitable statistical parameters of crystal size are selected to complete the measurement of three-dimensional size of the crystal.

Images