An Automatic Generation Method for Interactive 2D and 3D Medical Image Registration Parameters

Abstract

The present invention relates to a method for automatically generating registration parameters for interactive 2D and 3D medical images, which solves the defects of poor search accuracy and low search efficiency of initial registration parameters compared with the prior art. The present invention includes the following steps: loading three-dimensional image data; implementing three-dimensional model rendering and reconstruction on a window, performing two-dimensional mapping on the three-dimensional model, and aligning and displaying the image after real-time mapping with the X-ray image to be registered; loading the two-dimensional image to be registered and drawing it on the window in the form of a two-dimensional texture; dragging the mouse to align the three-dimensional model with the two-dimensional image to obtain appropriate registration parameters. The present invention comprehensively considers the difference in the image quality of the three-dimensional volume data drawn on the window after maximum intensity projection and the image quality of the X-ray image to be registered, can complete the alignment of the two images with high quality, can obtain a set of appropriate initial registration parameters and print and output them on the console.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and in particular to an interactive method for automatically generating 2D and 3D medical image registration parameters. Background Art

[0002] With the progress of medical imaging devices, for the same patient, images containing accurate anatomical information such as CT and MRI can be acquired; at the same time, images containing functional information such as SPECT can also be acquired. However, diagnosing by observing different images requires spatial imagination and the subjective experience of doctors. Adopting the correct image registration method can accurately fuse various types of information into the same image, enabling doctors to more conveniently and precisely observe lesions and structures from various angles. At the same time, by registering dynamic images acquired at different times, the changes in lesions and organs can be quantitatively analyzed, making medical diagnosis, surgical planning, and radiotherapy planning more accurate and reliable. Therefore, registering different-modal and different-dimensional medical images has become a hot and frontier topic in current medical image informatics research.

[0003] In the application process of medical image registration processing, finding a set of optimal initial registration parameters can bring great convenience to the subsequent registration process. The initial registration of medical images, that is, the rough registration of medical images, is to roughly align the images to facilitate the next fine registration step. However, how to find a suitable set of initial registration parameters has always been a difficult problem. Therefore, it is of great significance to seek a suitable method to optimize the initial registration parameters with high quality and high efficiency.

[0004] Previous research on 2D and 3D medical image registration mainly divides the main registration methods into two categories: feature-based registration methods and intensity-based registration methods.

[0005] Feature-based registration method. Feature-based registration method first preprocesses the images to be registered, that is, the process of feature extraction. Then, the features extracted are used to complete the matching between the features of the two images. Since there are many types of features that can be used in images, various feature-based methods have emerged: (1) Point feature-based registration. Point features are one of the most commonly used image features in image registration, and are divided into two types: external feature points and internal feature points. (2) Line feature-based registration. Line segments are another feature that is easy to extract in images. The Hough transform is an effective method for extracting lines in images. The Hough transform can transform a curve or a line with a given shape in the original image into a point position in the transformed spatial domain. It makes all the points on the curve or line with a given shape in the original image converge to a certain point position in the transformed domain to form a peak. In this way, the problem of detecting lines or curves in the original image becomes the problem of finding peak points in the transformed space. Correctly establishing the corresponding relationship between the line segments extracted from the two images respectively remains the key point and difficulty of this method. Considering the slope of the line segment and the positional relationship of the endpoints comprehensively, a histogram of these information indicators can be constructed, and the matching of the line segments can be achieved by finding the aggregation bundle of the histogram. (3) Contour and curve feature-based registration. (4) Surface feature-based registration. The most typical algorithm is the "head and cap method", that is, a surface model called "head" is extracted from the image, and the point set extracted from the contour of another image is called: "cap". The point set of the "cap" is transformed onto the "head" using a rigid body transformation or a selective affine transformation, and then an optimization algorithm is used to make the mean square distance from each point of the "cap" to the surface of the "head" the smallest.

[0006] Gray-based registration method. It directly uses the gray information of the image for registration, thus avoiding the errors caused by segmentation. Therefore, it has the characteristics of high accuracy, strong robustness, and can achieve automatic registration without preprocessing. The registration methods based on gray processing mainly include: one category is to directly calculate representative elements such as ratios and directions through the gray levels of the image; the other category is to use all the gray information during the registration process. The first method is represented by the moment and principal axis method, and the second method is generally called voxel formality. (1) Moment and principal axis method: The moment and principal axis method means first calculating the centroids and principal axes of two images using the principle of the mass distribution of an object in classical mechanics, and then achieving the registration of the two images through transformations such as translation and rotation. Using this method, the image can be modeled as a point distribution in an elliptical region. Such a distribution can be described by the first-order and second-order moments of the positions of these points. This method is sensitive to data loss and requires the entire object to appear completely in both images. Generally speaking, the registration accuracy is poor, so currently it is more used for rough registration to initially align the two images and reduce the search steps of the subsequent main registration method. (2) Voxel similarity method: The voxel similarity method is a type of method that has been studied more currently. Since it uses all the gray information in the image, this method is generally relatively stable and can obtain quite accurate results. Another advantage of this method is that it is completely automatic and does not require special preprocessing. However, because this method requires a large amount of complex calculations, it has only been applied in practice in recent years.

[0007] Therefore, in the field of 2D and 3D medical image registration, the registration method based on gray information has become a popular research direction at present, and its accuracy has exceeded that of the feature-based registration method. During the registration process based on gray information, how to find a set of suitable initial registration parameters in a short time to facilitate subsequent fine registration has become a technical problem that urgently needs to be solved. Summary of the Invention

