A multi-modal image coordinated head lesion multi-level automatic positioning method

Abstract

The application provides a multi-modal image cooperative head lesion multi-level automatic positioning method, which fuses rigid initial registration and non-rigid accurate registration algorithm based on thin plate spline algorithm, optimizes the thin plate spline algorithm from the local and global levels through the Delaunay triangulation method, balances the registration accuracy, the calculation efficiency and the robustness by adopting a multi-target constraint processing mechanism, simultaneously replaces the traditional dense scheme by light-weight 7 key position points to reduce interference, realizes accurate positioning by combining head posture geometric analysis and initial alignment, and takes into account the theoretical reliability and engineering practicability, thereby providing strong technical support for multi-modal image navigation; the application effectively improves the accuracy of the lesion position cutting point, reduces the dependence on artificial experience in the clinic, and realizes the positioning function of the high-precision and low-cost image navigation system by MRI image reconstruction and dimension reduction, face key position point selection, head posture preliminary alignment and fine registration between virtual and real images.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of medical image navigation system optimization, and in particular to a multi-modal image collaborative head lesion multi-level automatic positioning method. BACKGROUND

[0002] Surgical clinical operation is usually limited by the fact that the lesion is embedded in the patient's soft and hard tissue, causing a visual field blind area. Magnetic resonance imaging technology can effectively display lesion location, size and shape information, and is often used for preoperative diagnosis. Computer-assisted surgical navigation technology, which combines medical imaging technology, computer science, positioning and tracking technology, and robotics, has a history of several decades. Surgeons can analyze patient anatomy and provide precise surgical path guidance during intraoperative treatment with the help of this technology. In image navigation robot treatment systems, given that magnetic resonance imaging (MRI) and computed tomography (CT) medical images are often presented in the form of tomographic slices, it is difficult to fully represent the patient's surface anatomy, so three-dimensional reconstruction of medical images is required. Precise alignment of the reconstructed three-dimensional virtual image with the patient's real intraoperative anatomy, i.e. patient virtual-real image registration, is a core step in achieving the precise treatment function of the system, and its performance directly affects the accuracy and safety of the operation, and must meet the engineering indicators of high registration accuracy, real-time, low clinical invasiveness, and controllable system cost. This is a key problem that needs to be solved.

[0003] The binocular infrared stereoscopic directional technology currently widely used in image navigation systems can achieve three-dimensional positioning through parallax calculation in theory, but has limitations in multiple dimensions in actual clinical scenarios. From the hardware perspective, binocular infrared cameras require high-sensitivity infrared sensors and narrow-band filter modules to avoid environmental light interference, and their hardware cost is about 3-5 times that of a monocular camera with the same parameters. In clinical operations, binocular infrared technology relies on active infrared markers attached to the patient's body surface to achieve feature matching, but the markers are easily displaced due to patient breathing and body position changes, and the attachment process can cause discomfort for the patient. More importantly, the operator needs to manually select the corresponding relationship between the markers and the anatomical landmark points in the MRI image during registration, which is affected by the operator's anatomical knowledge reserve differences, leading to subjective bias in the registration results. In contrast, using a monocular camera instead of a binocular infrared camera can achieve systematic optimization while ensuring performance. Monocular cameras can directly extract stable anatomical feature points from the patient's body surface without relying on artificial markers, and have higher accuracy. In addition, binocular cameras often have image missing phenomena in single view, while monocular cameras can easily make up for missing information through multi-frame image parallax synthesis technology, and have strong noise resistance. Combining the above advantages, monocular cameras can better replace binocular cameras to improve system performance.

[0004] There is a multi-modal image conflict problem in real positioning tracking process using medical images, and image registration technology is an effective means to solve this problem. Specifically, image registration is a kind of image processing technology that aligns different images of the same object obtained by different or same sensors and perspectives in the same scene. Geometric transformation such as rigid, affine, projection, perspective and global transformation is needed in the process. After registration, the images will be located in a common coordinate system. Because point features can represent the geometric invariance between images, the image registration task can be simplified to a feature point based registration problem. Point set registration can be divided into rigid registration and non-rigid registration. Among them, rigid registration usually includes rotation, translation and scaling, and the solution of deformation and occlusion belongs to the category of non-rigid registration. For human head lesion positioning, the following matters need attention: 1) The self-motion of the human head collected by optical images can be approximately regarded as rigid motion, and the human head position remains unchanged while there is no organizational deformation during MRI data collection, so the human head model reconstructed from the collected MRI slices can be approximately regarded as a rigid body. 2) Considering that the collection time of MRI slice images and monocular data is different, the human facial expression and muscle state cannot be consistent, so there is a certain degree of local non-rigid deformation between them. 3) At the same time, due to the inconsistency of the collection means and equipment of the two modal data, there is a large global non-rigid deformation between them.

[0005] Considering the above needs and comparing with some existing image registration algorithms, the thin plate spline (TPS) algorithm can realize global continuous smooth deformation by minimizing the bending energy functional, while strictly meeting the corresponding constraints of the landmark points. In solving the problems of B-spline control point dependence, Demons gradient driven limitation and MLS neighborhood radius sensitivity, TPS algorithm shows three core advantages of global deformation coordination, deformation rationality guaranteed by physical mechanism and landmark point error self-adaptive ability; and the deformation naturalness guaranteed by its energy minimization mechanism is highly consistent with the topological structure maintenance requirement.

[0006] However, the existing TPS-based medical image registration research has gaps in the aspects of object, deformation scene and modality combination, which needs targeted algorithm improvement and innovation.

[0007] In addition, the existing non-rigid registration method often fails to achieve accurate alignment when dealing with images with large initial pose differences, and the deformation model is easily disturbed by global misplacement and produces overcorrection phenomenon.

