Method and system for myocardial strain measurement based on group registration algorithm

Abstract

The application provides a kind of magnetic resonance image myocardial strain determination method and system based on group registration algorithm, comprising: step S1: left ventricular myocardium is segmented from diastolic end image of cardiac magnetic resonance movie image, and the image sequence is preprocessed according to the segmentation result, and the direction field is calculated;Step S2: using group registration algorithm based on global or local low rank difference measure to register image sequence, and estimate displacement field;Step S3: left ventricular myocardial strain related parameters are calculated using direction field and displacement field.The application uses group registration based method to eliminate cumulative error, so that myocardial strain can be more accurately determined, especially the strain in the later diastolic phase.The global or local low rank difference used in the method can more accurately measure the motion change between multiple movie images, and is more robust to artifacts.In addition, the joint use of spatial regularization term and cycle time regularization term further improves the accuracy of the method.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of medical image processing, in particular to a myocardial strain measurement method and system based on a group registration algorithm, and more particularly to a cardiac magnetic resonance cine image myocardial strain measurement method and system based on a group registration algorithm. BACKGROUND

[0002] Myocardial strain refers to the percentage change in length of local myocardium in a specified direction relative to the stress-free state (end diastole). It can be used to quantitatively describe the diastolic and systolic functions of the left ventricle and has the potential to identify subclinical left ventricular dysfunction.

[0003] Cardiac magnetic resonance cine imaging is currently the gold standard for non-invasive assessment of cardiac structure and function. It can reflect the anatomical structure at each stage of the cardiac cycle and each layer of the heart, and has the advantages of high resolution, good soft tissue contrast, and no ionizing radiation, so it is widely used in the measurement of clinical indicators such as stroke volume, ejection fraction, and myocardial strain. Cardiac cine imaging is usually based on two imaging sequences: balanced Steady State Free Precession (bSSFP) sequence and Spoiled Gradient Echo sequence. Cine images are usually a sequence of two-dimensional images, containing the entire cardiac cycle in time and the entire heart through multiple layers in space. Cardiac magnetic resonance cine imaging has different imaging angles, including short-axis view, 2-chamber long-axis view (perpendicular long-axis), 4-chamber long-axis view (horizontal long-axis), 3-chamber long-axis view, etc.

[0004] Cardiovascular Magnetic Resonance-Feature Tracking (CMR-FT) is a class of methods for measuring myocardial strain from conventional cine images of cardiac magnetic resonance. CMR-FT measures circumferential strain and radial strain using short-axis images, and longitudinal strain and radial strain using long-axis images (usually 2-chamber). CMR-FT is fast, non-invasive, and does not require the use of any contrast agent. Compared with tagging-Magnetic Resonance Imaging (tagging-MRI), CMR-FT does not require additional time-consuming sequences, and can be used for retrospective analysis of patients who have not undergone strain imaging. Currently, most CMR-FT methods are based on optical flow or sequentially pairwise registration. Optical flow-based methods are fast, but highly dependent on local structures in the images, such as speckles, which are rarely seen in the myocardial region of cine images. In contrast, pairwise registration-based methods enhance the robustness to local artifacts by exploiting contextual information. However, since its motion estimation is also frame-by-frame, the tracking error accumulates throughout the cardiac cycle, resulting in inaccurate myocardial strain measurements, especially in the diastolic phase, which is referred to as the drift effect.

[0005] Patent document CN103340628B discloses a heart real-time movie imaging image processing method and system. The method comprises: setting a reference point from a selected left ventricular center in a heart real-time movie imaging image, positioning a right ventricle according to the reference point; selecting a region from the heart real-time movie imaging image, and segmenting the region to obtain boundary points; fitting the boundary points to obtain a left ventricular center and a left ventricular area; recording a fitted left ventricular center changing with time to obtain a respiratory motion signal, and determining an end-expiratory period; calculating a cross-correlation coefficient between all frame images at the end-expiratory period, selecting two frame images with the smallest cross-correlation coefficient, and taking a left ventricular area obtained by fitting in the two frame images as a heart diastolic end period image, and taking a left ventricular area obtained by fitting in the two frame images as a heart systolic end period image. However, the invention does not use a group registration algorithm based on global or local low-rank difference measurement to register image sequences to estimate a displacement field, and uses a direction field and the displacement field to calculate left ventricular myocardial strain related parameters. SUMMARY

[0006] Aiming at the defects in the prior art, the present application aims to provide a magnetic resonance image myocardial strain measurement method and system based on a group registration algorithm.

[0007] According to the present application, a magnetic resonance image myocardial strain measurement method based on a group registration algorithm is provided, comprising:

[0008] Step S1: segmenting the left ventricular myocardium from the end-diastolic image of the cardiac magnetic resonance movie image, and pre-processing the image sequence and calculating the direction field according to the segmentation result;

[0009] Step S2: registering the image sequence using a group registration algorithm based on global or local low-rank difference measurement, and estimating the displacement field;

[0010] Step S3: calculating the left ventricular myocardial strain related parameters using the direction field and the displacement field.

[0011] Preferably, in the step S1:

[0012] The magnetic resonance movie image comprises two-dimensional cardiac movie images and three-dimensional cardiac movie images; and the acquisition sequence of the magnetic resonance movie image comprises a balanced steady-state free precession movie imaging sequence and a spoiled gradient echo movie imaging sequence.

[0013] Preferably, in the step S1:

[0014] The left ventricular myocardium is segmented from the end-diastolic image; the image is pre-processed according to the segmentation result, including cropping and gray scale normalization; the unit direction vector at each myocardial pixel is calculated according to the segmentation result to form a direction field; for the short-axis image, the circumferential direction field and the radial direction field are calculated; for the long-axis image, the longitudinal direction field and the radial direction field are calculated.

