Dynamic dr image motion correction method and system based on video stream

Abstract

The application provides a dynamic DR image motion correction method and system based on a video stream, and relates to the technical field of image correction.The method comprises the following steps: calculating the spatial correlation coefficient of a first topographic kernel and a second topographic kernel according to the local gray distribution difference of the two, and constructing a time-varying baseline connecting the first topographic kernel and the second topographic kernel based on the spatial correlation coefficient; when the length of the time-varying baseline in the current frame changes relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, moving the first topographic kernel until the change falls within the preset threshold to obtain an updated first topographic kernel.The application extracts multiple types of topographic kernels, combines spatial correlation and phase matching to complete dynamic updating of the kernel and transmission of the motion vector, relies on the kernel to construct a motion deformation grid to realize pixel-level image transformation, can effectively compensate for the motion offset in the dynamic DR image, suppresses the motion artifacts, and improves the overall clarity of the video stream image.

Description

Technical Field

[0001] This invention relates to the field of image correction technology, and in particular to a method and system for dynamic DR image motion correction based on video streams. Background Technology

[0002] During the acquisition and processing of Dynamic Digital Radiography (Dynamic DR) video sequences, the relative motions of different tissue regions within the subject at various frequencies and amplitudes (e.g., the superposition of low-frequency periodic displacement and high-frequency oscillatory motion) result in non-rigid and inconsistent motion patterns for different tissue structures within the sequence. These motion patterns may cause periodic positional shifts and local deformations in certain local image regions (such as high-density markers or image blocks with stable texture structures) between different frames, thus affecting subsequent image registration, motion estimation, and motion compensation tasks to some extent.

[0003] For example, in a series of dynamic DR chest imaging sequences, the sequence contains 30 consecutive frames, each with a resolution of 1024×1024 pixels and a grayscale depth of 12 bits. Between frames 10 and 15 of the sequence, two motion components with different characteristics can be observed. One component has a period of approximately 15 frames and a vertical displacement amplitude of approximately 8 pixels; the other component has a period of approximately 4 frames and a horizontal displacement amplitude of approximately 2 pixels. At this point, observing the region where the left diaphragm meets the apex of the heart (coordinates approximately pixels (320, 480) to (380, 520)) and the region where the right cardiophrenic angle overlaps with the pericardial fat pad (pixels (680, 400) to (720, 450)), it can be found that the pixel grayscale distribution gradient within these regions is relatively complex, and the relative motion direction and amplitude between different image blocks vary with the frame number.

[0004] When processing this sequence using existing motion estimation methods based on optical flow fields (such as the Horn-Schunck algorithm or the Lucas-Kanade method), invalid values ​​or holes sometimes appear in the optical flow field calculation within the aforementioned coordinate region, meaning that effective optical flow vectors cannot be obtained at certain pixel locations. Furthermore, if feature point registration-based methods (such as SIFT or SURF feature matching) are used, feature point mismatches or systematic drift along the texture direction may occur in the overlapping areas with high texture repetition. This makes it difficult to accurately distinguish between the overall displacement (coupled motion) dominated by low-frequency periodic displacement and the local relative slip (independent motion) dominated by high-frequency oscillatory motion. A potential consequence of these problems is that in subsequent motion-compensated image frames, the local texture at the intersection of the lower edge of the aortic arch and the left main bronchus (e.g., the edge structure near pixel (512, 256)) sometimes exhibits unnatural breaks or misalignments; that is, the originally continuous edge contours appear as stepped or torn artifacts after compensation. Based on observations, these artifacts are misaligned by about 2 to 3 pixels in some frames, causing some interference with the visual continuity of the video sequence and the accuracy of subsequent processing. Summary of the Invention

[0005] This invention provides a method and system for motion correction of dynamic DR images based on video streams, which effectively eliminates problems such as broken or misaligned textures at the edges of anatomical structures and artifacts in dynamic DR images.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0007] Firstly, a dynamic DR image motion correction method based on video streams, the method comprising:

[0008] Step 1: Obtain a sequence of multiple consecutive frames of images from the dynamic DR video stream, and extract the first morphology kernel, the second morphology kernel, and the third morphology kernel for each frame.

[0009] Step 2: Calculate the spatial correlation coefficient between the first and second morphological base kernels based on the difference in local grayscale distribution. Construct a time-varying baseline connecting the first and second morphological base kernels based on the spatial correlation coefficient. When the change in the length of the time-varying baseline in the current frame relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, move the first morphological base kernel until the change falls back to within the preset threshold to obtain the updated first morphological base kernel.

[0010] Step 3: Match the local phase consistency distribution of the third morphological base kernel with the updated local phase consistency distribution of the first morphological base kernel and the local phase consistency distribution of the second morphological base kernel respectively. Select the one with the larger matching degree as the coupling object. Establish a propagation link between the third morphological base kernel and the corresponding coupling object. Transmit the motion vector of the coupling object to the third morphological base kernel to obtain the predicted motion vector of the third morphological base kernel.

[0011] Step 4: Use the updated first topology kernel, second topology kernel, and third topology kernel carrying the predicted motion vector as control nodes to construct a motion deformation mesh; use the motion deformation mesh to perform pixel-by-pixel transformation on the current frame to obtain the motion-corrected image frame.

[0012] Secondly, a dynamic DR image motion correction system based on video streams includes:

[0013] The extraction module is used to acquire a series of consecutive multi-frame image sequences in a dynamic DR video stream, and extract the first morphology kernel, the second morphology kernel and the third morphology kernel for each frame image;

[0014] The update module is used to calculate the spatial correlation coefficient between the first morphology base kernel and the second morphology base kernel based on the difference in local grayscale distribution, and to construct a time-varying baseline connecting the first morphology base kernel and the second morphology base kernel based on the spatial correlation coefficient; when the change in the length of the time-varying baseline in the current frame relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, the first morphology base kernel is moved until the change falls back to within the preset threshold, and the updated first morphology base kernel is obtained.

[0015] The matching module is used to match the local phase consistency distribution of the third morphological base kernel with the updated local phase consistency distribution of the first morphological base kernel and the local phase consistency distribution of the second morphological base kernel, respectively. The one with the larger matching degree is selected as the coupling object. A propagation link is established between the third morphological base kernel and the corresponding coupling object, and the motion vector of the coupling object is transmitted to the third morphological base kernel to obtain the predicted motion vector of the third morphological base kernel.

[0016] The transformation module is used to construct a motion deformation mesh by using the updated first topography kernel, second topography kernel, and third topography kernel carrying the predicted motion vector as control nodes; and to perform pixel-by-pixel transformation on the current frame using the motion deformation mesh to obtain the motion-corrected image frame.

[0017] Thirdly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0018] The above-described solution of the present invention has at least the following beneficial effects:

[0019] This invention uses three types of dedicated morphological kernels accurately extracted from dynamic DR images as core reference nodes for motion correction. It establishes a stable motion benchmark based on the inherent imaging characteristics of the thoracic cavity's anatomical features, avoiding problems such as invalid motion vectors, feature mismatches, and motion trajectory drift that often occur in traditional algorithms in areas with complex overlapping textures and gray-level gradients within the thoracic cavity. By using fixed morphological kernels as the motion analysis carrier, the targeted nature of inter-frame motion information extraction is improved, enabling the differentiation between low-frequency overall coupled motion and high-frequency local independent motion of human thoracic tissues, overcoming the technical bottleneck of existing technologies' inability to distinguish multi-level composite motion. A spatial correlation coefficient is constructed based on the local gray-level distribution differences of the dual morphological kernels, and a dynamically iterative time-varying baseline mechanism is established, capable of adapting to small-amplitude nonlinear displacements and local deformations between consecutive frames in a dynamic DR video stream. Adaptive updates of the first morphological kernel are achieved through baseline length variation threshold constraints, real-time correction of benchmark node position deviations caused by inter-frame motion offsets, effectively overcoming the shortcomings of fixed feature nodes being unable to adapt to dynamic tissue motion and excessive cumulative motion errors in long-term frame sequences.

