Image processing method and device, electronic equipment and computer readable storage medium

Abstract

The application provides an image processing method and device, electronic equipment and a computer readable storage medium. The image correction method comprises the following steps: photographing an imaging target and obtaining an initial frame image and a continuous frame image; calculating initial three-dimensional feature points of the initial frame image; calculating continuous frame three-dimensional feature points of the continuous frame image; matching the initial three-dimensional feature points and the continuous frame three-dimensional feature points to obtain an RT matrix; and correcting the image by using the RT matrix to obtain a reconstructed image. According to the example embodiment of the application, the accuracy of the motion-corrected image is improved and the calculation complexity is reduced by using a simple but more accurate image recording device to obtain the motion information of the object.

Description

Technical Field

[0001] This application relates to the field of image processing, and more specifically, to an image processing method, apparatus, electronic device, and computer-readable storage medium for rigid body motion. Background Technology

[0002] Currently, when correcting images acquired from awake targets, correction is mainly performed on two types of motion: non-rigid motion and rigid motion. Non-rigid motion refers to smaller-scale movements, such as the movement of a mouse's torso, while rigid motion refers to larger-scale movements, such as the movement of a mouse's head.

[0003] For the correction of rigid motion, the following two methods are mainly used: motion correction method based on projection domain and motion correction method based on image domain.

[0004] Image-domain-based motion correction methods typically employ Motion-Incorporated Reconstruction (MIR) algorithms that incorporate motion information. The MIR algorithm, proposed by Feng Qiao and Tinsu Pan et al. in 2006, is a voxel-based motion correction algorithm. The MIR algorithm is based on two assumptions: first, that the imaging target only moves between adjacent data frames and remains stationary within each data frame; and second, that the imaging target consists of countless infinitesimally small spatial points, and that the concentration of radioactive material at these points does not change with their movement. In other words, during data acquisition, the concentration of radioactive material at all parts of the imaging target remains constant, and the concentration of radioactive material is only redistributed in space as the imaging target moves.

[0005] The process of image processing using an image domain-based motion correction method includes the following steps.

[0006] First, by using image registration techniques between CT images corresponding to each data frame on a PET (Positron Emission Tomography) image, the mapping relationship matrix M between image voxels is obtained. t This is also known as the motion matrix. Registration of CT images is achieved using a non-rigid transformation model based on cubic B-splines and a mean squared error similarity criterion.

[0007] Motion matrix M t The distribution of radioactive material concentration in the field of view (FOV) region satisfies formula (1).

[0008] X t =Mt X0 formula (1)

[0009] Among them, X t X represents the distribution of radioactive material concentration in the FOV region at time t, and X0 represents the distribution of radioactive material concentration in the FOV region at time 0.

[0010] Secondly, an imaging model of the PET system for stationary objects is constructed, as shown in formula (2).

[0011] E(Y)=GX Formula (2)

[0012] Where Y is a one-dimensional vector representing the projection data, E(Y) is the expected value of the projection data Y; X is a one-dimensional vector composed of image voxels; and G is the system response matrix. Furthermore, two-dimensional and three-dimensional images can also be converted into one-dimensional images for the above processing.

[0013] The system response matrix G can be calculated through actual measurement, Monte Carlo (MC) simulation, or mathematical analysis.

[0014] The calculation of the system response matrix G using actual measurement methods includes: measuring the same point source at different positions in the imaging field of the PET system, and then constructing the entire system response matrix G by processing the responses in the projection space. These responses include the geometric information of the system and the physical information of the detection process, and the responses are parameterized to correct the point source position and smooth the projection noise.

[0015] The calculation of the system response matrix G using the Monte Carlo simulation method includes: in space, simulating each voxel within the field of view based on its size and position, obtaining the ratio of its coincidence count to the total number of emitted photons, and obtaining the probability that this voxel is detected by the current response line (LOR), which is the system response matrix G.

[0016] Computational methods based on geometric models, also known as analytical methods, are categorized into many types based on the different geometric models used, such as line tracing and solid angle methods. Line tracing is a simple and effective method for calculating the system response matrix G. This method abstracts the response line as a real ray, using the length of this ray passing through each voxel as the probability that a photon emitted by that voxel will be detected by this response line, thus obtaining the system response matrix G. Solid angle models are more accurate and complex than line tracing. This method uses the size of the spatial vertices formed by each voxel and a pair of detectors as the probability that a photon emitted by that voxel will be detected by this pair of detectors, thus obtaining the system response matrix G. Compared to line tracing, solid angle models better reflect the actual PET imaging process.

[0017] Furthermore, based on the assumption that the imaging target remains stationary within each data frame, Equation (2) is modified into Equation (3).

[0018] E{Y t} = GM t X0 formula (3)

[0019] Since all data frames correspond one-to-one with the image sequence X0,…,X m-1 Projected data Y0,…,Y m-1 and motion matrices M0,…,M m-1 , where m is the number of data frames. Therefore, formula (2) can be modified into the discrete form shown in formula (4).

[0020] E(Y0) = GM0X0

[0021]

[0022] E(Y m-1 ) = GM m-1 X0 formula (4)

[0023] Where M0≡I, is the identity matrix. Let:

[0024]

[0025] G * =[(GM0)] T (GM1) T ,…,(GM m-1 ) T ] T Formula (6)

[0026] Formula (4) can be rewritten as formula (7), which is consistent with the form of formula (2) of the PET imaging model for stationary targets.