[0008] The purpose of the present invention is to solve the defects of poor search accuracy and low search efficiency of the initial registration parameters in the prior art, and provide an interactive automatic generation method for 2D and 3D medical image registration parameters to solve the above problems.

[0009] In order to achieve the above purpose, the technical solution of the present invention is as follows:

[0010] An interactive automatic generation method for 2D and 3D medical image registration parameters, including the following steps:

[0011] 11) Loading three-dimensional image data: Obtaining the three-dimensional image data of the interactive medical image to be registered;

[0012] 12) Implement 3D model rendering and reconstruction on the window, perform 2D mapping on the 3D model, and align and display the real-time mapped image with the X-ray image to be registered;

[0013] 13) Load the 2D image to be registered and draw it on the window in the form of a 2D texture: Read the corresponding 2D image, convert it into 2D texture data using the method of texture mapping, map the texture pixels in the texture space to the pixels in the window space, and draw the image on the rendering window;

[0014] 14) Drag the mouse to align the 3D model with the 2D image to obtain appropriate registration parameters: Use mouse interaction to drag the 3D model for rotation and translation operations, so that the 3D model can be aligned with the 2D image visually, and then automatically generate a set of registration parameters.

[0015] The implementation of 3D model rendering and reconstruction on the window, performing 2D mapping on the 3D model, and aligning and displaying the real-time mapped image with the X-ray image to be registered includes the following steps:

[0016] 21) Process the 3D image data into data suitable for the OpenGL rendering pipeline: For the input sequence of several 3D slice images, i.e., DICOM data, write a Python script using the pydicom library of Python to convert the DICOM data into a binary file suitable for the system, and then store the remaining DICOM data in the binary file using a C++ container;

[0017] 22) Load the data into the real-time rendering and drawing pipeline: Store the binary data into the rendering pipeline through the corresponding functions of OpenGL and a fixed rendering process;

[0018] 23) Drawing of the 2D window: Project the corresponding 3D CT image using the maximum intensity projection algorithm and display it on the window in the form of volume rendering;

[0019] 24) Design the rotation interaction method of the 3D model using the arcball algorithm: Adopt the idea of arcball, store the information of each mouse interaction model rotation into a quaternion, and control the model rotation by converting the quaternion into the form of Euler angles and rotation matrices;

[0020] 25) Set the interaction method for translating the 3D model with the mouse: Set the translation interaction method of the model through the window coordinate calculation of QT and mouse interaction;

[0021] 26) Use the interaction result to control the model for real-time rendering and display.

[0022] The processing of the 3D image data into data suitable for the OpenGL rendering pipeline includes the following steps:

[0023] 31) Use a Python script to read DICOM data in the release folder; open the cmd console and convert the DICOM data into a binary bin file in a specific format through commands entered in the cmd for subsequent use in the rendering process, and store the binary bin file in the corresponding folder;

[0024] 32) Use a for loop to iterate through all slice data, use the sprintf_s function in C++ to format and output the data path to a string, read a series of data under the path, and use the fopen_s function in C++ to open the binary file;

[0025] 33) Use the fread_s function in C++ to read the information in the binary bin file, which includes the position of the slice in the image sequence, the spacing of the slice image pixels along the X-axis direction, the spacing of the slice image pixels along the Y-axis direction, the width of the slice, and the height of the slice, store the remaining DICOM data in a container, and adjust the size of the container to prevent data overflow;

[0026] 34) After reading the data, sort all slice data based on the position of the slices;

[0027] 35) Load the DICOM-formatted data of the three-dimensional volume data through the information in the read binary file.

[0028] Loading the data into the pipeline of real-time rendering and drawing includes the following steps:

[0029] 41) The rendering pipeline calls a function to generate a three-dimensional texture object and binds the three-dimensional texture in the pipeline;

[0030] 42) Perform three-dimensional texture mapping. The parameters of the mapping function are the width and height of each two-dimensional slice image and the depth of the data, and texture filtering is performed on the mapped three-dimensional texture;

[0031] 43) Generate and set vertex data attributes and bind the vertex data in the rendering pipeline;

[0032] 44) Write OpenGL shaders, write vertex shaders. For each vertex Vertex sent to the GPU, vertex shading is performed once. Its function is to transform the three-dimensional coordinates of each vertex in the virtual space into two-dimensional coordinates displayed on the window and carry depth information for the z-buffer; write fragment shaders to calculate the color and other attributes of each pixel; compile and link the written shaders, and delete the shaders after completion;

[0033] 45) Pass the read 3D texture into the shader, and through texture sampling, perform texture drawing on the 3D model, and end the loading.

[0034] The drawing of the 2D window includes the following steps:

[0035] 51) Modify the ray casting algorithm based on the principle of maximum density projection to implement a shader with the maximum density projection function;

[0036] 52) Use the GetUniformLocation function of OpenGL to obtain the location labels of the maximum density value, texture loading, camera position parameters, and the model view projection transformation matrix MVP matrix in the maximum density projection shader;

[0037] 53) Load the position parameters of the volume data in space one by one through the corresponding location labels, and then load the 3D texture map;

[0038] 54) Call the maximum density projection shader to project the loaded 3D volume data, so as to draw the model after maximum density projection on the 2D window.