[0008] If non-rigid registration is directly implemented without rigid registration, the overall spatial relationship of the anatomical structure will be distorted due to the lack of global pose constraints, which not only cannot guarantee the registration accuracy, but also may introduce false error. At the same time, the non-rigid algorithm in the table has limitations in complex deformation modeling capability, the fusion degree of cross-modal features is not good, and convergence lag problem is easy to appear in large-scale displacement scene, which is difficult to adapt to the actual application requirements.

[0009] The mathematical modeling method has the advantages of quantitative analysis of complex motion, systematic tracing of error sources, and accurate description of the coupling relationship of multiple parameters, which provides a reliable technical path for solving the above problems. Therefore, the head motion posture can be analyzed in detail by the mathematical modeling method, the rigid registration is completed to realize the preliminary alignment of the global pose, and the reliable spatial reference is provided for the subsequent non-rigid registration, and then the local deformation is corrected in detail, so as to meet the strict requirements of the task on the registration accuracy. SUMMARY

[0010] In view of the problems in the prior art, the application provides a multi-modal image cooperative head lesion multi-level automatic positioning method, which aims to improve the accuracy of lesion position cutting point, reduce the dependence on manual experience in clinical practice, and realize the positioning function of high-precision and low-cost image navigation system.

[0011] A multi-modal image cooperative head lesion multi-level automatic positioning method, comprising the following steps:

[0012] Step 1: Collecting head MRI slice images containing lesions and head optical images, and based on the head MRI slice images, constructing an MRI three-dimensional head surface model of the outer skin of the head and a three-dimensional lesion model of the lesion site;

[0013] Step 2: Adjusting the MRI three-dimensional head surface model to a preset posture sequence by Euler angle transformation, and for each posture, cutting the corresponding MRI two-dimensional body surface image at a specified step, and grouping these MRI two-dimensional body surface images into an MRI two-dimensional body surface image dataset after dimension reduction of the MRI three-dimensional head surface model;

[0014] Step 3: Manually labeling the positions of 7 key position points in the head MRI image dataset and the head optical image and extracting the coordinates of each key position point; the 7 key position points are: left outer canthus, left inner canthus, right outer canthus, right inner canthus, nose tip, left corner of the mouth and right corner of the mouth;

[0015] Step 4: Spatially solving the head posture of the image in the yaw, lateral flexion and pitch directions in the MRI two-dimensional body surface image dataset based on a geometric position estimation method, combining the corresponding rotation angle of the MRI three-dimensional head surface model, and establishing a displacement distance-angle analysis model;

[0016] Step 5: Initial alignment of the optical image pair to the MRI two-dimensional body surface image by using the offset distance-angle analytical model;

[0017] Step 6: Subsite registration of multi-modal images and automatic positioning of head lesions in optical images, specifically: establishing a lesion point position map relative to key position points by Delaunay triangulation method, and improving the thin plate spline function by using the relative position relationship to realize the position mapping of the lesion point in the MRI two-dimensional body surface image to the target lesion point in the optical image.

[0018] Further, the 7 key position points are respectively the tip of the nose, the left outer canthus, the left inner canthus, the right inner canthus, the right outer canthus, the left corner of the mouth and the right corner of the mouth, and the coordinates of the tip of the nose, the left outer canthus, the left inner canthus, the right inner canthus, the right outer canthus, the left corner of the mouth and the right corner of the mouth are respectively 、 、 、 、 、 、 ;

[0019] Step 4 is specifically:

[0020] Step 4.1: Introducing the position of the root of the nose above the midpoint of the line connecting the inner corners of the two eyes as an auxiliary reference , taking the frontal view as the standard to calculate the offset degree in each direction;

[0021] Step 4.2: Lateral flexion angle estimation;

[0022] The specific calculation of the lateral flexion degree is as follows:

[0023]

[0024] ;

[0025] Wherein, is the inclination angle of the line passing through the root of the nose and the tip of the nose;

[0026] Step 4.3: Calculation of the yaw degree:

[0027] The calculation process is as follows:

[0028]

[0029]

[0030]

[0031] Wherein, the ratio of the actual distance between the left inner canthus and the right inner canthus in the yaw state to the corresponding distance in the reference state ; , , X-coordinate of the tip of the nose, left inner corner of the eye, right inner corner of the eye in the initial state, respectively;

[0032] Step 4.4: Determine the degree of pitch:

[0033] Quantitative analysis of head pitch angle

[0034]

[0035]

[0036]

[0037] wherein, , is the original item weight, and the difference between the vertical coordinates of the root of the nose and the tip of the nose in the reference state is calculated and the real-time state difference ; , respectively represent the tip of the nose coordinate and the vertical coordinate of the root of the nose in the initial state, represent the normalized dynamic offset.

[0038] Further, step 5 is specifically:

[0039] Step 5.1: According to the formula in steps 4.1-4.4, the head posture in the head optical image is spatially solved to obtain the posture parameters in the yaw, lateral bending and pitch directions, and the offset degree in each direction is quantified;

[0040] Step 5.2: Combine the offset distance-angle analysis model established in step 4 to perform spatial transformation and image interception on the MRI three-dimensional head surface model, realize the preliminary posture alignment between the head optical image and the MRI three-dimensional head surface model, and reduce the average error of the coordinates of the corresponding key position points between the MRI two-dimensional body surface image and the head optical image to tend to be stable.

[0041] Further, step 6 is specifically:

[0042] Step 6.1: Collect the coordinates of the key position points in the MRI two-dimensional body surface image to establish a source image point set , ; Collect the coordinates of the key position points in the head optical image to establish a target image point set , ;

[0043] Step 6.2: Process the MRI two-dimensional body surface image and optical image respectively. Perform Delaunay triangulation with the seven selected key locations as vertices to determine the positional relationship of the remaining points relative to the key locations: inside the nearest triangle, on the side of the nearest triangle, outside all Delaunay triangles and outside the circumcircle, and outside all Delaunay triangles and inside the circumcircle.

[0044] Step 6.3: Based on the lesion point Modify the deformation function in the thin plate spline function according to the different positional relationships with the Delaunay triangle. and constraint functions To perform local structural constraints;

