Heart T1 quantitative differential homeomorphic registration method based on differential model and physical constraint

Abstract

The invention discloses a heart T1 quantitative differential homeomorphic registration method based on a differential model and physical constraints. The method comprises the following steps: acquiring a heart T1 weighted image of a target object; the heart T1 weighted image is input into a trained motion correction network, a corrected image is obtained, the motion correction network comprises an image synthesis sub-network and a deformable registration sub-network, and the image synthesis sub-network and the deformable registration sub-network are both constructed based on a UNet architecture; the image synthesis sub-network is used for fitting a heart T1 image and a heart M0 image by using the registered heart T1 weighted image output by the deformable registration sub-network, and then synthesizing a heart T1 weighted image without relative motion by using the fitted heart T1 image, the heart M0 image and a signal attenuation model; the deformable registration sub-network is configured to perform intensity-based registration using an output image of the image synthesis sub-network and an original image. According to the invention, the image registration effect and efficiency can be improved at the same time.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, and more specifically, to a method for quantitative differential homeomorphic registration of cardiac T1 based on differential models and physical constraints. Background Technology

[0002] Myocardial longitudinal relaxation time (T1) serves as a marker for various myocardial pathophysiological symptoms such as myocardial inflammation, diffuse fibrosis, and infiltration, and is therefore useful for diagnosing related myocardial diseases. Quantifying myocardial T1 values ​​requires pixel-wise fitting of multiple cardiac T1-weighted images acquired at different reversal times within the same slice slice using a nonlinear curve fitting model. However, a patient's respiration and cardiac motion will cause varying degrees of misalignment in the spatial location and shape of the heart among the cardiac T1-weighted images within the same slice slice. These misalignments will lead to errors in quantifying myocardial T1 values. Typically, image registration techniques are used to align these cardiac T1-weighted images for more accurate quantification. However, these cardiac T1-weighted images are multi-contrast images with significant contrast differences. This makes the registration task particularly difficult; general intensity-based registration methods are ineffective, and general multimodal methods struggle to solve the registration problem. Furthermore, the ideal medical image registration mapping should be a differential homeomorphic mapping, i.e., the solved deformation field mapping is a continuously differentiable and reversible bijective. Such registration mapping ensures the consistency of the topological structure before and after image registration. Furthermore, incorporating a physical model during registration helps to obtain physically possible deformations along the longitudinal relaxation axis. However, currently, there is a lack of effective schemes for image registration that consider differential homeomorphic registration and physical constraints.

[0003] In recent years, deep learning-based image analysis algorithms have developed rapidly in the field of medical image processing. For example, VoxelMorphp is a deep learning-based medical image registration framework, laying the foundation for deep learning-based medical image registration research. Subsequently, a differential homeomorphic registration network based on a probabilistic differential model was proposed, which can achieve a registration field that more closely matches actual physical deformations. To reduce network training parameters without compromising network performance, thereby lowering training and testing time costs, Heinrich et al. proposed an efficient sparse deformable convolutional module, Obelisk, which can effectively replace certain layers in UNet, reducing network parameters while maintaining good network performance. Currently, for the problem of significant contrast variations in cardiac T1-weighted images, some methods treat it as a multimodal image registration problem for feature-based registration. Others first use sparse weighted representation or disentangled representation learning to transform multi-contrast cardiac T1-weighted images into a series of images with the same contrast for intensity-based registration. In addition, researchers have proposed a model-based cardiac T1 mapping motion correction method that incorporates a signal attenuation model into the registration network to address the distortion problem that may exist in general registration methods and does not conform to actual physics.

[0004] While the aforementioned methods for quantitative T1 registration of the heart achieve good performance, they either neglect the differential homeomorphic property that ideal medical image registration should possess, which may lead to inconsistent anatomical structures in the registered image compared to the original image; or they fail to consider physically possible deformations along the longitudinal relaxation axis, which may result in deformation fields that do not conform to physical reality. Both of these situations can lead to registration errors in cardiac T1-weighted images, thereby causing errors in the quantification of myocardial T1 values. Furthermore, current deep learning-based methods for quantitative T1 registration of the heart employ complex network models with numerous parameters, resulting in slow network training and testing speeds.

