In vivo temporal composite imaging method and system based on respiratory phase tracking

Abstract

This invention discloses a method and system for in vivo temporal composite imaging based on respiratory phase tracking, belonging to the field of medical imaging. The method includes the following steps: acquiring a digital image sequence of a living subject; calculating the quasi-characteristic curve of respiratory motion by analyzing statistical information between images in the image sequence; analyzing the quasi-characteristic curve in the time and frequency domains to obtain respiratory feature images, respiratory feature curves, respiratory frequency, and respiratory phase curves; dividing the respiratory cycle and phase intervals according to requirements; obtaining selected frame images through region image retrieval, displacement correction, and outlier image removal; and then fusing the images of the target phase interval to obtain an image with enhanced signal-to-noise ratio. This method can obtain target phase images with effectively suppressed respiratory motion artifacts, improving the signal-to-noise ratio and visualization level of the image.

Description

Technical Field

[0001] This invention relates to the field of medical imaging, and more specifically, to a method and system for in vivo time-lapse imaging based on respiratory phase tracking. Background Technology

[0002] With the continuous maturation of medical imaging technologies, such as X-ray, computed tomography, magnetic resonance imaging, and ultrasound imaging, these technologies are now widely used in clinical practice, making significant contributions to disease treatment and medical discovery and progress.

[0003] However, during in vivo imaging, interference caused by physiological movements such as respiration, heartbeat, and vascular pulsation can lead to distortion of tissue structures and severe image artifacts in the acquired images.

[0004] In related technologies, for example, Chinese patent document CN114521152A discloses a medical system with a respiratory monitoring system. This system uses an external respiratory monitoring system to collect current motion signals and synchronizes the imaging dataset with the measured motion signals based on the difference between the current signals and the expected motion signals. The drawback of this approach is the need for additional motion monitoring equipment and the requirement to maintain good synchronization between devices during imaging. Another example is Chinese patent document CN114795182A, which discloses a method and components for eliminating magnetic resonance imaging artifacts using navigation echo. Phase recognition is performed using a navigation echo sequence, and when it is determined that the patient is not in motion, a preset scanning imaging sequence is used for scanning imaging. The drawback of this approach is that imaging is only triggered at the end of respiratory motion, discarding other motion phase images, which leads to a decrease in the effective frame rate and a longer acquisition time. For example, Chinese patent document CN111681740A discloses a respiratory separation strain imaging method based on in vivo ultrasound images, which is an invention patent previously filed by the applicant. This invention solves the problems of determining the respiratory cycle and performing image registration separately during the exhalation and inhalation phases. However, this scheme only utilizes images from the "exhalation" and "inhalation" states, discarding images from other motion states (phases). Furthermore, threshold filtering further reduces the temporal density of the selected image sequence. Another example is Chinese patent document CN114240815A, which discloses a multi-threaded strain imaging method and device based on in vivo ultrasound images, also an invention patent previously filed by the applicant. This scheme utilizes the correlation characteristics of in vivo ultrasound images, combined with threshold setting, displacement correction, image registration, strain imaging, and multi-threaded strain fusion, to obtain accurate, high temporal density strain images. However, it only relies on multi-phase motion suppression based on the "exhalation" and "inhalation" states, resulting in the loss of respiratory motion phases other than "exhalation" and "inhalation," and lacks further image fusion, denoising, and enhancement. In summary, how to effectively suppress the interference of physiological movements such as respiration on imaging and obtain live imaging that meets real-time and multi-phase tracking requirements without the aid of additional monitoring equipment is a problem that needs to be solved by existing technologies. Summary of the Invention

[0005] 1. The technical problem that the invention aims to solve

[0006] The purpose of this invention is to overcome the shortcomings of the prior art and provide a live time composite imaging method and system based on respiratory phase tracking. By obtaining the respiratory motion phase of the image and combining the methods of target phase interval retrieval, displacement correction and multi-image fusion, a target phase image with suppressed motion artifacts and improved signal-to-noise ratio and an image set of multiple respiratory phases are obtained.

[0007] 2. Technical Solution

[0008] The objective of this invention is achieved through the following technical solutions.

[0009] The in vivo temporal composite imaging method based on respiratory phase tracking includes the following steps:

[0010] Acquire digital image sequences of living objects;

[0011] Analyze the statistical information between images in the image sequence to calculate the quasi-characteristic curve of respiratory motion;

[0012] By analyzing the quasi-characteristic curves in the time and frequency domains, respiratory characteristic images, respiratory characteristic curves, respiratory frequency, and respiratory phase curves are obtained.

[0013] Divide the respiratory cycle and phase intervals;

[0014] Region image retrieval and displacement correction: In the image sequence of the target phase interval, select one frame image in the initial period according to the requirements, and select the region of interest as the reference image according to the target area position. Retrieve and perform displacement correction in the images of subsequent periods to obtain the selected frame image sequence.

[0015] Outlier image removal: Based on the similarity between the selected frame images, outlier images in the selected frame image sequence are identified and removed to obtain an updated selected frame image sequence;

[0016] Image fusion in the target phase range: Multi-image fusion is performed on selected frames of the target phase range to obtain an image with enhanced signal-to-noise ratio.

[0017] Furthermore, in the step of acquiring digital image sequences of a living subject, for a living subject in respiratory motion, image sequence I of its target tissue is continuously acquired for at least two respiratory cycles.

[0018] Furthermore, the specific steps for calculating the quasi-characteristic curve of respiratory motion are as follows: take the m-th image from the N-frame image sequence I as the template frame I. ref Using all images as target frames I tar , by I ref with I tar The quasi-characteristic curve γ was calculated. m (n), where m = 1, 2, ..., M, n = 1, 2, ..., N, N and M are natural numbers, M frames of images cover at least one complete respiratory cycle, and a total of M quasi-feature curves are obtained.