[0045] Step 6.4: Introduce the lesion point relative centroid coordinates of control points It is used to replace the absolute coordinates of the traditional fixed coordinate system, thereby measuring the distribution relationship between lesion points and control points and completing global constraints.

[0046] Further:

[0047] Step 6.2 specifically involves:

[0048] Step 6.2.1: Build a KD tree based on the coordinates of the lesion point. Perform the first partition of the root node according to the x-coordinate, perform the second partition based on the y-coordinate, and then loop according to the x-coordinate until all partitions are completed and find the control point closest to the lesion point.

[0049] Step 6.2.2: Extract the triangle containing the nearest control point, through the lesion point. The nearest triangle is obtained by taking the perpendicular distance from the sides of each triangle; let the nearest triangle in the MRI two-dimensional body surface image be... The vertex coordinates are respectively .

[0050] Further:

[0051] The lesion point within the nearest triangle will be recorded as the lesion point. The lesion point closest to the side of the triangle is denoted as the lesion point. Lesions located outside the Delaunay triangle and outside its circumcircle are denoted as lesion points. All lesions outside the Delaunay triangle but inside its circumcircle are denoted as lesion points. ;

[0052] Step 6.3 specifically involves:

[0053] Step 6.3.1: When the lesion point When inside the nearest triangle,

[0054] In the barycentric coordinates are defined as and satisfy the following equation:

[0055]

[0056] where, are all calculated by area ratio, ;

[0057]

[0058] where, are the control points except the nearest triangle vertex; are the control point weights, to strengthen the influence of the nearest control point, the parameter value is adjusted according to experience, so that is slightly larger than ; , is a local spatial structure constraint function; represents a hyperparameter, represents a radial basis function, represents a triangle vertex, represents a local spatial structure constraint function;

[0059] Step 6.3.2: When the lesion point is on the nearest triangle edge, increase the influence of the common edge vertex on , and constrain the position relationship between the common edge vertex and to ensure that the position relationship is relatively unchanged before and after transformation:

[0060] There are and , that is, located on the common edge , for example, the vertices are

[0061]

[0062] where, are the control points except the common edge vertex, is the transformation ratio, and the calculation formula is:

[0063]

[0064] where, represents the longitudinal coordinate of the lesion point, , both represent the longitudinal coordinate of the common edge vertex, , both represent the longitudinal coordinate of the common edge The x-coordinate of the vertex;

[0065] Step 6.3.3: When the lesion point When the lesion point is outside all circumcircles of the nearest triangle, there is no need to perform a secondary subdivision of the Delaunay triangle. The barycentric coordinates of the nearest triangle are extended to generalized barycentric coordinates. Using the generalized barycentric coordinates as constraints, the relative positions of the generalized barycentric coordinates remain unchanged before and after the transformation.

[0066]

[0067] in, These are the control points other than the nearest triangle vertex. Represents the number of nearest triangles;

[0068] Step 6.3.4: When the lesion point When the lesion is outside the nearest triangle and inside any circumcircle, a second Delaunay triangle subdivision must be performed first, that is, the existing lesion point... Together with the control points, a Delaunay triangle is established, representing the lesion point. The position of the vertex is determined by the positions of the other vertices that together form the Delaunay triangle, and is jointly optimized using the side length ratio and angle difference, introducing an adaptive weight allocation mechanism:

[0069]

[0070] in, The remaining control points are those that, together with the other vertices forming the Delaunay triangle, are excluded. This is the side length ratio constraint factor. For angle constraint factors, For the transformation of the ray and The angle formed For the transformed ray and The angle formed.

[0071] Further:

[0072] Step 6.4 specifically involves:

[0073] The global constraints are:

[0074]

[0075] in, This represents the centroid coordinates of the lesion point relative to the control point. Represents the coordinates of the seven control points;

[0076] The centroid of the control point in the target image is used in the calculation Instead of the control point to calculate;

[0077]

[0078] Therefore, the global constraint term Can be defined as:

[0079]

[0080] So far, the registration constraint equation is modified as:

[0081] .

[0082] The beneficial effects of the present application: fusion of rigid initial registration and non-rigid accurate registration algorithm based on thin plate spline algorithm, optimization of thin plate spline algorithm from local and global level by Delaunay triangulation method, balance of registration accuracy, calculation efficiency and robustness by using multi-objective constraint processing mechanism, at the same time, light weight 7 key position points are used to replace the traditional dense scheme to reduce interference, combined with head posture geometric analysis and initial alignment to realize accurate positioning, taking into account the theoretical reliability and engineering practicability, providing strong technical support for multi-modal image navigation; the present application can effectively improve the accuracy of lesion position cutting point by MRI image reconstruction and dimension reduction, facial key position point selection, head posture preliminary alignment and fine registration between virtual and real images, reduce the dependence on artificial experience in clinic, realize the positioning function of high-precision, low-cost image navigation system. BRIEF DESCRIPTION OF DRAWINGS

[0083] Figure 1 It is a head lesion multi-level automatic positioning method flow chart;

[0084] Figure 2 It is a MRI 3D surface model reconstruction flow chart;

[0085] Figure 3 It is a lateral flexion curve in the offset distance-angle analytical model;

[0086] Figure 4 It is a yaw curve in the offset distance-angle analytical model;

[0087] Figure 5 It is a pitch curve in the offset distance-angle analytical model;

[0088] Figure 6 It is a MRI two-dimensional surface image and optical image initial registration flow chart;

[0089] Figure 7 It is a head looking straight without offset initial state schematic diagram after initial registration result;

[0090] Figure 8 Fig. 6 is a combined posture schematic diagram of right lateral flexion, right yaw and head-up after the initial registration result;

[0091] Figure 9 Fig. 7 is a combined posture schematic diagram of left lateral flexion, left yaw and head-down after the initial registration result;

[0092] Figure 10 Fig. 8 is a position distribution diagram of lesion points and control points;

