WO2023082306A1

WO2023082306A1 - Image processing method and apparatus, and electronic device and computer-readable storage medium

Info

Publication number: WO2023082306A1
Application number: PCT/CN2021/131635
Authority: WO
Inventors: 李鑫宇; 肖鹏; 兰嘉友; 董超群; 谢庆国
Original assignee: 苏州瑞派宁科技有限公司
Priority date: 2021-11-12
Filing date: 2021-11-19
Publication date: 2023-05-19
Also published as: CN114170146A

Abstract

Provided in the present application are an image processing method and apparatus, and an electronic device and a computer-readable storage medium. The image correction method comprises: photographing an imaging target and acquiring an initial image frame and continuous image frames; calculating initial three-dimensional feature points of the initial image frame; calculating continuous frame three-dimensional feature points of the continuous image frames; matching the initial three-dimensional feature points and the continuous frame three-dimensional feature points, so as to obtain an RT matrix; and correcting an image by using the RT matrix, so as to obtain a reconstructed image. According to the example embodiments of the present application, motion information of an object is acquired by using an image recording device which is simple but has higher precision, thereby improving the accuracy of a motion corrected image, and reducing the calculation complexity.

Description

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

technical field

The present application relates to the field of image processing, and in particular, relates to an image processing method, device, electronic device and computer-readable storage medium for rigid body motion.

Background technique

At present, when correcting images collected by awake imaging targets, the following two types of motion are mainly corrected: non-rigid motion and rigid motion. Among them, the non-rigid motion mainly refers to the motion of a small range, for example, the movement of the mouse trunk, while the rigid motion mainly refers to the motion of a large range, such as the motion of the head of the mouse.

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

Motion correction methods based on the image domain usually use the Motion Incorporated Reconstruction (MIR) algorithm combined with motion information. The MIR algorithm was proposed by Feng Qiao and Tinsu Pan et al. in 2006. It is a voxel-based motion correction algorithm. The MIR algorithm is based on two assumptions: one is that the imaging target only moves between adjacent data frames, and the imaging target remains stationary within each data frame; the other is that the imaging target is composed of countless spatial points, and the The size is infinitesimally small, and the concentration of radioactive substances in these spatial points does not change with the movement of these points, that is, during the data acquisition process, the concentration of radioactive substances in each part of the imaging target is constant, and the concentration of radioactive substances only increases with the space of the imaging target. Movement and redistribution in space.

The image processing process using the motion correction method based on the image domain includes the following steps.

First, by using the image registration technology between the CT images corresponding to each data frame on the PET (Positron Emission Tomography, PET) image, the mapping relationship matrix M _t between image voxels is obtained, also called the motion matrix. Among them, the registration of CT images is realized by using a non-rigid transformation model based on cubic B-splines and a mean square error similarity criterion.

The motion matrix M _t , and the distribution of the concentration of radioactive substances in the field of view (Field of View, FOV for short) area satisfy the formula (1).

X _t = M _t X ₀ formula (1)

Wherein, X _t represents the distribution of the concentration of radioactive substances in the FOV area at time t, and X ₀ represents the distribution of the concentration of radioactive substances in the FOV area at time 0.

Secondly, construct the imaging model of the PET system for the stationary object, as shown in formula (2).

E(Y)＝GX formula (2)

Among them, 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; G is the system response matrix. In addition, two-dimensional and three-dimensional images can also be converted to one-dimensional images for the above-mentioned processing.

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

Calculating the system response matrix G using the actual measurement method includes: measuring the same point source at different positions of the PET system imaging field of view, 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 physics of the process is detected, and the response is parameterized to correct for point source locations and smooth projection noise.

Calculating the system response matrix G using the Monte Carlo simulation method includes: in space, according to the size of the voxel and its position in the space, simulate each voxel in the field of view, and obtain the coincidence count and the total number of emitted photons. Ratio, the probability that this voxel is detected by the current LineOfResponse (LOR) is obtained, which is the system response matrix G.

Calculation methods based on geometric models, that is, mathematical analysis methods, are divided into many types due to different geometric models, such as line tracing, solid angle and other methods. Line tracing is a simple and effective method for calculating the system response matrix G. This method abstracts the response line as an actual ray, and takes the length of the ray passing through each voxel as the probability that the photon emitted by the voxel is detected by the response line, so as to obtain the system response matrix G. The solid corner model is a more accurate and complex model than the line tracing method. In this method, the calculated vertical angle formed by each voxel in the space and a pair of detectors is taken as the probability that the photons emitted by the voxel are detected by the pair of detectors, so as to obtain the system response matrix G. Compared with the line tracing method, the solid angle model is more in line with the actual PET imaging process.

Again, based on the assumption that the imaging target remains stationary within each data frame, formula (2) is modified into formula (3).

E{Y _t }＝GM _t X ₀ formula (3)