[0027] E(Y * ) = G * X0 formula (7)

[0028] Since the obtained imaging model of the PET system for the imaging target is consistent with the imaging model of the PET system for the static target, the 4D projection data Y in formulas (5) and (6) can be used. * and the system response matrix G that incorporates motion information from each data frame * Replace the system response matrix and projection data matrix in Equation (2). Therefore, the voxel-based iterative reconstruction algorithm MLEM, which is applied to Equation (2) for PET imaging models of stationary targets, is also applicable to Equation (7).

[0029] The iterative update equation of the MLEM algorithm based on equation (2) is shown in equation (8).

[0030]

[0031] in, Let j be the value of voxel j in the image obtained in the k-th iteration. Y represents the value of voxel j in the image obtained in the (k+1)th iteration. i G is the count of coincidence events on the i-th response line (LOR) that was actually measured. ij denoted by , M represents the probability that a gamma photon emitted by voxel j is detected by the scintillation crystals at both ends of the i-th LOR, N represents the total number of response lines (LORs), and j0 represents the total number of voxels.

[0032] Y * G * Substituting into the MLEM iterative update equation, we can obtain the update equation based on formula (7), as shown in formula (9).

[0033]

[0034] Among them, X k S is the estimate of the reference position image obtained in the k-th iteration, and S is the sensitivity image, which can be calculated by formula (10).

[0035]

[0036] For ease of understanding, the following operation rules are adopted in the following matrix operations: If matrix A and matrix B have the same number of rows and columns, then AB and A / B represent the multiplication and division of corresponding elements, respectively.

[0037] As can be seen from formula (9), the data from all data frames are involved in the reconstruction of the reference image.

[0038] The drawback of existing image correction techniques is the difficulty in obtaining the mapping matrix M. t In such cases, it's not feasible to use only one type of medical imaging equipment; an additional 4D CT scanner is required for CT scanning, and the object's motion must be consistent across both devices. This not only increases imaging costs but also limits the algorithm to handling regular, small-range movements, such as breathing and heartbeats. Furthermore, it's difficult to guarantee that the target's motion is completely consistent across both devices, and image registration techniques typically involve optimizing numerous parameters. Therefore, methods for obtaining voxel mapping information through image registration are only suitable for small-range motion correction; the accuracy of the obtained motion information cannot be guaranteed when the motion range is large. Additionally, registration using CT images requires significant computational resources, resulting in low accuracy. Summary of the Invention

[0039] This application proposes an image processing method, apparatus, electronic device, and computer-readable storage medium for rigid body motion, in order to solve the problems of high acquisition cost, computational complexity, and low accuracy when acquiring motion information of moving imaging targets.

[0040] According to one aspect of this application, an image processing method is proposed, the image processing method comprising: capturing an imaging target and acquiring an initial frame image and consecutive frame images; calculating initial three-dimensional feature points of the initial frame image; calculating consecutive frame three-dimensional feature points of the consecutive frame images; matching the initial three-dimensional feature points and the consecutive frame three-dimensional feature points to obtain a RT matrix; and using the RT matrix to correct the image to obtain a reconstructed image.

[0041] According to some embodiments of this application, an image recording device is used to capture an imaging target and acquire initial frame images and continuous frame images. The image recording device includes any device that can record planar images or 3D images of the imaging target, such as a camera, DV, camcorder, 3D scanner, etc.

[0042] According to some embodiments of this application, a pair of image recording devices of the same model are placed parallel to each other on both sides of the imaging target to form a binocular vision system for capturing the imaging target.

[0043] According to some embodiments of this application, the motion of the imaging target is recorded in the form of video using the image recording device. The video includes an initial frame image at the initial moment and consecutive frame images at non-initial moments, wherein the initial frame image and the consecutive frame images each include a pair of images captured by the image recording device at the same moment.

[0044] According to some embodiments of this application, the step of calculating the initial three-dimensional feature points of the initial frame image includes: extracting two-dimensional feature points of the initial frame image using a scale-invariant feature transform method; removing the background image of the two-dimensional feature points of the initial frame image; matching the two-dimensional feature points of the initial frame image after removing the background; and triangulating the matching information to obtain the initial three-dimensional feature points of the initial frame image.

[0045] According to some embodiments of this application, calculating the continuous frame three-dimensional feature points of the continuous frame image includes: extracting two-dimensional feature points of the continuous frame image using a scale-invariant feature transform method; removing the background image of the two-dimensional feature points of the continuous frame image; matching the two-dimensional feature points of the continuous frame image after removing the background; and triangulating the matching information to obtain the initial three-dimensional feature points of the continuous frame image.

[0046] According to some embodiments of this application, the step of triangulation matching information includes: calculating the depth Z of any point P on the imaging target in the coordinates of the image recording device using formula (7-1); and calculating the x-coordinate X and y-coordinate Y of any point P on the imaging target in the coordinates of the camera using formula (7-2).

[0047]

[0048]

[0049] Where T is the optical distance between the two image recording devices, f is the focal length of the image recording device, x1 and xr are the abscissa distances of the intersections of the projections of point P on the image frames of the two image recording devices and the optical axis image, respectively; y1 is the ordinate distance of the intersections of the projections of point P on the image frame of one image recording device and the optical axis image.

[0050] According to some embodiments of this application, the step of obtaining the RT matrix includes: minimizing the error between the initial three-dimensional feature points and the three-dimensional feature points of the consecutive frames by using a singular value decomposition method to obtain the RT matrix.