[0093] Figure 11 Fig. 9 is an external position of a circumscribed circle of a Delaunay triangle;

[0094] Figure 12 Fig. 10 is a comparison diagram of secondary subdivision of a Delaunay triangle;

[0095] Figure 13 Fig. 11 is a comparison diagram of position distribution after transformation of key position points;

[0096] Figure 14 Fig. 12 is a registration deviation in an initial state;

[0097] Figure 15 Fig. 13 is a registration deviation in state 1;

[0098] Figure 16 Fig. 14 is a registration deviation in state 2;

[0099] Figure 17 Fig. 15 is a comparison of lesion point position distribution after transformation;

[0100] Figure 18 Fig. 16 is an evaluation index result of different deformation algorithms. DETAILED DESCRIPTION

[0101] The application will be described in detail below with reference to the drawings. The embodiments of the application are described in detail below, and examples of the embodiments are shown in the drawings, wherein the same or similar reference numerals represent the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the drawings are exemplary and are only used to explain the application, and should not be interpreted as a limitation on the application. The left, middle, right, up, down and other orientation terms in the embodiments of the application are only relative concepts or are referenced to the normal use state of the product, and should not be considered as limiting.

[0102] A multi-modal image collaborative head lesion multi-level automatic positioning method, as shown in Figure 1 includes the following steps:

[0103] Step 1: Collecting head MRI slice images and head optical images containing lesions, as shown in Figure 2As shown, by noise filtering, threshold segmentation, retaining the largest island, contour splicing, image interpolation, volume rendering, surface smoothing and other operations, and based on the head MRI slice image, the MRI three-dimensional head surface model of the head lateral skin and the three-dimensional lesion model of the lesion site are constructed;

[0104] Step 2: Adjust the MRI three-dimensional head surface model to the preset posture sequence by Euler angle transformation, and intercept the corresponding MRI two-dimensional body surface image at the specified step for each posture, and group these MRI two-dimensional body surface images to form the MRI two-dimensional body surface image dataset after dimension reduction of the MRI three-dimensional head surface model; in the angle transformation parameter setting, the single rotation angle needs to be selected as a moderate value: too small angle will significantly increase the experimental workload and reduce the efficiency; too large angle may cause the posture sampling interval to be too wide, which affects the accuracy of subsequent posture analysis, so the optimal single rotation angle range needs to be determined through pre-experiment;

[0105] Step 3: Manually mark the positions of 7 key position points in the head MRI image dataset and the head optical image and extract the coordinates of each key position point; the 7 key position points are: left outer canthus, left inner canthus, right outer canthus, right inner canthus, nose tip, left corner of the mouth and right corner of the mouth; the coordinates of the 7 key position points are shown in the following table:

[0106] Table 1: Key position points and their coordinates

[0107]

[0108] Step 4: Based on the geometric position estimation method, the head posture of the image in the yaw, lateral flexion and pitch directions in the MRI two-dimensional body surface image dataset is spatially calculated, and the corresponding MRI three-dimensional head surface model rotation angle is combined to establish a shift distance-angle analysis model; specifically:

[0109] Step 4.1: Considering that facial expressions may change unconsciously, the change range of the corners of the mouth is relatively large, therefore, the posture estimation does not consider the coordinate positions of the corners of the mouth; the nose root position above the midpoint of the line connecting the inner corners of the eyes is introduced as an auxiliary reference, and the front face view is taken as the standard to calculate the offset degree in each direction;

[0110] Step 4.2: As shown, the lateral flexion angle is estimated; Figure 3

[0111] The neutral position in facial anatomy is taken as the reference for calculating the left and right lateral flexion offsets of the multi-modal data, at this time the coronal plane of the face is perpendicular to the ground, and this state is defined as the zero reference position of rotation offset, i.e. the offset degree is 0°; in this standard posture, a Cartesian plane rectangular coordinate system is established, and the inclination angle of the straight line passing through the nose root and the nose tip is calculated ​, which is normalized to evaluate the degree of rotation offset; however, it is worth noting that due to the tilt of the upper body of the human body, the head will be offset to be parallel to the ground, at which time the absolute value of the tilt angle is , the left and right outer canthus are used to distinguish the left and right directions of the head; therefore, the specific calculation of the degree of lateral flexion is as follows:

[0112]

[0113] ;

[0114] wherein, is the tilt angle of the line passing through the root of the nose and the tip of the nose;

[0115] Step 4.3: As shown in Figure 4 , the yaw degree is calculated as:

[0116] The quantitative determination of the yaw state is achieved by the feature point distance ratio method; however, due to the possible large difference in the left and right horizontal offset amplitude, when the full profile situation occurs, the key position points on the other side will be unable to be obtained due to complete occlusion, resulting in the failure of traditional posture calculation methods based on key position points, and therefore most current researches do not consider such extreme offset scenarios; to solve this problem, we expand the application range to the visible angle of the outer canthus; specifically, two-dimensional images of the target in the reference position (front view) and in the yaw state are collected, the left and right inner canthus feature points in the images are extracted, and the ratio of the actual distance between the left and right inner canthi in the yaw state to the corresponding distance in the reference state is calculated ; this ratio can be used as a quantitative index to establish a yaw motion determination model based on the ratio coefficient, and the numerical size is used to accurately evaluate the horizontal motion state of the target, and the change in the horizontal coordinate value of the tip of the nose is used to measure the left and right motion directions, thereby achieving non-contact measurement of two-dimensional space yaw displacement; the calculation process is as follows:

[0117]

[0118]

[0119]

[0120] wherein, the ratio of the actual distance between the left and right inner canthi in the yaw state to the corresponding distance in the reference state is calculated ; , , the horizontal coordinates of the tip of the nose, the left inner corner of the eye, and the right inner corner of the eye in the initial state, respectively

[0121] Step 4.4: As shown inFigure 5 The pitch degree determination is shown as follows:

[0122] For the quantitative evaluation of the head pitch posture, a calculation model based on feature point vector analysis is constructed. By synchronously collecting two-dimensional images of the target in the reference orthoview and the pitch state, the difference in the sagittal coordinates of the root of the nose and the tip of the nose in the reference state is calculated and the difference in the real-time state , and a posture calculation equation is constructed using the ratio relationship to convert the coordinate components of the two vectors into pitch degree values. This method realizes quantitative analysis of the head pitch angle by comparing the changes in the feature vectors in different states:

[0123]

[0124]

[0125]

[0126] In view of the inherent defects in the above formula in the process of quantifying head lowering, a normalized dynamic offset is introduced , a dynamic normalization model based on real-time calculation of the root coordinates is constructed to accurately map the geometric characteristics of head space rotation; the dimensionless displacement ratio is used to replace the traditional absolute distance measurement, effectively decoupling the many-to-one mapping relationship between distance and angle in the nonlinear deformation process, and eliminating the quantitative ambiguity of the original model in the large-angle head-down scenario from the theoretical level; therefore, is optimized as:

[0127]

[0128] wherein , is the weight of the original term, The smaller the correction effect is better; the difference in the sagittal coordinates of the root of the nose and the tip of the nose in the reference state is calculated and the difference in the real-time state ; , represent the coordinates of the tip of the nose and the longitudinal coordinates of the root of the nose in the initial state, represents the normalized dynamic offset;

[0129] Step 5: Use the offset distance-angle analysis model to initially align the optical image pair with the MRI two-dimensional surface image. In the biomechanical analysis of head movement, if the muscle synergy and joint biomechanical effects are not considered, the yaw, lateral flexion and pitch of the head can be approximately regarded as independent rigid body movements; as shown in Figure 6 , specifically:

[0130] Step 5.1: Perform spatial calculation of head posture in head optical images according to the formulas in steps 4.1-4.4 to obtain the posture parameters in the three directions of yaw, flexion and pitch, and quantify the degree of offset in each direction; while the quantification of flexion motion directly uses angle values, without the need to construct an analytical relationship between it and the offset distance.

[0131] Step 5.2: Combining the offset distance-angle analytical model established in Step 4, perform spatial transformation (single rotation in three directions) and image cropping on the MRI three-dimensional head surface model to achieve preliminary pose alignment between the head optical image and the MRI three-dimensional head surface model, so that the average error of the coordinates of the corresponding key points between the MRI two-dimensional body surface image and the head optical image is reduced to a stable value.

[0132] Step 6: Multimodal image segmentation and automatic localization of head lesions in optical images. This involves establishing a position map of the lesion relative to key locations using Delaunay triangulation, and improving the thin-plate spline function (TPS) based on relative positional relationships to achieve positional mapping from lesion points in the MRI 2D surface image to target lesion points in the optical image. While pose uniformity is effective in reducing global deformation caused by rigid head movements, the complexity of human head movements makes it difficult to achieve complete alignment of different modal data using pose uniformity alone. Affine transformation can compensate for the shortcomings of this technique in global registration, and the thin-plate spline function can decompose into an affine transformation component to handle rigid body motion. Given the significant differences in color and texture between MRI 2D surface images and optical images, and the diverse facial morphologies during data acquisition, simple rigid transformations are difficult to implement. The nonlinear part of the thin-plate spline function can solve this problem. Specifically:

[0133] Step 6.1: Collect the coordinates of key points in the 2D MRI surface images to establish the source image point set. , Collect the coordinates of key points in the head optical images to establish a target image point set. , Given that the treatment areas vary for different diseases and the number of lesions is not unique, for ease of analysis, we will take the location of a single lesion as an example, and its location in the MRI two-dimensional body surface image and the head optical image are as follows: , ;

[0134] Step 6.2: Process the MRI 2D surface image and optical image separately, performing Delaunay triangulation with the seven selected key locations as vertices, such as... Figure 10 and Figure 11The position relationship of the rest of the points relative to the key position points is determined as shown: inside the nearest triangle, on the nearest triangle edge, outside all Delaunay triangles and outside the circumscribed circle, outside all Delaunay triangles and inside the circumscribed circle; specifically:

[0135] Step 6.2.1: Establish a KD tree according to the coordinates of the lesion points, complete the first division of the root node according to the x coordinate, perform secondary division based on the y coordinate, and then cycle with the x coordinate until all divisions are completed to find the control points closest to the lesion points;

[0136] Step 6.2.2: Extract the triangle in which the nearest control point is located, and obtain the nearest triangle through the perpendicular distance of the lesion point to each triangle edge; let the nearest triangle in the MRI two-dimensional surface image be , and the vertex coordinates be

[0137] Step 6.3: According to the different position relationships of the lesion point and the Delaunay triangle, the influence of the key position points is also different, therefore, the deformation function and the constraint function in the thin plate spline function are modified for local structure constraint; specifically:

[0138] The original TPS algorithm formula is:

[0139]

[0140]

[0141]

[0142]

[0143] wherein, , is a linear matrix, , is a data fitting error, is a bending energy, is a radial basis function; due to the complex facial deformation, a second-order term is introduced for calculation; represents a two-dimensional Euclidean space, represents the weight coefficient corresponding to the key position point coordinate , represents the source image key position point coordinate, represents the lesion point coordinate, represents the independent variable in the radial basis function;

[0144]

[0145] the lesion point within the nearest triangle is recorded as the lesion point the lesion point on the edge of the nearest triangle is recorded as the lesion point the lesion point outside all the Delaunay triangles and outside the circumcircle is recorded as the lesion point the lesion point outside all the Delaunay triangles and inside the circumcircle is recorded as the lesion point ;