[0005] In summary, further improvements are needed to the existing technology to provide a quantitative registration method for cardiac T1 that simultaneously considers differential homeomorphic properties and physical deformation constraints. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of the prior art and provide a quantitative differential homeomorphic registration method for cardiac T1 based on differential models and physical constraints. This method includes the following steps:

[0007] Obtain the T1-weighted image of the heart of the target object;

[0008] The T1-weighted cardiac image is input into a trained motion correction network to obtain a corrected image. The motion correction network includes an image synthesis subnetwork and a deformable registration subnetwork.

[0009] The image synthesis subnetwork and the deformable registration subnetwork are both built on the UNet architecture. The image synthesis subnetwork is used to fit a heart T1 map and a heart M0 map using the registered heart T1 weighted image output from the deformable registration subnetwork. Then, the fitted heart T1 map, heart M0 map and signal attenuation model are used to synthesize a heart T1 weighted image without relative motion. The deformable registration subnetwork is used to perform intensity-based registration using the output image of the image synthesis subnetwork and the original image.

[0010] Compared with the prior art, the advantage of this invention is that it proposes a cardiac T1 quantitative differential homeomorphic registration method based on deep learning and physical models, which is of great significance for improving the registration accuracy of cardiac T1 weighted images and the accuracy of myocardial T1 value quantification.

[0011] Other features and advantages of the invention will become clear from the following detailed description of exemplary embodiments of the invention with reference to the accompanying drawings. Attached Figure Description

[0012] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments of the invention and, together with their description, serve to explain the principles of the invention.

[0013] Figure 1 This is a flowchart of a method for quantitative differential homeomorphic registration of cardiac T1 based on differential models and physical constraints according to an embodiment of the present invention;

[0014] Figure 2 This is a schematic diagram of the overall process of a motion correction network according to an embodiment of the present invention;

[0015] Figure 3 This is a schematic diagram of the overall process of a cardiac T1 quantitative differential homeomorphic registration method based on a differential model and physical constraints according to an embodiment of the present invention. Detailed Implementation

[0016] Various exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be noted that, unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps set forth in these embodiments do not limit the scope of the invention.

[0017] The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the invention or its application or use.

[0018] Techniques, methods, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, such techniques, methods, and equipment should be considered part of the specification.

[0019] In all the examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0020] It should be noted that similar labels and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be discussed further in subsequent figures.

[0021] In summary, this invention proposes a motion correction network incorporating a physical model, comprising an image synthesis subnetwork and a deformable registration subnetwork. Cardiac T1 data is simultaneously input into both subnetworks for "motion-free" image synthesis based on a signal attenuation model and intensity-based image registration, respectively. The image synthesis subnetwork relies on the deformable registration subnetwork to fit optimized cardiac T1 and M0 images with gradually reduced motion artifacts, and then uses the signal attenuation model to synthesize a series of "motion-free" weighted cardiac T1 images from the cardiac T1 and M0 images. Simultaneously, the deformable registration subnetwork relies on the series of "motion-free" weighted cardiac T1 images synthesized by the image synthesis subnetwork to perform intensity-based registration, indirectly addressing the problem of significant contrast variations. Simultaneous training of both subnetworks continuously optimizes their performance, ultimately yielding cardiac T1 images free of motion artifacts, which can be directly used for myocardial T1 value assessment.

[0022] Specifically, see Figure 1 As shown, the provided method for quantitative differential homeomorphic registration of cardiac T1 based on differential models and physical constraints includes the following steps:

[0023] Step S110: Construct a dataset containing multiple samples of cardiac T1-mapped imaging and annotations of the inner and outer myocardial contours.

[0024] In the example experiments, the cardiac T1 mapping dataset T1Dataset210 from the Harvard Dataverse was used as the training and testing dataset for the network. T1Dataset210 was imaged using a 1.5T Philips Achieva system and a 32-channel cardiac coil, containing cardiac T1 mapping images of 210 patients (134 males, aged 57 ± 14 years) with known or suspected cardiovascular disease. Each patient dataset consisted of five short axial slices covering the left ventricle from the base to the top. The inner and outer myocardial contours of all images (N = 11550) in this dataset were manually annotated for results analysis.