[0015] Preferably, in the step S2:

[0016] The displacement field of the image sequence is estimated using a group registration algorithm based on global or local low-rank difference measurement, and the group registration algorithm is composed of a motion model, an objective function and an optimization strategy:

[0017] The motion model is used to describe the motion mode of each organ in the image, and determines the flexibility and smoothness of the generated displacement field, including a non-parametric model or a parametric model;

[0018] The objective function is composed of a difference measurement and a regularization term, wherein the difference measurement describes the registration degree of the image under different displacement fields, and the regularization term restricts the range of variation of the displacement field itself;

[0019] Optimization strategy: find the optimal displacement field by minimizing the objective function; the optimization strategy of the group registration algorithm is based on algorithms including gradient descent method, conjugate gradient method and LBFGS algorithm, and adopts a multi-resolution optimization framework from coarse to fine.

[0020] Preferably, the objective function comprises:

[0021] The adopted difference measure is based on global low-rankness or local low-rankness of the registered movie image sequence.

[0022] The difference measure based on global low-rankness is defined as follows: flatten the two-dimensional images in a movie image sequence into column vectors and arrange them into a Casorati matrix F in the time direction, and take the nuclear norm of the matrix as the difference measure of the image sequence:

[0023]

[0024] Where, ||F|| denotes the nuclear norm of the F matrix, which is defined as the sum of all non-zero singular values of the F matrix, σ * k represents the kth non-zero singular value of F, and rank(F) represents the rank of F. k

[0025] The difference measure based on local low-rankness is defined as follows: extract multiple image blocks from a movie image sequence at equal intervals along the spatial direction, and construct Casorati matrices based on the spatio-temporal signals carried by each image block, calculate the nuclear norm of each Casorati matrix, and then take the sum of the nuclear norms of all Casorati matrices as the difference measure of the image sequence:

[0026]

[0027] Where, B m m represents the Casorati matrix composed of the mth image block, and M represents the number of image blocks.

[0028] Preferably, the objective function comprises:

[0029] The adopted regularization term includes a spatial regularization term and a temporal regularization term; the spatial regularization term is a diffusion smoothing term based on first-order difference or a thin-plate bending energy based on second-order difference; the temporal regularization term is the square sum of first-order or second-order difference of the displacement field along the time dimension; the difference in the temporal regularization term can be a cyclic difference, so that the deformation of the last image and the deformation of the first image remain continuous; the strengths of the temporal and spatial regularization are controlled by two regularization coefficients respectively.

[0030] Preferably, the optimization strategy comprises:

[0031] ​The group registration algorithm iteratively solves a minimization problem to obtain an optimal displacement field, and the algorithm includes gradient descent method, conjugate gradient method and LBFGS algorithm; the group registration algorithm adopts a multi-resolution optimization framework from coarse to fine, first registers the images with lower resolution, and then initializes the registration of the images with higher resolution by using the obtained control point grid; under this framework, the group registration algorithm generates global deformation at low resolution to align the target to a preset degree, and generates local deformation at high resolution to make adjustment.

[0032] Preferably, in the step S3:

[0033] According to the obtained displacement field, the position of each pixel point in the center myocardial region of the end-diastolic image in each frame of image is tracked; based on the group registration method, the pixel points to be tracked are first mapped to an average coordinate system, and then mapped to the remaining frame images;

[0034] In combination with the calculated direction field, the strain of each pixel in the myocardial region is calculated according to the definition of myocardial strain, and a myocardial strain map is obtained; for the short-axis image, a circumferential strain map and a radial strain map are calculated; for the long-axis image, a longitudinal strain map and a radial strain map are calculated; the strain in the whole left ventricular myocardial region is averaged to obtain a curve of global myocardial strain changing with time; the global myocardial strain is differentiated to obtain a curve of myocardial strain rate changing with time.

[0035] According to the group registration algorithm-based magnetic resonance image myocardial strain measurement system provided by the application, comprising:

[0036] Module M1: The left ventricular myocardium is segmented from the end-diastolic image of the cardiac magnetic resonance movie image, and the image sequence is preprocessed and the direction field is calculated according to the segmentation result;

[0037] Module M2: The image sequence is registered by using the group registration algorithm based on global or local low-rank difference measurement, and the displacement field is estimated;

[0038] Module M3: The left ventricular myocardial strain related parameters are calculated by using the direction field and the displacement field.

[0039] Preferably, in the module M1:

[0040] The magnetic resonance movie image includes two-dimensional cardiac movie image and three-dimensional cardiac movie image; the acquisition sequence of the magnetic resonance movie image includes balanced steady-state free precession movie imaging sequence and spoiled gradient echo movie imaging sequence;

[0041] Segmenting left ventricular myocardium from end-diastolic images; according to the segmentation results, pre-processing the images, including cropping and gray scale normalization; according to the segmentation results, calculating the unit direction vector at each myocardial pixel to form a direction field; for short-axis images, calculating the circumferential direction field and the radial direction field; for long-axis images, calculating the longitudinal direction field and the radial direction field;

[0042] In the module M2:

[0043] Estimating the displacement field of the image sequence using a group registration algorithm based on global or local low-rank difference measure, the group registration algorithm consisting of a motion model, an objective function, and an optimization strategy:

[0044] Motion model: used to describe the motion mode of each organ in the image, determining the flexibility and smoothness of the generated displacement field, including non-parametric model or parametric model;

[0045] Objective function: composed of difference measure and regularization term, where the difference measure describes the registration degree of the image under different displacement fields, and the regularization term constrains the range of variation of the displacement field itself;

