Head focus multi-level automatic positioning method based on multi-modal image cooperation

Abstract

The invention provides a multi-modal image collaborative multi-level automatic positioning method for a head focus, and the method comprises the steps: fusing rigid initial registration with a non-rigid precise registration algorithm based on a thin-plate spline algorithm, and optimizing the thin-plate spline algorithm from local and global levels through a Delaunay triangulation method. A multi-target constraint processing mechanism is adopted to balance registration precision, calculation efficiency and robustness, meanwhile, seven lightweight key position point characterization replaces a traditional dense scheme to reduce interference, head posture geometric analysis and initial alignment are combined to achieve accurate positioning, theoretical reliability and engineering practicability are both considered, and the method is suitable for large-scale popularization and application. Powerful technical support is provided for multi-mode image navigation; through MRI image reconstruction and dimension reduction, face key position point selection, head posture preliminary alignment and fine registration between virtual and real images, the accuracy of a focus position entry point is effectively improved, the dependence on human experience in clinic is reduced, and a high-precision and low-cost image navigation system positioning function is realized.

Description

Technical Field

[0001] This invention relates to the field of medical image navigation system optimization, specifically to a multi-modal image-based multi-level automatic localization method for head lesions. Background Technology

[0002] Surgical procedures are often limited by blind spots caused by lesions embedded in the patient's soft and hard tissues. Magnetic resonance imaging (MRI) can effectively display information such as the location, size, and shape of lesions and is often used for preoperative diagnosis. Computer-aided surgical navigation technology, which integrates medical imaging, computer science, positioning and tracking technology, and robotics, has been developing for decades. Surgeons can use this technology to analyze the patient's anatomical structure during surgery and provide precise surgical path guidance. In image-guided robotic treatment systems, since medical images such as MRI and CT scans are mostly presented in the form of tomographic slices, it is difficult to fully represent the patient's surface anatomy. Therefore, three-dimensional reconstruction of medical images is necessary. Precisely aligning the reconstructed three-dimensional virtual image with the patient's actual intraoperative anatomical structure—that is, registering the patient's virtual and real images—is the core of achieving the system's precise treatment function. Its performance directly affects surgical accuracy and safety, and must simultaneously meet engineering indicators such as high registration accuracy, real-time performance, low clinical invasiveness, and controllable system cost. This is a critical problem that urgently needs to be solved.

[0003] While the binocular infrared stereotactic technology widely used in current image navigation systems can theoretically achieve 3D positioning through parallax calculation, it suffers from multiple limitations in practical clinical scenarios. From a hardware perspective, binocular infrared cameras require high-sensitivity infrared sensors and narrowband filtering modules to avoid ambient light interference, with hardware costs approximately 3-5 times that of a monocular camera with similar parameters. In clinical practice, binocular infrared technology relies on active infrared markers attached to the patient's skin for feature matching. However, these markers are prone to displacement due to patient breathing and positional changes, and the attachment process can cause discomfort. More importantly, the registration process requires the operator to manually select the correspondence between the marker and anatomical landmarks in the MRI image. This process is affected by variations in the operator's anatomical knowledge, leading to subjective biases in the registration results. In contrast, using a monocular camera instead of a binocular infrared camera can achieve systematic optimization while maintaining performance. A monocular camera can directly extract stable anatomical feature points from the patient's skin, eliminating the need for manual markers and offering higher accuracy. Furthermore, binocular cameras often exhibit image loss in single-view situations, while monocular cameras can more easily compensate for missing information through multi-frame image parallax synthesis technology, and have strong noise resistance. Combining these advantages, monocular cameras can effectively replace binocular cameras to improve system performance.