[0019] Furthermore, the specific steps for analyzing quasi-feature curves are as follows: after transforming all quasi-feature curves from the time domain to the frequency domain, record the quasi-feature curve γ corresponding to the m-th frame image. m The maximum spectral amplitude of (n) is A m If A pFor all A m The maximum value in, i.e., A p =max(A1,A2,…,A) M If the p-th frame image is denoted as the breathing feature image, then the quasi-feature curve γ is... p (n) is denoted as the respiratory characteristic curve, where p = 1, 2, ..., M;

[0020] In γ p In the spectrum of (n), A p The frequency at which it occurs is the respiratory rate f. res Remove γ using a low-pass filter p (n) greater than f res The frequency components are used to obtain the filtered respiratory characteristic curve γ'. p (n), will γ' p (n) Subtract its mean amplitude to obtain the zero-mean curve γ' p,0 (n);

[0021] With γ' p,0 Using (n) and its Hilbert transform as the real and imaginary parts respectively, the phase curves distributed in the interval [-π, π] are calculated.

[0022] Furthermore, in the steps of dividing the respiratory cycle and phase intervals, based on the phase curve... The extreme points are used to divide all images into Ψ respiratory cycles, where Ψ is a natural number. The maximum or minimum point in the phase curve is found and used as the boundary of the respiratory cycle.

[0023] In the phase curve The image sequence I is divided into Ψ respiratory cycles by dividing the image sequence into Ψ respiratory cycles by traversing the image sequence completely in this way, with each increase in phase amplitude from -π to π.

[0024] Within the phase amplitude range [-π, π], all images are divided into Φ consecutive phase intervals as needed to obtain the phase intervals belonging to the target. Image sequences I with different periods c a,c ;

[0025] Where Φ is a natural number, a = 1, 2, ..., Φ, and c = 1, 2, ..., Ψ.

[0026] Furthermore, the specific steps for region image retrieval and displacement correction are as follows: [The text abruptly shifts to a different topic] ...image sequence I within the target phase interval... a,c In the process, a frame of image is selected from the initial cycle, and a region of interest is selected from it as a reference image based on the target area location. This process is then applied to the image sequence I within the target phase interval and the first respiratory cycle. a,cImages were selected from the first quartile to the last quartile in ascending order of time coordinates, and then those with the largest respiratory amplitude γ were chosen. p The image is used as the initial image;

[0027] The region of interest (ROI) with dimensions A×L in the initial image, where A is the axial direction and L is the lateral direction, and the top-left corner coordinates are (z0, x0), is taken as the template frame R. ref Insert image set R a,c .

[0028] Furthermore, the specific steps for region image retrieval and displacement correction are as follows: In I a,c In the next period of the in-phase interval, take the value of R for each frame image. ref Regions of the same position and size are used as detection frames R tar In the neighboring translation detection frame, the range is axial -A / 3 to A / 3 and lateral -L / 3 to L / 3, so that R tar With R ref The two-dimensional cross-correlation coefficient is maximized, and the translated image is updated to R. tar ;

[0029] Compare all R tar With R ref The cross-correlation coefficients, and the R value corresponding to the maximum value. tar Update to template frame R ref Insert set R a,c At the end, using the updated R ref The same operation is performed on the next cycle of the same phase interval, and this process is repeated until all cycles have completed the retrieval and displacement correction. This process is repeated for all phase intervals to obtain an image set R of Φ phases. a,c , where a = 1, 2, ..., Φ.

[0030] Furthermore, the specific steps for removing outlier images are as follows: based on the similarity between the selected images, determine and remove R-values. a,c The outlier image in the image is used to obtain the updated frame selection image R. a,c For the target phase interval Based on all selected frame images R a,c Calculate the cross-correlation coefficient Γ between adjacent frames a,c The selected frame image R of the c-th cycle a,c Compared with the frame selection image R of the previous cycle a,c-1 The cross-correlation coefficient between adjacent frames in this period is calculated. By traversing all periods in this way, the cross-correlation coefficient Γ between adjacent frames can be obtained. a,c Where c = 1, 2, ..., Ψ;

[0031] Based on the median of the cross-correlation coefficients of historical neighboring frames Γmed Set threshold Γ th =Γ med -0.15, where the historical neighbor frame cross-correlation number of the c-th period represents the set of cross-correlation numbers between adjacent frames with a period range of 1 to (c-1);

[0032] The cross-correlation coefficient between the selected frame in the current period and the selected frame in the previous period is lower than Γ. th If the selected frame of that phase interval within the current period is discarded, then the selected frame of that phase interval is discarded.

[0033] Furthermore, the specific steps of image fusion are as follows: for the target phase interval Selected frame image R a,c Multi-image fusion is performed to obtain an image R with enhanced signal-to-noise ratio. a , consisting of all phase intervals The resulting enhanced image R a That is, the image set R that constitutes multiple breathing phases, where a = 1, 2, ..., Φ.

[0034] The imaging system based on the above-described in vivo temporal composite imaging method using respiratory phase tracking,

[0035] Includes an acquisition module for acquiring digital image sequences of live objects;

[0036] The calculation module is used to analyze the statistical information between images in the image sequence and calculate the quasi-characteristic curve of respiratory motion;

[0037] The quasi-feature curve module analyzes the time and frequency domains to obtain respiratory feature images, respiratory feature curves, respiratory rate, and respiratory phase curves.