[0046] Optimization strategy: find the optimal displacement field by minimizing the objective function; the optimization strategy of the group registration algorithm is based on algorithms including gradient descent method, conjugate gradient method and LBFGS algorithm, and adopts a multi-resolution optimization framework from coarse to fine;

[0047] The objective function includes:

[0048] The difference measure adopted is based on the global low-rank or local low-rank of the registered movie image sequence.

[0049] The difference measure based on global low-rank is defined as follows: flatten the two-dimensional images in a movie image sequence into column vectors, and arrange them into a Casorati matrix F in the time direction, and take the kernel norm of the matrix as the difference measure of the image sequence:

[0050]

[0051] Where, ||F|| * represents the kernel norm of the F matrix, which is defined as the sum of all non-zero singular values of the F matrix, σ k is the kth non-zero singular value of F, and rank(F) is the rank of F;

[0052] The difference measurement based on local low-rankness is defined as follows: multiple image blocks are extracted from a movie image sequence along the spatial direction at equal intervals, and Casorati matrices are constructed based on the space-time signals carried by the image blocks respectively, the nuclear norm of each Casorati matrix is calculated, and the sum of the nuclear norms of all Casorati matrices is taken as the difference measurement of the image sequence:

[0053]

[0054] Wherein, B m is a Casorati matrix composed of the mth image block, and M is the number of image blocks;

[0055] The objective function comprises:

[0056] The regularization term adopted comprises a spatial regularization term and a temporal regularization term; the spatial regularization term is a diffusion smoothing term based on first-order difference or a thin-plate bending energy based on second-order difference; the temporal regularization term is the square sum of first-order or second-order difference of the displacement field along the time dimension; the difference in the temporal regularization term can be a cyclic difference, so that the deformation of the last image and the deformation of the first image remain continuous; the strengths of the temporal and spatial regularization are controlled by two regularization coefficients respectively;

[0057] The optimization strategy comprises:

[0058] The group registration algorithm iteratively solves a minimization problem to obtain an optimal displacement field, and the algorithm used comprises gradient descent method, conjugate gradient method and LBFGS algorithm; the group registration algorithm adopts a multi-resolution optimization framework from coarse to fine, and first registers lower resolution images, and then initializes higher resolution registration by using the obtained control point grid; under this framework, the group registration algorithm first generates a global deformation at low resolution to align the target to a preset degree, and generates a local deformation at high resolution for adjustment;

[0059] In the module M3:

[0060] According to the obtained displacement field, the position of each pixel point in the center muscle region of the end-diastolic image in each of the other frames is tracked; based on the group registration method, the pixel points to be tracked are first mapped to an average coordinate system, and then mapped to the remaining frames;

[0061] In combination with the calculated direction field, the strain of each pixel in the myocardial region is calculated according to the definition of myocardial strain, and a myocardial strain map is obtained; for short-axis images, a circumferential strain map and a radial strain map are calculated; for long-axis images, a longitudinal strain map and a radial strain map are calculated; the strain in the entire left ventricular myocardial region is averaged to obtain a curve of global myocardial strain change over time; the global myocardial strain is differentiated to obtain a curve of myocardial strain rate change over time.

[0062] Compared with the prior art, the present application has the following beneficial effects:

[0063] 1. The present application uses group registration to estimate myocardial motion, thereby eliminating cumulative errors and more accurately determining the diastolic myocardial strain in the later timing;

[0064] 2. The present application uses global or local low-rank difference to more accurately measure the motion change between multiple movie images and is more robust to artifacts;

[0065] 3. The present application further improves the accuracy of the method by using the combination of spatial regularization and cyclic time regularization. BRIEF DESCRIPTION OF DRAWINGS

[0066] Other features, objects and advantages of the present application will become more apparent from the following detailed description of non-limiting embodiments with reference to the attached drawings:

[0067] Figure 1 Flow chart of a myocardial strain measurement method based on a group registration algorithm of an embodiment of the present application;

[0068] Figure 2 Schematic diagram of the circumferential and radial myocardial strain curve measurement results of an embodiment of the present application and the prior art;

[0069] Figure 3 Schematic diagram of the end-systolic circumferential and radial myocardial strain map measurement results of an embodiment of the present application and the prior art. DETAILED DESCRIPTION

[0070] The present application will be described in detail below with reference to specific embodiments. The following embodiments will help those skilled in the art to further understand the present application, but do not limit the present application in any form. It should be noted that those skilled in the art can make several changes and improvements without departing from the concept of the present application. These all belong to the protection scope of the present application.

[0071] Example 1:

[0072] In view of the limitations of the prior art, the present application proposes a myocardial strain measurement method based on a group registration algorithm of a cardiac magnetic resonance movie image. The method uses group registration to estimate myocardial motion, thereby eliminating the drift effect commonly seen in optical flow or pair registration methods. At the same time, the global or local low-rank difference measurement and the joint use of spatial and cyclic time regularization further improve the accuracy of the method.

[0073] The application discloses a cardiac magnetic resonance cine image myocardial strain measurement method based on a group registration algorithm. The method comprises the following steps: (S1) segmenting left ventricular myocardium from an end-diastolic image of a cine image, and pre-processing and calculating a direction field according to a segmentation result; (S2) registering the image sequence by using a group registration algorithm based on global or local low-rank difference measurement, and estimating a displacement field; (S3) calculating left ventricular myocardial strain related parameters by using the direction field and the displacement field, including but not limited to a myocardial strain map, a myocardial strain curve and a myocardial strain rate curve. Compared with existing methods based on optical flow or sequential pair registration, the group registration based method eliminates cumulative errors, so that the myocardial strain can be more accurately measured, especially the diastolic strain in the later time sequence. The global or local low-rank difference used in the method can more accurately measure the motion change between multiple cine images, and is more robust to artifacts. In addition, the joint use of a spatial regularization term and a cycle time regularization term further improves the accuracy of the method.