[0004] In the process of using medical images for real-world localization and tracking, there is a problem of multimodal image conflict. Image registration technology is an effective means to solve this problem. Specifically, image registration is an image processing technique that aligns different images of the same object acquired from different or the same sensors and perspectives in the same scene. The process requires geometric transformations, such as rigid, affine, projection, perspective, and global geometric transformations. After registration, the images will be located in a common coordinate system. Since point features can represent the geometric invariance between images, the image registration task can be simplified into a feature point-based registration problem. Point set registration can be divided into two categories: rigid registration and non-rigid registration. Rigid registration usually includes rotation, translation, and scaling, while the resolution of problems such as deformation and occlusion falls under the category of non-rigid registration. For localization of lesions in the human head, the following precautions should be taken: 1) The autonomous movement of the human head acquired by optical image acquisition can be approximated as rigid body movement. During MRI data acquisition, the position of the human head remains unchanged and there is no tissue deformation. Therefore, the human head model reconstructed from the acquired MRI slices can be approximated as a rigid body. 2) Considering the different acquisition times of MRI slice images and monocular data, facial expressions and muscle states cannot remain consistent. Therefore, there is a certain degree of local non-rigid deformation between the two. 3) At the same time, due to the inconsistency in the acquisition methods and equipment of the two modal data, there is significant global non-rigid deformation between them.

[0005] Considering the above requirements, compared with some existing image registration algorithms, the Thin Plate Spline (TPS) algorithm constructs a deformation field by minimizing the bending energy functional, achieving continuous and smooth deformation across the entire domain while strictly satisfying the marker point correspondence constraints. In addressing the control point dependency of B-splines, the gradient-driven limitations of Demons, and the neighborhood radius sensitivity of MLS, the TPS algorithm exhibits three core advantages: global deformation coordination, deformation rationality guaranteed by physical mechanisms, and adaptive capability for marker point errors. Furthermore, the deformation naturalness guaranteed by its energy minimization mechanism is highly consistent with the requirements of maintaining topological structure.

[0006] However, existing TPS-based medical image registration research falls short of the actual needs of this project in terms of processing objects, deformation scenarios, and modal combinations, and urgently requires targeted algorithm improvements and innovations.

[0007] In addition, existing non-rigid registration methods often fail to achieve accurate alignment when dealing with images with large initial pose differences, and their deformation models are susceptible to global misalignment interference, resulting in overcorrection.

[0008] If rigid registration is skipped and non-rigid registration is performed directly, the lack of global pose constraints will cause distortion of the overall spatial relationship of the anatomical structure. This not only fails to guarantee registration accuracy but may also introduce artifact errors. At the same time, in-table non-rigid algorithms have limitations in complex deformation modeling capabilities, poor cross-modal feature fusion, and are prone to convergence lag in large-scale displacement scenarios, making them difficult to adapt to practical application requirements.

[0009] Mathematical modeling methods, with their ability to quantitatively analyze complex motions, systematically trace the sources of errors, and accurately characterize multi-parameter coupling relationships, provide a reliable technical path to solving the aforementioned challenges. Therefore, a detailed analysis of head motion posture can be performed using mathematical modeling methods to achieve rigid registration and preliminary alignment of the global pose, providing a reliable spatial reference for subsequent non-rigid registration. Only then can fine-tuning of local deformations be performed to meet the project's stringent requirements for registration accuracy. Summary of the Invention

[0010] To address the problems in existing technologies, this invention provides a multi-modal image-based multi-level automatic localization method for head lesions, aiming to improve the accuracy of lesion location entry points, reduce reliance on human experience in clinical practice, and achieve high-precision, low-cost image navigation system localization functionality.

[0011] A multi-modal image-based automated localization method for head lesions includes the following steps:

[0012] 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.

[0013] 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.

[0014] 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;

[0015] 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.

[0016] Step 5: Initially align the optical image with the MRI 2D body surface image using the offset distance-angle analytical model;

[0017] 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.

[0018] Further, 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 ;

[0019] Step 4 specifically involves:

[0020] 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;

[0021] Step 4.2: Estimate the lateral flexion angle;

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

[0023]

[0024] ;

[0025] in, The angle of inclination of the straight line between the root of the nose and the tip of the nose;

[0026] Step 4.3: Yaw Calculation:

[0027] The calculation process is as follows:

[0028]

[0029]

[0030]

[0031] Among them, the actual distance between the left and right inner canthi under yaw conditions is calculated. Distance corresponding to the baseline 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;

[0032] Step 4.4: Determine the pitch level:

[0033] Quantitative analysis of head pitch angle

[0034]

[0035]

[0036]

[0037] 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.