[0051] According to some embodiments of this application, the error is minimized by formula (9-1):

[0052]

[0053] Where R and T represent the rotation matrix and translation matrix, respectively, J represents the error, and p i q i The initial set of three-dimensional feature points is P = {p1, ..., p2}. n} and the set of three-dimensional feature points of the consecutive frames Q = {q1, ..., q n The i-th element in}, where i and n are both natural numbers, i ≤ n.

[0054] According to some embodiments of this application, the Singular Value Decomposition (SVD) method is performed using formula (10-1):

[0055] [U,S,V]=SVD(H) Formula (10-1)

[0056] Where U, S, and V represent the first unitary matrix, diagonal matrix, and second unitary matrix after singular value decomposition, respectively; SVD represents the singular value decomposition of the physical quantity within the brackets; and H is the covariance matrix between the re-centered point sets, which is determined by formula (10-2):

[0057]

[0058] Where, p i 'and q i ' represents the three-dimensional feature points after re-centering, p i '=p i -μ P q i '=q i -μ Q .

[0059] According to some embodiments of this application, the rotation matrix R and the translation matrix T are determined by formulas (11-1) and (11-2), respectively:

[0060] R = VU T Formula (11-1)

[0061] T = -Rμ P +μ Q Formula (11-2)

[0062] Where, μ P and μ Q These are the center positions of sets P and Q, respectively.

[0063] According to some embodiments of this application, the images include CT images, MR images, PET images, PET-CT images, PET-MR images, and CT-MR images.

[0064] According to some embodiments of this application, the step of correcting the PET image using the RT matrix to obtain a reconstructed image includes: performing a PET scan on the imaging target to obtain sine wave data; and correcting the sine wave data using the RT matrix to reconstruct the image.

[0065] According to some embodiments of this application, before correcting the image using the RT matrix, the image processing method further includes correcting the voxel values ​​of the image using formula (14-1):

[0066] X t =W t RTX0 formula (14-1)

[0067] Among them, X t Let Xt represent the concentration distribution of radioactive material in the FOV region at time t, X0 represent the initial concentration distribution of radioactive material in the FOV region, and Wt represent the concentration distribution of radioactive material in the FOV region at time t. t The weights for the interpolation of all voxel positions at time t.

[0068] According to some embodiments of this application, the steps for obtaining the RT matrix include: setting a plurality of radioactive point sources on the imaging target; identifying the radioactive point sources by reconstructing the images within the acquisition time frame; obtaining the pose of the imaging target based on the surface formed by three non-collinear radioactive point sources, wherein the pose includes the position of the imaging target and the rotation angle of the imaging target around the three-dimensional coordinate axis; and obtaining the RT matrix by calculating the motion between consecutive frames.

[0069] According to one aspect of this application, an image processing apparatus is proposed, the image processing apparatus comprising: a frame image acquisition unit for capturing an imaging target and acquiring an initial frame image and consecutive frame images; an initial three-dimensional feature point calculation unit for calculating initial three-dimensional feature points of the initial frame image; a consecutive frame three-dimensional feature point calculation unit for calculating consecutive frame three-dimensional feature points of the consecutive frame images; an RT matrix calculation unit for matching the initial three-dimensional feature points and the consecutive frame three-dimensional feature points to obtain an RT matrix; and an image reconstruction unit for correcting the image using the RT matrix to obtain a reconstructed image.

[0070] According to one aspect of this application, an electronic device is provided, comprising one or more processors; a storage device for storing a computer program; and, when the computer program is executed by the one or more processors, causing the one or more processors to perform the method as described above.

[0071] According to one aspect of this application, a computer-readable storage medium is provided having program instructions stored thereon, which, when executed, implement the method as described in any of the preceding claims.

[0072] According to some example embodiments of this application, a simple yet highly accurate image recording device is used to acquire motion information of an object while it is moving, ensuring the accuracy of motion information acquisition and improving the accuracy of motion-corrected images. At the same time, it avoids performing a large number of complex image registration calculations, reducing computational complexity. Moreover, the provided algorithm is not only applicable to small-scale motion correction but also to larger-scale motion correction, significantly reducing costs. Attached Figure Description

[0073] To more clearly illustrate the technical solutions in the embodiments of this application, the accompanying drawings used in the description of the embodiments will be briefly introduced below.

[0074] Figure 1 A flowchart of an image processing method according to an example embodiment of this application is shown.

[0075] Figure 2 A triangulation schematic diagram according to an example embodiment of this application is shown.

[0076] Figure 3 A schematic diagram of a voxel mapping process in a two-dimensional graphic according to an example embodiment of this application is shown.

[0077] Figure 4 An image correction apparatus for rigid body motion is shown according to an example embodiment of this application.

[0078] Figure 5a The correction results are shown when the object is stationary.

[0079] Figure 5b The reconstruction results before motion correction are shown.

[0080] Figure 5c The reconstruction result after motion correction of an object is shown according to an example embodiment of this application.

[0081] Figure 6 A block diagram of another electronic device for rigid body motion according to an embodiment of this application is shown. Detailed Implementation

[0082] Exemplary embodiments will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limited to the embodiments set forth herein; rather, they are provided so that this application will be thorough and complete, and will fully convey the concept of the exemplary embodiments to those skilled in the art. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted.