[0074] According to the application, a magnetic resonance image myocardial strain measurement method based on a group registration algorithm is provided, which comprises the following steps: Figures 1-3 , including:

[0075] Step S1: segmenting left ventricular myocardium from an end-diastolic image of a cardiac magnetic resonance cine image, and pre-processing and calculating a direction field according to a segmentation result;

[0076] Specifically, in the step S1:

[0077] The magnetic resonance cine image comprises a two-dimensional cardiac cine image and a three-dimensional cardiac cine image; and the acquisition sequence of the magnetic resonance cine image comprises a bSSFP (balanced steady-state free precession) cine imaging sequence and a spoiled GRE (spoiled gradient echo) cine imaging sequence.

[0078] Specifically, in the step S1:

[0079] Segmenting left ventricular myocardium from the end-diastolic image; pre-processing the image according to the segmentation result, including cropping and gray scale normalization; calculating a unit direction vector at each myocardial pixel according to the segmentation result to form a direction field; for a short-axis image, calculating a circumferential direction field and a radial direction field; for a long-axis image, calculating a longitudinal direction field and a radial direction field.

[0080] Step S2: registering the image sequence by using a group registration algorithm based on global or local low-rank difference measurement, and estimating a displacement field;

[0081] Specifically, in the step S2:

[0082] The displacement field of the image sequence is estimated using a group registration algorithm based on a global or local low-rank dissimilarity measure, which consists of a motion model, an objective function and an optimization strategy:

[0083] The motion model is used to describe the motion manner of each organ in the image, and determines the flexibility and smoothness of the generated displacement field, including a non-parametric model or a parametric model;

[0084] The objective function is composed of a dissimilarity measure and a regularization term, wherein the dissimilarity measure describes the registration degree of the image under different displacement fields, and the regularization term restricts the range of variation of the displacement field itself;

[0085] The optimization strategy is to find the optimal displacement field by minimizing the objective function; the optimization strategy of the group registration algorithm is based on the gradient descent method, the conjugate gradient method and the LBFGS algorithm, and a multi-resolution optimization framework from coarse to fine is adopted.

[0086] Specifically, the objective function includes:

[0087] The adopted dissimilarity measure is based on the global low-rank or local low-rank of the registered movie image sequence.

[0088] The dissimilarity measure based on global low-rank (Globally Low-Rank, GLR) is defined as follows: the two-dimensional images in a movie image sequence are flattened into column vectors and arranged into a Casorati matrix F in the time direction, and the nuclear norm of the matrix is taken as the dissimilarity measure of the image sequence:

[0089]

[0090] Wherein, ||F|| * represents the nuclear norm of the F matrix, which is defined as the sum of all non-zero singular values of the F matrix, σ k is the kth non-zero singular value of F, and rank(F) is the rank of F;

[0091] The dissimilarity measure based on local low-rank is defined as follows: multiple image blocks are extracted at equal intervals along the spatial direction from a movie image sequence, and Casorati matrices are constructed based on the space-time signals carried by each image block, the nuclear norm of each Casorati matrix is calculated, and the sum of the nuclear norms of all Casorati matrices is taken as the dissimilarity measure of the image sequence:

[0092]

[0093] Wherein, B m is the Casorati matrix composed of the mth image block, and M is the number of image blocks.

[0094] The objective function comprises:

[0095] The regular term adopted comprises a spatial regular term and a temporal regular term; the spatial regular term is a diffusion smoother based on first-order difference or a bending energy of a thin plate based on second-order difference; the temporal regular term is a square sum of first-order or second-order difference of displacement field along time dimension; optionally, the difference in the temporal regular term can be a cyclic difference, so that the deformation of the last frame image and the deformation of the first frame image remain continuous. The strength of the temporal and spatial regularization is controlled by two regularization coefficients respectively.

[0096] Specifically, the optimization strategy comprises:

[0097] The group registration algorithm iteratively solves a minimization problem to obtain an optimal displacement field, and the algorithm used includes gradient descent method, conjugate gradient method and LBFGS algorithm; the group registration algorithm adopts a multi-resolution optimization framework from coarse to fine, first registers lower resolution images, and then initializes the registration of higher resolution images using the obtained control point grid; under this framework, the group registration algorithm first generates global deformation at low resolution to align the target to a preset degree, and then generates local deformation at high resolution for adjustment.

[0098] Step S3: calculating left ventricular myocardial strain related parameters using the direction field and the displacement field.

[0099] Specifically, in the step S3:

[0100] According to the obtained displacement field, the position of each pixel point in the myocardial region of the end-diastolic image in each of the other frames of images is tracked; based on the group registration method, the pixel points to be tracked are first mapped to an average coordinate system, and then mapped to the remaining frames of images;

[0101] In combination with the calculated direction field, the strain at each pixel in the myocardial region is calculated according to the definition of myocardial strain, and a myocardial strain map is obtained; for short-axis images, a circumferential strain map and a radial strain map are calculated; for long-axis images, a longitudinal strain map and a radial strain map are calculated; the strain in the entire left ventricular myocardial region is averaged to obtain a curve of global myocardial strain change over time; the global myocardial strain is differentiated to obtain a curve of myocardial strain rate change over time.

[0102] Embodiment 2:

[0103] Embodiment 2 is a preferred example of Embodiment 1, which more specifically illustrates the present application.