[0038] Further, step 5 specifically involves:

[0039] 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;

[0040] 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.

[0041] Further, step 6 specifically involves:

[0042] 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. , ;

[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 In order 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] exist The centroid coordinates in the middle are defined as And satisfy the following formula:

[0055]

[0056] in, All are calculated from the area ratio. ;

[0057]

[0058] 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;

[0059] 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:

[0060] To exist and ,Right now Located on the public side Taking the above as an example, the vertices are respectively

[0061]

[0062] in, These are the control points excluding the vertices on the common edges. To change the scale, the calculation formula is:

[0063]

[0064] 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;

[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 excluding 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 centroids of control points in the target image are used in the calculation. Calculations are performed in place of control points;

[0077]

[0078] Therefore, global constraint terms It can be defined as:

[0079]

[0080] Therefore, the registration constraint equations are now corrected to:

[0081] .

[0082] The beneficial effects of this invention are as follows: It integrates rigid initial registration with a non-rigid precise registration algorithm based on thin-plate spline algorithm, optimizes the thin-plate spline algorithm from local and global levels through Delaunay triangulation, and adopts a multi-objective constraint processing mechanism to balance registration accuracy, computational efficiency and robustness. At the same time, it replaces the traditional dense scheme with a lightweight representation of 7 key location points to reduce interference, and achieves accurate positioning by combining head posture geometric analysis and initial alignment. It takes into account both theoretical reliability and engineering practicality, and provides strong technical support for multimodal image navigation. This invention effectively improves the accuracy of lesion location entry points through MRI image reconstruction and dimensionality reduction, selection of key facial location points, initial alignment of head posture and fine registration between virtual and real images, reduces the reliance on human experience in clinical practice, and realizes the positioning function of a high-precision and low-cost image navigation system. Attached Figure Description

[0083] Figure 1 A schematic diagram of a multi-level automatic localization method for head lesions;

[0084] Figure 2 Flowchart for MRI 3D body surface model reconstruction;

[0085] Figure 3 The diagram shows the lateral buckling curve in the offset distance-angle analytical model.

[0086] Figure 4 This is a diagram of the yaw curve in the offset distance-angle analytical model.

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

[0088] Figure 6 Flowchart for the initial registration of MRI two-dimensional body surface images and optical images;

[0089] Figure 7 This is a schematic diagram of the initial state without offset when the head is looking straight ahead after the initial registration results.

[0090] Figure 8 A schematic diagram of the combined attitude of right yaw, right flex, and pitch after the initial registration results;

[0091] Figure 9 This is a schematic diagram of the combined attitude of left yaw, left flexion, and nose-down after the initial registration results.

[0092] Figure 10 A diagram showing the distribution of lesion points and control points;

[0093] Figure 11 This refers to the external position of the circumcircle of the Delaunay triangle;

[0094] Figure 12 Comparison diagram of the second subdivision of the Delaunay triangle;

[0095] Figure 13 A comparison diagram of the transformed positional distribution of key location points;

[0096] Figure 14 This represents the initial state registration deviation;

[0097] Figure 15 Registration bias for state 1;

[0098] Figure 16 Registration bias for state 2;

[0099] Figure 17 Comparison of the distribution of lesion sites after location change;

[0100] Figure 18 The evaluation index results are for different deformation algorithms. Detailed Implementation

[0101] The present invention will now be described in detail with reference to the accompanying drawings. Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention. The directional terms such as left, center, right, top, and bottom in the embodiments of the present invention are only relative concepts or referenced to the normal use state of the product, and should not be considered restrictive.

[0102] A multi-modal image-based, multi-level automatic localization method for head lesions, such as... Figure 1 As shown, it includes the following steps:

[0103] Step 1: Acquire head MRI slice images and head optical images containing the lesion, such as... Figure 2As shown, through operations such as noise filtering, threshold segmentation, preservation of the largest island, contour stitching, image interpolation, volume rendering, and surface smoothing, and based on head MRI slice images, a three-dimensional MRI head surface model of the outer skin of the head and a three-dimensional lesion model of the lesion site are constructed.

[0104] Step 2: The MRI 3D head surface model is successively adjusted to a preset pose sequence using Euler angle transformation. For each pose, corresponding MRI 2D body surface images are extracted at a specified step size. These MRI 2D body surface images are combined into a dataset of MRI 2D body surface images after dimensionality reduction of the MRI 3D head surface model. In the setting of angle transformation parameters, a moderate value should be selected for the single rotation angle: too small an angle will significantly increase the experimental workload and reduce efficiency; too large an angle may affect the accuracy of subsequent pose analysis due to the wide pose sampling interval. Therefore, the optimal range of single rotation angles needs to be determined through pre-experimentation.

[0105] Step 3: Manually label the positions of 7 key points in the head MRI image dataset and head optical images, 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 the nose, left corner of the mouth, and right corner of the mouth; the coordinates of the 7 key points are shown in the table below:

[0106] Table 1 Key Location Points and Their Coordinates

[0107]

[0108] Step 4: Based on the geometric position estimation method, spatial calculation of head posture in the three directions of yaw, flexion, and pitch is performed on the images in the MRI two-dimensional body surface image dataset. Combined with the corresponding rotation angle of the MRI three-dimensional head surface model, an offset distance-angle analytical model is established; specifically:

[0109] Step 4.1: Considering that facial expressions are easily altered unconsciously, and the corners of the mouth vary relatively greatly, the pose estimation does not consider the coordinate positions of the corners of the mouth on both sides; instead, it introduces the position of the root of the nose above the midpoint of the line connecting the inner corners of both eyes. As an auxiliary reference, the degree of offset in each direction is calculated using the frontal view as the standard;

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

[0111] Using the anatomically neutral facial position as the baseline for calculating left and right lateral flexion offsets in multimodal data, where the coronal plane of the face is perpendicular to the ground, this state is defined as the zero reference position for rotational offset, i.e., the degree of offset is 0°. Under this standard pose, a Cartesian coordinate system is established to calculate the tilt angle of the line connecting the nasal root and the nasal tip. This is normalized to assess the degree of rotational offset; however, it is worth noting that since the tilt of the upper body causes the head to shift to be parallel to the ground, the absolute value of the tilt angle is... The left and right lateral flexion directions of the head are distinguished by using the vertical coordinates of the left and right outer canthi; therefore, the specific calculation of the degree of lateral flexion is as follows:

[0112]

[0113] ;

[0114] in, The angle of inclination of the straight line between the root of the nose and the tip of the nose;

[0115] Step 4.3: As Figure 4 As shown, the yaw rate is calculated:

[0116] The quantification of yaw state is achieved through the feature point distance ratio method. However, since the horizontal displacement of the head can vary significantly, in the case of a completely profile view, key position points on the other side will be completely obscured and cannot be obtained, causing traditional pose calculation methods based on key position points to fail. Therefore, most current research has not considered such extreme displacement scenarios. To solve this problem, we expand the applicable scope to the outer canthus view. Specifically, we acquire two-dimensional images of the target in the reference position (front view) and in the yaw state, respectively, extract the feature points of the left and right inner canthi in the images, and calculate the actual distance between the left and right inner canthi in the yaw state. Distance corresponding to the reference state ratio This ratio can be used as a quantitative indicator to establish a yaw motion determination model based on a proportional coefficient. The numerical value allows for accurate assessment of the target's horizontal motion state, while the change in the nose tip's horizontal coordinate value measures the left and right motion direction, enabling non-contact measurement of yaw displacement in two-dimensional space. The calculation process is shown below:

[0117]

[0118]

[0119]

[0120] 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 their initial states are respectively...

[0121] Step 4.4: As Figure 5 As shown, the pitch degree is determined:

[0122] For the quantitative assessment of head pitch posture, a computational model based on feature point vector analysis is constructed. By simultaneously acquiring two-dimensional images of the target in a baseline frontal view and in a pitch state, the difference in the ordinates of the nasal root and nasal tip under key feature point coordinates is calculated using these coordinates. Difference from real-time status This method constructs an attitude calculation equation using a ratio relationship, transforming the coordinate components of two vectors into pitch values. By comparing the changes in feature vectors under different states, it achieves quantitative analysis of head pitch angle.

[0123]

[0124]

[0125]

[0126] Considering the above formula To address the inherent flaws in quantifying the head-down process, a normalized dynamic offset is introduced. By constructing a dynamic normalized model based on real-time calculation of the nasal root coordinates, the geometric characteristics of head spatial rotation are accurately mapped. A dimensionless displacement ratio is used instead of the traditional absolute distance metric, effectively decoupling the many-to-one mapping relationship between distance and angle during nonlinear deformation, theoretically eliminating the quantification ambiguity of the original model in large-angle head-down scenarios. Therefore... Optimized to:

[0127]

[0128] in, , For the weights of the original terms, The smaller the value, the better the correction effect; calculate the difference in ordinates between 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. Indicates the normalized dynamic offset;

[0129] Step 5: Initially align the optical image with the MRI 2D body surface image using the offset distance-angle analytical model; In the biomechanical analysis of head movement, if muscle synergy and joint biomechanical effects are temporarily disregarded, the three spatial motion modes of head—yaw, lateral flexion, and pitch—can be approximated as independent rigid body movements; such as Figure 6 As shown, 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 11As shown, this determines the positional relationships of all points except the key point relative to the key point: 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; specifically:

[0135] 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.

[0136] 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

[0137] Step 6.3: Based on the lesion point The influence of key points varies depending on the different positional relationships with the Delaunay triangle; therefore, the deformation function in the thin-plate spline function needs to be modified. and constraint functions This is to perform local structural constraints; specifically:

[0138] The original TPS algorithm formula is:

[0139]

[0140]

[0141]

[0142]

[0143] in, , It is a linear matrix. , For data fitting error, For bending energy, It is a radial basis function; due to the complex facial deformation, Introduce second-order terms for calculation; Representing two-dimensional Euclidean space, Represents the coordinates of key locations. The corresponding weighting coefficients, Indicates the coordinates of key points in the source image. Indicates the coordinates of the lesion point. This represents the independent variable in the radial basis function;

[0144]

[0145] 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. ;

[0146] Only the source image points were constrained by the deformation function The distance difference between the transformed point and the target image point, in the context of determining the specific location of the lesion point on the optical image in this paper, neglecting the correlation between the control point and the global structure, is prone to local excessive distortion, affecting the lesion point localization accuracy, and thus failing to form an effective constraint; for example, Figure 12 As shown, to ensure that the spatial structure controlled by the key location points remains unchanged before and after the transformation, the Delaunay triangulation method is used to establish the location map of the key location points. This method can ensure that the triangles formed do not overlap, can cover the entire plane, and that each point is not inside the circumcircle of the triangle that does not contain that point.

[0147] Step 6.3.1: When the lesion point When within the nearest triangle, increase the number of nearest triangle vertex pairs. The degree of influence is considered, and secondly, the position of the center of gravity coordinates is used as a constraint to ensure that the position of the center of gravity coordinates remains relatively unchanged before and after the transformation.

[0148] exist The centroid coordinates in the middle are defined as And satisfy the following formula:

[0149]

[0150] in, All are calculated from the area ratio. ;

[0151]

[0152] 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; This indicates hyperparameters (manually adjusted by the user). Represents radial basis functions. Represents the vertices of a triangle. Represents the local spatial structure constraint function;

[0153] Step 6.3.2: When the lesion point When on the nearest triangle (possibly more than one) 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:

[0154] To exist and ,Right now Located on the public side Taking the above as an example, the vertices are respectively ;

[0155]

[0156] in, These are the control points excluding the vertices on the common edges. To change the scale, the calculation formula is:

[0157]

[0158] 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;

[0159] 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.

[0160]

[0161] in, These are the control points excluding the nearest triangle vertex. This represents the number of nearest triangles; when n=2, ; Indicates subscript The set, , , 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] Therefore, global constraint terms It can be defined as:

[0172]

[0173] Building upon the initial rigid registration, a geometric topological attribute fusion strategy is further introduced to construct a multi-level local and global constraint system for lesion point location changes, and optimization iterations are carried out based on thin-plate spline functions. At the local structural constraint level, the spatial positional relationship between the lesion point and its nearest neighbor triangular unit is accurately determined, and deformation functions are corrected for different geometric configurations, while local spatial structural constraint conditions are constructed. At the global structural constraint level, a constraint model is established based on the joint mapping relationship between the centroid coordinates of the lesion point relative to the control point and the centroid coordinates of the control point, and a surface torsion energy constraint function is introduced to form a composite constraint system. Thus, the registration constraint equation is corrected as follows:

[0174] .

[0175] Example verification and result analysis

[0176] Given that single-person data is representative of typical samples in terms of methodology, subsequent analyses use single-person data as standardized research objects. This analytical framework can be directly extended to other sample analysis scenarios through parameterized adjustments, thereby ensuring the universality and scalability of the research method.

[0177] The offset distance-angle analytical model was used to register three different head states in the optical image. The initial state was a frontal view without any angular offset. State 1 was a combined movement state of right flexion, right yaw, and head tilt; State 2 was a combined movement state of left flexion, left yaw, and head tilt. After the optical image was initially aligned with the MRI 2D body surface image, a comparison of the coordinate distribution of key points between the MRI 2D body surface image and the head optical image was shown, as follows. Figures 7 to 9 As shown; in order to provide a unified reference framework for the coordinate data of different samples and facilitate comparison, the coordinates of each key location point are normalized using the nose tip coordinate as the standard.