[0083] The described features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. Numerous specific details are provided in the following description to give a full understanding of embodiments of this disclosure. However, those skilled in the art will recognize that the technical solutions of this disclosure can be practiced without one or more of these specific details, or other methods, components, materials, apparatus, or operations may be employed. In these cases, well-known structures, methods, apparatuses, implementations, materials, or operations will not be shown or described in detail.

[0084] The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily have to be performed in the described order. For example, some operations / steps can be broken down, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0085] The terms "first," "second," etc., in the specification, claims, and accompanying drawings of this application are used to distinguish different objects, not to describe a specific order. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or apparatus that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or apparatuses.

[0086] The specific embodiments according to this application will now be described in detail with reference to the accompanying drawings.

[0087] Figure 1 A flowchart illustrating an image processing method according to an example embodiment of this application is shown. Referring below... Figure 1 A detailed description of an image processing method according to an example embodiment of this application will be provided.

[0088] According to some example embodiments of this application, Figure 1 The example embodiment shown corrects motion artifacts generated during PET scanning of a lucid imaging target. Since the instrument's final image is only focused on the head of the imaging target, the motion of the moving head of the imaging target can be simplified as rigid body motion.

[0089] In step S101, the imaging target is captured and the initial frame image and consecutive frame images are obtained.

[0090] According to some example embodiments of this application, the motion of the imaging target is recorded in video form by placing a pair of image recording devices parallel to each other on both sides of the imaging target. These image recording devices include any device capable of recording planar or 3D images of the imaging target, such as cameras, digital video cameras, camcorders, 3D scanners, etc. The video captured by the pair of image recording devices includes an initial frame image at the initial moment and consecutive frame images at non-initial moments. The initial frame image and consecutive frame images each include images captured by the pair of image recording devices at the same moment. Those skilled in the art should understand that in this embodiment, the image recording devices are placed parallel to each other to facilitate a unified coordinate system during subsequent calculations and processing, reducing computational complexity. In fact, the image recording devices can be set at any angle and in any number; appropriate coordinate system conversions can be performed during subsequent data processing.

[0091] According to some embodiments, the initial frame image is the frame image at the beginning of a video captured by an image recording device, and the images at other times in the video are consecutive frame images.

[0092] According to some embodiments, in step S101, a pair of identical image recording devices, such as cameras, are placed parallel to each other on the left and right sides of the imaging target to form a binocular vision system, which records the motion information of the object in the form of video.

[0093] In step S103, the initial three-dimensional feature points of the initial frame image are calculated.

[0094] According to some example embodiments of this application, the initial three-dimensional feature points of the initial frame image are calculated through the following steps:

[0095] First, the scale-invariant feature transform (SIFT) method is used to extract two-dimensional feature points from the initial frame image.

[0096] Then, the background image of the two-dimensional feature points of the initial frame image is removed using the principle of digital subtraction angiography (DSA) to ensure that the extracted two-dimensional feature points are all on the imaging target image in the initial frame image.

[0097] Finally, the two-dimensional feature points in the initial frame images captured by each of the pair of image recording devices are matched and the matching information is triangulated to obtain the initial three-dimensional feature points of the initial frame images.

[0098] In step S105, the three-dimensional feature points of the consecutive frames of the consecutive frame images are calculated.

[0099] In step S105, the method for calculating the three-dimensional feature points of consecutive frames is the same as the method for calculating the three-dimensional feature points of the initial frame in step S103, and will not be described again here.

[0100] The triangulation principle in steps S103 and S105 will be explained below. Triangulation refers to the process of recovering the three-dimensional position information of an object's feature points based on their two-dimensional projection positions on different image recording devices. For example... Figure 2 The diagram shown is a triangulation principle diagram according to an example embodiment of this application.

[0101] like Figure 2 As shown, let the coordinates of the intersection point of the camera's optical axis and the image be (cx, cy), the distance between the optical centers of the two cameras be D, and the focal length be f. The projections of a feature point P in three-dimensional space onto the left and right images are Pl and Pr, respectively. The horizontal coordinate distances between the feature point P and the intersection point of the optical axis and the image are xl and xr, respectively, and the vertical coordinate distances are yl and yr.

[0102] Using similar triangles, the depth Z of point P in the coordinates of the image recording device can be obtained by formula (11).

[0103]

[0104] The x-coordinate X and y-coordinate Y of point P are calculated using formula (12).

[0105]

[0106] Therefore, the set of three-dimensional feature points for each frame of the image can be calculated using formulas (11) and (12).

[0107] In step S107, the initial three-dimensional feature points and the three-dimensional feature points of the continuous frames are matched to obtain the RT matrix.

[0108] According to some example embodiments of this application, the steps for calculating the RT matrix include:

[0109] According to steps S103 and S105, the initial three-dimensional feature point set P = {p1, ..., p2} of the initial frame image is obtained respectively. n The set of three-dimensional feature points in consecutive frames Q = {q1, ..., q} n}

[0110] The error J between the initial 3D feature points and the 3D feature points of consecutive frames is minimized by the Singular Value Decomposition (SVD) method, as shown in Equation (13).

[0111]

[0112] Where J is the error between the initial 3D feature points and the 3D feature points of consecutive frames, R is the rotation matrix, and Rp i Represents the relationship between rotation matrix R and the i-th point p in point set P. i Multiplication, where T is the translation matrix, p i q i Let i and n represent the i-th three-dimensional feature points in sets P and Q, respectively, where i is a natural number less than n, and n is the number of three-dimensional feature points in sets P and Q.