[0104] The application also provides a magnetic resonance image myocardial strain measurement system based on a group registration algorithm, which can be realized by performing the process steps of the magnetic resonance image myocardial strain measurement method based on the group registration algorithm, i.e., the magnetic resonance image myocardial strain measurement method based on the group registration algorithm can be understood by those skilled in the art as a preferred embodiment of the magnetic resonance image myocardial strain measurement system based on the group registration algorithm.

[0105] According to the application, a magnetic resonance image myocardial strain measurement system based on a group registration algorithm is provided, comprising:

[0106] Module M1: segmenting left ventricular myocardium from an end-diastolic image of a cardiac magnetic resonance cine image, and pre-processing the image sequence and calculating a direction field according to the segmentation result;

[0107] Specifically, in the module M1:

[0108] The magnetic resonance cine image comprises a two-dimensional cardiac cine image and a three-dimensional cardiac cine image; and the acquisition sequence of the magnetic resonance cine image comprises a balanced steady-state free precession cine imaging sequence and a spoiled gradient echo cine imaging sequence.

[0109] Segmenting left ventricular myocardium from an end-diastolic image; pre-processing the image according to the segmentation result, including cropping and gray scale normalization; calculating a unit direction vector at each myocardial pixel according to the segmentation result to form a direction field; for a short-axis image, calculating a circumferential direction field and a radial direction field; for a long-axis image, calculating a longitudinal direction field and a radial direction field;

[0110] Module M2: registering the image sequence using a group registration algorithm based on a global or local low-rank difference measure to estimate a displacement field;

[0111] In the module M2:

[0112] The group registration algorithm based on a global or local low-rank difference measure is used to estimate the displacement field of the image sequence, and the group registration algorithm is composed of a motion model, an objective function, and an optimization strategy:

[0113] Motion model: used to describe the motion mode of each organ in the image, and determines the flexibility and smoothness of the generated displacement field, including a non-parametric model or a parametric model;

[0114] Objective function: composed of a difference measure and a regularization term, wherein the difference measure describes the registration degree of the image under different displacement fields, and the regularization term restricts the range of variation of the displacement field itself;

[0115] Optimization strategy: find the optimal displacement field by minimizing the objective function; the optimization strategy of the group registration algorithm is based on algorithms including gradient descent method, conjugate gradient method and LBFGS algorithm, and adopts a multi-resolution optimization framework from coarse to fine;

[0116] The objective function includes:

[0117] The adopted difference measure is based on global low-rankness or local low-rankness of the registered movie image sequence.

[0118] The difference measure based on global low-rankness is defined as follows: flatten the two-dimensional images in a movie image sequence into column vectors and arrange them into a Casorati matrix F in the time direction, and take the nuclear norm of the matrix as the difference measure of the image sequence:

[0119]

[0120] Where, ||F|| * represents the nuclear norm of the F matrix, which is defined as the sum of all non-zero singular values of the F matrix, σ k is the kth non-zero singular value of F, and rank(F) is the rank of F;

[0121] The difference measure based on local low-rankness is defined as follows: extract multiple image blocks from a movie image sequence at equal intervals along the spatial direction, and construct Casorati matrices based on the space-time signals carried by each image block, calculate the nuclear norm of each Casorati matrix, and then take the sum of the nuclear norms of all Casorati matrices as the difference measure of the image sequence:

[0122]

[0123] Where, B m is the Casorati matrix composed of the mth image block, and M is the number of image blocks;

[0124] The objective function includes:

[0125] The adopted regularization term includes a spatial regularization term and a temporal regularization term; the spatial regularization term is a diffusion smoothing term based on first-order difference or a thin-plate bending energy based on second-order difference; the temporal regularization term is the square sum of first-order or second-order difference of the displacement field along the time dimension; the difference in the temporal regularization term can be a cyclic difference, so that the deformation of the last image and the deformation of the first image remain continuous; the strengths of the temporal and spatial regularization are controlled by two regularization coefficients respectively;

[0126] The optimization strategy includes:

[0127] The group registration algorithm solves the minimum problem iteratively to obtain an optimal displacement field, and the algorithm includes gradient descent method, conjugate gradient method and LBFGS algorithm; the group registration algorithm adopts a multi-resolution optimization framework from coarse to fine, and first registers the images with lower resolution, and then initializes the registration of the images with higher resolution by using the obtained control point grid; under this framework, the group registration algorithm generates global deformation at low resolution to align the target to a preset degree, and generates local deformation at high resolution to adjust;

[0128] Module M3: calculating left ventricular myocardial strain related parameters by using the direction field and the displacement field.

[0129] In the module M3:

[0130] According to the obtained displacement field, the position of each pixel point in the myocardial region of the end-diastolic image in each frame of image is tracked; the pixel point to be tracked is mapped into an average coordinate system and then mapped into the remaining frames of image based on the group registration method;

[0131] In combination with the calculated direction field, the strain of each pixel in the myocardial region is calculated according to the definition of myocardial strain, and a myocardial strain map is obtained; for the short-axis image, a circumferential strain map and a radial strain map are calculated; for the long-axis image, a longitudinal strain map and a radial strain map are calculated; the strain in the whole left ventricular myocardial region is averaged to obtain a curve of global myocardial strain changing with time; the global myocardial strain is differentiated to obtain a curve of myocardial strain rate changing with time.

[0132] Embodiment 3:

[0133] Embodiment 3 is a preferred example of Embodiment 1, and is used to more specifically illustrate the present application.

[0134] In this embodiment, the FFD model is used as the motion model of the group registration algorithm; the global or local low-rank difference is used as the difference measurement between the cardiac magnetic resonance movie images; the thin-plate bending energy is used as the spatial regularization term; the second-order cyclic time regularization term is used; and the gradient projection algorithm is used to minimize the objective function to obtain the optimal displacement field.