[0038] The module for dividing respiratory cycles and phase intervals is used to divide the image into Ψ respiratory cycles and Φ phase intervals;

[0039] The region image retrieval and displacement correction module is used to retrieve and perform displacement correction in the image to obtain a selected frame image sequence;

[0040] The outlier removal module is used to identify and remove outlier images from the selected frame image sequence based on the similarity between the selected images, so as to obtain the updated selected frame images.

[0041] The image fusion module is used to perform multi-image fusion on selected frames of the target phase interval to obtain an image with enhanced signal-to-noise ratio, and the enhanced images obtained from all phase intervals constitute a multi-breathing phase image set.

[0042] 3. Beneficial effects

[0043] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0044] The present invention relates to a live in vivo time-lapse imaging method and system based on respiratory phase tracking. Based on phase amplitude, the continuously acquired image sequence is divided into different periods and phase intervals. By performing region image retrieval, displacement correction and multi-image fusion, the target phase image with good motion suppression can be obtained, which improves the noise suppression effect and signal-to-noise ratio. During imaging, the target tissue images of multiple respiratory phases can be tracked in real time and dynamically. The parallel data processing and computing structure has extremely high real-time application potential. Attached Figure Description

[0045] Figure 1 This is a schematic flowchart of a live-cell temporal composite imaging method based on respiratory phase tracking in one embodiment of the present invention.

[0046] Figure 2 This is a schematic diagram of simultaneous acquisition of ultrasound images during thermal ablation of living tissue in another embodiment of the present invention.

[0047] Figure 3 This is a schematic diagram of the respiratory feature curve in the expiratory state of a respiratory feature image in another embodiment of the present invention.

[0048] Figure 4 This is a schematic diagram of the respiratory feature curve in the inhalation state in another embodiment of the present invention.

[0049] Figure 5 This is a schematic diagram of respiratory phase curves of respiratory feature images in the expiratory state, as shown in another embodiment of the present invention.

[0050] Figure 6 This is a schematic diagram of respiratory phase curves of respiratory feature images in the inspiratory state, as shown in another embodiment of the present invention.

[0051] Figure 7 This is a schematic diagram of respiratory feature curves during two respiratory cycles when the respiratory feature image is in the expiratory state, according to another embodiment of the present invention.

[0052] Figure 8 The respiratory motion phase interval in another embodiment of the present invention A schematic diagram of the cross-correlation coefficient curves between neighboring frames;

[0053] Figure 9 The respiratory phase interval in another embodiment of the present invention A schematic diagram of B-mode ultrasound after multi-image fusion;

[0054] The labels in the diagram are as follows: 01, B-mode ultrasound imaging probe; 02, ultrasound imaging system; 03, microwave ablation needle; 04, microwave heater; 05, computer. Detailed Implementation

[0055] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0056] Example 1

[0057] Combination Figure 1 The present invention provides a live-cell temporal composite imaging method based on respiratory phase tracking, comprising the following steps:

[0058] Acquiring digital image sequences of live objects

[0059] For a living subject undergoing respiratory motion, continuously acquire image sequences (Sequence I) of the target tissue area over at least two respiratory cycles. It is important to note that the living subject does not need to hold its breath during image acquisition; spontaneous breathing or mechanically assisted breathing is permitted. Spontaneous movement and movement caused by external forces should be avoided during image acquisition. The acquired image data can be images directly generated by the imaging system or images that have undergone further processing (such as downsampling, filtering, and compression).

[0060] In a specific technical solution, imaging methods such as ultrasound imaging, magnetic resonance imaging, and computed tomography can be selected according to the needs of diagnosis or treatment. Based on the location, structure, and size of the target tissue, imaging modes, imaging probes, and imaging systems with different imaging performance can be selected. Specifically, image data should be acquired continuously and at equal intervals.

[0061] Application scenarios include, but are not limited to, conventional imaging and image-guided diagnosis and treatment. For applications such as thermotherapy, elastography, temperature imaging, or flow imaging, target tissue is heated or stressed while acquiring image data. Thermotherapy includes, but is not limited to, techniques such as radiofrequency ablation, microwave ablation, laser ablation, and focused ultrasound ablation; elastography includes, but is not limited to, elastography methods utilizing quasi-static pressure, acoustic radiation force, and shear force; temperature imaging includes, but is not limited to, thermal strain imaging, thermal imaging of backscattered energy changes, and thermal imaging of frequency shifts.

[0062] Calculate the quasi-characteristic curve of respiratory motion

[0063] The quasi-characteristic curve of respiratory motion is calculated by analyzing the statistical information between images in image sequence I. Specifically, the m-th frame of image sequence I, consisting of N frames, is used as template frame I. ref Using all images as target frames I tar , by I ref with I tar The quasi-characteristic curve γ was calculated. m(n), where m = 1, 2, ..., M, n = 1, 2, ..., N, and N and M are natural numbers. M frames of images cover at least one complete respiratory cycle, resulting in a total of M quasi-feature curves. It is worth noting that quasi-feature curves can be based on features between images, or on the amplitude and gradient information of each image. For example, the features of each image can be independently represented using indicators such as mean, standard deviation, energy, contrast, and entropy. Taking normalized cross-correlation as an example, the quasi-feature curve γ is calculated. m The process for (n) is as follows:

[0064]

[0065] Where COV represents the calculation of the covariance between the template frame and the target frame, and σ represents the standard deviation of the image.

[0066] Specifically, methods for analyzing statistical information between images include, but are not limited to, cross-correlation, mutual information, gray-level histograms, or gray-level co-occurrence matrices; methods for calculating cross-correlation between images include, but are not limited to, normalized cross-correlation algorithms, mean absolute difference algorithms, absolute error sum algorithms, and error sum of squares algorithms.