[0113] Using formula (14) to average the coordinates of each point in set P and set Q, the center position μ of set P and set Q can be obtained. P and μ Q .

[0114]

[0115] To eliminate the influence of the translation matrix T on the solution of the rotation matrix R, the point set is re-centered using formula (15) to obtain new point sets P′ and Q′, i.e., the coordinates of each point minus the center position. At this time, the center position coincides with the origin of the coordinate system. These new point sets P′ and Q′ respectively include the re-centered three-dimensional feature points p i 'and q i ', p i '=p i -μ P q i '=q i -μ Q .

[0116]

[0117] The covariance matrix H between the re-centered point sets is solved using formula (16).

[0118]

[0119] Singular value decomposition is performed on the covariance matrix H using formula (17).

[0120] [U,S,V]=SVD(H) Formula (17)

[0121] Where U, S, and V represent the first unitary matrix, the diagonal matrix, and the second unitary matrix after singular value decomposition, respectively. SVD represents the singular value decomposition of the physical quantity in parentheses. For example, if the covariance matrix H is an m×n matrix, then the first unitary matrix U after singular value decomposition is an m×m unitary matrix, the diagonal matrix S is an m×n diagonal matrix with the elements on the diagonal being the singular values ​​of the covariance matrix H, and the second unitary matrix V is an n×n unitary matrix.

[0122] The rotation matrix R is calculated based on the decomposition result using formula (18).

[0123] R = VU T Formula (18)

[0124] The translation matrix T is calculated using the rotation matrix R obtained from formula (18) via formula (19).

[0125] T = -Rμ P +μ Q Formula (19)

[0126] The R and T matrices obtained through formulas (18) and (19) represent the six-axis motion information of the rigid body motion of an object. The rigid body motion of an object can be decomposed into translational and rotational motions around three axes in space, namely the x, y, and z axes. The resulting motions of the object around these three axes can form an RT matrix. The RT matrix is ​​also called the rotation-translation (RT) matrix, used to describe the six-axis motion information of an object. The six axes are translational and rotational motions in the x, y, and z axes.

[0127] According to some embodiments, the R and T matrices obtained by formulas (18) and (19) are combined to obtain the RT matrix.

[0128] According to some embodiments of this application, the steps for obtaining the RT matrix include: setting a plurality of radioactive point sources on the imaging target; identifying the radioactive point sources by reconstructing the images within the acquisition time frame; obtaining the pose of the imaging target based on the surface composed of three non-collinear radioactive point sources, wherein the pose includes the position of the imaging target and the rotation angle of the imaging target around the three-dimensional coordinate axis; and obtaining the RT matrix by calculating the motion between consecutive frames.

[0129] In step S109, the image is corrected using the RT matrix to obtain a reconstructed image.

[0130] According to some example embodiments of this application, the step of correcting PET images using an RT matrix includes:

[0131] A PET scan of the live target is performed to obtain sinusoidal data. The sinusoidal data is divided into time frames, and the data in the current time frame is corrected to the initial position using formula (9) using the RT matrix in step S107. That is, the mapping relationship matrix M in formula (9) is replaced by the RT matrix. t This allows us to obtain a corrected and reconstructed image.

[0132] According to some embodiments of this application, before performing step S109, the voxel values ​​of the image need to be corrected using formula (20).

[0133] X t =W t RTX0 formula(20)

[0134] Among them, X t Let Xt represent the concentration distribution of radioactive material in the FOV region at time t, and X0 represent the concentration distribution of radioactive material in the FOV region at time 0. Wt t The weights for the interpolation of all voxel positions at time t.

[0135] The following details how to correct the voxel values ​​of PET images using formula (20).

[0136] To replace the complex image registration process and obtain the mapping relationship between image voxels at different times, the center point of each voxel can be transformed into spatial coordinates based on motion information to obtain its corresponding voxel position. Since the point obtained after spatial coordinate transformation is not necessarily the voxel center of the image, in order to improve the positioning accuracy, the voxel position at time t is transformed back to the voxel position corresponding to the initial time based on motion information. Then, the voxel value of the transformed voxel position is calculated by interpolation algorithm using the voxel values ​​around the transformed voxel position, as shown in formula (21). i Activity value f(r) at the location i ,0) is considered to be related to the spatial point r i The weighted sum of the voxel values ​​of all neighboring voxels.

[0137]

[0138] in, Represents a spatial point r i The set of subscripts of nearby voxels, ω ij For weights. x j (0) corresponds to X0 in formula (20), f(r) i ,0) and X in formula (21) t Correspondingly, the matrix representation of formula (21) is shown in formula (20).

[0139] According to some embodiments, based on spatial point r i The ratio of the area or volume of overlap between the voxel centered at its position and its neighboring voxels to the total area or volume of the voxels is used as the weight value ω. ij .

[0140] Similarly, the activity values ​​f(r0,0), f(r1,0), ..., f(r) of continuous spatial points can be obtained through formula (21). n-1 ,0).

[0141] Figure 3 A schematic diagram illustrating a voxel mapping process in a two-dimensional graphic according to an example embodiment of this application is shown. Figure 3 As shown in the figure, the gray dots represent a voxel center A at time t that is mapped to its corresponding original position B at time 0. The voxel value of the original position B can be calculated by interpolating the voxels near the position of B at time 0.