[0039] The method for designing the rotation interaction mode of the 3D model using the arcball algorithm includes the following steps:

[0040] 61) Map the 2D window coordinates after mouse interaction. Imagine a unit hemisphere located at the center of the window and adjust the range of the 2D window coordinates to the interval [-1....1].

[0041] Formula: pt.x = (pt.x * AdjustWidth) - 1.0f, pt.y = 1.0f - (pt.y * AdjustHeight); where pt is the defined 3D coordinate, AdjustWidth is the scaling factor of the width, and AdjustHeight is the scaling factor of the length; AdjustWidth = 1.0f / ((NewWidth - 1.0f) * 0.5f), NewWidth and NewHeight are the width and height of the 2D window; AdjustHeight = 1.0f / ((NewHeight - 1.0f) * 0.5f);

[0042] 62) The coordinates of two points in the window are normalized to two points on the hemispherical surface through the mapping formula. If the two-dimensional window coordinates are not on the unit hemisphere centered at the window center, the two-dimensional coordinates are scaled to the hemisphere, and the scaling factor is set to norm = 1.0 / FuncSqrt(length), where length = (pt.x * pt.x) + (pt.y * pt.y); the two-dimensional coordinates are mapped and converted into two points in the hemispherical space. If the two-dimensional coordinates are on the unit hemisphere, the Z-direction coordinate pt.z is calculated according to the X-direction axis coordinate pt.x and the Y-direction coordinate pt.y, and the calculation formula is pt.z = FuncSqrt(1.0f - length), where length = (pt.x * pt.x) + (pt.y * pt.y);

[0043] 63) Set that when the left mouse button is pressed, a starting point is generated, and the coordinate value at the time of pressing is mapped into three-dimensional coordinates and stored in vector form through the previous coordinates; when the mouse is released, an ending point is generated, and the coordinate value at the time of release is also mapped into three-dimensional coordinates and stored in vector form to obtain the direction vector of the starting point and the ending point during the rotation interaction process;

[0044] 64) Set a combined quaternion q = [v, w] = [x, y, z, w] which consists of two parts. One is the scalar w, which is equal to cosθ / 2, where θ is the rotation angle, and the other is the vector v, which is equal to sinθ / 2 times the vector along the rotation axis; the result of the operation of two quaternions is the result of their rotational combination, so the rotational combination operation is represented by quaternion cross product; the cross product of two rotation vectors records the direction of the rotation axis, and the dot product of two vectors records the rotation angle;

[0045] 65) Call the corresponding function to save the previously calculated cross product and dot product in the form of a quaternion. The two rotation direction vectors are OP1 and OP2. First, find the inner product s of the two vectors: s = OP1 · OP2, then find the outer product v of the two vectors: v = OP1 × OP2. Denote the quaternion q = [s, v], and normalize it. At this time, q is the rotation quaternion;

[0046] 66) Convert the quaternion into a rotation matrix, and the formula is TM = QuatToMatrix(q). The TM matrix is the corresponding rotation matrix, and the rotation matrix acts on the model transformation matrix rendered by OpenGL to control the model to perform rotation operations;

[0047] 67) Convert the quaternion into Euler angles, and the formula is Rotate = QuatToeulerAngles(q). Rotate is the corresponding Euler angle, which can display the rotation angles of the model around the three coordinate axes during the real-time interaction process;

[0048] 68) Perform a rotation interaction operation on the model through the mouse to obtain the quaternion for each rotation interaction; continuously store and display the rotation information through quaternion multiplication.

[0049] The described interactive method for setting the mouse to pan the 3D model includes the following steps:

[0050] 71) When the right mouse button is pressed, record the current 2D window coordinates, i.e., the starting point. When the mouse stops sliding, record the 2D window coordinates at this time, i.e., the ending point;

[0051] 72) Call the QPointF function of QT to calculate the differences in the X-axis direction and Y-axis direction between the two 2D window coordinates, thereby determining the translation of the model in the X-axis direction and the translation in the Y-axis direction;

[0052] 73) Call the mouse wheel mechanism of QT to control the translation transformation effect of the model on the Z-axis. When the mouse wheel slides up, the model moves in the positive Z-axis direction. When the mouse wheel slides down, the model moves in the negative Z-axis direction;

[0053] 74) Substitute the translation parameters of the model in the X-axis direction, Y-axis direction, and Z-axis direction into the translation matrix, and substitute the calculated translation matrix into the model transformation matrix, thereby controlling the translation operation of the model in real time.

[0054] The described method of using the interaction result to control the model for real-time rendering and display includes the following steps:

[0055] 81) Convert the quaternion that saves the rotation information into a rotation matrix, and convert the vector that saves the translation information into a translation matrix;

[0056] 82) Set the model transformation matrix as the product of the translation matrix, rotation matrix, and scaling matrix; the model transformation matrix acts on the model to control the transformation of the model in the world coordinate system, and converts the object from the model coordinate system to the world coordinate system;

[0057] 83) Set the view matrix to convert the object from the world coordinate system to the view coordinate system;

[0058] 84) Set the projection matrix to perform a projection transformation on the object model, and convert the object from the view coordinate system to the clip coordinate system;

[0059] 85) Set the viewport transformation to convert the object from the clip coordinate system to the window coordinate system, thereby rendering the drawing result on the window in real time.

[0060] Beneficial effects