[0146] only the source image point is constrained by the deformation function The distance difference between the point obtained after transformation and the target image point, in the application background of determining the specific position of the lesion point on the optical image, ignoring the association of the control point with the global structure, local over-distortion is easy to occur, which affects the positioning accuracy of the lesion point, so as to form an effective constraint; as shown in Figure 12 To ensure that the spatial structure controlled by the key position points before and after transformation is unchanged, the position map of the key position points is established by borrowing the Delaunay triangle division method which can ensure that the triangles formed are mutually non-overlapping, can cover the entire plane, and each point is not in the circumcircle of the triangle not including the point;

[0147] Step 6.3.1: When the lesion point is within the nearest triangle, the influence degree of the nearest triangle vertex pair is increased, and then the barycentric coordinates are used as constraints to ensure that the positions of the barycentric coordinates before and after transformation are relatively unchanged;

[0148] The barycentric coordinates in are defined as and satisfy the following formula:

[0149]

[0150] Wherein, are all obtained by area ratio calculation, ;

[0151]

[0152] Wherein, is the control point other than the nearest triangle vertex; is the control point weight, in order to strengthen the influence degree of the nearest control point, the parameter value is adjusted according to experience, so that is slightly larger than ; , is a local spatial structure constraint function; represents a hyperparameter (manually adjusted by man), denotes a radial basis function, denotes a triangle vertex, denotes a local spatial structure constraint function;

[0153] Step 6.3.2: When the lesion point is on the edge of the nearest triangle (possibly more than one), increase the influence degree of the common edge vertex pair , and ensure the relative invariance of the position relationship before and after the transformation by taking the position relationship of the common edge vertex and as a constraint:

[0154] Take the existence of and , i.e. is located on the common edge , as an example, the vertices are ;

[0155]

[0156] wherein, is the rest of the control points except the common edge vertex, is the transformation scale, and the calculation formula is:

[0157]

[0158] wherein, denotes the longitudinal coordinate of the lesion point, , both denote the longitudinal coordinate of the common edge vertex, , both denote the horizontal coordinate of the common edge vertex;

[0159] Step 6.3.3: When the lesion point is outside all the circumcircles of the nearest triangle, there is no need to subdivide the Delaunay triangle twice, and the lesion point is expanded to the generalized barycentric coordinate of the nearest triangle, and the generalized barycentric coordinate is taken as a constraint to ensure the relative invariance of the position of the generalized barycentric coordinate before and after the transformation:

[0160]