[0142] According to some embodiments, the RT matrix can be obtained using point labeling technology. Specific steps include: attaching small, lightweight radioactive point sources to the head of an experimental subject (e.g., a mouse); reconstructing images within the acquisition time frame; identifying highlighted points in the images as point sources; assuming the point source positions are stationary relative to the rat's head; and obtaining the rat's head pose based on the surface formed by three non-collinear point sources. The head pose includes the object's position and three rotation angles around the x, y, and z axes, used to represent the object's orientation. The motion matrix RT is obtained by calculating the motion between consecutive frames.

[0143] According to some example embodiments of this application, a simple yet more accurate external camera is used to acquire motion information of a moving object while it is in motion, ensuring the accuracy of motion information acquisition and improving the accuracy of motion-corrected images. At the same time, it avoids performing a large number of complex image registration calculations, reducing computational complexity. Moreover, the provided algorithm is not only applicable to small-scale motion correction but also to larger-scale motion correction, significantly reducing costs.

[0144] Figure 4 An image processing apparatus according to an example embodiment of this application is shown. For example... Figure 4 The image processing device shown includes a frame image acquisition unit 401, an initial three-dimensional feature point calculation unit 403, a continuous frame three-dimensional feature point calculation unit 405, an RT matrix calculation unit 407, and an image reconstruction unit 409.

[0145] The frame image acquisition unit 401 uses a camera to capture images of the imaging target and acquires initial frame images and consecutive frame images. Camera settings and image acquisition can be achieved using... Figure 1 The methods described in the embodiments will not be repeated here.

[0146] The initial 3D feature point calculation unit 403 is used to calculate the initial 3D feature points of the initial frame image. The calculation of the initial 3D feature points can utilize... Figure 1 The method described in step S103 of the embodiment will not be repeated here.

[0147] The continuous frame 3D feature point calculation unit 405 is used to calculate the continuous frame 3D feature points of a continuous frame image. The calculation of continuous frame 3D feature points can utilize... Figure 1 The method described in step S105 of the embodiment will not be repeated here.

[0148] The RT matrix calculation unit 407 is used to match initial 3D feature points with 3D feature points from consecutive frames to obtain the RT matrix. The RT matrix can be calculated using... Figure 1 The method described in step S107 of the embodiment will not be repeated here.

[0149] Image reconstruction unit 409 is used to correct PET images using the RT matrix to obtain reconstructed images. Image reconstruction can utilize... Figure 1 The method described in step S109 of the embodiment will not be repeated here.

[0150] Figure 5a The correction results are shown when the object is stationary. Figure 5b The reconstruction results before motion correction are shown. Figure 5c The image shows the reconstruction result of an object after motion correction according to an example embodiment of this application. Figures 5a-5c It can be seen that by using the image processing method provided in this application, artifacts caused by object motion can be basically eliminated, thereby improving the quality of the reconstructed image.

[0151] Those skilled in the art should note that the above embodiments correct PET images using the RT matrix. In fact, the above method can also be applied to motion correction of other medical images such as CT images and MR images, or to motion correction of one or more types of images including PET-CT images, PET-MR images, and CT-MR images. In the correction process of these medical images, only the corresponding M... t This can be achieved using a matrix, which is easily accomplished by those skilled in the art based on the inspiration of this invention, and will not be elaborated further here.

[0152] Figure 6 A block diagram of another electronic device for rigid body motion according to an embodiment of this application is shown. Figure 6 The electronic device shown is merely an example and should not impose any limitation on the functionality and scope of use of the embodiments of this application.

[0153] like Figure 6 As shown, this electronic device is presented in the form of a general-purpose computing device. The components of this electronic device may include, but are not limited to: at least one processor 210, at least one memory 220, a bus 230 connecting different system components (including memory 220 and processor 210), a display unit 240, etc. The memory 220 stores program code, which can be executed by the processor 210, causing the processor 210 to perform the methods described in this specification according to the various exemplary embodiments of this application. For example, the processor 210 can perform, for example... Figure 1 The method shown.

[0154] The memory 220 may include a readable medium in the form of volatile memory cells, such as random access memory (RAM) 2201 and / or cache memory 2202, and may further include read-only memory (ROM) 2203.

[0155] The memory 220 may also include a program / utility 2204 having a set (at least one) of program modules 2205, including but not limited to: an operating system, one or more application programs, other program modules, and program data, each or some combination of these examples may include an implementation of a network environment.

[0156] Bus 230 can represent one or more of several types of bus structures, including a memory cell bus or memory cell controller, a peripheral bus, a graphics acceleration port, a processing unit, or a local bus using any of the various bus structures.

[0157] The electronic device can also communicate with one or more external devices 300 (e.g., keyboard, pointing device, Bluetooth device, etc.), one or more devices that enable a user to interact with the image correction device, and / or any device that enables the electronic device to communicate with one or more other computing devices (e.g., router, modem, etc.). This communication can be performed via input / output (I / O) interface 250. Furthermore, the electronic device can communicate with one or more networks (e.g., local area network (LAN), wide area network (WAN), and / or public networks, such as the Internet) via network adapter 260. Network adapter 260 can communicate with other modules of the electronic device via bus 230. It should be understood that, although not shown in the figures, other hardware and / or software modules can be used in conjunction with the electronic device, including but not limited to: microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data backup storage systems.