[0020] A local phase-consistency distribution matching mechanism is employed to select the optimal coupling object. Motion propagation links are established based on high-matching feature associations to achieve motion vector transfer and deduction. Phase-consistency features exhibit strong robustness to image grayscale changes, local deformations, and texture interference. Compared to traditional grayscale matching and gradient matching methods, this effectively improves the accuracy of motion association matching between different morphological bases, avoiding motion vector distortion caused by erroneous coupling. The predicted motion vectors for the target region are deduced, realistically restoring the non-rigid and non-uniform motion characteristics of different tissue regions within the thoracic cavity. A motion deformation mesh is constructed, and pixel-by-pixel fine-grained motion correction is achieved based on the deformation mapping relationship covering the entire mesh domain. The deformation mesh constructed through multi-node linkage can adapt to the complex non-rigid deformation characteristics of the entire thoracic cavity, overcoming the limitation of traditional local motion correction in not being able to take into account global image deformation. This effectively eliminates problems such as anatomical structure edge texture breaks, misalignments, and artifacts in dynamic DR images, improving the imaging continuity between frames and the overall image quality of dynamic DR video streams. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the dynamic DR image motion correction method based on video stream provided in an embodiment of the present invention.

[0022] Figure 2 This is a schematic diagram of a dynamic DR image motion correction system based on video stream provided in an embodiment of the present invention. Detailed Implementation

[0023] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0024] like Figure 1 As shown, embodiments of the present invention propose a dynamic DR image motion correction method based on video streams, the method comprising the following steps:

[0025] Step 1: Obtain a sequence of multiple consecutive frames of images from the dynamic DR video stream, and extract the first morphology kernel, the second morphology kernel, and the third morphology kernel for each frame.

[0026] Step 2: Calculate the spatial correlation coefficient between the first and second morphological base kernels based on the difference in local grayscale distribution. Construct a time-varying baseline connecting the first and second morphological base kernels based on the spatial correlation coefficient. When the change in the length of the time-varying baseline in the current frame relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, move the first morphological base kernel until the change falls back to within the preset threshold to obtain the updated first morphological base kernel.

[0027] Step 3: Match the local phase consistency distribution of the third morphological base kernel with the updated local phase consistency distribution of the first morphological base kernel and the local phase consistency distribution of the second morphological base kernel respectively. Select the one with the larger matching degree as the coupling object. Establish a propagation link between the third morphological base kernel and the corresponding coupling object. Transmit the motion vector of the coupling object to the third morphological base kernel to obtain the predicted motion vector of the third morphological base kernel.

[0028] Step 4: Use the updated first topology kernel, second topology kernel, and third topology kernel carrying the predicted motion vector as control nodes to construct a motion deformation mesh; use the motion deformation mesh to perform pixel-by-pixel transformation on the current frame to obtain the motion-corrected image frame.

[0029] In this embodiment of the invention, three types of chest cavity-specific topographic kernels are extracted as stable motion benchmarks, avoiding vector failure, mismatch, and trajectory drift problems that occur in complex textures and grayscale regions in traditional algorithms. This effectively distinguishes between the overall coupled motion of the chest cavity and its local independent motion. Simultaneously, a time-varying baseline is constructed using the grayscale difference of dual kernels, and the benchmark nodes are adaptively updated to eliminate the accumulated error of fixed feature nodes and adapt to small nonlinear deformations between frames. A motion propagation link is established based on more robust phase consistency matching to deduce the motion vectors of each region, realistically restoring the non-rigid and non-uniform motion characteristics of the chest cavity. A global motion deformation mesh is constructed based on multiple kernels, achieving pixel-by-pixel fine-grained correction, overcoming the limitations of traditional local correction, effectively eliminating image edge breaks and misalignment artifacts, and improving the inter-frame continuity and imaging quality of dynamic DR video.

[0030] In a preferred embodiment of the present invention, step 1 includes:

[0031] Step 100: Obtain a continuous multi-frame image sequence from the dynamic DR video stream, and extract a first morphological kernel, a second morphological kernel, and a third morphological kernel for each frame. The first morphological kernel is located on the gradient ridge line at the junction of the top of the left diaphragm and the apex of the heart; the second morphological kernel is located on the gray-level extreme value band at the overlap of the pericardial fat pad and the diaphragm in the right cardiophrenic angle region; and the third morphological kernel is located at the phase-consistent peak point at the intersection of the lower edge of the aortic arch and the left main bronchus. Specifically, it includes:

[0032] The system continuously acquires real-time video streams output by dynamic DR equipment, and extracts 8 to 16 consecutive image sequences from the continuous video frames as a group of time-series processing units. This method processes standard DR chest cavity anteroposterior images with a resolution of 1024×1024 pixels and 12-bit grayscale depth. For each image frame, three types of specific morphological kernels are extracted. All kernels are uniformly calculated using a 32×32 pixel local window, ensuring consistent calculation scale across the entire sequence of frames. All three types of morphological nuclei are selected from exclusive regions with fixed anatomical positions, unique physiological movement characteristics, and extremely strong imaging stability in the human thoracic cavity. Specifically, the first morphological nucleus is located on the gradient ridge line at the junction of the top of the left diaphragm and the apex of the heart within the pixel coordinates (320,480) to (380,520) of the standard DR thoracic image. This region is dominated by low-frequency respiratory movements and is minimally affected by high-frequency disturbances of the heartbeat. The junction of the diaphragm and the apex of the heart will undergo stable and periodic vertical and overall displacement with respiration. It can stably follow the low-frequency periodic coupled movements of the thoracic cavity to undergo synchronous deformation and displacement, and can reflect the overall low-frequency periodic movement pattern of the left thoracic cavity tissues.

[0033] The second morphological core is located on the extreme grayscale band overlapping the pericardial fat pad and diaphragm in the right cardiophrenic angle region within the pixel coordinates (680, 400) to (720, 450) of a standard DR chest image. This fat pad exhibits stable grayscale and clear boundaries, serving as a stable reference for right-sided chest cavity movement and aiding in matching the correlation characteristics of left and right chest cavity movements. The third morphological core is located at the peak point of the aortic arch crossing the left main bronchus, near the pixel coordinates (512, 256) of a standard DR chest image. This anatomical location is minimally affected by overall respiratory displacement, primarily exhibiting high-frequency, small-amplitude local micro-deformations due to heartbeat and slight tracheal vasodilation and vasoconstriction. The phase-consistency peak characteristics in this region are extremely sensitive to minute pixel displacements and local texture deformations, capturing subtle inter-frame structural changes and characterizing high-frequency, subtle deformation movements in the chest cavity. These three differentiated morphological cores correspond to low-frequency overall coupled movement, stable reference movement, and high-frequency independent local movement of the chest cavity, respectively, comprehensively covering the composite motion patterns of dynamic DR images.

[0034] This embodiment, relying on the fixed anatomical features of the human thoracic cavity, selectively chooses three types of morphological kernels with different locations and functions, abandoning the traditional algorithm's method of randomly selecting points across the entire domain and blindly extracting features. The three types of kernels correspond to the overall motion of the thoracic cavity, the gray-level extreme value region, and the phase feature region, respectively. They have high feature recognition and strong inter-frame stability, and can adapt to the composite motion characteristics of low-frequency overall coupled motion and high-frequency local independent motion superimposed in dynamic DR thoracic images. This avoids the problems of invalid features, redundant features, and feature drift that are prone to occur in traditional feature extraction methods from the source.

[0035] In a preferred embodiment of the present invention, step 2 includes:

[0036] Step 200a: Extract the gray values ​​of all pixels within the local window surrounding the first morphology base kernel to form a first gray-level distribution sequence; extract the gray values ​​of all pixels within the corresponding window surrounding the second morphology base kernel to form a second gray-level distribution sequence; calculate the cumulative distribution function of the first gray-level distribution sequence and the second gray-level distribution sequence respectively, and calculate the L1 norm distance between the two cumulative distribution functions, specifically including:

[0037] Based on the accurate localization of the three types of morphological base kernels completed in step 100, a local rectangular calculation window with a fixed size of 32×32 pixels is extracted from the center pixels of the first and second morphological base kernels, respectively. The total number of pixels in a single window is fixed at 1024 pixels, and the calculation scale and shape of the two windows are completely consistent. The dynamic DR image processed by this method is 12-bit grayscale, and the grayscale values ​​of all pixels are integers, with the global value range fixed at [0, 4095]. For the two local windows, a uniform pixel traversal order from top to bottom, from left to right, and row-first-column is adopted to read the original grayscale integer value of the corresponding coordinate position in the image pixel matrix pixel by pixel. The original grayscale values ​​of all 1024 pixels in the first morphological base kernel window are completely acquired and arranged in the pixel traversal order to generate the first grayscale distribution sequence; the same acquisition method is used to acquire the original grayscale values ​​of all 1024 pixels in the second morphological base kernel window to generate the second grayscale distribution sequence.