Since all data frames correspond to image sequences X ₀ ,…,X _m-1 , projection data Y ₀ ,…,Y _m-1 and motion matrices M ₀ ,…,M _m-1 , where m is the Quantity, therefore, Equation (2) can be modified to the discrete form shown in Equation (4).

Where M ₀ ≡I is the identity matrix. make:

G ^* ＝[(GM ₀ ) ^T ,(GM ₁ ) ^T ,...,(GM _m-1 ) ^T ] ^T formula (6)

Equation (4) can be rewritten as Equation (7), which is consistent with the form of Equation (2) of PET imaging model for stationary objects.

E(Y ^* )＝G ^* _X0 formula (7)

Since the obtained imaging model of the PET system for the imaging target is in the same form as the imaging model of the PET system for the static target, the 4D projection data Y ^* in formula (5) and formula (6) and the y* of each data frame can be used The system response matrix G ^* of motion information replaces the system response matrix and projection data matrix in formula (2). Therefore, the voxel-based iterative reconstruction algorithm MLEM, which is applied to the PET imaging model formula (2) of a stationary object, is also applicable to formula (7).

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

The value of voxel j, Y _i is the coincidence event count on the i-th line of response (LOR) actually measured, _Gij represents the probability that the gamma photon emitted by voxel j is detected by the scintillation crystals at both ends of the i-th LOR , M represents the total number of lines of response (LOR), N represents the sum of voxels, and j ₀ is the voxel of the inner layer.

Substituting Y ^* and G ^* into the MLEM iterative update equation can get the update equation based on formula (7), as shown in formula (9).

Among them, X ^k is the estimation of the reference position image obtained in the kth iteration, and S is the sensitivity image, which can be calculated by formula (10).

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

It can be seen from the formula (9) that the data of all data frames are involved in the reconstruction of the reference image.

The disadvantage of the existing technology for image correction is that when obtaining the mapping relationship matrix M _t , only one medical imaging device can be used, and an additional 4D CT device is required for CT scanning, and the movement of the object in the two devices needs to be consistent . It not only increases the cost of imaging, but also limits the algorithm to only target regular small-scale movements, such as breathing, heartbeat and other movements. It is also difficult to ensure that the movement of the imaging target in the two devices is completely consistent, and image registration Usually involves an optimization process with a large number of parameters. Therefore, the method of obtaining mapping information between voxels through image registration technology is only suitable for small-scale motion correction, and the accuracy of the obtained motion information cannot be guaranteed when the motion range is large. In addition, when using CT images for registration, a large amount of calculation is required, and the accuracy is not high.

Contents of the invention

This application proposes an image processing method, device, electronic device and computer-readable storage medium for rigid body motion, to solve the problems of high acquisition cost, complex calculation and low accuracy when acquiring motion information of an active imaging target .

According to one aspect of the present application, an image processing method is proposed, the image processing method includes: photographing an imaging target and acquiring an initial frame image and continuous frame images; calculating initial three-dimensional feature points of the initial frame images; calculating the continuous Continuous frame 3D feature points of frame images; matching the initial 3D feature points with the continuous frame 3D feature points to obtain an RT matrix; using the RT matrix to correct the image to obtain a reconstructed image.

According to some embodiments of the present application, the image recording device is used to shoot the imaging target and obtain the initial frame image and continuous frame images. The image recording device includes any device that can record the imaging target plane image or 3D image, such as a camera, DV, video camera, 3D scanner etc.

According to some embodiments of the present application, a pair of image recording devices of the same model are placed in parallel on both sides of the imaging target to form a binocular vision system to photograph the imaging target.

According to some embodiments of the present application, the image recording device is used to record the motion of the imaging target in the form of a video, and the video includes an initial frame image at an initial moment and continuous frame images at a non-initial moment, wherein the initial The frame images and the continuous frame images respectively include images captured by a pair of image recording devices at the same time.

According to some embodiments of the present application, the step of calculating the initial three-dimensional feature points of the initial frame image includes: using a scale-invariant feature transformation method to extract two-dimensional feature points of the initial frame image; removing two-dimensional feature points of the initial frame image The background image of the three-dimensional feature points; matching the two-dimensional feature points of the initial frame image after background removal; triangulating the matching information, thereby obtaining the initial three-dimensional feature points of the initial frame image.

According to some embodiments of the present application, the calculating the continuous frame 3D feature points of the continuous frame images includes: using a scale invariant feature transformation method to extract the continuous frame image 2D feature points; removing the continuous frame images The background image of the two-dimensional feature points; matching the two-dimensional feature points of the continuous frame images after the background is removed; and triangulating the matching information, so as to obtain the initial three-dimensional feature points of the continuous frame images.

According to some embodiments of the present application, the step of triangulating the matching information includes: calculating the depth Z of any point P on the imaging target under the coordinates of the image recording device by formula (7-1); by formula (7-2) Calculate the abscissa X and ordinate Y of any point P on the imaging target under the camera coordinates;