[0158] From the above description of the embodiments, those skilled in the art will readily understand that the exemplary embodiments described herein can be implemented by software or by combining software with necessary hardware. The technical solutions according to the embodiments of this application can be embodied in the form of a software product, which can be stored in a computer-readable storage medium (such as a CD-ROM, USB flash drive, or external hard drive) or on a network, including several computer program instructions to cause a computing device (such as a personal computer, server, or network device, etc.) to execute the methods described above according to the embodiments of this application.

[0159] Software products may employ any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example,, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: electrical connections with one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0160] Computer-readable storage media may include data signals propagated in baseband or as part of a carrier wave, carrying readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. A readable storage medium may also be any readable medium other than a readable storage medium that can transmit, propagate, or transfer a program for use by or in connection with an instruction execution system, apparatus, or device. The program code contained on the readable storage medium may be transmitted using any suitable medium, including but not limited to wireless, wired, optical fiber, RF, etc., or any suitable combination thereof.

[0161] Program code for performing the operations of this application can be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java and C++, and conventional procedural programming languages ​​such as C or similar languages. The program code can execute entirely on the user's computing device, partially on the user's computing device, as a standalone software package, partially on the user's computing device and partially on a remote computing device, or entirely on a remote computing device or server. In cases involving remote computing devices, the remote computing device can be connected to the user's computing device via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computing device (e.g., via the Internet using an Internet service provider).

[0162] The aforementioned computer-readable medium carries one or more program instructions that, when executed by a device, cause the computer-readable medium to perform the aforementioned functions.

[0163] Those skilled in the art will understand that the above modules can be distributed in the device as described in the embodiments, or they can be modified accordingly to be uniquely different from one or more devices in this embodiment. The multiple modules of the above embodiments can be combined into one module, or a single module can be further divided into multiple sub-modules.

[0164] From the above description of the embodiments, those skilled in the art will readily understand that the exemplary embodiments described herein can be implemented by software or by combining software with necessary hardware. The technical solutions according to the embodiments of this application can be embodied in the form of a software product, which can be stored in a computer-readable storage medium (such as a CD-ROM, USB flash drive, or external hard drive) or on a network, including several computer program instructions to cause a computing device (such as a personal computer, server, or network device, etc.) to execute the methods described above according to the embodiments of this application.

[0165] Software products may employ any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example,, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: electrical connections with one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0166] Computer-readable storage media may include data signals propagated in baseband or as part of a carrier wave, carrying readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. A readable storage medium may also be any readable medium other than a readable storage medium that can transmit, propagate, or transfer a program for use by or in connection with an instruction execution system, apparatus, or device. The program code contained on the readable storage medium may be transmitted using any suitable medium, including but not limited to wireless, wired, optical fiber, RF, etc., or any suitable combination thereof.

[0167] Program code for performing the operations of this application can be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java and C++, and conventional procedural programming languages ​​such as C or similar languages. The program code can execute entirely on the user's computing device, partially on the user's computing device, as a standalone software package, partially on the user's computing device and partially on a remote computing device, or entirely on a remote computing device or server. In cases involving remote computing devices, the remote computing device can be connected to the user's computing device via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computing device (e.g., via the Internet using an Internet service provider).

[0168] The aforementioned computer-readable medium carries one or more program instructions that, when executed by a device, cause the computer-readable medium to perform the aforementioned functions.

[0169] Those skilled in the art will understand that the above modules can be distributed in the device as described in the embodiments, or they can be modified accordingly to be uniquely different from one or more devices in this embodiment. The multiple modules of the above embodiments can be combined into one module, or a single module can be further divided into multiple sub-modules.

[0170] According to some example embodiments of this application, a simple yet more accurate external camera is used to acquire motion information of a moving object while it is in motion, ensuring the accuracy of motion information acquisition and improving the accuracy of motion-corrected images. At the same time, it avoids performing a large number of complex image registration calculations, reducing computational complexity. Moreover, the provided algorithm is not only applicable to small-scale motion correction but also to larger-scale motion correction, significantly reducing costs.

[0171] While this application provides the operational steps of the methods described in the above embodiments or flowcharts, the methods may include more or fewer operational steps based on conventional or non-inventive methods. For steps where there is no logically necessary causal relationship, the execution order of these steps is not limited to the execution order provided in the embodiments of this application.

[0172] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

[0173] The embodiments of this application have been described in detail above. Specific examples have been used to illustrate the principles and implementation methods of this application. The descriptions of the embodiments above are only for the purpose of helping to understand the method and core ideas of this application. Furthermore, any changes or modifications made by those skilled in the art based on the ideas of this application, and on the specific implementation methods and application scope of this application, are all within the scope of protection of this application. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. An image processing method, characterized by, The image processing method comprises: photographing an imaging target and obtaining an initial frame image and a continuous frame image; calculating initial three-dimensional feature points of the initial frame image; calculating continuous frame three-dimensional feature points of the continuous frame image; matching the initial three-dimensional feature points and the continuous frame three-dimensional feature points to obtain an RT matrix; According to the motion information, the center point of each voxel of the image is spatially transformed to obtain a corresponding voxel position, and the voxel value of the transformed voxel position is calculated by using the voxel values around the transformed voxel position through an interpolation algorithm to correct the voxel value of the image, and the formula of the interpolation algorithm used is X t = W t RTX0, wherein X t represents the distribution of the concentration of the radioactive substance in the FOV region at time t, X0 represents the distribution of the concentration of the radioactive substance in the FOV region at the initial time, W t is the weight of the voxel value around the transformed voxel position. correcting an image by using the RT matrix to obtain a reconstructed image.