[0038] Based on two sets of gray-level distribution sequences, gray-level frequency statistics are performed to obtain window gray-level statistics, i.e., traversing the entire gray-level range. ∈[0,4095], respectively count the gray values ​​in the first window and the second window that are equal to the current gray level. The total number of pixels is used to obtain the grayscale frequency statistics function of the first window. With the gray-level frequency statistics function of the second window Based on the gray-level frequency statistics, gray-level cumulative distribution functions are constructed for the two windows respectively. The cumulative distribution function is used to characterize the cumulative proportion of pixels with gray values ​​less than or equal to the current gray level within the window, and can reflect the overall gray-level distribution pattern of the local area, that is: ;

[0039] in, This represents the total number of pixels in the window, with a fixed value of 1024. The first morphological base window at gray level The cumulative distribution function value at that location; This represents the cumulative distribution function value of the second morphology kernel window at gray level g; grayscale levels for the first window Corresponding pixel frequency; For the second window grayscale level The corresponding pixel frequency.

[0040] To quantify the difference in overall grayscale distribution between the local regions of the first and second morphological base kernels, global difference quantization calculations were performed on the two sets of cumulative distribution functions. The resulting L1 norm distance, used to characterize the difference in grayscale distribution between the two regions, is: ;

[0041] In the formula, A quantitative indicator for the difference in grayscale distribution. The larger the value, the more significant the difference in the overall gray-scale distribution patterns of the local regions of the two base nuclei; The smaller the value, the higher the consistency of gray-level distribution between the two regions, and the stronger the stability of the regional imaging features.

[0042] Step 201a: Calculate the absolute cosine value between the principal local gradient direction angle within the window containing the first morphological base kernel and the principal local gradient direction angle within the window containing the second morphological base kernel; multiply the L1 norm distance by the absolute cosine value, take the reciprocal, and normalize the result to obtain the spatial correlation coefficient, specifically including:

[0043] Sobel gradient calculations were performed on the 32×32 pixel local windows corresponding to the first and second topology kernels, respectively, and the lateral gradient of each pixel within the window was calculated. With longitudinal gradient Based on the gradient data of all 1024 pixels within the window, the overall principal gradient direction angle of the window is statistically fitted. The principal gradient direction angle is used to characterize the overall dominant trend of local window grayscale changes, and the specific calculation formula is as follows: ;

[0044] in, This is the mean of the horizontal gradient of all pixels within the window. The mean of the vertical gradient of all pixels within the window can be used to suppress single-point gradient noise interference and obtain stable main direction features of the region. The gradient principal direction angle of the current window, with its value constrained within a certain range. The principal direction angles of the gradient at the first topographic base window are calculated using this formula. Principal direction angle of gradient in the second morphology base kernel window .

[0045] The gradient direction similarity between the two windows is calculated by solving for the absolute cosine of the angle between the two principal direction angles. This cosine is used to quantify the consistency of the gray-scale spatial variation trends in the two regions. The calculation formula is as follows: ;

[0046] in, The gradient direction similarity coefficient has a value range of [value range missing]. The closer the value is to 1, the stronger the consistency of the gray-level spatial change direction between the two local windows and the more similar the regional motion deformation trend. The L1 norm distance of the gray-level distribution calculated in step 200a is... Similarity coefficient with gradient direction Multiplying these yields a composite feature quantity that combines the differences in grayscale amplitude and structural orientation. The composite difference feature can simultaneously characterize the differences in gray-level distribution and gradient structure in the local region of the dual-kernel matrix. To eliminate the influence of the L1 norm distance dimension and achieve feature inverse mapping (the smaller the difference, the higher the correlation), the inverse of the composite feature is taken to obtain the initial correlation coefficient Rraw. The maximum-minimum normalization method is used to perform global normalization on the initial correlation coefficient, linearly mapping the values ​​to the standard interval. The final spatial correlation coefficient after normalization is obtained. The closer the coefficient is to 1, the stronger the consistency of gray-level distribution, gradient structure, and spatial deformation in the local regions of the first morphological base core and the second morphological base core.

[0047] Calculate the difference between the principal direction angles of the two sets of gradients, and solve for the absolute value of the cosine of the corresponding angle. This value is used to quantify the consistency of the direction of gray-level spatial change between two local regions. The closer the value is to 1, the closer the gray-level gradient change trends are between the two regions. The gray-level distribution L1 norm distance calculated in step 200a is then used as the quantization factor. Multiplying the result by the absolute value of the cosine of the current gradient direction yields a composite feature quantity that integrates the differences in grayscale distribution and gradient direction. This composite feature quantity can simultaneously characterize the differences in grayscale amplitude and structural change direction. Taking the reciprocal of this composite feature quantity eliminates the numerical bias caused by the L1 distance dimension, preventing a single grayscale difference index from dominating the evaluation result. The result after taking the reciprocal is normalized globally using a maximum-minimum normalization algorithm, linearly mapping the value to the standard interval [0,1]. Finally, the spatial correlation coefficient between the first morphological base kernel and the second morphological base kernel is obtained. The closer the coefficient is to 1, the stronger the consistency of grayscale distribution, gradient structure, and spatial deformation in the local regions of the first and second morphological base kernels.

[0048] Step 202a: Starting from the coordinates of the first morphological base kernel on the image plane and ending at the coordinates of the second morphological base kernel, a directed line segment is drawn as the initial representation of the time-varying baseline. The spatial correlation coefficient is then assigned as the weight of the directed line segment. Specifically, this includes:

[0049] Using the Cartesian coordinate system of DR images as a unified reference, the real-time two-dimensional coordinates of the first topographic kernel are read. ( Using this as the starting endpoint of the time-varying baseline, the two-dimensional coordinates of the second topographic kernel, which exhibits highly stable position and no significant deformation drift, are read. ( The two-point connected directed line segment is constructed based on the coordinates of two sets of feature reference points as the time-varying baseline for inter-frame dynamic tracking. The spatial shape of this baseline can synchronously follow the low-frequency respiratory motion of the first morphological base to produce small expansion and contraction and offset, matching the low-frequency periodic overall motion characteristics of the left side of the thoracic cavity.

[0050] The spatial correlation coefficient R, obtained by normalization in step 201a and with a fixed value range of [0,1], is directly assigned as the exclusive feature weight of the time-varying baseline of the current frame. This weight carries the multi-dimensional correlation features of the consistency of local gray-level distribution and gradient deformation direction of the dual-base kernel. When R approaches 1, it indicates strong motion coupling of the dual-base kernel and high reliability of the baseline topology, which can be used as a reliable motion benchmark. When R approaches 0, it indicates weak feature correlation of the dual-base kernel and poor baseline reference. By binding the baseline with quantized weights, traditional geometric line segments are upgraded into dynamic benchmarks that integrate structural features and motion correlation.

[0051] Step 200b: Calculate the ratio of the Euclidean length of the time-varying baseline in the current frame to the Euclidean length of the time-varying baseline in the previous frame. If the ratio is greater than a preset upper threshold or less than a preset lower threshold, it is determined to exceed the preset threshold. Using the current position of the first topographic base core as the center, and taking the straight line perpendicular to the direction of the line connecting the second topographic base core on the corresponding time-varying baseline as a reference, move the first topographic base core gradually along the normal direction of the corresponding reference line, with a step size of the current frame's time-varying baseline length multiplied by a preset scaling factor, where the preset scaling factor is less than 1. After each movement, recalculate the ratio until the ratio falls within a preset tolerance range, stop moving, and record the current position of the first topographic base core as the updated first topographic base core. Specifically, this includes:

[0052] The coordinates of the two ends of the time-varying baseline corresponding to the previous frame and the current frame are extracted respectively, and the length of the previous frame baseline is calculated using the Euclidean distance formula. Current frame baseline length Further, the ratio of inter-frame baseline length changes is calculated using the following formula: For the characteristics of low-frequency respiratory motion in the diaphragm of patients with respiratory distress syndrome (DR), the preset baseline length ratio convergence range is [0.85, 1.15], that is, the lower threshold is 0.85 and the upper threshold is 1.15. This threshold range is suitable for the amplitude of low-frequency deformation of the thoracic cavity under stable breathing conditions, and can effectively distinguish between normal physiological small-amplitude baseline fluctuations and abnormal baseline deviations.

[0053] A threshold determination is performed on the calculated ratio k. If ∈[0.85,1.15], the baseline deformation of the current frame is determined to be a normal low-frequency small-amplitude motion fluctuation of the thoracic cavity, and the position of the first morphological base is accurate and without offset, requiring no correction; if >1.15 or A value <0.85 indicates that the inter-frame baseline deformation is out of tolerance, meaning that the first topographic base core has drifted due to non-rigid deformation of the thoracic cavity, slight equipment vibration, or other disturbances, rendering the original reference coordinates invalid. An adaptive iterative correction process must be initiated. The correction direction is determined based on the baseline topology. Using the current time-varying baseline line as the axial reference, a normal reference direction perpendicular to the baseline is constructed as the sole iterative movement direction of the base core. This ensures that the correction direction conforms to the low-frequency motion displacement pattern of the left thoracic cavity tissue, avoiding irregular offset correction. Specifically, a baseline axial vector is constructed using the reference points at both ends of the current frame's time-varying baseline, and the coordinates of the current frame's baseline starting point are taken. ( ), endpoint coordinates ( Construct baseline axial vector This vector represents the axial extension direction of the baseline from the right cardiophrenic angle to the left apex of the heart. Based on the orthogonality property of planar vectors, two sets of normal vectors perpendicular to the baseline are derived, and the calculation formula is as follows: ;

[0054] The low-frequency respiratory motion in this anatomical region only has effective vertical periodic displacement perpendicular to the baseline axis. The lateral displacement parallel to the baseline is all invalid noise offset. Therefore, only the normal vector is retained as the correction dimension. By matching the displacement projection, a single normal direction that fits the real physiological displacement trend is selected and locked as the only iterative movement direction of the first morphological base core, thus eliminating disordered offset correction.

[0055] To avoid overfitting and position jumps caused by large step size correction, a fixed iteration scaling factor of 0.08 is preset, and the formula for calculating the correction step size in a single iteration is as follows: The base coordinates are corrected iteratively using a small-step fine-tuning method. After each coordinate iteration update, the inter-frame baseline length ratio is immediately recalculated. Continuously perform loop verification and correction until... Once the numerical values ​​fall back to the convergence interval [0.85, 1.15], the kernel position is determined to have returned to the optimal matching state, and the baseline topology is restored to stability, at which point the iteration process is terminated immediately. The kernel coordinates after iteration convergence are recorded as the updated first topographic kernel, completing the dynamic calibration of the low-frequency motion core benchmark. This calibration process relies on the dual-kernel correlation features and quantization threshold constraints, and is data-driven and adaptively iterative throughout, effectively eliminating the cumulative benchmark error of long video frames.

[0056] This embodiment integrates two dimensions of parameters—grayscale distribution difference and gradient direction difference—to calculate the spatial correlation coefficient of the dual-core system and construct a weighted time-varying baseline, overcoming the limitations of traditional single-feature correlation calculation. Grayscale differences are quantified by the L1 norm distance of the cumulative grayscale distribution, and structural differences are characterized by the cosine of the gradient direction. The correlation coefficient after multi-dimensional fusion accurately reflects the inter-frame correlation state of the two core cores. The time-varying baseline dynamically adapts to changes in the basic position between frames, freeing the motion reference from fixed coordinates and improving the comprehensiveness and accuracy of motion correlation analysis. Through baseline length ratio threshold determination and iterative movement correction mechanisms, the first morphological core is adaptively and dynamically updated, effectively solving the shortcomings of traditional fixed feature nodes in adapting to subtle dynamic deformations of human tissue and excessive cumulative errors in long-sequence frame motion. Through refined step-size iterative correction, the reference offset problem caused by small-amplitude nonlinear displacements between frames can be corrected in real time, ensuring the positional accuracy of the core reference node and continuously maintaining the stability of the motion reference.

[0057] In a preferred embodiment of the present invention, step 3 includes:

[0058] Step 300: Calculate the peak values ​​of the cross-correlation surfaces of the phase consistency direction spectra between the third morphological base kernel and the updated first morphological base kernel, and between the third morphological base kernel and the second morphological base kernel, respectively. Use the corresponding peak values ​​of the cross-correlation surfaces as the first and second matching degrees, respectively. Specifically, this includes:

[0059] Based on the updated high-precision first topography kernel obtained in step 200b, a fixed-size 32×32 pixel local feature window is uniformly used to extract local image regions of the updated first topography kernel, the fixed second topography kernel, and the high-frequency response third topography kernel. A multi-scale Log-Gabor filter with fixed parameters is used to extract phase features. The frequency domain expression of the filter is: ;

[0060] in, This is the Log-Gabor filter response function in the frequency domain; The frequency domain independent variable represents any frequency component in the image's frequency domain space, with a fixed center frequency. Bandwidth coefficient A four-scale, six-directional traversal filter was used, with a scale sequence of 0.8, 1.2, 1.6, and 2.0, and an angle sequence of 0, . 6、 / 3、 / 2、2 / 3、5 / 6.

[0061] Convolve each window image with all filter kernels to obtain the complex filter response. ,in pixel coordinates Complex filter response at point; The pixel filtering amplitude represents the saliency of local structural features at the current location; It is the imaginary unit of complex numbers; The original phase component of the pixel represents the structural orientation information of local texture and edges, and is obtained through the four-quadrant arctangent formula. Solve for the phase components of each pixel, where The value of the imaginary part of the complex filter response; The real part of the complex filter response is obtained by solving. The value range is [-π, π]. Combining multi-scale, multi-directional amplitude and phase data, the phase consistency formula is used... Phase consistency values ​​are calculated pixel by pixel, generating 32×32-dimensional phase consistency direction spectrum matrices for the three types of base kernels. For pixels The phase consistency value of the position takes a range of [0,1]. The total number of filter groups is obtained by fixing the four scales and six directions. =24; The filter group number; For the first The pixel amplitude corresponding to the group filter; For the first The pixel phase components corresponding to the group filter; These are the cosine and sine values ​​of the phase angle, respectively.

[0062] Two independent matching combinations are constructed: the third morphology kernel and the updated first morphology kernel, and the third morphology kernel and the second morphology kernel. A two-dimensional normalized cross-correlation algorithm is used to complete the feature matching calculation. The calculation formula is as follows: ;

[0063] Using the baseline feature matrix as a template, iterate through all integer pixel offset coordinates. This yields discrete NCC sampling data points across the entire domain, where... Offset coordinates The normalized cross-correlation coefficient at point [-1, 1] takes values ​​in the range [-1, 1]. The reference phase-consistent direction spectrum matrix; The phase consistency direction spectrum matrix to be matched; The mean value of all pixels in the reference matrix; The mean value of all pixels in the matrix to be matched; These are the pixel coordinates within the matrix. Using the baseline feature matrix as a fixed template, iterate through all feasible integer pixel offset coordinates. The NCC sampling data points are obtained by calculating the discrete distribution of the entire domain point by point.

[0064] To eliminate data discontinuities and numerical abrupt changes in discrete sampling, a bilinear interpolation algorithm is used to complete the global continuous fitting and compensation. For any sub-pixel to be fitted point, arbitrarily selected sub-pixel coordinates are used. Locate the four nearest integer discrete sampling points, namely the top left sampling point. upper right sampling point Lower left sampling point Sampling point in the lower right corner The four points correspond to the known NCC values ​​as follows: Calculate the lateral offset weights of sub-pixel points relative to the integer grid. With vertical offset weight The weight values ​​range from [0,1]. The formula for calculating the fitted values ​​of two horizontal sampling points on the same vertical axis is: ;

[0065] Based on the two sets of horizontal interpolation results, a quadratic vertical interpolation is performed to obtain the continuous NCC values ​​of the final sub-pixel positions. The calculation formula is as follows: .

[0066] Point-by-point interpolation is performed across all sub-pixel intervals in the entire domain to eliminate discrete data breaks and abrupt changes, generating a smooth, continuous, local distortion-free, and sub-pixel-precision two-dimensional cross-correlation surface. The absolute value of the surface peak is extracted as the maximum structural matching coupling strength. The surface peak corresponding to the third morphological base kernel and the updated first morphological base kernel is defined as the first matching degree, used to characterize the structural coupling and motion correlation between the high-frequency micro-deformation region of the thoracic cavity and the low-frequency overall motion reference region. The surface peak corresponding to the third morphological base kernel and the second morphological base kernel is defined as the second matching degree, used to characterize the structural coupling and motion correlation between the high-frequency micro-deformation region of the thoracic cavity and the fixed stable reference region. The values ​​of the first and second matching degrees are in the range of [-1, 1]. The closer the value is to 1, the higher the structural matching degree and the stronger the motion coupling consistency. The value is close to 0, indicating that the two regions are structurally independent and have no obvious motion correlation. The value is close to -1, indicating that the structural characteristics are significantly different and the motion laws are opposite.

[0067] Step 301: Compare the first matching degree with the second matching degree, and select the morphology kernel corresponding to the larger value as the coupling object; construct a straight-line propagation link with the third morphology kernel and the coupling object as the two endpoints, whereby the motion vector at any point on the straight-line propagation link is obtained by linear interpolation of the motion vectors at the two endpoints; propagate the motion vector of the coupling object between the current frame and the previous frame along the straight-line propagation link to the endpoint where the third morphology kernel is located to obtain the predicted motion vector of the third morphology kernel, specifically including:

[0068] Based on the first and second matching degrees obtained in step 300, the optimal coupling screening and high-frequency motion vector derivation are completed. Let the first matching degree be... The peak value of the structure matching between the third morphological core and the updated first morphological core is , and the second matching degree is . The peak value corresponding to the structural matching between the third and second morphological base nuclei is in the range of [-1, 1]. The optimal coupling object is determined through numerical comparison; when the following condition is met... When the updated first topological base kernel is selected as the coupling object of the third topological base kernel; when the condition is met... When the second morphological base kernel is selected as the coupling object of the third morphological base kernel, a larger matching degree value indicates a higher similarity in the phase structure of the two sets of features and a more consistent inter-frame deformation motion coupling law. Using this as the coupling object can ensure the accuracy and rationality of motion vector propagation.

[0069] Let the current frame two-dimensional coordinates of the third topographic kernel be... The current frame two-dimensional coordinates of the optimally coupled object are: ,by and A linear propagation link is constructed for two fixed endpoints. This link constrains the local deformation vector to exhibit a linear and continuous distribution, conforming to the smooth deformation characteristics of thoracic soft tissue. The inter-frame two-dimensional motion vector of the optimally coupled object is defined as the reference propagation vector. The reference motion vector is obtained by calculating the coordinate offset between the current frame and the previous frame, using the following formula: ;

[0070] In the formula, , These are the horizontal and vertical components of the motion vector of the optimally coupled object, respectively; , The current frame coordinates of the optimally coupled object; , The coordinates of the optimally coupled object in the previous frame.

[0071] Based on the principle of two-point linear interpolation, the global propagation derivation of the reference motion vector along the propagation link is completed. Link normalized distance weighting coefficients are defined. The value range is [0,1], and the weight calculation formula is: ;

[0072] In the formula, Let be the coordinates of any point on the propagation link, the denominator be the total Euclidean length of the propagation link, and the numerator be the Euclidean distance from the current point to the coupled object. =0 corresponds to the endpoint of the optimally coupled object, which fully inherits the reference motion vector; =1 corresponds to the endpoint of the third topographic base kernel, which is the position of the predicted vector to be solved. Since the link motion vector follows a linear, uniform, and gradual change law, the predicted motion vector at the third topographic base kernel is obtained by linear interpolation of the vectors at both ends. The motion vector of the third topographic base kernel itself in the previous frame is initially zero-constrained, and is derived linearly only based on the coupled reference vector. The final formula for calculating the predicted motion vector is: ;

[0073] In the formula, , These are the lateral and longitudinal components of the predicted motion vector of the third morphological base core, respectively. This calculation method can smoothly transmit the reference low-frequency motion vector along a straight link to the high-frequency micro-deformation region, preserving the overall low-frequency motion correlation of the thoracic cavity while adapting to the small-amplitude high-frequency deformation characteristics of the region where the third morphological base core is located, ultimately obtaining the two-dimensional predicted motion vector.

[0074] This embodiment employs phase consistency features, which offer stronger anti-interference capabilities, to achieve matching coupling. Compared to traditional grayscale matching and gradient matching methods, it is unaffected by image grayscale fluctuations, minor local deformations, and texture overlap, resulting in higher matching accuracy. By quantifying the matching degree to select the optimal coupling object and constructing a dedicated motion propagation link, it achieves precise transmission of motion vectors. This effectively distinguishes between the overall coupled motion of the thoracic cavity and independent local motions, deducing the true motion state of local areas. It avoids problems such as motion vector distortion and motion estimation bias caused by incorrect coupling, realistically restoring the non-rigid and non-uniform motion characteristics of thoracic tissues.

[0075] In a preferred embodiment of the present invention, step 4 includes:

[0076] Step 400a: Using the updated first topography kernel, second topography kernel, and third topography kernel carrying the predicted motion vector as the initial control node set, four auxiliary nodes are inserted at equal intervals on the boundary of the planar region covered by the initial control node set, forming a planar point set containing seven nodes, specifically including:

[0077] The updated first morphological kernel obtained from the iterative correction in step 200b, the second morphological kernel with fixed coordinates throughout the process, and the third morphological kernel obtained from step 301 carrying high-frequency predicted motion vectors are used together as the initial core control nodes. These three core nodes correspond to the low-frequency overall motion reference, the static stability reference, and the high-frequency micro-deformation motion region of the thoracic cavity, respectively, and can comprehensively characterize the deformation characteristics of different regions of the thoracic cavity. The two-dimensional image coordinates of the three core nodes are recorded and denoted as the coordinates of the updated first morphological kernel. Second morphological base coordinates Third morphological base coordinates Extract the maximum and minimum x and y coordinates, and the minimum x coordinate of the three sets of node coordinates respectively. Maximum x-coordinate Minimum y-coordinate Maximum ordinate The outer rectangular region surrounding all core nodes is constructed using the coordinate extrema. This rectangular region is the effective motion deformation coverage area of ​​the thoracic cavity, which can completely encompass the range of low-frequency overall deformation and local high-frequency micro-deformation.

[0078] Perform equidistant auxiliary node insertion operations on the four boundaries of the circumscribed rectangle to evenly supplement the four boundary auxiliary nodes and improve the node area coverage. The four boundaries are the top boundary, bottom boundary, left boundary, and right boundary of the rectangle. The coordinates of the equidistant auxiliary nodes on each boundary are calculated as follows: the coordinates of the equidistant auxiliary nodes on the top boundary are... The coordinates of the equidistant auxiliary nodes at the lower boundary are: The coordinates of the equidistant auxiliary nodes on the left boundary are The coordinates of the equidistant auxiliary nodes on the right boundary are: Four auxiliary nodes are evenly distributed along the boundary of the deformation region, which can constrain the topology of the mesh edge and avoid problems such as boundary mesh distortion and edge deformation fitting distortion. The above three core control nodes are integrated with the four equidistant boundary auxiliary nodes to form a complete planar node set containing seven nodes. The central region of this node set focuses on the core deformation feature points of the thoracic cavity, and the boundary region achieves full coverage. The nodes are evenly distributed and scale-adapted.

[0079] Step 401a: Traverse the triangle combinations formed by any three points in the planar point set, calculate the circumcircle of each triangle, and select the circumcircle that does not contain any nodes in the planar point set other than the three vertices of the corresponding triangle as candidate circles. Select the candidate circle with the largest radius, and add the center coordinates of the circle with the largest radius as a new control node to the planar point set. Specifically, this includes:

[0080] Traverse the initial seven-node planar node set and combine any three non-collinear nodes to generate all valid triangular elements. Perform circumcircle calculation for each group of triangular elements. Let the coordinates of the three vertices of the triangle be... The center of the circumcircle is calculated using the general formula for solving the circumcircle of a triangle. With the radius of the circumcircle The formula for calculating the coordinates of the circumcircle's center is: ;

[0081] The formula for calculating the radius of the circumcircle is: ;

[0082] In the formula, Let the area be the triangle. .

[0083] After solving for the parameters of all circumcircles of the triangles, each circumcircle is checked for emptiness. All remaining nodes in the current planar node set, excluding the three vertices of the current triangle, are traversed. The Euclidean distance between each remaining node and the center of the circumcircle is calculated. If the distance from all remaining nodes to the center is greater than the radius of the current circumcircle, it is determined that there are no other nodes interfering within the circumcircle, and this circumcircle is marked as a candidate empty circle. The radius values ​​of all candidate empty circles are traversed, and the optimal candidate empty circle with the largest radius is selected. The region containing this circle is the sparsest region of the current mesh and the region with the weakest deformation fitting accuracy. The center coordinates of the optimal empty circle are used as a new encrypted control node and added to the original seven-node planar node set, ultimately forming an encrypted planar node set containing eight sets of control nodes, effectively improving the subsequent mesh's adaptability to local small deformations and edge deformations.

[0084] Step 402a: Obtain all control nodes of the motion deformation mesh and the edges of all triangular elements in the motion deformation mesh, forming an initial planar graph. Calculate the Euclidean distance between the two endpoints of each edge in the initial planar graph as the weight of the corresponding edge. Sort all edges in the initial planar graph in ascending order of weight, establish an empty edge set, and starting from the edge with the smallest weight, sequentially determine whether the current edge forms a cycle with any existing edge in the empty edge set. If the two control nodes connected by the current edge do not belong to the same connected component in the current connected component of the empty edge set, then add the current edge to the empty edge set; otherwise, skip to the next step. Pass through the current edge; repeat this process until the number of edges in the empty edge set is one less than the total number of control nodes, obtaining the minimum spanning tree; use each edge in the minimum spanning tree as a constraint edge, and perform an edge flipping operation on the triangular elements in the motion deformable mesh, that is, for each internal edge shared by two triangles in the motion deformable mesh, calculate the minimum interior angle of the two triangles before flipping and the minimum interior angle of the two triangles after flipping. If the minimum interior angle after flipping is greater than the minimum interior angle before flipping, then perform the flipping; traverse all shared internal edges until there are no more flippable edges, obtaining the optimized motion deformable mesh, specifically including:

[0085] Integrate all control nodes (the encrypted eight-node set) and the connecting edges between nodes to generate an initial planar graph containing all mesh vertices and cell partitioning edges, completing the initial global triangular mesh partitioning. Calculate the weights of all connecting edges in the initial planar graph, using the Euclidean distance between the two endpoints of the edge as the edge weight. Sort all mesh edges in ascending order of weight, and initialize an empty edge set to store constraint edges for the minimum spanning tree. Initialize a disjoint-set data structure for determining node connectivity. The disjoint-set data structure uses standardized initialization rules, with the initialization formula as follows: In the formula, Number any control node. For nodes The parent node. Initially, the eight control nodes are each an independent, single-connected branch, with each node's root node being itself, and the nodes are not connected to each other. Configure a disjoint-set data structure root node lookup function with path compression; the complete iterative calculation formula is: ;

[0086] In the formula, For nodes The unique root node of the connected component. When the node's parent node is itself, this node is the root node of the current branch; when the node's parent node is not itself, recursively search upwards for the root node and update the path synchronously, compressing the node hierarchy and improving the efficiency of connectivity determination. Traverse the sorted grid edges, call the root node search function to solve for the root nodes of the two ends. If the two nodes do not belong to the same connected component, it means that adding the edge will not form a topological cycle, so the edge is added to the empty edge set, and the connected components of the two nodes are merged; if the two nodes belong to the same connected component, the edge is skipped directly to avoid cycle generation. Continue iterative filtering until the number of edges inside the empty edge set meets the requirement. ,in The current total number of control nodes is 8. The minimum number of edges in the spanning tree is 7. All edges within the final empty edge set are defined as mesh fixed constraint edges, which lock the overall mesh topology framework.

[0087] Under the constraint of minimum spanning tree topology, an iterative edge flipping optimization operation is performed for all unconstrained triangles sharing internal edges within the mesh. All unconstrained shared internal edges are traversed; each shared edge is shared by two adjacent triangle elements, forming a quadrilateral topology. All interior angles of the two original triangles before flipping are calculated, and the minimum interior angle of each triangle is extracted and summed to obtain the total minimum interior angle before flipping. Perform edge flipping and topological reconstruction, delete existing shared edges, connect opposite vertices of quadrilaterals to generate new shared edges, and reconstruct two new sets of triangular units; solve for all interior angles of the two new sets of triangles after flipping, extract the minimum interior angle of each triangle, and sum them to obtain the total minimum interior angle of the entire flipped structure. ;like This demonstrates that after flipping, the acute angle of the triangular unit decreases, its shape becomes more regular, the distortion is reduced, and the interpolation accuracy is improved, thus preserving the topological result of this edge flipping; if This proves that flipping cannot optimize cell morphology and may even exacerbate mesh distortion and increase interpolation error. Therefore, the flipping should be immediately canceled to restore the original mesh topology.

[0088] After a full round of traversal and judgment of all unconstrained shared edges, all internal shared edges satisfy the condition that the sum of the minimum interior angles after flipping is less than or equal to the sum of the minimum interior angles before flipping. There is no internal edge that can improve mesh regularity or reduce element distortion by flipping. At this point, all triangular elements have reached the globally optimal form under the current topology, and there is no room for further optimization. The iteration loop is then terminated immediately. The final optimized motion deformation mesh has all triangular elements with uniform shape, no elongated distortion, no acute angle distortion, and a stable topological structure. It also takes into account both global topological constraints and local element smoothness, and can match the composite motion characteristics of low-frequency large-scale rigid deformation and high-frequency local non-rigid micro-deformation of the thoracic cavity.

[0089] Step 400b: Using the motion deformation mesh as the transformation substrate, for each pixel in the current frame, locate the triangular unit in the motion deformation mesh that the corresponding pixel falls into; extract the motion vectors of the three vertices of the triangular unit, and calculate the target position coordinates of the corresponding pixel using the thin-plate spline interpolation function, specifically including:

[0090] A global traversal is performed on all pixels of the current frame image. For any pixel in the original image, a triangular cell containment detection is performed to determine the unique triangular mesh cell to which the pixel belongs, thus achieving the local deformation constraint region localization. Specifically, let the coordinates of the pixel to be detected be... The coordinates of the three vertices of a certain triangular unit are as follows: By constructing three sets of side vectors corresponding to the triangular unit of the pixel, and solving the system of equations, the three sets of centroid weight coefficients of the pixel relative to the triangular unit are obtained. The calculation formula is: ;

[0091] If both conditions are met If the current pixel falls completely inside the triangle unit, the triangle unit is locked as the local deformation constraint unit corresponding to the pixel, and the traversal ends; if the condition is not met, the traversal continues to the next triangle unit until a unique belonging unit is matched, ensuring that each pixel corresponds to only one set of triangle constraint units, with no repetition and no omission.

[0092] After matching the triangular units to which a pixel belongs, the original coordinates and inter-frame motion vectors of the three vertex control nodes of that triangular unit are extracted. The planar coordinates of the eight sets of grid control nodes are obtained by densifying and expanding the nodes in step 401a, and are the fixed spatial coordinates for the current frame. The inter-frame two-dimensional motion vectors of all control nodes are calculated in step 301. A thin-plate spline interpolation model is used to fit the continuous, non-rigid deformation of the thoracic soft tissue without abrupt changes, adapting to the composite motion law of low-frequency overall deformation during respiration and high-frequency micro-deformation during heartbeat. The complete two-way coordinate mapping transformation formula is: ;

[0093] In the formula, The coordinates of the two-dimensional target after pixel deformation correction; is the global linear deformation coefficient, used to fit the overall low-frequency rigid displacement, rotation and scaling deformation of the thoracic cavity; For the first Nonlinear weight coefficients corresponding to each grid control node; This represents the total number of all grid control nodes ultimately determined in step 401a; For the first The original coordinates of each grid control node; For the current pixel and the first The Euclidean distance between the control nodes.

[0094] To uniquely solve for all global linear deformation coefficients and nonlinear weighting coefficients, a complete set of thin-plate spline interpolation constraint equations is constructed based on the original coordinates and the deformed target coordinates of eight sets of control nodes. The complete formula is as follows: ;

[0095] in, , They are the first The horizontal (X-direction) and vertical (Y-direction) nonlinear radial basis function weight coefficients corresponding to each control node are uniquely solved using a matrix solution method. Six global linear deformation coefficients and 16 node nonlinear weight coefficients are fixed in the current frame. The global linear deformation coefficients have globally uniform values ​​to ensure the consistency of low-frequency motion deformation across the entire image, while the nonlinear weight coefficients have node-specific, independent values ​​to ensure the fitting accuracy of local high-frequency micro-deformations. All the solved fixed coefficients are substituted into the thin-plate spline interpolation mapping formula, and the coordinates of each original image pixel in the current frame are calculated point-by-point to obtain the deformation target position coordinates corresponding to each pixel. After the global pixel-by-pixel calculation is completed, a continuous, smooth, one-to-one mapping relationship is established between the original pixel coordinates and the correction target coordinates. The resulting deformation coordinate field has no pixel jumps, no local distortions, and no topological errors, realizing a composite non-rigid deformation mapping that couples low-frequency overall respiratory motion with high-frequency heartbeat micro-motion in dynamic DR chest images.

[0096] Step 401b: Transfer the original grayscale value of the corresponding pixel in the current frame to the pixel position at the target location coordinates; after completing the grayscale value transfer of all pixels in the current frame, detect the pixel positions in the image that have not yet been assigned values, and use inverse distance weighted interpolation to fill the grayscale values ​​of the pixel positions that have not yet been assigned values, to obtain the complete image frame after motion correction, specifically including:

[0097] A global pixel grayscale forward migration operation is performed. The pixel grayscale value at the original coordinates (x, y) of the current frame is retrieved from the corresponding deformation correction target coordinates (x', y') calculated in step 400b. The pixel grayscale value corresponding to the original coordinate position is then mapped one-to-one to the deformed target coordinate position, completing the grayscale migration of pixels in the main image region. This forward mapping process achieves distortion-free deformation correction in the main image region. Due to coordinate remapping caused by pixel deformation offset, some blank pixel positions without assigned values ​​will appear at image edges and grid cell gaps. These blank positions lack effective grayscale information and need to be filled in using an interpolation algorithm.

[0098] For all blank pixels generated after deformation correction, an inverse distance-weighted interpolation algorithm is used to fill the grayscale. This algorithm follows the spatial grayscale correlation law of near-dense and far-sparse areas, adapting to the local grayscale continuity characteristics of the image. Let the coordinates of the blank pixels to be filled be... Select all valid pixels within a preset 8-neighborhood that have completed grayscale assignment. The preset neighborhood is fixed as a 3×3 pixel rectangle centered on the blank pixel. The number of valid neighboring pixels is denoted as . , No. The effective pixel coordinates are Grayscale value First, calculate the Euclidean distance between the pixel to be determined and each of its effective neighboring pixels. The interpolation weights for a single pixel are calculated based on distance. The weight calculation formula is as follows: ;

[0099] In the formula, For the first The interpolation weights for each neighboring effective pixel are calculated as follows: the closer the pixels are spatially, the larger the weight value, and the greater its contribution to the grayscale fitting of blank pixels; the farther the pixels are spatially, the smaller the weight value, which can effectively weaken the grayscale interference from distant irrelevant pixels and ensure the realism and accuracy of local grayscale fitting. The final formula for calculating the interpolated grayscale value of blank pixels is: ;

[0100] In the formula, To obtain the final fitted grayscale value for the blank pixels to be filled, weighted normalization is applied to avoid grayscale deviations caused by differences in the number of neighboring pixels. All blank pixels in the image are traversed one by one, and inverse distance-weighted grayscale filling is performed point by point to eliminate image gaps, black borders, and tortuosity. Finally, a corrected image frame is generated that is globally complete, with continuous grayscale, smooth deformation, and no geometric distortion. This eliminates the composite motion artifacts caused by the low-frequency movement of chest breathing and the high-frequency micro-deformation of the heartbeat in dynamic DR video streams, completing motion correction for the entire frame.

[0101] This embodiment achieves refined construction and optimization of the motion deformation mesh through multiple rounds of node expansion, optimal node selection, minimum spanning tree constraints, and edge flipping optimization. Compared to traditional fixed and sparse mesh structures, the optimized mesh has a uniform node distribution, regular unit structure, and no distorted mesh units. It can fully adapt to the complex non-rigid deformation characteristics of the entire thoracic cavity, overcoming the limitation of traditional local correction meshes that cannot take into account global deformation. Combining the advantages of high-precision deformation mapping of thin-plate spline interpolation and hole filling of inverse distance weighted interpolation, a dual-layer correction mechanism is formed. Thin-plate spline interpolation achieves refined displacement correction for each pixel, and inverse distance weighted interpolation fills in image grayscale holes, solving problems such as image edge breakage, texture misalignment, local artifacts, and pixel holes that are prone to occur in traditional correction algorithms. This effectively improves the continuity and consistency of inter-frame imaging in dynamic DR video streams and optimizes the overall image quality.

[0102] like Figure 2 As shown, embodiments of the present invention also provide a dynamic DR image motion correction system based on video streams, comprising:

[0103] The extraction module is used to acquire a series of consecutive multi-frame image sequences in a dynamic DR video stream, and extract the first morphology kernel, the second morphology kernel and the third morphology kernel for each frame image;

[0104] The update module is used to calculate the spatial correlation coefficient between the first morphology base kernel and the second morphology base kernel based on the difference in local grayscale distribution, and to construct a time-varying baseline connecting the first morphology base kernel and the second morphology base kernel based on the spatial correlation coefficient; when the change in the length of the time-varying baseline in the current frame relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, the first morphology base kernel is moved until the change falls back to within the preset threshold, and the updated first morphology base kernel is obtained.

[0105] The matching module is used to match the local phase consistency distribution of the third morphological base kernel with the updated local phase consistency distribution of the first morphological base kernel and the local phase consistency distribution of the second morphological base kernel, respectively. The one with the larger matching degree is selected as the coupling object. A propagation link is established between the third morphological base kernel and the corresponding coupling object, and the motion vector of the coupling object is transmitted to the third morphological base kernel to obtain the predicted motion vector of the third morphological base kernel.

[0106] The transformation module is used to construct a motion deformation mesh by using the updated first topography kernel, second topography kernel, and third topography kernel carrying the predicted motion vector as control nodes; and to perform pixel-by-pixel transformation on the current frame using the motion deformation mesh to obtain the motion-corrected image frame.

[0107] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0108] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0109] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for motion correction of dynamic DR images based on video streams, characterized in that, The method includes: Step 1: Obtain a sequence of multiple consecutive frames of images from the dynamic DR video stream, and extract the first morphology kernel, the second morphology kernel, and the third morphology kernel for each frame. Step 2: Calculate the spatial correlation coefficient between the first and second morphological base kernels based on the difference in local grayscale distribution. Construct a time-varying baseline connecting the first and second morphological base kernels based on the spatial correlation coefficient. When the change in the length of the time-varying baseline in the current frame relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, move the first morphological base kernel until the change falls back to within the preset threshold to obtain the updated first morphological base kernel. Step 3: Match the local phase consistency distribution of the third morphological base kernel with the updated local phase consistency distribution of the first morphological base kernel and the local phase consistency distribution of the second morphological base kernel respectively. Select the one with the larger matching degree as the coupling object. Establish a propagation link between the third morphological base kernel and the corresponding coupling object. Transmit the motion vector of the coupling object to the third morphological base kernel to obtain the predicted motion vector of the third morphological base kernel. Step 4: Use the updated first topology kernel, second topology kernel, and third topology kernel carrying the predicted motion vector as control nodes to construct a motion deformation mesh; use the motion deformation mesh to perform pixel-by-pixel transformation on the current frame to obtain the motion-corrected image frame.

2. The dynamic DR image motion correction method based on video stream according to claim 1, characterized in that, The first morphological base is located on the gradient ridge line at the junction of the top of the left diaphragm and the apex of the heart. The second morphological base is located on the gray extreme value band at the overlap of the pericardial fat pad and the diaphragm in the right cardiophrenic angle region. The third morphological base is located at the phase-consistent peak point at the intersection of the lower edge of the aortic arch and the left main bronchus.

3. The dynamic DR image motion correction method based on video stream according to claim 2, characterized in that, Based on the difference in local grayscale distribution between the first and second morphological base kernels, the spatial correlation coefficient between them is calculated, and a time-varying baseline connecting the first and second morphological base kernels is constructed based on the spatial correlation coefficient, including: Extract the gray values ​​of all pixels within the local window surrounding the first morphology base kernel to form a first gray-level distribution sequence, and extract the gray values ​​of all pixels within the corresponding window surrounding the second morphology base kernel to form a second gray-level distribution sequence. Calculate the cumulative distribution function of the first gray-level distribution sequence and the second gray-level distribution sequence respectively, and calculate the L1 norm distance between the two cumulative distribution functions; Calculate the absolute cosine value between the principal direction angle of the local gradient within the window containing the first morphological base kernel and the principal direction angle of the local gradient within the window containing the second morphological base kernel; multiply the L1 norm distance by the absolute cosine value, take the reciprocal, and normalize the result to obtain the spatial correlation coefficient; Starting from the coordinates of the first morphological base kernel on the image plane and ending at the coordinates of the second morphological base kernel, a directed line segment is drawn as the initial representation of the time-varying baseline, and the spatial correlation coefficient is assigned as the weight of the directed line segment.

4. The dynamic DR image motion correction method based on video stream according to claim 3, characterized in that, When the change in the length of the time-varying baseline in the current frame relative to the corresponding baseline length in the previous frame exceeds a preset threshold, the first topography kernel is moved until the change falls back below the preset threshold, resulting in an updated first topography kernel, including: Calculate the ratio of the Euclidean length of the time-varying baseline in the current frame to the Euclidean length of the time-varying baseline in the previous frame. If the ratio is greater than a preset upper threshold or less than a preset lower threshold, it is determined to exceed the preset threshold. Using the current position of the first topographic base core as the center and the straight line perpendicular to the direction of the line connecting the second topographic base core on the corresponding time-varying baseline as a reference, move the first topographic base core gradually along the normal direction of the corresponding reference line with a step size of the current frame's time-varying baseline length multiplied by a preset scaling factor, where the preset scaling factor is less than 1. After each move, recalculate the ratio until the ratio falls within a preset tolerance range, stop moving, and record the position of the first topographic base core at this time as the updated first topographic base core.

5. The dynamic DR image motion correction method based on video stream according to claim 4, characterized in that, Step 3 includes: Calculate the peak values ​​of the cross-correlation surfaces of the phase consistency direction spectra between the third morphological base and the updated first morphological base, and between the third morphological base and the second morphological base, respectively, and take the corresponding peak values ​​of the cross-correlation surfaces as the first matching degree and the second matching degree, respectively. Compare the first matching degree with the second matching degree, and select the morphological base kernel corresponding to the larger value as the coupling object; construct a straight-line propagation link with the third morphological base kernel and the coupling object as the two endpoints, and obtain the motion vector of any point on the straight-line propagation link by linear interpolation of the motion vectors of the two endpoints; The motion vector of the coupled object between the current frame and the previous frame is propagated along a linear propagation link to the endpoint where the third topography base kernel is located, thus obtaining the predicted motion vector of the third topography base kernel.

6. The dynamic DR image motion correction method based on video stream according to claim 5, characterized in that, Using the updated first topographic kernel, second topographic kernel, and third topographic kernel carrying the predicted motion vector as control nodes, a motion deformation mesh is constructed, including: Using the updated first topographic kernel, second topographic kernel, and third topographic kernel carrying the predicted motion vector as the initial control node set, four auxiliary nodes are inserted at equal intervals on the boundary of the planar region covered by the initial control node set to form a planar point set containing seven nodes. Traverse any three points in the planar point set to form a triangle combination, calculate the circumcircle of each triangle, and take the circumcircle of the planar point set that does not contain any nodes in the planar point set other than the three vertices of the corresponding triangle as candidate circles. Select the candidate circle with the largest radius, and add the center coordinates of the circle with the largest radius as a new control node to the planar point set. Using all control nodes as vertices, a motion deformation mesh is obtained.

7. The dynamic DR image motion correction method based on video stream according to claim 6, characterized in that, Using all control nodes as vertices, the motion deformation mesh is obtained, including: Obtain all control nodes of the motion deformation mesh and the edges of all triangular elements in the motion deformation mesh, and form an initial planar graph from all control nodes and edges; Calculate the Euclidean distance between the two endpoints of each edge in the initial planar graph as the weight of the corresponding edge; sort all edges in the initial planar graph in ascending order of weight to establish an empty edge set; starting from the edge with the smallest weight, check whether the current edge forms a cycle with any existing edge in the empty edge set; if the two control nodes connected by the current edge are not in the same connected component in the current connected component of the empty edge set, add the current edge to the empty edge set; otherwise, skip the current edge; repeat this process until the number of edges in the empty edge set is the total number of control nodes minus 1, thus obtaining the minimum spanning tree; Each edge in the minimum spanning tree is used as a constraint edge. An edge flipping operation is performed on the triangular elements in the motion deformable mesh. That is, for each internal edge in the motion deformable mesh that is shared by two triangles, the minimum interior angle of the two triangles before flipping and the minimum interior angle of the two triangles after flipping are calculated. If the minimum interior angle after flipping is greater than the minimum interior angle before flipping, then the flipping is performed. All shared internal edges are traversed until there are no more edges that can be flipped, and the optimized motion deformable mesh is obtained.

8. The dynamic DR image motion correction method based on video stream according to claim 7, characterized in that, The current frame is transformed pixel-by-pixel using a motion-deformed mesh to obtain a motion-corrected image frame, including: Using a motion deformation mesh as the transformation substrate, for each pixel in the current frame, the triangular unit in the motion deformation mesh into which the corresponding pixel falls is located; the motion vectors of the three vertices of the triangular unit are extracted, and the target position coordinates of the corresponding pixel are calculated by using a thin plate spline interpolation function; The original grayscale value of the corresponding pixel in the current frame is transferred to the pixel position at the target position coordinates. After the grayscale value transfer of all pixels in the current frame is completed, the pixel positions in the image that have not yet been assigned a value are detected. The unassigned pixel positions are filled with grayscale using inverse distance weighted interpolation to obtain the complete image frame after motion correction.

9. A dynamic DR image motion correction system based on video stream, wherein the system implements the method as described in any one of claims 1 to 8, characterized in that, include: The extraction module is used to acquire a series of consecutive multi-frame image sequences in a dynamic DR video stream, and extract the first morphology kernel, the second morphology kernel and the third morphology kernel for each frame image; The update module is used to calculate the spatial correlation coefficient between the first morphology base kernel and the second morphology base kernel based on the difference in local grayscale distribution, and to construct a time-varying baseline connecting the first morphology base kernel and the second morphology base kernel based on the spatial correlation coefficient; when the change in the length of the time-varying baseline in the current frame relative to the length of the corresponding baseline in the previous frame exceeds a preset threshold, the first morphology base kernel is moved until the change falls back to within the preset threshold, and the updated first morphology base kernel is obtained. The matching module is used to match the local phase consistency distribution of the third morphological base kernel with the updated local phase consistency distribution of the first morphological base kernel and the local phase consistency distribution of the second morphological base kernel, respectively. The one with the larger matching degree is selected as the coupling object. A propagation link is established between the third morphological base kernel and the corresponding coupling object, and the motion vector of the coupling object is transmitted to the third morphological base kernel to obtain the predicted motion vector of the third morphological base kernel. The transformation module is used to construct a motion deformation mesh by using the updated first topography kernel, second topography kernel, and third topography kernel carrying the predicted motion vector as control nodes; and to perform pixel-by-pixel transformation on the current frame using the motion deformation mesh to obtain the motion-corrected image frame.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 8.

Images