[0135] As Figure 1As shown, the application discloses a cardiac magnetic resonance cine image myocardial strain measurement method based on a group registration algorithm. An embodiment of the application is introduced below, but is not as a limitation of the application. The embodiment adopts a free-form deformation model (FFD) as a motion model of the group registration algorithm; adopts local low-rank difference as a difference measurement between cardiac magnetic resonance cine images; adopts a bending energy of a thin plate as a spatial regularization term; adopts a second-order cyclic time regularization term; and adopts a gradient projection algorithm to minimize an objective function to obtain an optimal displacement field. Specifically, the method comprises the following steps: (S1) determining an end diastole image by automatic / semi-automatic software, and segmenting left ventricular myocardium from the end diastole image. According to the segmentation result, each frame of image is cropped to 128*128 with a left ventricular center as a center, and the gray scale of each frame of image is linearly normalized to [0, 1]. Then, according to the segmentation result, a unit direction vector at each myocardial pixel is calculated to form a direction field. For a short-axis image, a circumferential direction field and a radial direction field are calculated. The circumferential direction is parallel to an endocardium and an epicardium, and the radial direction is perpendicular to the circumferential direction. For a long-axis image, a longitudinal direction field and a radial direction field are calculated. The longitudinal direction is parallel to the endocardium and the epicardium, and the radial direction is perpendicular to the longitudinal direction.

[0136] (S2) using a group registration algorithm based on global or local low-rank difference measurement to estimate a displacement field of the image sequence. The group registration algorithm comprises three parts, i.e., a motion model, an objective function and an optimization strategy, which are described below.

[0137] (S21) motion model: the group registration algorithm adopts an FFD model as a motion model. The FFD model controls the displacement of pixels by moving control points and B-spline interpolation, so that the displacement field generated is flexible and smooth, and is suitable for describing local morphological changes of soft tissues such as myocardium. In addition, in order to register all images to an average coordinate system, it is required that the mean value of the displacement field in the time dimension is zero. For the FFD model, since there is a linear relationship between the displacement field and the displacement of the control point grid, the constraint can be converted into the mean value of the control point grid in the time dimension being zero.(S22) objective function: the objective function of the group registration algorithm is composed of a difference measurement and a regularization term, which are introduced below.

[0138] The difference measure adopted by the algorithm is based on the locally low-rankness (LLR) between the registered images. Since the cardiac cine images often contain irregular temporal signal variations such as blood flow artifacts, they do not necessarily obey the globally low-rankness (GLR) even after alignment. Therefore, the LLR difference measure is more robust to artifacts than the GLR difference measure. Specifically, a Casorati matrix is constructed using a number of image patches extracted from the images at equal intervals, and its nuclear norm is computed and summed up as the LLR difference measure between the images, i.e.

[0139]

[0140] where Bm is the Casorati matrix composed of the m-th image patch, and M is the number of image patches. m

[0141] To make the generated deformation consistent with the real situation of the tissue, the thin-plate-spline energy is adopted as the spatial regularization term to encourage the deformation to be spatially smooth. Similarly, to ensure the temporal smoothness of the deformation, a second-order temporal regularization term is also introduced. Since the actual images are discrete in both spatial and temporal dimensions, the partial derivatives in the regularization terms are all computed using difference approximations. For a periodic image sequence such as the cardiac cine images, the difference in the temporal regularization term can be further replaced by a cyclic difference so that the deformation of the last image also remains continuous with that of the first image. The strengths of the temporal and spatial regularizations are controlled by two regularization coefficients, respectively.

[0142] (S23) Optimization strategy: Combining the above difference measure, the temporal and spatial regularization terms, and the constraints on the control point mesh, the optimization problem to be solved by the group registration algorithm is obtained.

[0143] The gradient projection algorithm is adopted to iteratively solve the optimization problem to obtain the optimal displacement field. In the case where the initial displacements of the control points are taken as zeros, the optimization algorithm is equivalent to centering the gradient of the control point mesh in the temporal dimension before each update of the gradient descent algorithm. The group registration algorithm also adopts a coarse-to-fine multi-resolution optimization framework, i.e., registering the images at a lower resolution first, and then initializing the registration of the images at a higher resolution using the obtained control point mesh.

[0144] ​(S3) To measure the myocardial strain, first calculate the displacement field from the end-diastole image to other frames based on the displacement field obtained in S2. The groupwise registration based method will first map each pixel in the end-diastole image to an average coordinate system, and then map to other frames, thus eliminating the drift effect. Then, calculate the Green-Lagrange strain tensor based on the definition. Combine with the direction field calculated in S1, calculate the strain at each pixel in the myocardial region, and obtain the myocardial strain map. For short-axis images, calculate the circumferential strain map and radial strain map: for long-axis images, calculate the longitudinal strain map and radial strain map

[0145] Take the average of the strain in the entire left ventricular myocardial region to obtain the curve of global myocardial strain over time. Then, take the derivative (difference) of the global myocardial strain to obtain the curve of myocardial strain rate over time.

[0146] Figure 2 An embodiment of the present application and the end-systolic circumferential and radial myocardial strain maps measured by the existing method are shown. The existing method includes the Farneback optical flow based method and the pairwise registration based method used by the commercial software Segment Medviso. The arrows in the figure indicate the infarct myocardium located in the septal region. In all measured radial strain maps, the strain in the septal region is significantly lower, which reflects the decrease in the activity of the infarct myocardium. However, compared with the Farneback and pairwise methods, the radial strain map measured by the groupwise method more accurately describes the infarct region. In addition, the circumferential strain map measured by the groupwise method is also smoother. (Circumferential Strain: CS; Radial Strain: RS)