[0161] wherein, is the rest of the control points except the nearest triangle vertex, represents the number of the nearest triangle; when n=2, ; denotes a set with subscript , , , These represent the coordinates of the barycenter. Indicates subscript The set ( and (representing the centroid coordinates of different nearest triangles);

[0162] Step 6.3.4: When the lesion point When the lesion is outside the nearest triangle and inside any circumcircle, the Delaunay triangle partitioning principle is not satisfied. Therefore, a second Delaunay triangle partitioning is required, i.e., the existing lesion point... Together with the control points, a Delaunay triangle is established, representing the lesion point. The position of the vertex is determined by the positions of the other vertices that together form the Delaunay triangle, and is jointly optimized using the side length ratio and angle difference, introducing an adaptive weight allocation mechanism:

[0163]

[0164] in, The remaining control points are those that, together with the other vertices forming the Delaunay triangle, are excluded. This is the side length ratio constraint factor. For angle constraint factors, For the transformation of the ray and The angle formed For the transformed ray and The angle formed;

[0165] Step 6.4: Introduce the lesion point relative centroid coordinates of control points This is used to replace the absolute coordinates of the traditional fixed coordinate system, thereby measuring the distribution relationship between lesion points and control points and completing global constraints; specifically:

[0166] The global constraints are:

[0167]

[0168] in, This represents the centroid coordinates of the lesion point relative to the control point. This represents the coordinates of the seven control points (to distinguish them from the subscript i);

[0169] To highlight relative positional relationships, ensure constraint consistency, and conform to spatial transformation logic, the centroids of control points in the target image are used in the calculation. Calculations are performed in place of control points;

[0170]

[0171] Thus, the global constraint term may be defined as:

[0172]

[0173] On the basis of the rigid preliminary registration in the early stage, a geometric topology attribute fusion strategy is further introduced, a local and global multi-level constraint system is constructed for the lesion point position transformation, and optimization iteration is carried out based on the thin plate spline function; at the local structure constraint level, the spatial position relationship between the lesion point and the nearest neighbor triangular element is accurately distinguished, the deformation function is modified for different geometric configurations, and a local spatial structure constraint condition is constructed; at the global structure constraint level, a constraint model is established based on the joint mapping relationship between the centroid coordinates of the lesion point and the control point, and a composite constraint system is formed by introducing a surface distortion energy constraint function; thus, the registration constraint equation is modified as:

[0174] .

[0175] Example verification and result analysis

[0176] Since the single-person data has typical sample representativeness in methodology, the subsequent analysis uses single-instance data as a standardized research object. This analysis framework can be directly extended to the analysis of other samples by parameter adjustment, thereby ensuring the universality and scalability of the research method.

[0177] The offset distance-angle analytical model is used to register three different states of the head in the optical image. The initial state is the frontal view without any angular offset. State 1 is a combined motion state of right flexion, right yaw, and head raising. State 2 is a combined motion state of left flexion, left yaw, and head lowering. After the optical image pair is roughly aligned with the MRI two-dimensional surface image, the coordinate distribution comparison of the key position points of the MRI two-dimensional surface image and the head optical image is shown in FIG. 8. In order to have a unified reference framework for the coordinate data of different samples for easy comparison, the coordinates of each key position point are normalized with the coordinate of the nose tip as the standard; Figures 7 to 9

[0178] From Figures 7 to 9 ​It can be seen that after initial registration, there are still different degrees of deviation in each control point, so further accurate deformation correction is needed; step (6) is used for accurate registration. In order to verify the registration effect of the method, B-spline, Demons, MLS and TPS original algorithms are selected as comparison methods to process three groups of experimental data. First, the coordinate position distribution of the key position points after deformation and the key position points of the target image is analyzed for intuitive verification, and the results are shown in Fig. 13. At the same time, the registration error of each key position point is quantitatively calculated, and the related results are shown in Figs. 14-16.

[0179] From the results of Figs. 14-16, it can be seen that in the three groups of experimental data, the registration accuracy of the method proposed in this paper is better, and it is more suitable for actual application requirements; in order to further accurately evaluate the registration deviation between key position points, this study uses mean absolute error (MAE) for quantitative analysis, and the related results are shown in Table 2.

[0180]

[0181] wherein, represents the coordinates of the transformed MRI key position points.

[0182] Table 2 MAE calculation results

[0183]

[0184] From the data in Table 2, the algorithm proposed in this application has the smallest single-point registration deviation and overall registration deviation of all key position points, with an accuracy improvement of about 97.45% compared to the original algorithm. However, the core demand of this study is not limited to optimizing the registration performance of key position points, but the starting point is to accurately quantify the position change of the lesion points. Therefore, this application has carried out position transformation on all pixel points in the whole image, and the spatial distribution results after transformation are shown in Fig. 17. In this figure, the first column corresponds to the original MRI image, the last column corresponds to the optical target image, and the middle five columns are the processing effect comparison of different algorithms.

[0185] The goal of this study is to make the position distribution of each point in the transformed image as close as possible to the target image; from Figure 17It can be seen that in the initial state, the transformation effect of the MLS algorithm is poor, and the rest of the algorithms have similar performance and high similarity with the target image; in state 1, only the algorithm proposed in this paper has high approximation degree with the target image; in state 2, the MLS algorithm and the method proposed in this paper have the best effect; in order to quantitatively evaluate the similarity between the image processed by the deformation algorithm and the target image, this paper uses mutual information (MI), root mean square error (RMSE) and cross correlation coefficient (CC) as performance indicators, and the calculation results are as follows Figure 18 .

[0186]

[0187]

[0188]

[0189] wherein, represents the edge probability distribution of the image, is the total number of pixels of the image; is the pixel value of the RGB image, is the pixel value of the transformed image of the MRI two-dimensional surface image, is the joint probability, is the edge probability, is the pixel position coordinate, , is the pixel mean value of the corresponding image in the first .

[0190] The MI, RMSE and CC quantitative results of the three groups of data are as follows Figure 18 It can be seen that the registration performance of the algorithm proposed in this application is significantly better than that of the comparative methods. Specifically, the MI values of B-spline, Demons, MLS, TPS (original) and the algorithm in this paper are 0.426, 0.387, 0.393, 0.324 and 0.904 respectively, and the CC values are 0.545, 0.522, 0.557, 0.430 and 0.954 respectively, indicating that the algorithm proposed in this paper has the highest similarity with the target image; the RMSE values are 6.649, 6.773, 6.479, 7.129 and 4.131 respectively, further confirming that the registration accuracy of the algorithm proposed in this paper is the best, fully verifying its performance advantage.

[0191] The above shows and describes the basic principles, main features and advantages of the present application. Those skilled in the art should understand that the present application is not limited to the above-mentioned embodiments, and the above-mentioned embodiments and descriptions in the specification are only to illustrate the principles of the present application. Without departing from the spirit and scope of the present application, various changes and improvements can be made to the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A multi-modal image-based, multi-level automatic localization method for head lesions, characterized in that, Includes the following steps: Step 1: Acquire head MRI slice images and head optical images containing lesions, and based on the head MRI slice images, construct a three-dimensional MRI head surface model of the outer skin of the head and a three-dimensional lesion model of the lesion site. Step 2: The MRI three-dimensional head surface model is successively adjusted to the preset pose sequence through Euler angle transformation, and the corresponding MRI two-dimensional body surface image is cropped for each pose according to the specified step size. These MRI two-dimensional body surface images are combined into an MRI two-dimensional body surface image dataset after the MRI three-dimensional head surface model is reduced in dimension. Step 3: Manually label the positions of 7 key points in the head MRI image dataset and head optical image, and extract the coordinates of each key point; the 7 key points are: left lateral canthus, left medial canthus, right lateral canthus, right medial canthus, tip of nose, left corner of mouth, and right corner of mouth; Step 4: Based on the geometric position estimation method, the head posture in the three directions of yaw, flexion and pitch of the central images in the MRI two-dimensional body surface image dataset is spatially calculated. Combined with the corresponding rotation angle of the MRI three-dimensional head surface model, an offset distance-angle analytical model is established. Step 5: Initially align the optical image with the MRI 2D body surface image using the offset distance-angle analytical model; Step 6: Multimodal image segmentation registration and automatic localization of head lesions in optical images. Specifically, the position map of the lesion point relative to key location points is established by using the Delaunay triangulation method, and the thin plate spline function is improved by using the relative position relationship to realize the position mapping from the lesion point in the MRI two-dimensional body surface image to the target lesion point in the optical image.

2. The multi-modal image-based automatic localization method for head lesions according to claim 1, characterized in that: The seven key locations are the tip of the nose, the left outer canthus, the left inner canthus, the right inner canthus, the right outer canthus, the left corner of the mouth, and the right corner of the mouth. The coordinates of the tip of the nose, the left outer canthus, the left inner canthus, the right inner canthus, the right outer canthus, the left corner of the mouth, and the right corner of the mouth are as follows: , , , , , and ; Step 4 is as follows: Step 4.1: Introduce the area above the midpoint of the line connecting the inner corners of both eyes, at the root of the nose. As an auxiliary reference, the frontal view is used as the standard to calculate the degree of offset in each direction; Step 4.2: Estimate the lateral flexion angle; The specific calculation of the degree of lateral flexion is as follows: ； ； in, The angle of inclination of the straight line between the root of the nose and the tip of the nose; Step 4.3: Yaw Calculation: The calculation process is as follows: ； ； ； Among them, the actual distance between the left and right inner canthi under yaw conditions is calculated. Distance corresponding to the reference state ratio ; , , The horizontal coordinates of the tip of the nose, the left inner corner of the eye, and the right inner corner of the eye in the initial state; Step 4.4: Determine the pitch level: Quantitative analysis of head pitch angle ； ； ； in, , As the weights of the original terms, calculate the difference in the ordinates of the bridge of the nose and the tip of the nose under the baseline condition. Difference from real-time status ; , These represent the initial coordinates of the nasal tip and the ordinate of the nasal root, respectively. This represents the normalized dynamic offset.

3. The multi-modal image-based automatic localization method for head lesions according to claim 2, characterized in that: Step 5 specifically involves: Step 5.1: Perform spatial calculation of the head posture in the head optical image according to the formulas in Steps 4.1-4.4, obtain the posture parameters in the three directions of yaw, lateral flexion and pitch, and quantify the degree of offset in each direction; Step 5.2: Combining the offset distance-angle analytical model established in Step 4, perform spatial transformation and image cropping on the MRI three-dimensional head surface model to achieve preliminary pose alignment between the head optical image and the MRI three-dimensional head surface model, thereby reducing the average error of the coordinates of corresponding key points between the MRI two-dimensional body surface image and the head optical image to a stable level.

4. The multi-modal image-based automatic localization method for head lesions according to claim 2, characterized in that: Step 6 specifically involves: Step 6.1: Collect the coordinates of key points in the 2D MRI surface images to establish the source image point set. , ; Collect the coordinates of key points in the head optical images to establish a target image point set. , ; Step 6.2: Process the MRI two-dimensional body surface image and optical image respectively. Perform Delaunay triangulation with the seven selected key locations as vertices to determine the positional relationship of the remaining points relative to the key locations: inside the nearest triangle, on the side of the nearest triangle, outside all Delaunay triangles and outside the circumcircle, and outside all Delaunay triangles and inside the circumcircle. Step 6.3: Based on the lesion point Modify the deformation function in the thin plate spline function according to the different positional relationships with the Delaunay triangle. and constraint functions To perform local structural constraints; Step 6.4: Introduce the lesion point relative centroid coordinates of control points It is used to replace the absolute coordinates of the traditional fixed coordinate system, thereby measuring the distribution relationship between lesion points and control points and completing global constraints.

5. The multi-modal image-based automatic localization method for head lesions according to claim 4, characterized in that: Step 6.2 specifically involves: Step 6.2.1: Build a KD tree based on the coordinates of the lesion point. Perform the first partition of the root node according to the x-coordinate, perform the second partition based on the y-coordinate, and then loop according to the x-coordinate until all partitions are completed and find the control point closest to the lesion point. Step 6.2.2: Extract the triangle containing the nearest control point, through the lesion point. The nearest triangle is obtained by taking the perpendicular distance from the sides of each triangle; let the nearest triangle in the MRI two-dimensional body surface image be... The vertex coordinates are respectively .

6. The multi-modal image-based automatic localization method for head lesions according to claim 5, characterized in that: The lesion point within the nearest triangle will be recorded as the lesion point. The lesion point closest to the side of the triangle is denoted as the lesion point. Lesions located outside the Delaunay triangle and outside its circumcircle are denoted as lesion points. All lesions outside the Delaunay triangle but inside its circumcircle are denoted as lesion points. ; Step 6.3 specifically involves: Step 6.3.1: When the lesion point When inside the nearest triangle, exist The centroid coordinates in the middle are defined as And satisfy the following formula: ； in, All are calculated from the area ratio. ; ； in, These are the control points other than the nearest triangle vertex; To assign control point weights and enhance the influence of the nearest control point, parameter values ​​are adjusted based on experience to achieve the desired effect. Slightly larger ; , For local spatial structure constraint functions; Indicates hyperparameters, Represents radial basis functions. Represents the vertices of a triangle. Represents the local spatial structure constraint function; Step 6.3.2: When the lesion point When on the nearest triangle edge, increase the number of vertices on the common edge. The degree of influence, and based on the common edge vertex and Positional relationships are used as constraints to ensure that the relative positional relationships remain unchanged before and after the transformation: To exist and ,Right now Located on the public side Taking the above as an example, the vertices are respectively ; ； in, These are the control points excluding the vertices on the common edges. To change the scale, the calculation formula is: ； in, Represents the vertical coordinate of the lesion point. , Both represent common edges Vertex ordinate, , Both represent common edges The x-coordinate of the vertex; Step 6.3.3: When the lesion point When the lesion point is outside all circumcircles of the nearest triangle, there is no need to perform a secondary subdivision of the Delaunay triangle. The barycentric coordinates of the nearest triangle are extended to generalized barycentric coordinates. Using the generalized barycentric coordinates as constraints, the relative positions of the generalized barycentric coordinates remain unchanged before and after the transformation. ； in, These are the control points other than the nearest triangle vertex. Represents the number of nearest triangles; Step 6.3.4: When the lesion point When the lesion is outside the nearest triangle and inside any circumcircle, a second Delaunay triangle subdivision must be performed first, that is, the existing lesion point... Together with the control points, a Delaunay triangle is established, representing the lesion point. The position of the vertex is determined by the positions of the other vertices that together form the Delaunay triangle, and is jointly optimized using the side length ratio and angle difference, introducing an adaptive weight allocation mechanism: ； in, The remaining control points are those that, together with the other vertices forming the Delaunay triangle, are excluded. This is the side length ratio constraint factor. For angle constraint factors, For the transformation of the ray and The angle formed For the transformed ray and The angle formed.

7. The multi-modal image-based automatic localization method for head lesions according to claim 3, characterized in that: Step 6.4 specifically involves: The global constraints are: ； in, This represents the centroid coordinates of the lesion point relative to the control point. Represents the coordinates of the seven control points; The centroids of control points in the target image are used in the calculation. Calculations are performed in place of control points; ； Therefore, global constraint terms It can be defined as: ； Therefore, the registration constraint equations are now corrected to: 。

Images