Respectively, the abscissa distance between the projection of point P on the image frames of the two image recording devices and the intersection point of the optical axis image; y1 is the intersection point of the projection of point P on the image frame of one of the image recording devices and the optical axis image The ordinate distance of .

According to some embodiments of the present application, the step of obtaining the RT matrix includes: minimizing errors between the initial 3D feature points and the 3D feature points of consecutive frames by using a singular value decomposition method to obtain the RT matrix.

According to some embodiments of the present application, the error is minimized by formula (9-1):

Among them, R and T represent the rotation matrix and translation matrix respectively, J represents the error, p _i and q _i represent the set of initial three-dimensional feature points P={p ₁ ,...,p _n } and the continuous frame The ith element in the set Q={q ₁ ,...,q _n } of three-dimensional feature points, i and n are both natural numbers, i≤n.

According to some embodiments of the present application, the singular value decomposition (Singular Value Decomposition, referred to as SVD) method is carried out by formula (10-1):

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

Among them, U, S, and V respectively represent the first unitary matrix, diagonal matrix and second unitary matrix after singular value decomposition, SVD represents the singular value decomposition of the physical quantities in brackets, and H represents the re-centered point set. The covariance matrix between , the covariance matrix H is determined by the formula (10-2):

Wherein, p _i ' and q _i ' respectively denote the re-centered three-dimensional feature points, p _i '=pi _-μ _P , q _i '=q _i -μ _Q .

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

R = V U ^T formula (11-1)

Among them, μ _P and μ _Q are the central positions of set P and set Q respectively.

According to some embodiments of the present application, the images include CT images, MR images, PET images, PET-CT images, PET-MR images, and CT-MR images.

According to some embodiments of the present application, using the RT matrix to correct the PET image to obtain a reconstructed image includes: performing a PET scan on the imaging target to obtain sinogram data; using the RT matrix to correct the sinusoid image data to reconstruct the image.

According to some embodiments of the present application, before using the RT matrix to correct the image, the image processing method further includes correcting the voxel value of the image by formula (14-1):

X _t =W _t RTX ₀ Formula (14-1)

Among them, X _t represents the distribution of radioactive material concentration in the FOV area at time t, X ₀ represents the distribution of radioactive material concentration in the FOV area at the initial time, and W _t is the interpolation weight of all voxel positions at time t.

According to some embodiments of the present application, the step of obtaining the RT matrix includes: setting several radioactive point sources on the imaging target, and identifying the radioactive point sources by reconstructing images within the acquisition time frame, according to three The plane composed of the non-collinear radioactive point sources obtains the pose of the imaging target, and the pose includes the position of the imaging target and the rotation angle of the imaging target around the three-dimensional coordinate axis, by calculating The movement between gets the RT matrix.

According to one aspect of the present application, an image processing device is proposed. The image processing device includes a frame image acquisition unit for capturing an imaging target and acquiring initial frame images and continuous frame images; an initial three-dimensional feature point calculation unit for calculating The initial three-dimensional feature point of the initial frame image; the continuous frame three-dimensional feature point calculation unit, used to calculate the continuous frame three-dimensional feature point of the continuous frame image; the RT matrix calculation unit, used to match the initial three-dimensional feature point and the The three-dimensional feature points of the continuous frames are used to obtain an RT matrix; the image reconstruction unit is used to correct the image by using the RT matrix to obtain a reconstructed image.

According to one aspect of the present application, an electronic device is proposed, including 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 as in any one of the preceding.

According to one aspect of the present application, a computer-readable storage medium is provided, on which program instructions are stored, and when the program instructions are executed, the method described in any one of the preceding items is implemented.

According to some exemplary embodiments of the present application, a simple but higher-precision image recording device is used to obtain motion information of an object while it is moving, which ensures the accuracy of motion information acquisition and improves the accuracy of motion-corrected images. At the same time, a large number of complicated image registration calculations are avoided, and the complexity of calculations is reduced. Moreover, the provided algorithm is not only applicable to small-range motion correction, but also applicable to larger-range motion correction, which greatly reduces the cost.

Description of drawings

In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following briefly introduces the drawings that need to be used in the description of the embodiments.

Fig. 1 shows a flowchart of an image processing method according to an example embodiment of the present application.

Fig. 2 shows a schematic diagram of triangulation according to an example embodiment of the present application.

Fig. 3 shows a schematic diagram of a mapping process of voxels on a two-dimensional graph according to an exemplary embodiment of the present application.

Fig. 4 shows an image correction device for rigid body motion according to an exemplary embodiment of the present application.

Figure 5a shows the correction results when the object is stationary.

Figure 5b shows the reconstruction results before object motion correction.

Fig. 5c shows a reconstruction result after object motion correction according to an exemplary embodiment of the present application.

Fig. 6 shows a block diagram of another electronic device for rigid body motion according to an embodiment of the present application.

Detailed ways