[0178] from Figures 7 to 9As can be seen, after initial registration, there are still different degrees of deviation among the control points, so further precise deformation correction is required. Precise registration is performed using step (6). To verify the registration effect of the proposed method, the original algorithms B-spline, Demons, MLS, and TPS are selected as comparison methods, and the three sets of experimental data are processed. 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 Figure 13. At the same time, the registration error of each key position point is quantified and calculated, and the relevant results are shown in Figures 14 to 16.

[0179] As can be seen from the results in Figures 14 to 16, the registration accuracy of the proposed method is better in all three sets of experimental data and is more in line with the actual application requirements. In order to further accurately evaluate the registration deviation between key location points, this study uses mean absolute error (MAE) for quantitative analysis, and the relevant results are shown in Table 2.

[0180]

[0181] in, This represents the coordinates of key MRI locations after transformation.

[0182] Table 2 MAE Calculation Results

[0183]

[0184] As shown in Table 2, the algorithm proposed in this application exhibits the smallest registration deviation at both single points and overall registration deviation across all key locations, achieving an accuracy improvement of approximately 97.45% compared to the original algorithm. However, the core objective of this research is not merely to optimize the registration performance of key locations; its starting point lies in accurately quantifying the positional changes of lesion points. To this end, this application performed positional transformations on all pixels in the entire image, and the resulting spatial distribution is shown in Figure 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 represent a comparison of the processing effects of different algorithms.