[0061] An automatic generation method for interactive 2D and 3D medical image registration parameters according to the present invention comprehensively considers, compared with the prior art, the difference in the image quality of the three-dimensional volume data drawn on the window after maximum intensity projection and the X-ray image quality to be registered, can complete the alignment of the two images with relatively high quality, can obtain a set of appropriate initial registration parameters and print and output them on the console.

[0062] The three-dimensional CT data is drawn on a two-dimensional window through maximum intensity projection, and the X-ray image to be registered is drawn on the two-dimensional window in the form of a two-dimensional texture. By interacting with the three-dimensional model with the mouse, the model is aligned with the X-ray image under visual intuitive observation. On the one hand, the search accuracy of the initial registration parameters is improved, and on the other hand, the search efficiency of the initial registration parameters is effectively improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is the sequence diagram of the method of the present invention;

[0064] Figure 2 is the existing two-dimensional hip X-ray image;

[0065] Figure 3 is the existing corresponding two-dimensional mask image;

[0066] Figure 4 is the rendering display diagram on the two-dimensional window after loading the three-dimensional CT data;

[0067] Figure 5 is the display diagram after the three-dimensional model is interacted with translation and rotation;

[0068] Figure 6a is the diagram of the two-dimensional X-ray photo and the three-dimensional model being simultaneously rendered and displayed on the window;

[0069] Figure 6b is the alignment display diagram of the two-dimensional X-ray photo and the three-dimensional model on the window after interacting with the mouse. DETAILED DESCRIPTION OF THE INVENTION

[0070] To have a further understanding and recognition of the structural features and achieved effects of the present invention, the following is a detailed description with reference to the preferred embodiments and the accompanying drawings:

[0071] As Figure 1 shown, an automatic generation method for interactive 2D and 3D medical image registration parameters according to the present invention includes the following steps:

[0072] The first step is to load three-dimensional image data: load the three-dimensional image data of the interactive medical image to be registered.

[0073] Step 2: Implement 3D model rendering and reconstruction on the window, perform 2D mapping on the 3D model, and align and display the real-time mapped image with the X-ray image to be registered.

[0074] Here, the maximum intensity projection is used to render and display the model on the window in a short time to achieve interactive operation of the 3D model; through manual operation, the rotation and translation parameters of the 3D model are obtained in real time and printed out, facilitating the subsequent registration and alignment process. Through interactive operation, the rotation of the 3D model is subjectively controlled, which is convenient and fast. Controlling the pose of the 3D model to achieve alignment with the X-ray image is conducive to quickly obtaining a set of suitable initial registration parameters in the subsequent process. At the same time, the difficulty of obtaining registration parameters through interactive registration lies in the rendering and interaction of the 3D model. To display the 3D model and the 2D image on the 2D window simultaneously, it is necessary to call OpenGL for 3D model rendering, combine the idea of arcball, use quaternion rotation to control the rotation transformation of the model, and at the same time, a series of rendering pipelines and shaders need to be called, and the interactive parameters need to be displayed accurately in real time. The specific steps are as follows:

[0075] (1) Process the 3D image data into data suitable for the OpenGL rendering pipeline: For the input sequence of several 3D slice images, i.e., dicom data, write a python script using the pydicom library of python to convert the dicom data into a binary file suitable for the system, and then use the container of C++ to store the remaining dicom data in the binary file.

[0076] A1) Use a python script to read the dicom data in the release folder; open the cmd console, and convert the dicom data into a binary bin file in a specific format through commands entered in the cmd for subsequent rendering processes, and store the binary bin file in the corresponding folder;

[0077] A2) Use a for loop to traverse all slice data, use the sprintf_s function of C++ to format and output the data path to a string, read a series of data under the path, and use the fopen_s function of C++ to open the binary file;

[0078] A3) Use the fread_s function of C++ to read the information in the binary bin file, including the position of the slice in the image sequence, the spacing of the slice image pixels along the X-axis direction, the spacing of the slice image pixels along the Y-axis direction, the width of the slice, and the height of the slice, store the remaining dicom data in a container, and adjust the size of the container to prevent data overflow;

[0079] A4) After reading the data, sort all slice data based on the position of the slice;

[0080] A5) Load the DICOM-formatted data of the three-dimensional volume data by reading the binary file information.

[0081] (2) Load the data into the pipeline of real-time rendering and drawing: Store the binary data into the rendering pipeline through the corresponding functions of OpenGL and the fixed rendering process.

[0082] B1) The rendering pipeline calls a function to generate a three-dimensional texture object and binds the three-dimensional texture in the pipeline;

[0083] B2) Perform three-dimensional texture mapping. The parameters of the mapping function are the width and height of each two-dimensional slice image and the depth of the data, and texture filtering is performed on the mapped three-dimensional texture;

[0084] B3) Generate and set the vertex data attributes and bind the vertex data in the rendering pipeline;

[0085] B4) Write the shaders of OpenGL, write the vertex shader. For each vertex Vertex sent to the GPU, vertex shading is performed once. Its function is to transform the three-dimensional coordinates of each vertex in the virtual space into two-dimensional coordinates displayed on the window and carry the depth information for the z-buffer; write the fragment shader to calculate the color and other attributes of each pixel; compile and link the written shaders, and delete the shaders after completion;

[0086] B5) Pass the read three-dimensional texture into the shader, and perform texture drawing of the three-dimensional model through texture sampling to end the loading.