Example embodiments will now be described more fully with reference to the accompanying drawings. Example embodiments may, however, be embodied in many forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided so that this application will be thorough and complete, and will fully convey the concept of example embodiments to those skilled in the art. The same reference numerals denote the same or similar parts in the drawings, and thus their repeated descriptions will be omitted.

The described features, structures, or characteristics may be combined in any suitable manner in one or more embodiments. In the following description, numerous specific details are provided in order to give a thorough understanding of embodiments of the present disclosure. However, those skilled in the art will appreciate that the technical solutions of the present disclosure may be practiced without one or more of these specific details, or other methods, components, materials, devices or operations may be used. In these instances, well-known structures, methods, devices, implementations, materials, or operations are not shown or described in detail.

The flow charts shown in the drawings are only exemplary illustrations, and do not necessarily include all contents and operations/steps, nor must they be performed in the order described. For example, some operations/steps can be decomposed, and some operations/steps can be combined or partly combined, so the actual order of execution may be changed according to the actual situation.

The terms "first", "second" and the like in the description and claims of the present application and the above drawings are used to distinguish different objects, rather than to describe a specific order. Furthermore, the terms "include" and "have", as well as any variations thereof, are intended to cover a non-exclusive inclusion. For example, a process, method, system, product or device comprising a series of steps or units is not limited to the listed steps or units, but optionally also includes unlisted steps or units, or optionally further includes For other steps or units inherent in these processes, methods, products or devices.

Specific embodiments according to the present application will be described in detail below with reference to the accompanying drawings.

Fig. 1 shows a flowchart of an image processing method according to an example embodiment of the present application. An image processing method according to an exemplary embodiment of the present application will be described in detail below with reference to FIG. 1 .

According to some exemplary embodiments of the present application, the exemplary embodiment shown in FIG. 1 corrects motion artifacts generated during PET scanning of awake imaging targets. Since the final imaging of the instrument is only aimed at the head of the imaging target, the movement of the head of the active imaging target can be simplified as a rigid body motion.

In step S101, an imaging target is photographed and initial frame images and continuous frame images are acquired.

According to some exemplary embodiments of the present application, the movement of the imaging target is recorded in the form of video by using a pair of image recording devices placed in parallel on both sides of the imaging target. The image recording device includes a plane image or a 3D image capable of recording the imaging target Any device, such as camera, DV, camcorder, 3D scanner, etc. The video captured by a pair of image recording devices includes an initial frame image at an initial moment and continuous frame images at a non-initial moment, wherein the initial frame image and continuous frame images respectively include images captured by a 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 in parallel for the purpose of unifying the coordinate system in the post-calculation process and reducing the computational complexity. In fact, the image recording equipment can be set at any angle and in any number, and the corresponding coordinate system conversion can be performed during data processing in the later stage.

According to some embodiments, the initial frame image is the frame image at the beginning moment of the video captured by the image recording device, and the images at other moments in the video are continuous frame images.

According to some embodiments, in step S101, a pair of image recording devices of the same model, such as cameras, are placed in parallel on the left and right sides of the imaging target to form a binocular vision system to record information about object movement in the form of video.

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

According to some exemplary embodiments of the present application, the initial three-dimensional feature points of the initial frame image are calculated through the following steps:

First, use the scale-invariant feature transform (SIFT for short) method to extract the two-dimensional feature points of the initial frame image from the initial frame image.

Then, use the digital subtraction angiography (Digital subtraction angiography, referred to as DSA) principle to remove the background image of the two-dimensional feature points of the initial frame image to ensure that the extracted two-dimensional feature points are all on the imaging target image in the initial frame image .

Finally, match the two-dimensional feature points in the initial frame images captured by a pair of image recording devices, and triangulate the matching information, so as to obtain the initial three-dimensional feature points of the initial frame images.

In step S105, three-dimensional feature points of consecutive frames of images are calculated.

In step S105, the method for calculating the three-dimensional feature points of the continuous 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 repeated here.

The principle of triangulation in step S103 and step S105 will be described below. The so-called triangulation refers to recovering the three-dimensional position information of the feature points of the object according to the two-dimensional projection positions of different image recording devices through the feature points of the object. As shown in FIG. 2 , it is a schematic diagram of triangulation according to an exemplary embodiment of the present application.

As shown in Figure 2, let the coordinates of the intersection point between the camera optical axis and the image be (cx, cy), the distance between the optical centers of the two cameras is D, the focal length is f, and a feature point P in the three-dimensional space is in the left and right images The projections on are respectively Pl and Pr, the abscissa distances between them and the intersection point of the optical axis image are xl and xr respectively, and the ordinate distances are yl and yr.

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

In step S107, the initial 3D feature points are matched with the 3D feature points of consecutive frames to obtain an RT matrix.

According to some example embodiments of the present application, the step of calculating the RT matrix includes:

According to step S103 and step S105, the set of initial 3D feature points P={p ₁ ,...,p _n } of the initial frame image and the set of 3D feature points of consecutive frames Q={q ₁ ,...,q _n } are respectively obtained.

The error J of 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 formula (13).