[0147] Figure 3An example of the present application is shown with the circumferential and radial myocardial strain curves determined by the existing methods. The existing methods include the Farneback optical flow based method and the pairwise registration based method used by the commercial software Segment Medviso. All methods successfully captured the reduction of myocardial strain in both MINF and DCM subjects. Since the cardiac magnetic resonance cine images are periodic, the myocardial strain in the last frame should be close to zero. It can be seen that the diastolic strains determined by the Farneback and pairwise methods have a significant bias as indicated by the arrows in the figure. In contrast, the two groupwise methods eliminate the accumulated error and determine more accurate diastolic strains. (Normal Subject, NOR; Myocardial Infarction, MINF; Dilated Cardiomyopathy, DCM; End-Diastole, ED; End-Systole, ES; Global Circumferential Strain, GCS; Global Radial Strain, GRS)

[0148] Those skilled in the art will appreciate that, in addition to implementing the system, apparatus and each module thereof provided by the present application in the form of pure computer readable program code, the system, apparatus and each module thereof provided by the present application can also be implemented in the form of logic gates, switches, application specific integrated circuits, programmable logic controllers and embedded microcontrollers, etc. by logically programming the method steps. Therefore, the system, apparatus and each module thereof provided by the present application can be considered as a hardware component, and the modules included therein for implementing various programs can also be considered as structures in the hardware component; the modules for implementing various functions can also be considered as both software programs for implementing methods and structures in the hardware component.

[0149] The specific embodiments of the present application are described above. It needs to be understood that the present application is not limited to the specific embodiments described above, and various changes or modifications can be made by those skilled in the art within the scope of the claims, which does not affect the essential content of the present application. The embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily without conflict.

Claims

1. A method of magnetic resonance image myocardial strain measurement based on a group registration algorithm, characterized by, Comprising: Step S1: segmenting left ventricular myocardium from end-diastolic image of cardiac magnetic resonance cine images, and pre-processing image sequence and calculating direction field according to segmentation result; Step S2: registering image sequence using group registration algorithm based on global or local low-rank difference measure, and estimating displacement field; Step S3: calculating left ventricular myocardial strain-related parameters using direction field and displacement field; In the step S2: Estimating displacement field of image sequence using group registration algorithm based on global or local low-rank difference measure, and group registration algorithm is composed of motion model, objective function and optimization strategy: Motion model: used for describing motion mode of each organ in image, and determining flexibility and smoothness of generated displacement field, including non-parametric model or parametric model; Objective function: composed of difference measure and regularization term, wherein difference measure describes registration degree of image under different displacement fields, and regularization term restricts range of displacement field itself that can be changed; Optimization strategy: finding optimal displacement field by minimizing objective function; optimization strategy of group registration algorithm is based on algorithms including gradient descent method, conjugate gradient method and LBFGS algorithm, and adopting multi-resolution optimization framework from coarse to fine; The objective function includes: Difference measure adopted is based on global low-rank or local low-rank of registered cine image sequence; The difference measure based on global low-rankness is defined as follows: flatten the two-dimensional images in a movie image sequence into column vectors and arrange them into a Casorati matrix in the temporal direction Take the nuclear norm of the matrix as the difference measure of the image sequence: wherein represent the nuclear norm of a matrix, defined as the sum of all non-zero singular values of a matrix, is the th non-zero singular value of is the rank of Difference measure based on local low-rank is defined as follows: multiple image blocks are extracted along spatial direction at equal intervals from a cine image sequence, and Casorati matrix is constructed based on space-time signal carried by each image block, core norm of each Casorati matrix is calculated, and sum of core norms of all Casorati matrices is taken as difference measure of the image sequence: wherein is a Casorati matrix composed of the first is the number of image blocks.

2. The method of claim 1, wherein the method further comprises: In the step S1: Magnetic resonance cine images include two-dimensional cardiac cine images and three-dimensional cardiac cine images; and acquisition sequence of magnetic resonance cine images includes balanced steady-state free precession cine imaging sequence and spoiled gradient echo cine imaging sequence.

3. The method of claim 1, wherein the method further comprises: determining a plurality of strain values for each of the plurality of cardiac phases; and determining a plurality of strain curves for each of the plurality of cardiac phases. In the step S1: Segmenting left ventricular myocardium from end-diastolic image; According to segmentation result, pre-processing image, including cropping and gray scale normalization; according to segmentation result, calculating unit direction vector at each myocardial pixel, forming direction field; for short-axis image, calculating circumferential direction field and radial direction field; for long-axis image, calculating longitudinal direction field and radial direction field.

4. The magnetic resonance image myocardial strain measurement method based on group registration algorithm according to claim 1, wherein: The objective function includes: Regularization term adopted includes spatial regularization term and temporal regularization term; spatial regularization term is diffusion smoothing term based on first-order difference or thin plate spline energy based on second-order difference; temporal regularization term is square sum of first-order or second-order difference of displacement field along time dimension; difference in temporal regularization term can be cyclic difference, so that deformation of last image and deformation of first image remain continuous; intensity of temporal and spatial regularization is controlled by two regularization coefficients respectively.

5. The magnetic resonance image myocardial strain measurement method based on group registration algorithm according to claim 1, wherein: The optimization strategy includes: The group registration algorithm iteratively solves a minimization problem to obtain an optimal displacement field, and the algorithm includes gradient descent, conjugate gradient and LBFGS algorithm; the group registration algorithm adopts a coarse-to-fine multi-resolution optimization framework, first registers the images with lower resolution, and then initializes the registration of the images with higher resolution by using the obtained control point grid; under this framework, the group registration algorithm first generates a global deformation at low resolution to align the target to a preset degree, and then generates a local deformation at high resolution to make adjustments.