[0087] (3) Drawing of the two-dimensional window: Project the corresponding three-dimensional CT image using the maximum intensity projection algorithm and display it in the form of volume rendering on the window.

[0088] Visualize the structures with high gray values in the volume data using the maximum intensity projection; first determine the positions of the light source points in space and the volume data, and then emit virtual rays outward from the light source points. The intersection points of the rays with the plane determine the positions of the pixels in the MIP image; emit a series of virtual rays from the back to the front in the slice direction of the volume data and project them onto a two-dimensional plane. When the rays pass through the volume data, equidistant sampling is performed on the rays, and the maximum value of the attributes among the sampling points is taken as the output of the ray. The scalar value of the position coordinates of each sampling point can be calculated through interpolation; the color value of the window pixel corresponding to the ray can be obtained through color mapping of the output value, and the final result is drawn on the projected two-dimensional plane.

[0089] C1) Modify the ray casting algorithm based on the principle of maximum intensity projection to implement a shader with the maximum intensity projection function;

[0090] C2) Use the GetUniformLocation function of OpenGL to obtain the location labels of the maximum density value in the maximum density projection shader, the location label of texture loading, the location label of the camera position parameter, and the location label of the model view projection transformation matrix MVP matrix;

[0091] C3) Load the position parameters of the volume data in space one by one through the corresponding location labels, and then load the three-dimensional texture map;

[0092] C4) Call the maximum density projection shader to project the loaded three-dimensional volume data, so as to draw the model after maximum density projection on the two-dimensional window.

[0093] The maximum density projection can accurately simulate the original volume data in a short time and reconstruct the three-dimensional model through volume rendering; by drawing the three-dimensional model on the two-dimensional window through the maximum density projection, the time-consuming of registration projection can be greatly reduced, and the practicability of interactive registration can be enhanced. Due to the defect that it is difficult to interact with the three-dimensional volume data generated by calling the maximum density projection in the medical image processing toolkits ITK and VTK, the present invention first converts the data type, and through the corresponding functions of OpenGL, performs the maximum density projection on the volume data and then displays it on the window. During the process, complex rendering pipeline binding and preprocessing, as well as shader writing are required, but the speed of generating the rendered model during the calling process is very fast.

[0094] (4) Design the rotation interaction method of the three-dimensional model using the arcball algorithm: Adopt the idea of arcball, store the information of each mouse interaction model rotation into a quaternion, and convert it into the form of Euler angles and rotation matrix through the quaternion to control the model rotation.

[0095] D1) Map the two-dimensional window coordinates after mouse interaction, imagine a unit hemisphere located at the center of the window, and adjust the range of the two-dimensional window coordinates to the interval [-1....1];

[0096] Formula: pt.x = (pt.x * AdjustWidth) - 1.0f, pt.y = 1.0f - (pt.y * AdjustHeight); where pt is the defined three-dimensional coordinate, AdjustWidth is the scaling factor of the width, and AdjustHeight is the scaling factor of the length; AdjustWidth = 1.0f / ((NewWidth - 1.0f) * 0.5f), NewWidth and NewHeight are the width and height of the two-dimensional window; AdjustHeight = 1.0f / ((NewHeight - 1.0f) * 0.5f);

[0097] D2) Normalize the coordinates of two points in the window to two points on the hemisphere using the mapping formula. If the two-dimensional window coordinates are not on the unit hemisphere centered at the window center, scale the two-dimensional coordinates to the hemisphere, and set the scaling factor to norm = 1.0 / FuncSqrt(length), where length = (pt.x * pt.x) + (pt.y * pt.y); Map the two-dimensional coordinates to two points in the hemisphere space. If the two-dimensional coordinates are on the unit hemisphere, calculate the Z-direction coordinate pt.z based on the X-direction axis coordinate pt.x and the Y-direction coordinate pt.y, and the calculation formula is pt.z = FuncSqrt(1.0f - length), where length = (pt.x * pt.x) + (pt.y * pt.y);

[0098] D3) Set that when the left mouse button is pressed, a starting point is generated, and the coordinate value at the time of pressing is mapped to three-dimensional coordinates and stored in vector form through the previous coordinates; when the mouse is released, an ending point is generated, and the coordinate value at the time of release is also mapped to three-dimensional coordinates and stored in vector form through the coordinates, obtaining the direction vectors of the starting point and the ending point during the rotation interaction process;

[0099] D4) Set a combined quaternion q = [v, w] = [x, y, z, w], which consists of two parts. One is the scalar w, which is equal to cosθ / 2, where θ is the rotation angle, and the other is the vector v, which is equal to sinθ / 2 times the vector along the rotation axis; The result of the operation of two quaternions is the result of their rotation combination, and thus the rotation combination operation is represented by quaternion cross product; The cross product of two rotation vectors records the direction of the rotation axis, and the dot product of two vectors records the rotation angle;

[0100] D5) Call the corresponding function to save the previously calculated cross product and dot product in the form of a quaternion. The two rotation direction vectors are OP1 and OP2. First, calculate the inner product s of the two vectors: s = OP1 · OP2, then calculate the outer product v of the two vectors: v = OP1 × OP2. Denote the quaternion q = [s, v], and normalize it. At this time, q is the rotation quaternion;

[0101] D6) Convert the quaternion to a rotation matrix, and the formula is TM = QuatToMatrix(q). The TM matrix is the corresponding rotation matrix, and the rotation matrix acts on the model transformation matrix rendered by OpenGL, thereby controlling the model to perform rotation operations;