Among them, J is the error between the initial three-dimensional feature point and the three-dimensional feature point of the continuous frame, R is the rotation matrix, Rp _i represents the multiplication of the rotation matrix R and the i-th point p _i in the point set P, T is the translation matrix, p _i , q _i represent the i-th three-dimensional feature point in the sets P and Q respectively, i is a natural number less than n, and n is the number of three-dimensional feature points in the sets P and Q.

Use the formula (14) to average the coordinates of each point in the set P and the set Q to obtain the center position of the set P and the set Q

Obtain new point sets P' and Q', that is, the coordinates of each point minus the center position. At this time, the center position coincides with the origin of the coordinates. 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 .

Singular value decomposition is performed on the covariance matrix H by formula (17).

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

Among them, U, S, and V respectively represent the first unitary matrix, diagonal matrix and second unitary matrix after singular value decomposition, and SVD represents the singular value decomposition of the physical quantities in brackets. For example, if the covariance matrix H is m ×n order matrix, then the first unitary matrix U after singular value decomposition is a unitary matrix of order m×m, the diagonal matrix S is a diagonal matrix of order m×n, and the elements on the diagonal are covariance matrices The singular value of H, the second unitary matrix V is a unitary matrix of order n×n.

Through the formula (18), the rotation matrix R is calculated according to the decomposed result.

R = V U ^T formula (18)

Through formula (19), the translation matrix T is calculated using the rotation matrix R obtained from formula (18).

The R and T matrices obtained by formula (18) and formula (19) represent the six-axis motion information of the rigid body motion of the object. The rigid body motion of an object can be decomposed into translational and rotational motions around three axes in space, namely x, y, and z axes, and the obtained motion of the object around the three axes can form an RT matrix. The RT matrix is also called a rotation-translation (Rotation-Translation, RT for short) matrix, and is used to describe six-axis motion information of an object. The six axes are translational and rotational movements in the directions of x, y and z axes.

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

According to some embodiments of the present application, the step of obtaining the RT matrix includes: setting several radioactive point sources on the imaging target, and identifying the radioactive point sources by reconstructing images within the acquisition time frame, according to three non-collinear radioactive point sources The pose of the imaging target is obtained from the surface composed of sources, the pose includes the position of the imaging target and the rotation angle of the imaging target around the three-dimensional coordinate axis, and the RT matrix is obtained by calculating the motion between consecutive frames.

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

According to some exemplary embodiments of the present application, the step of correcting the PET image using the RT matrix includes:

Perform PET scans on living targets to obtain sinogram data. Divide the sinogram data according to time frames, and use the RT matrix in step S107 to correct the data in the current time frame to the initial position through formula (9). That is, the RT matrix is used to replace the mapping relationship matrix M _t in formula (9), so as to obtain a corrected reconstructed image.

According to some embodiments of the present application, before step S109 is performed, the voxel value of the image needs to be corrected by formula (20).

X _t =W _t RTX ₀ formula (20)

Among them, X _t represents the distribution of radioactive material concentration in the FOV area at time t, X ₀ represents the distribution of radioactive material concentration in the FOV area at time 0, and W _t is the interpolation weight of all voxel positions at time t.

How to correct the voxel value of the PET image by formula (20) will be described in detail below.

In order to replace the complex image registration process to obtain the mapping relationship between image voxels at different times, the spatial coordinate transformation of the center point of each voxel can be performed according to the 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 corresponding voxel position at the initial time according to the motion information, and the transformed voxel position is used The surrounding voxel values are calculated by an interpolation algorithm to calculate the voxel value of the transformed voxel position, as shown in formula (21). The activity value f( _ri ,0) located at the spatial point r _i is regarded as the weighted sum of the voxel values of all voxels adjacent to the spatial point r _i .

in,

Indicates the subscript set of voxels near the spatial point r _i , and ω _ij is the weight. x _j (0) corresponds to X ₀ in formula (20), and f(r _i ,0) corresponds to X _t in formula (21). The matrix representation of formula (21) is shown in formula (20).

According to some embodiments, the weight value ω _ij is taken as the ratio of the overlapping area or volume between the voxel centered on the position of the spatial point r _i and its nearby voxels to the total area or volume of the voxel.

Similarly, the activity values f(r ₀ ,0),f(r ₁ ,0),…,f(r _n-1 ,0) of continuous space points can be obtained by formula (21).

Fig. 3 shows a schematic diagram of a mapping process of voxels on a two-dimensional graph according to an exemplary embodiment of the present application. As shown in Figure 3, the gray dot in the figure indicates that a voxel center A at time t is mapped to its corresponding original position B at time 0, and the voxel value of the original position B can be compared to B near the position at time 0 The interpolation calculation of voxels is obtained.