[0067] Calculate respiratory feature images, respiratory feature curves, respiratory rate, and respiratory phase curves.

[0068] Quasi-characteristic curves are analyzed in the time and frequency domains to obtain respiratory feature images, respiratory feature curves, respiratory rates, and respiratory phase curves corresponding to the image sequences. Specifically, after transforming all quasi-feature curves from the time domain to the frequency domain, the quasi-feature curve γ corresponding to the m-th frame image is recorded. m The maximum spectral amplitude of (n) is A m If A p For all A m The maximum value in, i.e., A p =max(A1,A2,…,A) M If the p-th frame image is denoted as the breathing feature image, then the quasi-feature curve γ is... p (n) is denoted as the respiratory characteristic curve, where p = 1, 2, ..., M.

[0069] The time-domain signal γ m(n) Converting to the frequency domain requires finer spectral resolution. The resolution of the converted spectrum is less than 0.1 Hz, which can be achieved using spectral refinement methods such as linear frequency modulation z-transform and refined fast Fourier transform. It is worth noting that M-frame images covering at least one complete respiratory cycle ensure a stable expiratory state during respiratory motion at point p in the respiratory feature image. If the M-frame images cover less than one complete respiratory cycle, but at least 1 second of breathing time, the respiratory feature image may be in a substable inspiratory state. Substable means that the inspiratory state has faster respiratory motion and fluctuations compared to the expiratory state; therefore, the inspiratory state is substable compared to the expiratory state. The respiratory feature curve γ... p (n) and phase curve The methods may change accordingly, but the methods in this invention remain effective.

[0070] In γ p In the spectrum of (n), A p The frequency at which it occurs is the respiratory rate f. res Design a low-pass filter and remove γ p (n) greater than f res The frequency components are used to obtain the filtered respiratory characteristic curve γ'. p (n). γ' p (n) Subtract its mean amplitude to obtain the zero-mean curve γ' p,0 (n); with γ' p,0 Using (n) and its Hilbert transform as the real and imaginary parts respectively, the phase curves distributed in the interval [-π, π] are calculated. By γ p (n) yields γ' p The process (n) can also employ other smoothing filters, including mean filtering, median filtering, local regression, Savitzky-Golay filtering, etc. Here, based on γ'... p,0 The phase curve is obtained by taking the arctangent of the ratio of the aforementioned imaginary and real parts (n). The separation and extraction of respiratory phases can also be achieved by combining methods such as matched filtering and phase correction.

[0071] Dividing the respiratory cycle and phase intervals

[0072] According to the phase curve The extreme points are used to divide all images into Ψ respiratory cycles, where Ψ is a natural number. Specifically, the maximum or minimum points in the phase curve are found and used as the boundaries of the respiratory cycles.

[0073] In the phase curve The image sequence I is divided into Ψ respiratory cycles by repeatedly traversing the image sequence in this way, where the phase amplitude increases from -π to π. It is worth noting that the respiratory phase changes periodically between [-π, π], therefore two adjacent phases can be selected. The point can be used as the boundary of a respiratory cycle; alternatively, it can be directly determined from the respiratory characteristic curve γ'. p The extreme points of (n) are used to divide the respiratory cycle.

[0074] The physiological motion process of the obtained respiratory cycle is related to the motion state of the respiratory feature image p, as shown by the phase curve. Taking the respiratory cycle formed by two adjacent maxima as an example, if the respiratory feature image p is in the exhalation state, then each respiratory cycle is a process of inhalation to exhalation and then back to inhalation; if the respiratory feature image p is in the inhalation state, then each respiratory cycle is a process of exhalation to inhalation and then back to exhalation.

[0075] Within the phase amplitude range [-π, π], all images are divided into Φ consecutive phase intervals as needed to obtain the phase intervals belonging to the target. Image sequences I with different periods c a,c Where Φ is a natural number, a = 1, 2, ..., Φ, and c = 1, 2, ..., Ψ. It is worth noting that when dividing the phase interval, it can be based on phase or time, or multiple phase intervals can be set manually; specifically, it can be a division with equal phase or time intervals, or a division with unequal phase or time intervals. Since multi-phase tracking has a parallel processing structure, increasing the number of phases tracked simultaneously does not significantly increase the computational burden.

[0076] Region image retrieval and displacement correction

[0077] Image sequence I in the target phase interval a,c In this process, a frame from the initial cycle is selected based on requirements, and a region of interest is chosen from it as a reference image based on the target area location; specifically, the image sequence I within the target phase interval and the first respiratory cycle is selected. a,c Images were selected from the first quartile to the last quartile in ascending order of time coordinates, and then those with the largest respiratory amplitude γ were chosen. p The image is used as the initial image. The selection of the first respiratory cycle image can also be based on factors such as time coordinates, phase amplitude, and the amplitude of the respiratory characteristic curve. The target phase interval is used as the initial image. For example, the region of interest with size A (axial) × L (horizontal) and top-left corner coordinates (z0, x0) in the initial image is taken as the template frame R. ref Insert an image set R representing the target phase interval and different periods.a,c .

[0078] Shift correction is performed on images in subsequent periods of the initial period to obtain the selected frame image sequence R. a,c Specifically, in I a,c In the next period of the in-phase interval, take the value of R for each frame image. ref Regions of the same position and size are used as detection frames R tar To correct for the offset caused by motion, in the neighboring translation detection frame (range: axial -A / 3~A / 3, lateral -L / 3~L / 3), R is made... tar With R ref The two-dimensional cross-correlation coefficient is the largest, R tar Update to the translated image; compare all R images. tar With R ref The cross-correlation coefficient of R ref Update to R corresponding to the maximum value. tar And put it into set R a,c At the end. Using the updated R ref The same operation is performed in the next cycle of the same phase interval, and so on, until all cycles have completed the search and displacement correction.