6. The method of claim 1, wherein the method further comprises: determining a plurality of strain values for each of the plurality of cardiac phases; and determining a plurality of strain curves for each of the plurality of cardiac phases. In the step S3: According to the obtained displacement field, the position of each pixel point in the center myocardial region of the end-diastolic image in each frame of image is tracked; based on the group registration method, the pixel point to be tracked is first mapped into an average coordinate system, and then mapped into each frame of image; Combined with the calculated direction field, the strain at each pixel in the myocardial region is calculated according to the definition of myocardial strain, and a myocardial strain map is obtained; for the short-axis image, a circumferential strain map and a radial strain map are calculated; for the long-axis image, a longitudinal strain map and a radial strain map are calculated; the average of the strain in the whole left ventricular myocardial region is taken to obtain a curve of the global myocardial strain changing with time; the derivative of the global myocardial strain is taken to obtain a curve of the myocardial strain rate changing with time.

7. A magnetic resonance image myocardial strain measurement system based on a group registration algorithm, characterized by, Comprise: Module M1: segmenting the left ventricular myocardium from the end-diastolic image of the cardiac magnetic resonance movie image, and pre-processing the image sequence and calculating the direction field according to the segmentation result; Module M2: using a group registration algorithm based on global or local low-rank difference measure to register the image sequence and estimate the displacement field; Module M3: calculating the left ventricular myocardial strain related parameters by using the direction field and the displacement field; In the module M2: The group registration algorithm based on global or local low-rank difference measure is used to estimate the displacement field of the image sequence, and the group registration algorithm is composed of a motion model, an objective function and an optimization strategy: The motion model is used to describe the motion mode of each organ in the image, and determines the flexibility and smoothness of the generated displacement field, including a non-parametric model or a parametric model; The objective function is composed of a difference measure and a regularization term, wherein the difference measure describes the registration degree of the image under different displacement fields, and the regularization term restricts the range of change of the displacement field itself; The optimization strategy finds the optimal displacement field by minimizing the objective function; the optimization strategy of the group registration algorithm is based on algorithms including gradient descent, conjugate gradient and LBFGS algorithm, and adopts a coarse-to-fine multi-resolution optimization framework; The objective function comprises: The difference measure adopted is based on the global low-rank or local low-rank of the registered movie image sequence; The difference measure based on global low-rankness is defined as follows: flatten the two-dimensional images in a movie image sequence into column vectors and arrange them into a Casorati matrix in the temporal direction Take the nuclear norm of the matrix as the difference measure of the image sequence: wherein represent the nuclear norm of a matrix, defined as the sum of all non-zero singular values of a matrix, is the rthnon-zero singular value of is the rank of ​ The difference measure based on local low-rank is defined as follows: a plurality of image blocks are extracted from a movie image sequence along the spatial direction at equal intervals, and Casorati matrices are constructed based on the space-time signals carried by each image block, the kernel norm of each Casorati matrix is calculated, and the sum of the kernel norms of all Casorati matrices is taken as the difference measure of the image sequence: wherein is a Casorati matrix composed of the first image blocks, is the number of image blocks.

8. The magnetic resonance image myocardial strain measurement system based on the group registration algorithm according to claim 7, characterized in that: In the module M1: The magnetic resonance cine images include two-dimensional cardiac cine images and three-dimensional cardiac cine images; the acquisition sequence of the magnetic resonance cine images includes a balanced steady-state free precession cine imaging sequence and a spoiled gradient echo cine imaging sequence; Segment the left ventricular myocardium from the end-diastolic image; According to the segmentation result, pre-process the images, including cropping and gray scale normalization; According to the segmentation result, calculate the unit direction vector at each myocardial pixel to form a direction field; For the short-axis images, calculate the circumferential direction field and the radial direction field; For the long-axis images, calculate the longitudinal direction field and the radial direction field; The objective function includes: The regular term adopted includes a spatial regular term and a temporal regular term; the spatial regular term is a diffusion smoothing term based on first-order difference or a thin-plate bending energy based on second-order difference; the temporal regular term is the square sum of the first or second difference of the displacement field along the time dimension; the difference in the temporal regular term can be a cyclic difference, so that the deformation of the last image and the deformation of the first image remain continuous; the strength of the temporal and spatial regularization is controlled by two regularization coefficients respectively; The optimization strategy includes: The group registration algorithm iteratively solves the minimization problem to obtain the optimal displacement field; the algorithm includes gradient descent method, conjugate gradient method and LBFGS algorithm; the group registration algorithm adopts a multi-resolution optimization framework from coarse to fine, first registers the lower resolution images, and then initializes the registration of higher resolution images with the obtained control point grid; under this framework, the group registration algorithm first generates global deformation at low resolution to align the target to a predetermined extent, and then generates local deformation at high resolution for adjustment; In the module M3: According to the obtained displacement field, track the position of each pixel point in the end-diastolic image myocardial region in each frame of image; based on the group registration method, first map the pixel points to be tracked to an average coordinate system, and then map them to the remaining frames of images; Combined with the calculated direction field, calculate the strain of each pixel in the myocardial region according to the definition of myocardial strain to obtain the myocardial strain map; for the short-axis images, calculate the circumferential strain map and the radial strain map; for the long-axis images, calculate the longitudinal strain map and the radial strain map; average the strain in the whole left ventricular myocardial region to obtain the curve of global myocardial strain change over time; take the derivative difference of the global myocardial strain to obtain the curve of myocardial strain rate change over time.

Images