According to some embodiments, the RT matrix may be acquired using point labeling techniques. The specific steps include: paste a small and light radioactive point source on the head of the experimental subject (such as a mouse), reconstruct the image within the acquisition time frame, identify the highlighted point in the image, that is, it is considered a point source, and the default point source position Keep still relative to the head of the rat. According to the surface composed of three non-collinear point sources, the head pose of the rat can be obtained. The head pose includes the position of the object and three rotation angles of the object around the xyz three axes, which are used to represent the orientation of the object. The motion matrix RT can be obtained by calculating the motion between consecutive frames.

According to some exemplary embodiments of the present application, a simple but higher-precision external camera is used to obtain motion information of a moving object while it is moving, thereby ensuring the accuracy of motion information acquisition and improving the accuracy of motion-corrected images. At the same time, a large number of complicated image registration calculations are avoided, and the complexity of calculations is reduced. Moreover, the provided algorithm is not only applicable to small-range motion correction, but also applicable to larger-range motion correction, which greatly reduces the cost.

Fig. 4 shows an image processing device according to an exemplary embodiment of the present application. The image processing device shown in FIG. 4 includes a frame image acquisition unit 401 , an initial 3D feature point calculation unit 403 , a continuous frame 3D feature point calculation unit 405 , an RT matrix calculation unit 407 and an image reconstruction unit 409 .

The frame image acquiring unit 401 uses a camera to capture an imaging target and acquire initial frame images and continuous frame images. The method described in the embodiment in FIG. 1 may be used for camera setting and image acquisition, and details are not repeated here.

The initial three-dimensional feature point calculating unit 403 is used to calculate the initial three-dimensional feature point of the initial frame image. The calculation of the initial three-dimensional feature points can use the method described in step S103 in the embodiment of FIG. 1 , which will not be repeated here.

The continuous frame three-dimensional feature point calculation unit 405 is used for calculating continuous frame three-dimensional feature points of the continuous frame images. The calculation of the three-dimensional feature points of consecutive frames can use the method described in step S105 in the embodiment of FIG. 1 , which will not be repeated here.

The RT matrix calculation unit 407 is used to match the initial 3D feature points with the 3D feature points of consecutive frames to obtain an RT matrix. The calculation of the RT matrix can use the method described in step S107 in the embodiment of FIG. 1 , which will not be repeated here.

The image reconstruction unit 409 is configured to use the RT matrix to correct the PET image to obtain a reconstructed image. The method described in step S109 in the embodiment of FIG. 1 can be used for image reconstruction, which will not be repeated here.

Figure 5a shows the correction results when the object is stationary. Fig. 5b shows a reconstruction result before object motion correction, and Fig. 5c shows a reconstruction result after object motion correction according to an exemplary embodiment of the present application. It can be seen from FIGS. 5 a to 5 c that by using the image processing method provided by the present application, artifacts caused by object motion can be basically eliminated, thereby improving the quality of the reconstructed image.

It should be noted by those skilled in the art that the PET image is corrected by the RT matrix in the above embodiment. In fact, the above method can also be applied to motion correction of other medical images such as CT images and MR images, or PET-CT images. , PET-MR images, CT-MR images in one or more images of motion correction, in the correction process of these medical images only need to replace the corresponding M _t matrix by RT matrix, which belongs to those skilled in the art What can be easily realized according to the teaching of the present invention will not be repeated here.

Fig. 6 shows a block diagram of another electronic device for rigid body motion according to an embodiment of the present application. The electronic device shown in FIG. 6 is only an example, and should not limit the functions and scope of use of this embodiment of the present application.

As shown in FIG. 6, the electronic device takes the form of a general-purpose computing device. Components of the 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 the memory 220 and the processor 210), a display unit 240, and the like. Wherein, the memory 220 stores program codes, and the program codes can be executed by the processor 210, so that the processor 210 executes the methods described in this specification according to various exemplary embodiments of the present application. For example, the processor 210 may execute the method as shown in FIG. 1 .

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

Memory 220 may also include programs/utilities 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, which Each or some combination of the examples may include the implementation of a network environment.

Bus 230 may represent one or more of several types of bus structures, including a memory cell bus or memory cell controller, a peripheral bus, an accelerated graphics port, a processing unit, or a local area using any of a variety of bus structures. bus.

The electronic device can also communicate with one or more external devices 300 (such as keyboards, pointing devices, bluetooth devices, etc.), and can also communicate with one or more devices that allow the user to interact with the image correction device, and/or communicate with the device that enables the user to interact with the image correction device The electronic device is capable of communicating with any device (eg, router, modem, etc.) that communicates with one or more other computing devices. Such communication may occur through input/output (I/O) interface 250 . Moreover, the electronic device can also communicate with one or more networks (eg, a local area network (LAN), a wide area network (WAN) and/or a public network such as the Internet) through the network adapter 260 . The network adapter 260 can communicate with other modules of the electronic device through the bus 230 . It should be understood that although not shown in the figures, other hardware and/or software modules may 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 system, etc.