[0079] By processing all phase intervals as described above, a total of Φ phase image sets R can be obtained. a,c (a = 1, 2, ..., Φ); It should be noted that known methods such as optimization based on amplitude invariance and optical flow can also be used to calculate the offset caused by motion.

[0080] Remove outlier images

[0081] Based on the similarity between selected frames, determine and remove the selected frame image sequence R. a,c The outlier images in the sequence are used to obtain the updated frame selection image sequence R. a,c Specifically, for the target phase interval Based on all selected frame images R a,c Calculate the cross-correlation coefficient Γ between adjacent frames a,c The selected frame image R of the c-th cycle a,c Compared with the frame selection image R of the previous cycle a,c-1 The cross-correlation coefficient between adjacent frames in this period is calculated. By traversing all periods in this way, the cross-correlation coefficient Γ between adjacent frames can be obtained. a,c Where c = 1, 2, ..., Ψ. The formula for calculating the cross-correlation coefficient can be a known method such as normalized cross-correlation or summation of absolute differences.

[0082] Based on the median of the cross-correlation coefficients of historical adjacent frames Γ med Set threshold Γ th =Γ med-0.15; Specifically, the historical neighbor frame cross-correlation coefficient of the c-th period represents the set of cross-correlation coefficients between adjacent frames in the period range of 1 to (c-1). If the cross-correlation coefficient between the selected frame in the current period and the selected frame in the previous period is lower than Γ... th If the selected frame in the current period is not selected, then the selected frame in that phase interval is discarded. It is worth noting that the threshold for judging outlier frames can be set according to the mean, median or upper quartile of the cross-correlation coefficients of historical neighboring frames, or the similarity of neighboring frame images can be calculated using methods such as mean square error, structural similarity, hash similarity and mutual information.

[0083] Image fusion of target phase range

[0084] For the target phase interval Selected frame image sequence R a,c Multi-image fusion is performed to obtain an image R with enhanced signal-to-noise ratio. a ; composed of all phase intervals The resulting enhanced image R a This refers to the image set R that constitutes the multi-breathing phase. It's worth noting that basic image fusion strategies include, but are not limited to, pixel-level, feature-level, and decision-level image fusion. Taking pixel-level image fusion methods as an example, these include, but are not limited to, well-known image fusion methods such as averaging, weighted averaging, pixel grayscale selection, and principal component analysis.

[0085] As the acquisition time increases, the number of images captured in the target phase gradually increases, thus significantly enhancing the motion suppression effect after multi-image fusion.

[0086] The in vivo temporal composite imaging system based on respiratory phase tracking of the present invention includes:

[0087] Data Acquisition Module:

[0088] For living subjects in respiratory motion, image sequences of their target tissues are continuously acquired over at least two respiratory cycles.

[0089] The acquired image data can be images directly generated by the imaging system or images that have undergone further processing, such as downsampling, filtering, and compression. Based on the location, structure, and size of the target tissue, imaging modes, probes, and systems with different imaging performance should be selected. During acquisition, image data should be acquired continuously and at equal intervals.

[0090] Calculation module:

[0091] For the image sequence acquired by the acquisition module, the statistical information between the images is analyzed to calculate the quasi-characteristic curve of respiratory motion. A portion of the images in the image sequence is taken as template frames, which must cover at least one complete respiratory cycle. All images are used as target frames, and the quasi-characteristic curve reflecting the periodicity and phase characteristics of respiratory motion is calculated from the template frames and the target frames.

[0092] Module for analyzing quasi-characteristic curves in the time and frequency domains:

[0093] The quasi-feature curves obtained by the calculation module are analyzed in the time and frequency domains to obtain respiratory feature images, respiratory feature curves, respiratory frequencies, and respiratory phase curves corresponding to the image sequences. Specifically, after converting the quasi-feature curves from the time domain to the frequency domain, the maximum spectral amplitude of each curve is recorded. The image with the maximum spectral amplitude among all curves is recorded as the respiratory feature image, its quasi-feature curve is recorded as the respiratory feature curve, and its frequency is the respiratory frequency. It should be noted that there may be other motion frequency components in this spectrum, such as heartbeats. By setting a preset frequency range and combining the amplitude relationship between spectral peaks, the frequencies of motions such as heartbeats can be extracted.

[0094] Design a low-pass filter and remove frequency components greater than the breathing frequency from the breathing characteristic curve to obtain the filtered breathing characteristic curve. Subtract its mean amplitude to obtain the zero-mean curve. Use the zero-mean curve and its Hilbert transform as the real and imaginary parts, respectively. Then, calculate the phase curve distributed between [-π, π] based on the arctangent of the ratio of the aforementioned imaginary and real parts formed by the zero-mean curve.

[0095] Module for dividing the respiratory cycle and phase intervals:

[0096] Based on the extreme points of the phase curve, all images are divided into several respiratory cycles. The maximum or minimum points in the phase curve are identified and used as the boundaries of the respiratory cycles. On the phase curve, images corresponding to each increase in phase amplitude from -π to π are grouped into the same respiratory cycle. This process is repeated to completely traverse the image sequence, dividing it into several respiratory cycles. The physiological movement process of the resulting respiratory cycle is related to the motion state of the respiratory feature image. Taking a respiratory cycle formed by two adjacent maximum points of the phase curve as an example, if the respiratory feature image is in the exhalation state, then each respiratory cycle is a process of inhalation to exhalation and then back to inhalation; if the respiratory feature image is in the inhalation state, then each respiratory cycle is a process of exhalation to inhalation and then back to exhalation.