[0025] Step S120: Construct a motion correction network, which includes an image synthesis subnetwork and a deformable registration subnetwork. Both the image synthesis subnetwork and the deformable registration subnetwork are built based on the UNet architecture, and the encoder-decoder of the set layer and the bottleneck layer are replaced with Obelisk modules.

[0026] like Figure 2 As shown, the cardiac T1 mapping motion correction network proposed in this invention comprises two parallel branches: (a) a UNet-based image synthesis subnetwork, used to fit a cardiac T1 map and a cardiac M0 map using the registered cardiac T1-weighted image output from the deformable registration subnetwork, and then using the fitted cardiac T1 map, cardiac M0 map, and signal attenuation model to synthesize a series of "motion-free" cardiac T1-weighted images. (b) a UNet-based deformable registration subnetwork, used to perform intensity-based registration using the "motion-free" cardiac T1-weighted images output from the image synthesis subnetwork and the original image. Therefore, the optimization of this network can be parameterized as Equation (1).

[0027]

[0028] Where T1 and M0 are the parameters of the exponential signal relaxation model, N is the number of all cardiac T1-weighted images within the same slice, and Φ i For the i-th deformation field, I i For the i-th original T1-weighted image of the heart, t i is the inversion time of the i-th T1-weighted cardiac image.

[0029] 1) Construction of image synthesis subnetwork

[0030] like Figure 2 As shown in branch a, the image synthesis sub-network architecture is similar to UNet, containing encoder and decoder parts and skip connections between them. To reduce model complexity, the third encoder-decoder layer and the bottleneck layer can be replaced with the Obelisk (Heinrich MP, Oktay O, Bouteldja N. OBELISK-Net: Fewer layers to solve 3D multi-organ segmentation with sparse deformable convolutions[J]. Medical imageanalysis,2019,54:1-9.) module, which improves model training and testing efficiency without affecting network performance. The input of the image synthesis sub-network is all N T1-weighted images of the heart within the same slice level (I i＝0,...,N-1The output is the predicted signal attenuation model parameters T1 and M0. After obtaining the cardiac T1 map and M0 map, the exponential signal attenuation model is used. and the inversion time t of each T1-weighted cardiac image i To synthesize the corresponding "non-relative motion" T1-weighted cardiac image.

[0031] 2) Construction of deformable registration subnetwork

[0032] like Figure 2 As shown in branch b, the deformable registration subnetwork architecture is similar to UNet, containing encoder and decoder parts and skip connections between them. Similarly, to reduce model complexity, the Obelisk module replaces the third encoder-decoder layer and the bottleneck layer in the original UNet structure. The input to the deformable registration subnetwork is also all N T1-weighted cardiac images (I1, I2, I3, I4) within the same slice level. i＝0,...,N-1 After fully extracting features, the network predicts the corresponding static deformation velocity fields for N T1-weighted images. Then, N differential homeomorphic integration layers are used to solve for the differentiable and invertible deformation fields from the N deformation velocity fields, achieving differential homeomorphic registration. The relationship between the deformation field and the static velocity field is expressed as:

[0033]

[0034] Where t is time, and the final deformation field Φ (t) It is obtained by integrating the static velocity field with respect to time [0,t]. In specific implementations, the integration of the velocity field is achieved using a scaling square layer.

[0035] Subsequently, the solved deformation field and the original image are input into the Spatial Transform Network (STN) for deformation operations (including resampling and interpolation) to obtain N registration result images.

[0036] Step S130: Train the motion correction network with the goal of minimizing the overall loss function.

[0037] The optimization objective of this invention is to minimize the difference between the original image and the "motionless" T1-weighted cardiac image computed using an image synthesis subnetwork. The image similarity loss Lsim can be expressed as:

[0038]