Through the description of the above implementations, those skilled in the art can easily understand that the example implementations described here can be implemented by software, or can be implemented by combining software with necessary hardware. The technical solutions according to the embodiments of the present application can be embodied in the form of software products, which can be stored in a computer-readable storage medium (which can be CD-ROM, U disk, mobile hard disk, etc.) or on the network, including several The computer program instructions enable a computing device (which may be a personal computer, a server, or a network device, etc.) to execute the above method according to the embodiments of the present application.

A software product may utilize any combination of one or more readable media. The readable medium may be a readable signal medium or a readable storage medium. The readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, device, or device, or any combination thereof. More specific examples (non-exhaustive list) of readable storage media include: electrical connection with one or more conductors, portable disk, hard disk, 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 of the foregoing.

A computer readable storage medium may include a data signal carrying readable program code in baseband or as part of a carrier wave traveling as part of a data signal. Such propagated data signals may take many forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination of the foregoing. A readable storage medium may also be any readable medium other than a readable storage medium that can send, propagate or transport a program for use by or in conjunction with an instruction execution system, apparatus or device. The program code contained on the readable storage medium may be transmitted by any suitable medium, including but not limited to wireless, cable, optical cable, RF, etc., or any suitable combination of the above.

Program codes for performing the operations of the present application can be written in any combination of one or more programming languages, including object-oriented programming languages—such as Java, C++, etc., as well as conventional procedural programming Language—such as C or a similar programming language. The program code may execute entirely on the user's computing device, partly on the user's device, as a stand-alone software package, partly on the user's computing device and partly on a remote computing device, or entirely on the remote computing device or server to execute. In cases involving a remote computing device, the remote computing device may be connected to the user computing device through any kind of network, including a local area network (LAN) or a wide area network (WAN), or may be connected to an external computing device (for example, using an Internet service provider). business to connect via the Internet).

The above-mentioned computer-readable medium carries one or more program instructions, and when the above-mentioned one or more program instructions are executed by one device, the computer-readable medium can realize the above-mentioned functions.

Those skilled in the art can understand that the above modules can be distributed in the device according to the description of the embodiment, or can be changed correspondingly to one or more devices that are only different from the embodiment. Multiple modules in the above embodiments may be combined into one module, or one module may be further split into multiple sub-modules.

Through the description of the above implementations, those skilled in the art can easily understand that the example implementations described here can be implemented by software, or by combining software with necessary hardware. The technical solutions according to the embodiments of the present application can be embodied in the form of software products, which can be stored in a computer-readable storage medium (which can be CD-ROM, U disk, mobile hard disk, etc.) or on the network, including several The computer program instructions enable a computing device (which may be a personal computer, a server, or a network device, etc.) to execute the above method according to the embodiments of the present application.

Those skilled in the art can understand that the above-mentioned modules can be distributed in the device according to the description of the embodiment, and corresponding changes can also be made in one or more devices that are only different from the embodiment. Multiple modules in the above embodiments may be combined into one module, or one module may be further split into multiple sub-modules.

According to some exemplary embodiments of the present application, a simple but higher-precision external camera is used to obtain motion information of a moving object while it is moving, thereby ensuring the accuracy of motion information acquisition and improving the accuracy of motion-corrected images. At the same time, a large number of complex image registration calculations are avoided, and the complexity of calculations is reduced. Moreover, the provided algorithm is not only applicable to small-range motion correction, but also applicable to larger-range motion correction, which greatly reduces the cost.

Although the present application provides the operation steps of the method described in the above embodiments or flowcharts, more or fewer operation steps may be included in the method based on routine or no creative effort. In steps that do not logically have a necessary causal relationship, the execution order of these steps is not limited to the execution order provided in the embodiment of the present application.

Each embodiment in this specification is described in a progressive manner, the same and similar parts of each embodiment can be referred to each other, and each embodiment focuses on the differences from other embodiments.

The above is a detailed introduction to the embodiments of the present application. In this paper, specific examples are used to illustrate the principles and implementation methods of the present application. The descriptions of the above embodiments are only used to help understand the methods and core ideas of the present application. At the same time, changes or deformations made by those skilled in the art based on the ideas of the application, specific implementation methods and application scopes of the application all belong to the scope of protection of the application. To sum up, the contents of this specification should not be understood as limiting the application.

Claims

An image processing method, characterized in that the image processing method comprises:

Shoot the imaging target and obtain the initial frame image and continuous frame images;

calculating initial three-dimensional feature points of the initial frame image;

calculating the continuous frame three-dimensional feature points of the continuous frame images;

matching the initial three-dimensional feature points and the three-dimensional feature points of the continuous frames to obtain an RT matrix;

The image is corrected by using the RT matrix to obtain a reconstructed image.
The image processing method according to claim 1, wherein an image recording device is used to shoot the imaging target and acquire initial frame images and continuous frame images.
The image processing method according to claim 2, wherein a pair of image recording devices of the same model are placed in parallel on both sides of the imaging target to form a binocular vision system to photograph the imaging target.
The image processing method according to claim 2, wherein the image recording device is used to record the motion of the imaging target in the form of a video, the video including an initial frame image at an initial moment and continuous frames at a non-initial moment images, wherein the initial frame image and the continuous frame images respectively include images taken by a pair of the image recording devices at the same time.
The image processing method according to claim 1, wherein the step of calculating the 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 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 background removal;