[0097] Within the phase amplitude range [-π, π], all images are divided into several consecutive phase intervals as needed to obtain image sequences of different periods belonging to the target phase interval.

[0098] Region image retrieval and displacement correction module:

[0099] In the image sequence of the target phase interval after being divided by the respiratory cycle and phase interval modules, a frame of the initial cycle is selected as required, and the region of interest is selected as the reference image based on the target area position.

[0100] For the image sequence within the target phase interval and the first respiratory cycle, images between the first and last quartiles are selected as candidate images in ascending order of time coordinates. The image with the largest respiratory characteristic amplitude is then chosen as the initial image. The selection of the first respiratory cycle image can also be based on time coordinates, phase amplitude, and the amplitude of the respiratory characteristic curve.

[0101] In the images of subsequent cycles following the initial cycle, displacement correction is performed to obtain a sequence of selected frames. Specifically, in the next cycle of the in-phase interval in the image sequence, a region with the same position and size as the template frame is selected as the detection frame for each frame. To correct for the offset caused by motion, the detection frame is translated in the neighborhood to maximize the two-dimensional cross-correlation coefficient between the detection frame and the template frame. The detection frame is then updated with the translated image. The cross-correlation coefficients of all detection frames and the template frame are compared, and the template frame is updated with the detection frame corresponding to the maximum value and placed at the end of the image set. The same operation is performed in the next cycle of the same phase interval using the updated template frame, and this process is repeated until all cycles have completed the retrieval and displacement correction.

[0102] By processing all phase intervals as described above, a set of images of several phases can be obtained, which is the frame selection image sequence.

[0103] Outlier Removal Module:

[0104] Based on the similarity between the selected frame images, outlier images in the selected frame image sequence obtained by the region image retrieval and displacement correction modules are identified and removed to obtain an updated selected frame image sequence. Specifically, for the target phase interval, the cross-correlation coefficient Γ between adjacent frames is calculated based on all selected frame images. a,c The cross-correlation coefficient between adjacent frames in a given period is calculated by comparing the selected frame image from one period with the selected frame image from the previous period. This process is repeated for all periods to obtain the cross-correlation coefficient Γ between adjacent frames. a,c The formula for calculating the cross-correlation coefficient can be a known method such as normalized cross-correlation or summation of absolute differences.

[0105] Based on the median of the cross-correlation coefficients of historical adjacent frames Γ med Set threshold Γ th =Γ med -0.15; If the cross-correlation coefficient between the selected frame in the current period and the selected frame in the previous period is lower than Γ th If so, the selected frame of that phase interval within the current period is discarded.

[0106] Image fusion module:

[0107] The selected frame image sequence of the target phase interval obtained by the outlier removal module is subjected to multi-image fusion to obtain the signal-to-noise ratio enhanced image. The enhanced images obtained from all phase intervals constitute a multi-breathing phase image set.

[0108] The in vivo temporal composite imaging method based on respiratory phase tracking of the present invention can track target tissue images of multiple respiratory phases in real time and dynamically. The respiratory phases to be tracked can be set manually or at equal phase intervals. Since multi-phase tracking has a parallel processing structure, increasing the number of phases tracked at the same time will not significantly increase the computational burden. However, due to the frame rate limitation of the imaging device, the number of phases tracked at the same time cannot be infinitely large, and as the number of phases increases, the motion suppression effect of retrieval and displacement correction may decrease. The method of the present invention is essentially an in vivo image processing method for respiratory motion, and therefore is not limited to specific imaging methods, such as ultrasound imaging, computed tomography, and magnetic resonance imaging.

[0109] Example 2

[0110] This embodiment uses the method described in Embodiment 1 above, employing microwave ablation to thermally ablate live porcine tumor tissue while continuously acquiring B-mode digital ultrasound image signals of the ablation area. A fused multi-respiratory-phase live image is obtained based on a respiratory phase tracking-based in vivo time-lapse composite imaging method. The specific steps are as follows:

[0111] Implementing device such as Figure 2 As shown, the 128-element linear ultrasound probe 01 with a center frequency of 10.5MHz is aligned with the target area, and the imaging depth is adjusted to 4cm so that the target tissue appears on the imaging plane. The trolley-type ultrasound imaging system 02 controls the ultrasound probe to acquire images. The imaging mode adopts column-by-column scanning iso-depth focusing, with a sampling rate of 40MHz and a frame rate of 50Hz.

[0112] The microwave ablation needle 03 is inserted into the target tissue through the wound on the body surface, along the normal direction of the imaging plane, and connected to the microwave heater 04. During the procedure, the computer 05 simultaneously triggers ultrasound image acquisition and microwave ablation, acquiring a 20-second image sequence. Each frame of the digital ultrasound image is 2048 (axial) * 128 (horizontal) pixels in size.

[0113] The thermal ablation lasted for 20 seconds, acquiring N=1000 digital ultrasound images. Based on the possible respiratory cycle, images acquired from 0 to 5 seconds (i.e., M=250) were selected as template frames I. ref Using all N frames of images as target frame I tar Calculate I ref with Itar The two-dimensional normalized cross-correlation yields the quasi-characteristic curve γ of respiratory motion. m (n) (m=1,2,…,M,n=1,2,…,N).

[0114] Using linear frequency modulated z-transform, all characteristic curves γ m (n) The refined spectrum is obtained through transformation, with a frequency band range of 0.1–5 Hz and 1000 spectral points. The respiratory rate is determined as f by comparing the peak amplitudes of the M spectral curves. res =0.34Hz, respiratory feature image p=20 and corresponding respiratory feature curve γ p like Figure 3 .