[0102] D7) Convert the quaternion to Euler angles, and the formula is Rotate = QuatToeulerAngles(q). Rotate is the corresponding Euler angle, which can display the rotation angles of the model around the three coordinate axes during the real-time interaction process;

[0103] D8) By performing a rotation interaction operation on the model with the mouse, the quaternion of each rotation interaction can be obtained; the quaternions can be multiplied successively, so that the rotation information can be continuously stored and displayed.

[0104] Quaternion rotation can avoid gimbal lock. Only a 4D quaternion is needed to perform rotation around any vector passing through the origin, which is convenient and fast, and is more efficient than the rotation matrix in some implementations; moreover, quaternion rotation can provide smooth interpolation. It is more convenient to interact through quaternion rotation, and a set of rotation parameters can be provided in real time and accurately. Through the mouse interaction of QT, the coordinate information of the mouse point is obtained. The principle of arcball needs to be used to design the rotation interaction. Generally, Euler rotation is used for model rotation. This design uses quaternion rotation, which is a bit more complex than Euler rotation because there is one more dimension, making it more difficult to understand and less intuitive.

[0105] (5) Set the interaction method for the mouse to translate the 3D model: Calculate the window coordinates of QT and set the translation interaction method of the model through mouse interaction.

[0106] E1) When the right mouse button is pressed, record the current 2D window coordinates, that is, the starting point. When the mouse stops sliding, record the 2D window coordinates at this time, that is, the ending point;

[0107] E2) Call the QPointF function of QT to calculate the differences between the two 2D window coordinates in the X-axis direction and the Y-axis direction, so as to determine the translation of the model in the X-axis direction and the Y-axis direction;

[0108] E3) Call the mouse wheel mechanism of QT to control the translation transformation effect of the model on the Z-axis. When the mouse wheel slides up, the model moves in the positive direction of the Z-axis. When the mouse wheel slides down, the model moves in the negative direction of the Z-axis;

[0109] E4) Substitute the translation parameters of the model in the X-axis direction, Y-axis direction and Z-axis direction into the translation matrix, and substitute the calculated translation matrix into the model transformation matrix, so as to control the translation operation of the model in real time.

[0110] (6) Use the interaction results to control the model to perform real-time rendering and display.

[0111] F1) Convert the quaternion that saves the rotation information into a rotation matrix, and convert the vector that saves the translation information into a translation matrix;

[0112] F2) Set the model transformation matrix as the product of the translation matrix, rotation matrix and scaling matrix; the model transformation matrix acts on the model to control the transformation of the model in the world coordinate system, and transforms the object from the model coordinate system to the world coordinate system;

[0113] F3) Set the view matrix to transform the object from the world coordinate system to the view coordinate system;

[0114] F4) Set the projection matrix to perform a projection transformation on the object model and transform the object from the view coordinate system to the clip coordinate system;

[0115] F5) Set the viewport transformation to transform the object from the clip coordinate system to the window coordinate system, so as to render the drawing result on the window in real time.

[0116] In the third step, load the two-dimensional image to be registered and draw it on the window in the form of a two-dimensional texture. Read the corresponding two-dimensional image, convert it into two-dimensional texture data by using the texture mapping method, map the texture pixels in the texture space to the pixels in the window space, and draw the image on the rendering window.

[0117] In the fourth step, drag the mouse to align the 3D model with the 2D image to obtain appropriate registration parameters. Using mouse interaction, drag the 3D model for rotation and translation operations, so that the 3D model can be aligned with the 2D image under visual observation, and then automatically generate a set of registration parameters for printing output.

[0118] As Figure 2 shown, it is a two-dimensional X-ray hip joint photo. As Figure 3 shown, it is a two-dimensional mask photo, which acts on the two-dimensional X-ray photo. As Figure 4 shown, the image is a projection map rendered and displayed on the window after maximum intensity projection, which is an image drawn on the window after the projection transformation of the 3D hip joint CT slice sequence image. At this time, the model has not undergone rotation and translation interaction operations, and is the display result in the initial posture. Figure 5 shown, the image is the posture display of the model after certain rotation and translation operations. Interact with the model through the left mouse button, drag the model to perform rotation transformation, and make certain rotation transformations around the X-axis, Y-axis, and Z-axis respectively. Interact with the model through the right mouse button, drag the model to perform translation transformation, and make certain translation transformations along the X-axis and Y-axis respectively. Scroll the mouse wheel forward or backward to drag the model to perform scaling transformation, that is, make a certain translation transformation along the Z-axis; this posture is the posture display obtained by rotating the model 80 degrees around the X-axis, translating 26 and 57 unit values along the X-axis and Y-axis respectively, and magnifying 1450 unit values along the Z-axis.

[0119] As Figure 6aAs shown, the two-dimensional X-ray hip joint image is drawn on the window as a background image with two-dimensional texture, and the three-dimensional hip joint CT image is drawn on the window through three-dimensional texture mapping and maximum intensity projection as a foreground image. The mouse can be dragged to perform interactive operations on the foreground image to align it with the background image. This image shows the display results of the foreground image and the background image in the initial pose. As Figure 6b shown, drag the mouse to perform rotation and translation interactive operations on the foreground model to achieve a certain alignment effect with the background X-ray image visually. In the image, it is shown that the right contour of the foreground model is roughly aligned with the background X-ray image, that is, the translation and rotation parameters in this pose are a set of relatively suitable initial configuration parameters. After processing by the method of the present invention, it can be found that a set of suitable 2D and 3D initial registration parameters can be obtained and printed out in a short time.