Triangulating the matching information, so as to obtain the initial three-dimensional feature points of the initial frame image.
The image processing method according to claim 1, wherein the calculating the continuous frame three-dimensional feature points of the continuous frame images comprises:

Extracting two-dimensional feature points of the continuous frame images by using a scale-invariant feature transformation method;

removing the background image of the two-dimensional feature points of the continuous frame images;

matching the two-dimensional feature points of the continuous frame images after background removal;

Triangulating the matching information, so as to obtain the initial three-dimensional feature points of the continuous frame images.
The image processing method according to claim 5 or 6, wherein the step of triangulating matching information comprises:

Calculate the depth Z of any point on the imaging target under the coordinates of the image recording device by formula (7-1);

Calculate the abscissa X and the ordinate Y of any point on the imaging target under the coordinates of the image recording device by formula (7-2);

Wherein, T is the optical distance between the two image recording devices, f is the focal length of the image recording devices, x1 and xr are respectively the projection and optical axis of the point on the image frames of the two image recording devices The abscissa distance of the image intersection point; y1 is the ordinate distance between the projection of the point on the image frame of one of the image recording devices and the intersection point of the optical axis image.
The image processing method according to claim 1, wherein the step of obtaining the RT matrix comprises:

The errors between the initial three-dimensional feature points and the three-dimensional feature points of the consecutive frames are minimized by a singular value decomposition method to obtain the RT matrix.
The image processing method according to claim 8, wherein the error is minimized by formula (9-1):

Among them, R and T represent the rotation matrix and translation matrix respectively, J represents the error, p i and q i represent the set of initial three-dimensional feature points P={p 1 ,...,p n } and the continuous frame The ith element in the set Q={q 1 ,...,q n } of three-dimensional feature points, i and n are both natural numbers, i≤n.
The image processing method according to claim 9, wherein the singular value decomposition method is carried out by formula (10-1):

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

Among them, U, S, and V respectively represent the first unitary matrix, diagonal matrix and second unitary matrix after the singular value decomposition, SVD represents the singular value decomposition, and the covariance matrix H between the re-centered point sets is passed by the formula (10-2) Determine:

Among them, p i ' and q i ' represent the three-dimensional feature points after re-centering respectively.
The image processing method according to claim 10, wherein the rotation matrix R and the translation matrix T are respectively determined by formulas (11-1), (11-2):

R = V U T formula (11-1)

T＝-Rμ P +μ Q formula (11-2)

Among them, μ P and μ Q are the central positions of set P and set Q respectively.
The image processing method according to claim 1, wherein the images include CT images, MR images, PET images, PET-CT images, PET-MR images and CT-MR images.
The image processing method according to claim 1, wherein said correcting the PET image using said RT matrix to obtain a reconstructed image comprises:

performing a PET scan on the imaging target to obtain sinogram data;

Correcting the sinogram data by using the RT matrix to reconstruct an image.
The image processing method according to claim 1, characterized in that before using the RT matrix to correct the image, the image processing method further comprises using the formula (14-1) to correct the voxels of the image The value is corrected:

X t =W t RTX 0 Formula (14-1)

Among them, X t represents the distribution of radioactive material concentration in the FOV area at time t, X 0 represents the distribution of radioactive material concentration in the FOV area at the initial time, and W t is the interpolation weight of all voxel positions at time t.
The image processing method according to claim 1, wherein the step of obtaining the RT matrix comprises:

Setting several radioactive point sources on the imaging target, identifying the radioactive point sources by reconstructing images within the acquisition time frame, and obtaining the imaging target according to a surface composed of 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-dimensional coordinate axis, and the RT matrix is obtained by calculating the motion between consecutive frames.
An image processing device, characterized in that the image correction device includes:

A frame image acquisition unit, configured to photograph an imaging target and acquire initial frame images and continuous frame images;

an initial three-dimensional feature point calculation unit, configured to calculate the initial three-dimensional feature point of the initial frame image;

A continuous frame 3D feature point calculation unit, configured to calculate the continuous frame 3D feature points of the continuous frame images;

An RT matrix calculation unit, configured to match the initial three-dimensional feature points and the three-dimensional feature points of the continuous frames to obtain an RT matrix;

An image reconstruction unit, configured to use the RT matrix to correct the image to obtain a reconstructed image.
An electronic device, characterized in that it comprises:

one or more processors;

storage means for storing computer programs;

When the computer program is executed by the one or more processors, the one or more processors are made to implement the method according to any one of claims 1-15.
A computer-readable storage medium, wherein program instructions are stored thereon, and the method according to any one of claims 1-15 is implemented when the program instructions are executed.