[0039] Where N is the number of T1-weighted cardiac images within the same slice slice, MSE is the mean squared error, M0 and T1 are the signal attenuation model parameters predicted by the image synthesis subnetwork, and t i Let I be the inversion time of the i-th T1-weighted cardiac image. i For the i-th T1-weighted image of the heart, Φi is the i-th deformation field predicted by the deformable registration network. This indicates a deformation operation performed using STN. This represents the image showing the registration result.

[0040] In addition, the L2 norm is used to apply regularization constraints to the deformation field. Its loss function can be calculated as the sum of the squares of the L2 norms of the gradients of the N deformation fields at point p, see Equation (4).

[0041]

[0042] In summary, in one embodiment, the total network loss function of the motion correction network can be expressed as:

[0043] L total =L sim +λ·L smooth (5)

[0044] After the motion correction network is trained, optimized network parameters, such as weights and biases, can be obtained. The trained motion correction network can then be used for registration of acquired cardiac T1 images to improve image quality and provide clinical guidance. The overall process of this invention is described in [link to invention]. Figure 3 As shown.

[0045] In summary, compared with the prior art, the present invention has the following advantages:

[0046] 1) This invention adds a velocity integral layer to the cardiac T1 quantitative registration network to solve for differentiable and reversible deformation fields, which can effectively maintain the topology of the image before and after registration deformation, so as to obtain a real, reasonable and accurate registration result image without introducing erroneous deformation.

[0047] 2) The motion correction network of the present invention uses the Obelisk module to replace some encoder-decoder layers in UNet, which reduces the number of network parameters and effectively reduces the time cost of network training and image registration.

[0048] 3) Experimental verification shows that the present invention improves the image registration effect and at the same time increases the efficiency of image registration.

[0049] This invention can be a system, method, and / or computer program product. A computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for causing a processor to implement various aspects of the invention.

[0050] Computer-readable storage media can be tangible devices capable of holding and storing instructions for use by an instruction execution device. Computer-readable storage media can be, for example, but not limited to, electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, mechanical encoding devices, such as punch cards or recessed protrusions storing instructions thereon, and any suitable combination thereof. The computer-readable storage media used herein are not to be construed as transient signals themselves, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through waveguides or other transmission media (e.g., light pulses through fiber optic cables), or electrical signals transmitted through wires.

[0051] The computer-readable program instructions described herein can be downloaded from computer-readable storage media to various computing / processing devices, or downloaded via a network, such as the Internet, local area network, wide area network, and / or wireless network, to an external computer or external storage device. The network may include copper transmission cables, fiber optic transmission, wireless transmission, routers, firewalls, switches, gateway computers, and / or edge servers. A network adapter card or network interface in each computing / processing device receives the computer-readable program instructions from the network and forwards them to the computer-readable storage media in the respective computing / processing device.

[0052] The computer program instructions used to perform the operations of this invention may be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, state setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Smalltalk, C++, Python, etc., and conventional procedural programming languages ​​such as "C" or similar languages. The computer-readable program instructions may be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving a remote computer, the remote computer may be connected to the user's computer via any type of network—including a local area network (LAN) or a wide area network (WAN)—or may be connected to an external computer (e.g., via the Internet using an Internet service provider). In some embodiments, electronic circuitry, such as programmable logic circuitry, field-programmable gate arrays (FPGAs), or programmable logic arrays (PLAs), is personalized by utilizing state information from the computer-readable program instructions. This electronic circuitry can execute the computer-readable program instructions to implement various aspects of the invention.

[0053] Various aspects of the present invention are described herein with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It should be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer-readable program instructions.

[0054] These computer-readable program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus to produce a machine such that, when executed by the processor of the computer or other programmable data processing apparatus, they create means for implementing the functions / actions specified in one or more blocks of the flowchart and / or block diagram. These computer-readable program instructions can also be stored in a computer-readable storage medium that causes a computer, programmable data processing apparatus, and / or other device to operate in a particular manner; thus, the computer-readable medium storing the instructions comprises an article of manufacture that includes instructions for implementing aspects of the functions / actions specified in one or more blocks of the flowchart and / or block diagram.

[0055] Computer-readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the instructions executed on the computer, other programmable data processing apparatus, or other device to perform the functions / actions specified in one or more boxes of a flowchart and / or block diagram.