[0120] The above shows and describes the basic principle, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. What is described in the above embodiments and the specification is only the principle of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection required by the present invention is defined by the appended claims and their equivalents.

Claims

1. An automatic generation method for interactive 2D and 3D medical image registration parameters, characterized in that, Including the following steps: 11) Loading three-dimensional image data: Obtaining the three-dimensional image data of the interactive medical image to be registered; 12) Implementing three-dimensional model rendering and reconstruction on the window, performing two-dimensional mapping on the three-dimensional model, and aligning and displaying the real-time mapped image with the X-ray image to be registered; The implementation of three-dimensional model rendering and reconstruction on the window, performing two-dimensional mapping on the three-dimensional model, and aligning and displaying the real-time mapped image with the X-ray image to be registered includes the following steps: 121) Processing the three-dimensional image data into data suitable for the OpenGL rendering pipeline: For the input sequence of several three-dimensional slice images, i.e., DICOM data, writing a Python script using the pydicom library in Python to convert the DICOM data into a binary file suitable for the system, and then storing the remaining DICOM data in the binary file using a C++ container; 122) Loading the data into the real-time rendering and drawing pipeline: Storing the binary data into the rendering pipeline through the corresponding functions of OpenGL and a fixed rendering process; 123) Drawing of the two-dimensional window: Projecting the corresponding three-dimensional CT image using the maximum intensity projection algorithm and displaying it in the form of volume rendering on the window; 124) Designing the rotation interaction method of the three-dimensional model using the arcball algorithm: Adopting the idea of arcball, storing the information of each mouse interaction model rotation into a quaternion, and controlling the model rotation by converting the quaternion into the form of Euler angles and rotation matrices; The design of the rotation interaction method of the three-dimensional model using the arcball algorithm includes the following steps: 1241) Mapping the two-dimensional window coordinates after mouse interaction, imagining a unit hemisphere located at the center of the window, and adjusting the range of the two-dimensional window coordinates to the interval [-1....1]; Formula: pt.x = (pt.x * AdjustWidth) - 1.0f, pt.y = 1.0f - (pt.y * AdjustHeight); where pt is the defined three-dimensional coordinate, AdjustWidth is the scaling factor of the width, and AdjustHeight is the scaling factor of the length; AdjustWidth = 1.0f / ((NewWidth - 1.0f) * 0.5f), NewWidth and NewHeight are the width and height of the two-dimensional window; AdjustHeight = 1.0f / ((NewHeight - 1.0f) * 0.5f); 1242) Normalize the coordinates of two points in the window to two points on the hemisphere through the mapping formula. If the two-dimensional window coordinates are not on the unit hemisphere centered at the window center, then scale the two-dimensional coordinates to the hemisphere, and the scaling factor is set to norm = 1.0 / FuncSqrt(length), where length = (pt.x * pt.x) + (pt.y * pt.y); Map and convert the two-dimensional coordinates into two points in the hemisphere space. If the two-dimensional coordinates are on the unit hemisphere, then calculate the Z-direction coordinate pt.z based on the X-direction axis coordinate pt.x and the Y-direction coordinate pt.y, and the calculation formula is pt.z = FuncSqrt(1.0f - length), where length = (pt.x * pt.x) + (pt.y * pt.y); 1243) Set that when the left mouse button is pressed, a starting point is generated, and the coordinate value when pressed is mapped into three-dimensional coordinates and stored in vector form through the previous coordinates; when the mouse is released, an ending point is generated, and the coordinate value when released is also mapped into three-dimensional coordinates and stored in vector form to obtain the direction vector of the starting point and the ending point during the rotation interaction process; 1244) Set a combined quaternion q = [v, w] = [x, y, z, w] which consists of two parts. One is the scalar w, which is equal to cosθ / 2, where θ is the rotation angle, and the other is the vector v, which is equal to sinθ / 2 times the vector along the rotation axis; The result of the operation of two quaternions is the result of their rotation combination, and thus the rotation combination operation is represented by quaternion cross multiplication; The cross product of two rotation vectors records the direction of the rotation axis, and the dot product of two vectors records the rotation angle; 1245) Call the corresponding function to save the previously calculated cross product and dot product in the form of a quaternion. The two rotation direction vectors are OP1 and OP2. First, calculate the inner product s of the two vectors: s = OP1 · OP2, then calculate the outer product v of the two vectors: v = OP1 × OP2. Denote the quaternion q = [s, v], and normalize it. At this time, q is the rotation quaternion; 1246) Convert the quaternion into a rotation matrix, and the formula is TM = QuatToMatrix(q). The TM matrix is the corresponding rotation matrix, and the rotation matrix acts on the model transformation matrix rendered by OpenGL to control the model to perform rotation operations; 1247) Convert the quaternion into Euler angles, and the formula is Rotate = QuatToeulerAngles(q). Rotate is the corresponding Euler angle, which can display the rotation angles of the model around the three coordinate axes during the real-time interaction process; 1248) Perform rotation interaction operations on the model through the mouse to obtain the quaternion of each rotation interaction; Continuously store and display the rotation information through quaternion multiplication; 125) Set the interaction method for translating the 3D model with the mouse: Calculate the window coordinates of QT and set the translation interaction method of the model through mouse interaction; 126) Use the interaction result to control the model to perform real-time rendering and display; 13) Load the 2D image to be registered and draw it on the window in the form of a 2D texture: Read the corresponding 2D image, convert it into 2D texture data using texture mapping, map the texture pixels in the texture space to the pixels in the window space, and draw the image on the rendering window; 14) Drag the mouse to align the 3D model with the 2D image to obtain appropriate registration parameters: Use mouse interaction to drag the 3D model for rotation and translation operations, so that the 3D model can be visually aligned with the 2D image, and then automatically generate a set of registration parameters for printing output.