[0185] The goal of this study is to make the positional distribution of each point in the transformed image as close as possible to that of the target image; by Figure 17As can be seen, in the initial state, except for the MLS algorithm which has a poor transformation effect, the other algorithms perform similarly and all have a high similarity to the target image. In state 1, only the algorithm proposed in this paper has a high degree of approximation to the target image. In state 2, the MLS algorithm and the method proposed in this paper have the best performance. 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. The calculation results are shown as follows: Figure 18 .

[0186]

[0187]

[0188]

[0189] in, Represents the edge probability distribution of the image. This represents the total number of pixels in the image. The pixel values ​​of the RGB image. These are the pixel values ​​of the transformed MRI two-dimensional body surface image. For joint probability, Marginal probability, These are the pixel position coordinates. , For the corresponding image in the first The average pixel value.

[0190] The quantitative results of MI, RMSE, and CC in the three sets of data are as follows: Figure 18 The results show that the registration performance of the proposed algorithm is significantly better than that of the comparison methods. Specifically, the MI values ​​of B-spline, Demons, MLS, TPS (original) and the proposed algorithm 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 proposed algorithm has the highest similarity to the target image; the RMSE values ​​are 6.649, 6.773, 6.479, 7.129, and 4.131, respectively, further confirming that the proposed algorithm has the best registration accuracy and fully verifying its performance advantage.

[0191] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention 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 In order 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