2. The image processing method of claim 1, wherein, photographing an imaging target by using an image recording device and obtaining an initial frame image and a continuous frame image.

3. The image processing method of claim 2, wherein, A pair of image recording devices of the same model are placed in parallel on two sides of the imaging target to form a binocular vision system to photograph the imaging target.

4. The image processing method of claim 2, wherein, The image recording device records the movement of the imaging target in the form of a video, and the video comprises an initial frame image at an initial time and a continuous frame image at a non-initial time, wherein the initial frame image and the continuous frame image each comprise images photographed by a pair of image recording devices at the same time.

5. The image processing method of claim 1, wherein, The step of calculating initial three-dimensional feature points of the initial frame image comprises: extracting two-dimensional feature points of the initial frame image by using a scale-invariant feature transformation method; removing background images of the two-dimensional feature points of the initial frame image; matching the two-dimensional feature points of the initial frame image after removing the background; triangulating the matching information to obtain the initial three-dimensional feature points of the initial frame image.

6. The image processing method of claim 1, wherein, The step of calculating continuous frame three-dimensional feature points of the continuous frame image comprises: extracting two-dimensional feature points of the continuous frame image by using a scale-invariant feature transformation method; removing background images of the two-dimensional feature points of the continuous frame image; matching the two-dimensional feature points of the continuous frame image after removing the background; triangulating the matching information to obtain the initial three-dimensional feature points of the continuous frame image.

7. The image processing method of claim 5 or 6, characterized in that, The step of triangulating the matching information comprises: calculating the depth Z of any point on the imaging target in the image recording device coordinates by formula (7-1); calculating the horizontal coordinate X and the vertical coordinate Y of any point on the imaging target in the image recording device coordinates by formula (7-2); wherein T is the optical distance of the two image recording devices, f is the focal length of the image recording device, x1 and xr are the horizontal coordinate distances of the intersection of the projection of the point on the image frame of the two image recording devices and the optical axis image, respectively, and y1 is the vertical coordinate distance of the intersection of the projection of the point on the image frame of one of the image recording devices and the optical axis image.

8. The image processing method of claim 1, wherein, The step of obtaining the RT matrix comprises: minimizing the error of the initial three-dimensional feature points and the continuous frame three-dimensional feature points by a singular value decomposition method to obtain the RT matrix.

9. The image processing method of claim 8, wherein, The error is minimized by formula (9-1): wherein R, T represent a rotation matrix and a translation matrix respectively, J represents the error, p i , q i are the i-th elements in the set of initial three-dimensional feature points P = {p1, …, p n n} and the set of continuous frame three-dimensional feature points Q = {q1, …, q n n} respectively, i, n are natural numbers, i ≤ n.

10. The image processing method of claim 9, wherein, The singular value decomposition method is performed by formula (10-1): [U, S, V] = SVD(H) formula (10-1) wherein U, S, and V represent the first unitary matrix, the diagonal matrix, and the second unitary matrix after singular value decomposition, respectively, SVD represents singular value decomposition, and the covariance matrix H between the re-centered point sets is determined by formula (10-2): where p i and q i respectively represent the re-centered three-dimensional feature points.

11. The image processing method of claim 10, wherein, The rotation matrix R and the translation matrix T are determined by formulas (11-1), (11-2), respectively: R = VU T Equation (11-1) T = -Rμ P + μ Q Equation (11-2) where μ P and μ Q are the center positions of the set P, the set Q, respectively.

12. The image processing method of claim 1, wherein, The image includes CT images, MR images, PET images, PET-CT images, PET-MR images and CT-MR images.

13. The image processing method of claim 1, wherein, The PET image is corrected using the RT matrix to obtain a reconstructed image, including: PET scanning is performed on the imaging target to obtain sinogram data; The sinogram data is corrected using the RT matrix to obtain a reconstructed image.

14. The image processing method of claim 1, wherein, The step of obtaining the RT matrix includes: A plurality of radioactive point sources are arranged on the imaging target, the radioactive point sources are identified by reconstructing the images in the acquisition time frame, the pose of the imaging target is obtained according to the plane formed by three non-collinear radioactive point sources, the pose includes the position of the imaging target and the rotation angle of the imaging target around the three coordinate axes, and the RT matrix is obtained by calculating the motion between the continuous frames.

15. An image processing apparatus characterized by comprising: The image processing method according to any one of claims 1-14 is used to work, and the image processing device includes: A frame image acquisition unit for photographing an imaging target and acquiring an initial frame image and a continuous frame image; An initial three-dimensional feature point calculation unit for calculating initial three-dimensional feature points of the initial frame image; A continuous frame three-dimensional feature point calculation unit for calculating continuous frame three-dimensional feature points of the continuous frame image; An RT matrix calculation unit for matching the initial three-dimensional feature points and the continuous frame three-dimensional feature points to obtain an RT matrix; An image reconstruction unit for correcting an image using the RT matrix to obtain a reconstructed image.

16. An electronic device, comprising: It includes: One or more processors; A storage device for storing a computer program; When the computer program is executed by the one or more processors, the one or more processors implement the method of any one of claims 1-14.

17. A computer readable storage medium characterized in that, Program instructions are stored thereon, which are executed to implement the method of any one of claims 1-14.

Images