2. The automatic generation method for interactive 2D and 3D medical image registration parameters according to claim 1, characterized in that, The processing of the 3D image data into data suitable for the OpenGL rendering pipeline includes the following steps: 21) Use a python script to read the dicom data in the release folder; Open the cmd console, and convert the dicom data into a binary bin file in a specific format by entering commands in cmd for subsequent use in the rendering process, and store the binary bin file in the corresponding folder; 22) Use a for loop to iterate through all the slice data, use the sprintf_s function in C++ to format the output of the data path into a string, read the data under the path, and use the fopen_s function in C++ to open the binary file; 23) Use the fread_s function in C++ to read the information in the binary bin file, which includes the position of the slice in the image sequence, the spacing of the slice image pixels along the X-axis direction, the spacing of the slice image pixels along the Y-axis direction, the width of the slice, and the height of the slice. Store the remaining dicom data in a container and adjust the size of the container to prevent data overflow; 24) After reading the data, sort all the slice data based on the position of the slice; 25) Load the dicom-formatted data of the 3D volume data through the information in the read binary file.

3. The automatic generation method for interactive 2D and 3D medical image registration parameters according to claim 1, characterized in that, The loading of the data into the real-time rendering and drawing pipeline includes the following steps: 31) The rendering pipeline calls a function to generate a 3D texture object and binds the 3D texture in the pipeline; 32) Perform 3D texture mapping. The parameters of the mapping function are the width and height of each 2D slice image and the depth of the data, and perform texture filtering on the mapped 3D texture; 33) Generate and set the vertex data attributes, and bind the vertex data in the rendering pipeline; 34) Write OpenGL shaders, write vertex shaders. For each vertex Vertex sent to the GPU, perform vertex shading once. Its function is to transform the three-dimensional coordinates of each vertex in the virtual space into two-dimensional coordinates displayed on the window, and carry depth information for the z-buffer; Write fragment shaders to calculate the color and other attributes of each pixel; Compile and link the written shaders, and then delete the shaders after completion; 35) Pass the read 3D texture into the shader, and perform texture sampling to draw the texture of the 3D model, and end the loading.

4. A method for automatically generating registration parameters for interactive 2D and 3D medical images according to claim 1, characterized in that, The drawing of the 2D window includes the following steps: 41) A shader that modifies the ray casting algorithm based on the principle of maximum density projection to implement the maximum density projection function; 42) Use the OpenGL's GetUniformLocation function to obtain the location tags of the maximum density value, texture loading, camera position parameters, and the model view projection transformation matrix MVP matrix in the maximum density projection shader; 43) Load the position parameters of the volume data in space one by one through the corresponding location tags, and then load the 3D texture map; 44) Call the maximum density projection shader to project the loaded 3D volume data, and thus draw the model after maximum density projection on the 2D window.

5. A method for automatically generating registration parameters for interactive 2D and 3D medical images according to claim 1, characterized in that, The interactive method for setting the mouse to pan the 3D model includes the following steps: 51) When the right mouse button is pressed, record the current 2D window coordinates, that is, the starting point. When the mouse stops sliding, record the 2D window coordinates at this time, that is, the ending point; 52) Call the QT's QPointF function to calculate the differences between the two 2D window coordinates in the X-axis direction and the Y-axis direction, so as to determine the translation of the model in the X-axis direction and the Y-axis direction; 53) Call the mouse wheel mechanism of QT to control the translation transformation effect of the model on the Z-axis. When the mouse wheel slides up, the model moves in the positive direction of the Z-axis. When the mouse wheel slides down, the model moves in the negative direction of the Z-axis; 54) Substitute the translation parameters of the model in the X-axis direction, Y-axis direction, and Z-axis direction into the translation matrix, and substitute the calculated translation matrix into the model transformation matrix, so as to control the translation operation of the model in real time.

6. A method for automatically generating registration parameters for interactive 2D and 3D medical images according to claim 1, characterized in that, The method of using the interactive result to control the model for real-time rendering and display includes the following steps: 61) Convert the quaternion that saves the rotation information into a rotation matrix, and convert the vector that saves the translation information into a translation matrix; 62) Set the model transformation matrix as the product of the translation matrix, rotation matrix, and scaling matrix; the model transformation matrix acts on the model to control the transformation of the model in the world coordinate system, and transforms the object from the model coordinate system to the world coordinate system; 63) Set the view matrix to transform the object from the world coordinate system to the view coordinate system; 64) Set the projection matrix to perform a projection transformation on the object model, and transform the object from the view coordinate system to the clip coordinate system; 65) Set the viewport transformation to transform the object from the clip coordinate system to the window coordinate system, so as to render the drawing result on the window in real time.

Images