[0056] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of an instruction containing one or more executable instructions for implementing a specified logical function. In some alternative implementations, the functions marked in the blocks may occur in a different order than those marked in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, can be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions. It will be known to those skilled in the art that implementation in hardware, implementation in software, and implementation using a combination of software and hardware are equivalent.

[0057] The various embodiments of the present invention have been described above. These descriptions are exemplary and not exhaustive, and are not limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or technical improvements to the embodiments in the market, or to enable others skilled in the art to understand the embodiments disclosed herein. The scope of the invention is defined by the appended claims.

Claims

1. A method for quantitative differential homeomorphic registration of cardiac T1 based on differential models and physical constraints, comprising the following steps: Obtain the T1-weighted image of the heart of the target object; The T1-weighted cardiac image is input into a trained motion correction network to obtain a corrected image. The motion correction network includes an image synthesis subnetwork and a deformable registration subnetwork. The image synthesis subnetwork and the deformable registration subnetwork are both built on the UNet architecture. The image synthesis subnetwork is used to fit a heart T1 map and a heart M0 map using the registered heart T1 weighted image output from the deformable registration subnetwork. Then, the fitted heart T1 map, heart M0 map and signal attenuation model are used to synthesize a heart T1 weighted image without relative motion. The deformable registration subnetwork is used to perform intensity-based registration using the output image of the image synthesis subnetwork and the original image.

2. The method according to claim 1, characterized in that, The motion correction network is trained with the goal of minimizing the following total loss function: L total =L sim +λ·L smooth Among them, L total L represents the total loss value. sim L represents the image similarity loss function. smooth Let λ represent the deformation field smoothness loss function, where λ is a set weighting coefficient.

3. The method according to claim 2, characterized in that, The image similarity loss function is set as follows: Where N is the number of T1-weighted cardiac images within the same slice slice, MSE is the mean squared error, M0 and T1 are the signal attenuation model parameters predicted by the image synthesis subnetwork, and t i Let I be the inversion time of the i-th T1-weighted cardiac image. i For the i-th T1-weighted image of the heart, Φ i Let i be the deformation field predicted by the deformable registration network. This indicates a deformation operation performed using a spatial transformation network. This represents the image showing the registration result.

4. The method according to claim 2, characterized in that, The deformation field smoothness loss function is set as follows: in, This represents the gradient of the deformation field at point p. Let represent the sum of squares of the L2 norms of the gradients of the N deformation fields at point p.

5. The method according to claim 1, characterized in that, During the training of the motion correction network, the optimization formula for the deformable registration subnetwork is expressed as: Where T1 and Mo are the parameters of the exponential signal relaxation model, N is the number of all cardiac T1-weighted images within the same slice, and Φ i For the i-th deformation field, I i For the i-th original T1-weighted image of the heart, t i is the inversion time of the i-th T1-weighted cardiac image.

6. The method according to claim 1, characterized in that, The image synthesis subnetwork uses the Obelisk module to replace the third-layer encoder-decoder and bottleneck layer in the original UNet structure.

7. The method according to claim 1, characterized in that, The deformable registration subnetwork uses the Obelisk module to replace the third-layer encoder-decoder and bottleneck layer in the original UNet structure.

8. The method according to claim 1, characterized in that, The input to the deformable registration subnetwork is all N T1-weighted cardiac images within the same slice level. After feature extraction, it predicts the corresponding deformable static velocity fields of the N T1-weighted images. Then, it uses N differential homeomorphic integration layers to solve for the differentiable and invertible deformation fields from the N deformation velocity fields to achieve differential homeomorphic registration. The relationship between the deformation field and the static velocity field is expressed as: Where t is time, and Φ is the deformation field. (t) It is obtained by integrating the static velocity field over time [0,t].

9. A computer-readable storage medium having a computer program stored thereon, wherein, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 8.

10. A computer device comprising a memory and a processor, wherein a computer program capable of running on the processor is stored in the memory, characterized in that, When the processor executes the computer program, it implements the steps of the method according to any one of claims 1 to 8.

Images