[0115] respiratory characteristic curve γ p The filtered respiratory characteristic curve γ' was obtained by using an FIR low-pass filter. p The filter's cutoff frequency is 1.5*f res The order is 81. The characteristic curve γ' after zero-mean processing... p,0 Calculated respiratory phase curve like Figure 5 If the M-frame images are insufficient to cover a complete respiratory cycle, and image sequence I is acquired starting from the inspiratory state, the respiratory feature image may be in the inspiratory state (p=119), resulting in a different respiratory feature curve γ. p like Figure 4 The corresponding respiratory phase curve like Figure 6 .

[0116] By finding the maximum points of the respiratory phase curve, the respiratory cycle boundaries with time coordinates of 2.56, 5.34, 8.02, 11, 13.76, 16.72, and 19.8 seconds were obtained, thus dividing the image sequence into Ψ = 7 respiratory cycles. Setting equal phase amplitude intervals of π / 3, the image sequence was divided into Φ = 6 phase intervals. Figure 7 As shown, for ease of demonstration, the respiratory characteristic curve within the time range of 2.8 to 8.8 seconds is used as an example. After the aforementioned division, two respiratory cycles and six phase intervals are obtained.

[0117] During the first respiratory cycle, initial images of six phase intervals were selected, and a region of interest with an axial dimension of 10 mm and a lateral dimension of 12 mm was selected as the reference image R, centered on the microwave ablation needle. ref Within the corresponding phase interval of subsequent cycles, the selected frame image R is obtained through retrieval and shift correction. a,c (a=1,2,...,Φ; c=1,2,...,Ψ). like Figure 7As shown, on the characteristic curves of the two respiratory cycles with times of 2.8 to 8.8 seconds, the distribution of selected frames in six phase intervals is marked with symbols.

[0118] Calculate the selected frame image R a,c The cross-correlation coefficients of historical neighboring frames are calculated, and the value Γ is used as the basis for the cross-correlation coefficients. med Set threshold Γ th =Γ med -0.15, after removing outlier images, the updated frame selection image R is obtained. a,c .like Figure 8 The phase interval is shown. and Distribution of cross-correlation coefficients between historical neighboring frames.

[0119] Select frame image R a,c By taking images within the same phase interval but with different periods, and calculating the pixel mean of multiple images spatially and pixel-wise, a target phase image R with further improved signal-to-noise ratio is obtained. a ; Traversing phases Then, the image set R of multiple respiratory phases is obtained. For example... Figure 9 The phase interval is shown. and The fused image. As the image acquisition time increases, the number of images captured in the target phase also gradually increases, thus significantly enhancing the motion suppression effect after multi-image fusion.

[0120] This invention is a data-driven, prospective method for motion suppression that extracts respiratory rate and phase using only image information, reducing the cost and complexity of in vivo imaging systems. Compared to retrospective in vivo imaging strategies, the process of acquiring, analyzing, capturing, and fusing images in this invention meets the logic of real-time applications. In the related processing of motion suppression, gating technology, displacement correction algorithms, and multi-image fusion methods are combined to improve the final motion suppression effect. Multi-motion phase tracking provides a rich and reliable image set for diagnosis and treatment, which has extremely high application value for both research and clinical practice.

[0121] The present invention has been described in detail above with reference to specific exemplary embodiments. However, it should be understood that various modifications and variations can be made without departing from the scope of the invention as defined by the appended claims. The detailed description and drawings should be considered illustrative only and not restrictive, and any such modifications and variations shall fall within the scope of the invention described herein. Furthermore, the background art is intended to illustrate the current state of development and significance of the technology and is not intended to limit the present invention or the scope of application of the present application.

Claims

1. A live-cell temporal composite imaging method based on respiratory phase tracking, comprising the following steps: Acquire digital image sequences of living objects; Analyze the statistical information between images in the image sequence to calculate the quasi-characteristic curve of respiratory motion; By analyzing the quasi-characteristic curves in the time and frequency domains, respiratory characteristic images, respiratory characteristic curves, respiratory frequency, and respiratory phase curves are obtained. Divide the respiratory cycle and phase intervals; Region image retrieval and displacement correction: In the image sequence of the target phase interval, select one frame image in the initial period according to the requirements, and select the region of interest as the reference image according to the target area position. Retrieve and perform displacement correction in the images of subsequent periods to obtain the selected frame image sequence. Outlier image removal: Based on the similarity between the selected frame images, outlier images in the selected frame image sequence are identified and removed to obtain an updated selected frame image sequence; Image fusion in the target phase range: Multi-image fusion is performed on selected frames of the target phase range to obtain an image with enhanced signal-to-noise ratio.

2. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: In the step of acquiring digital image sequences of a living subject, for a living subject in respiratory motion, image sequence I of its target tissue is continuously acquired for at least two respiratory cycles.

3. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: The specific steps for calculating the quasi-characteristic curve of respiratory motion are as follows: take the m-th image from the N-frame image sequence I as the template frame I. ref Using all images as target frames I tar , by I ref with I tar The quasi-characteristic curve γ was calculated. m (n), where m = 1, 2, ..., M, n = 1, 2, ..., N, N and M are natural numbers, M frames of images cover at least one complete respiratory cycle, and a total of M quasi-feature curves are obtained.

4. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: The specific steps for analyzing quasi-characteristic curves are as follows: after transforming all quasi-characteristic curves from the time domain to the frequency domain, record the quasi-characteristic curve γ corresponding to the m-th frame image. m The maximum spectral amplitude of (n) is A m If A p For all A m The maximum value in, i.e., A p =max(A1,A2,…,A) M If the p-th frame image is denoted as the breathing feature image, then the quasi-feature curve γ is... p (n) is denoted as the respiratory characteristic curve, where p = 1, 2, ..., M; In γ p In the spectrum of (n), A p The frequency at which it occurs is the respiratory rate f. res Remove γ using a low-pass filter p (n) greater than f res The frequency components are used to obtain the filtered respiratory characteristic curve γ'. p (n), will γ' p (n) Subtract its mean amplitude to obtain the zero-mean curve γ' p,0 (n); With γ' p Using 0(n) and its Hilbert transform as the real and imaginary parts respectively, the phase curve φ(n) distributed in the interval [-π,π] is calculated.

5. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: In the step of dividing the respiratory cycle and phase interval, all images are divided into Ψ respiratory cycles according to the extreme points of the phase curve φ, where Ψ is a natural number. The maximum or minimum points in the phase curve are found and used as the boundaries of the respiratory cycle. On the phase curve φ, the image corresponding to each increase in phase amplitude from -π to π is divided into the same respiratory cycle. This process is repeated to completely traverse the image sequence I, dividing the image sequence into Ψ respiratory cycles. Within the phase amplitude range [-π, π], all images are divided into Φ consecutive phase intervals as needed, resulting in the target phase interval φ. a Image sequences with different periods c a,c ; Where Φ is a natural number, a = 1, 2, ..., Φ, and c = 1, 2, ..., Ψ.

6. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: The specific steps for region image retrieval and displacement correction are as follows: in the target phase interval, the image sequence I... a,c In the process, a frame of image is selected from the initial cycle, and a region of interest is selected from it as a reference image based on the target area location. This process is then applied to the image sequence I within the target phase interval and the first respiratory cycle. a,c Images were selected from the first quartile to the last quartile in ascending order of time coordinates, and then those with the largest respiratory amplitude γ were chosen. p The image is used as the initial image; The region of interest (ROI) with dimensions A×L in the initial image, where A is the axial direction and L is the lateral direction, and the top-left corner coordinates are (z0, x0), is taken as the template frame R. ref Insert image set R a,c .

7. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 6, characterized in that: The specific steps for region image retrieval and displacement correction are as follows: In I a,c In the next period of the in-phase interval, take the value of R for each frame of the image. ref Regions of the same position and size are used as detection frames R tar In the neighboring translation detection frame, the range is axial -A / 3 to A / 3 and lateral -L / 3 to L / 3, so that R tar With R ref The two-dimensional cross-correlation coefficient is maximized, and the translated image is updated to R. tar ; Compare all R tar With R ref The cross-correlation coefficients, and the R value corresponding to the maximum value. tar Update to template frame R ref Insert set R a,c At the end, using the updated R ref The same operation is performed on the next cycle of the same phase interval, and this process is repeated until all cycles have completed the retrieval and displacement correction. This process is repeated for all phase intervals to obtain an image set R of Φ phases. a,c , where a = 1, 2, ..., Φ.

8. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: The specific steps for removing outlier images are as follows: based on the similarity between the selected images, determine and remove R-values. a,c The outlier image in the image is used to obtain the updated frame selection image R. a,c For the target phase interval φ a Based on all selected frame images R a,c Calculate the cross-correlation coefficient Γ between adjacent frames a,c The selected frame image R of the c-th cycle a,c Compared with the frame selection image R of the previous cycle a,c-1 The cross-correlation coefficient between adjacent frames in this period is calculated. By traversing all periods in this way, the cross-correlation coefficient Γ between adjacent frames can be obtained. a,c Where c = 1, 2, ..., Ψ; Based on the median of the cross-correlation coefficients of historical neighboring frames Γ med Set threshold Γ th =Γ med -0.15, where the historical neighbor frame cross-correlation number of the c-th period represents the set of cross-correlation numbers between adjacent frames with a period range of 1 to (c-1); The cross-correlation coefficient between the selected frame in the current period and the selected frame in the previous period is lower than Γ. th If the selected frame of that phase interval within the current period is discarded, then the selected frame of that phase interval is discarded.

9. The in vivo temporal composite imaging method based on respiratory phase tracking according to claim 1, characterized in that: The specific steps of image fusion are as follows: for the target phase interval φ a Selected frame image R a,c Multi-image fusion is performed to obtain an image R with enhanced signal-to-noise ratio. a , from all phase intervals φ a The resulting enhanced image R a That is, the image set R that constitutes multiple breathing phases, where a = 1, 2, ..., Φ.

10. An imaging system based on the in vivo temporal composite imaging method based on respiratory phase tracking as described in any one of claims 1-9, characterized in that, Includes an acquisition module for acquiring digital image sequences of live objects; The calculation module is used to analyze the statistical information between images in the image sequence and calculate the quasi-characteristic curve of respiratory motion; The quasi-feature curve module analyzes the time and frequency domains to obtain respiratory feature images, respiratory feature curves, respiratory rate, and respiratory phase curves. The module for dividing respiratory cycles and phase intervals is used to divide the image into Ψ respiratory cycles and Φ phase intervals; The region image retrieval and displacement correction module is used to retrieve and perform displacement correction in the image to obtain a selected frame image sequence; The outlier removal module is used to identify and remove outlier images from the selected frame image sequence based on the similarity between the selected images, so as to obtain the updated selected frame images. The image fusion module is used to perform multi-image fusion on selected frames of the target phase interval to obtain an image with enhanced signal-to-noise ratio, and the enhanced images obtained from all phase intervals constitute a multi-breathing phase image set.

Images