Method and system for compressed sensing reconstruction of medical images

Abstract

The application provides a compressed sensing reconstruction method and system of medical images, and relates to the technical field of magnetic resonance imaging.The method comprises the following steps: acquiring an undersampled frequency domain signal collected by a radio frequency coil of a medical imaging device, extracting physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal; performing non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and the sampling density distribution, to generate a complex image space pixel matrix; calling a preset wavelet transform base to perform multi-layer coefficient decomposition on the complex image space pixel matrix, to acquire a high-frequency and low-frequency sparse transform coefficient set containing high-frequency details and low-frequency backgrounds.The application realizes the unity of tissue edge sharpness preservation and global phase error compensation while suppressing undersampling artifacts.

Description

Technical Field

[0001] This invention relates to the field of magnetic resonance imaging technology, and in particular to a method and system for compressed sensing reconstruction of medical images. Background Technology

[0002] Compressed sensing theory utilizes the sparsity of images in a specific transform domain to reconstruct spatial pixel information of images from sampled data below the requirements of the Nyquist sampling theorem. It acquires undersampled frequency domain signals through radio frequency coils, combines the sparse priors of images in the wavelet domain, spatial difference domain, or total variation domain, and uses nonlinear reconstruction algorithms to restore the image, thereby shortening the scanning time.

[0003] Existing wavelet domain sparse constraint-based reconstruction methods decompose complex images into low-frequency approximation coefficients and high-frequency detail coefficients under a discrete wavelet basis. Dealiasing reconstruction is achieved by truncating or weighting the amplitude of the sparse coefficients. However, such methods apply a uniform constraint strength to the entire image. Under low signal-to-noise ratio sampling conditions, this uniform processing method is difficult to accurately quantify the true gray-level drop between adjacent pixels, resulting in a decrease in the sharpness of the true tissue boundaries in the reconstruction results. Some aliasing artifacts are preserved as edge structures, thus producing Gibbs ring artifacts. At the same time, the undersampling process and the physical field inhomogeneity of the RF coil introduce spatially uneven phase errors. Existing methods mostly perform phase correction in local areas and lack continuous extrapolation and summarization of the phase shift trend of pixels across the entire domain. This makes it difficult to accurately evaluate the global phase error, resulting in phase cancellation during subsequent signal superposition, causing local signal loss or uneven band artifacts in the image. Summary of the Invention

[0004] This invention provides a compressed sensing reconstruction method and system for medical images, which achieves a balance between preserving tissue edge sharpness and compensating for global phase error while suppressing undersampling artifacts.

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

[0006] Firstly, a compressed sensing reconstruction method for medical images, the method comprising:

[0007] Step 1: Acquire the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal;

[0008] Step 2: Perform non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image spatial pixel matrix;

[0009] Step 3: Call the preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix to obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background.

[0010] Step 4: Based on the complex image spatial pixel matrix, calculate the grayscale spatial change rate between adjacent pixels point by point to generate an edge response distribution map that reflects the sharpness of tissue edges.

[0011] Step 5: Based on the complex image spatial pixel matrix, calculate the cumulative phase deviation of each pixel across the entire domain to generate a global phase compensation field that reflects the global error distribution.

[0012] Step 6: Perform geodesic distance-based regional expansion processing on the edge response distribution map to generate an expanded edge protection area;

[0013] Step 7: Perform spatial domain point-by-point weighted fusion calculation between the extended edge protection region and the global phase compensation field to generate a multi-constraint reconstruction parameter set;

[0014] Step 8: Based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis, generate a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients. Then, perform a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets according to the weighted coefficient matrix to obtain the target medical image.

[0015] Secondly, the compressed sensing reconstruction system for medical images includes:

[0016] The acquisition module is used to acquire the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal.

[0017] The conversion module is used to perform non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image space pixel matrix;

[0018] The decomposition module is used to call a preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix, and obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background.

[0019] The calculation module is used to calculate the gray-scale spatial change rate between adjacent pixels point by point based on the complex image spatial pixel matrix, and generate an edge response distribution map that reflects the sharpness of tissue edges; based on the complex image spatial pixel matrix, it calculates the cumulative phase deviation of each pixel point across the entire domain, and generates a global phase compensation field that reflects the global error distribution.

[0020] The reconstruction module is used to perform geodesic distance-based regional expansion processing on the edge response distribution map to generate an expanded edge protection area; and to perform spatial domain point-by-point weighted fusion operation on the expanded edge protection area and the global phase compensation field to generate a multi-constraint reconstruction parameter set.

[0021] The matching module is used to generate a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis. Then, it performs a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets based on the weighted coefficient matrix to obtain the target medical image.

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

[0023] By employing non-uniform Fourier transform to generate complex images, the true frequency distribution of undersampling is preserved, overcoming the problem that local correction cannot assess global errors. Morphological dilation is used to expand the edge protection region, ensuring smooth transitions in adjacent edge regions. Edge protection and phase compensation are fused point-by-point in the spatial domain and weighted coefficients are generated by wavelet subband matching. Weighted inverse transform is performed on high and low frequency coefficients, enabling the two types of constraints to work synergistically in the reconstruction process. This achieves joint optimization of edge sharpness preservation and global phase error compensation, effectively suppressing Gibbs ring artifacts and banded brightness unevenness, and improving the reconstruction fidelity of medical images under undersampling conditions. Attached Figure Description

[0024] Figure 1 This is a flowchart illustrating the compressed sensing reconstruction method for medical images provided in an embodiment of the present invention.

[0025] Figure 2 This is a schematic diagram of a compressed sensing reconstruction system for medical images provided in an embodiment of the present invention.

[0026] Figure 3 This is a schematic diagram of the process provided by an embodiment of the present invention, which generates a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis, and performs a weighted inverse transformation on the high-frequency and low-frequency sparse transform coefficient sets based on the weighted coefficient matrix to obtain the target medical image. Detailed Implementation

[0027] 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.

[0028] like Figure 1 As shown, an embodiment of the present invention proposes a compressed sensing reconstruction method for medical images, the method comprising the following steps:

[0029] Step 1: Acquire the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal;

[0030] Step 2: Perform non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image spatial pixel matrix;

[0031] Step 3: Call the preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix to obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background.

[0032] Step 4: Based on the complex image spatial pixel matrix, calculate the grayscale spatial change rate between adjacent pixels point by point to generate an edge response distribution map that reflects the sharpness of tissue edges.

[0033] Step 5: Based on the complex image spatial pixel matrix, calculate the cumulative phase deviation of each pixel across the entire domain to generate a global phase compensation field that reflects the global error distribution.

[0034] Step 6: Perform geodesic distance-based regional expansion processing on the edge response distribution map to generate an expanded edge protection area;

[0035] Step 7: Perform spatial domain point-by-point weighted fusion calculation between the extended edge protection region and the global phase compensation field to generate a multi-constraint reconstruction parameter set;

[0036] Step 8: Based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis, generate a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients. Then, perform a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets according to the weighted coefficient matrix to obtain the target medical image.

[0037] In this embodiment of the invention, by using non-uniform Fourier transform to generate complex images, the true frequency distribution of undersampling is preserved, overcoming the problem that local correction cannot assess global errors. Morphological dilation is used to expand the edge protection region, ensuring smooth transition of the adjacent edge region. After fusing edge protection and phase compensation point by point in the spatial domain, weighted coefficients are generated by wavelet subband matching. Weighted inverse transformation is performed on high and low frequency coefficients, so that the two types of constraints work together in the reconstruction process. This achieves joint optimization of edge sharpness preservation and global phase error compensation, effectively suppressing Gibbs ring artifacts and banded brightness unevenness, and improving the reconstruction fidelity of medical images under undersampling conditions.

[0038] In a preferred embodiment of the present invention, step 1 above, acquiring the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging device, and extracting the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal, includes:

[0039] Step 11: Acquire the raw undersampled frequency domain signal received by the radio frequency coil of the medical imaging equipment, analyze the radio frequency pulse timing control information carried by the raw undersampled frequency domain signal to obtain the acquisition sequence containing spatial coordinates and timestamps, specifically including:

[0040] The raw undersampled frequency domain signal received by the multi-channel radio frequency coil of the magnetic resonance imaging device in one scan cycle is acquired. The radio frequency pulse timing control information of this scan is extracted from the pulse sequence log file of the device console, including the emission time of the radio frequency excitation pulse, the start and end time of the application of the frequency coding gradient, the step number of the phase coding gradient, and the opening time and duration of the signal readout window.

[0041] Using the emission time of the radio frequency excitation pulse as the time origin, the signal readout time of each readout sampling point is offset relative to this origin to obtain the relative timestamp of that sampling point. Based on the above relative timestamp, the spatial position of the sampling point is correlated with the frequency encoding gradient and phase encoding gradient during the readout period. For the readout direction, the resonant frequency offset generated by the frequency encoding gradient in the spatial position is linearly related to the spatial coordinates. The spatial coordinates of the sampling point in the readout direction can be calculated from the measured frequency offset using the hydrogen proton gyromagnetic ratio. For the phase encoding direction, the phase encoding gradient is applied successively with a fixed amplitude, generating a cumulative phase proportional to the spatial coordinates in that direction within a fixed duration. The spatial coordinates of the sampling point in the phase encoding direction are determined accordingly.

[0042] The spatial coordinates of the two directions mentioned above are bound to the corresponding relative timestamps as a collection record item. All records are arranged in ascending order of timestamps to construct a collection sequence containing spatial coordinates and timestamps. Each element in the sequence records the time of occurrence of a sampling event and its spatial sampling position.

[0043] Step 12: Based on the spatial coordinates and timestamps contained in the acquisition sequence, arrange the spatial coordinates of each sampling point in ascending order of timestamps, calculate the spatial displacement vector between sampling points corresponding to adjacent timestamps, and obtain the physical acquisition trajectory parameters, specifically including:

[0044] The generated acquisition sequence is arranged in ascending order of timestamps. Along the sequence direction, the spatial coordinates of two adjacent records are extracted one by one. The coordinate differences between the two coordinates in the readout direction and the phase encoding direction are calculated respectively and combined into a two-dimensional spatial displacement vector. The magnitude of this vector is equal to the square root of the sum of the squares of the differences in the two directions, which represents the spatial distance between adjacent sampling points. The direction angle of this vector is equal to the arctangent value corresponding to the ratio of the differences in the two directions, which represents the direction of movement. This spatial displacement vector describes the magnitude and direction of movement of the sampling position of the RF coil in the imaging plane between two adjacent sampling times.

[0045] The above calculation is performed on all adjacent records in the sequence one by one, and each spatial displacement vector is bound to its start and end timestamps. The parameters are arranged in chronological order to form physical acquisition trajectory parameters, which completely describe the motion path of the RF coil sampling position along time during a single scan.

[0046] Step 13: Based on the spatial displacement vector contained in the physical acquisition trajectory parameters, divide the target imaging space into multiple uniform grids, calculate the number of sampling points falling within each uniform grid, and obtain the initial sampling density distribution, specifically including:

[0047] Based on the generated physical acquisition trajectory parameters, the spatial coordinates of all sampling points are extracted to form a complete set of spatial coordinates. The field of view of the target imaging space, the total length of the readout direction and the total length of the phase encoding direction are divided into several columns in the readout direction and several rows in the phase encoding direction according to the preset grid division density, forming a regular grid array in the imaging plane. The side length of the readout direction of each grid unit is equal to the total length of the readout direction divided by the number of columns, and the side length of the phase encoding direction is equal to the total length of the phase encoding direction divided by the number of rows.

[0048] For each grid cell in the above-mentioned regular grid array, it is sequentially determined whether the coordinates of each sampling point in the complete set of spatial coordinates fall within the boundary range of the grid cell: that is, the readout coordinates of the sampling point are not less than the readout start coordinates and less than the readout end coordinates of the grid, and the phase encoding coordinates are not less than the phase encoding start coordinates and less than the phase encoding end coordinates. The sampling points that meet the conditions are counted grid by grid to obtain the number of sampling points falling into each grid cell. The grid cells and their corresponding sampling point counts are arranged into a two-dimensional counting matrix according to the row and column positions, which is the initial sampling density distribution. The larger the value of each element in the matrix, the higher the frequency of the grid area being crossed or covered by the RF coil sampling trajectory.

[0049] Step 14: Based on the number of sampling points included in the initial sampling density distribution, calculate the proportion of the number of sampling points in each uniform grid to the total number of global sampling points to obtain the spatial sampling density distribution, specifically including:

[0050] Based on the generated initial sampling density distribution, the counts of all grid sampling points are summed row by row and column by column to obtain the total number of global sampling points for this scan. The proportion of the sampling point count of each grid cell to the total number of global sampling points is calculated for each grid cell. This proportion is a dimensionless value between 0 and 1; a value of zero indicates that the grid area is not covered by any sampling trajectory, while a larger value indicates a denser concentration of sampling points in the grid area. The proportions of each grid cell are arranged in rows and columns into a matrix of the same size as the initial sampling density distribution, which is the spatial sampling density distribution.

[0051] In this embodiment of the invention, the acquisition sequence containing spatial coordinates and timestamps is obtained by parsing the radio frequency pulse timing control information, which overcomes the technical problem of the traditional method that assumes the sampling position as a uniform grid point and ignores the actual trajectory deviation, resulting in a decrease in the accuracy of the subsequent Fourier transform.

[0052] In a preferred embodiment of the present invention, step 2 above, which involves performing a non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image spatial pixel matrix, includes:

[0053] Step 21: Based on the spatial frequency correspondence between the spatial coordinates contained in the physical acquisition trajectory parameters and the pixel coordinates of the target imaging space, calculate the spatial phase rotation factor corresponding to each sampling point to obtain the set of phase rotation factors, specifically including:

[0054] Based on the physical acquisition trajectory parameters, the k-space coordinates of each sampling point in the target imaging space are extracted. Let the k-th point be... The spatial coordinates of the k sampling points are: ,in , This represents the total number of sampling points in this scan.

[0055] Determine the size of the complex image space pixel matrix generated after the non-uniform Fourier transform of the target imaging space, i.e., the pixel resolution of the reconstructed image in the readout direction and the phase coding direction. Let the number of pixels in the readout direction be... The number of pixels in the phase encoding direction is Then the first in the image space Line number The spatial coordinates of the column pixels are: ;

[0056] in and These are the field-of-view lengths for the readout direction and the phase encoding direction, respectively. , .

[0057] For the k spatial sampling points With the first in image space Line number Column pixels Based on the kernel function relationship of the non-uniform Fourier transform, the spatial phase rotation factor of the sampling point with respect to that pixel is calculated. : ;

[0058] in The imaginary unit ( ), It is a complex number with a modulus of 1 and a phase angle of . This characterizes the phase deflection amount contributed by the frequency component carried by the sampling point in k-space to the signal at that pixel location in the image space.

[0059] The first The spatial phase rotation factors corresponding to each sampling point and all pixels in the image space are arranged according to their pixel row and column positions. The complex matrix is ​​bound to the k-space coordinates of the sampling point and stored. All of them are traversed. After sampling points, a set of phase rotation factors is formed. : ;

[0060] The set is organized by sampling point index, with each sampling point entry containing its k-space coordinates and a phase rotation factor matrix of the same size as the target reconstructed image.

[0061] Step 22: Based on the set of phase rotation factors and the undersampled frequency domain signal, perform point-by-point complex multiplication to obtain the modulation frequency domain signal corresponding to each sampling point, specifically including:

[0062] Based on the acquired undersampled frequency domain signal, extract the first... The original complex values ​​corresponding to each sampling point : ;

[0063] in For the real channel component, For the imaginary channel component, The imaginary unit ( ), which respectively correspond to the free induction attenuation signal components of the RF coil induced magnetization vector in two orthogonal directions.

[0064] Based on the generated set of phase rotation factors Take out the first Phase rotation factor matrix corresponding to each sampling point For each pixel position in the image space , to the original complex value Phase rotation factor corresponding to this position Performing complex multiplication yields the first... The modulation frequency domain signal components of each sampling point at the pixel location : ;

[0065] The real and imaginary parts of the complex multiplication are calculated using the following formulas: ; ;

[0066] by amplitude As a scaling factor, phase angle and The sum of the phase angles is the new phase angle, generating a complex value with the amplitude modulated and the phase rotated.

[0067] Traverse all For each of the sampling points, perform the pixel-by-pixel complex multiplication operation described above to obtain the modulation frequency domain signal matrix corresponding to each sampling point. : ;

[0068] Each matrix The size is the same as the target reconstructed image, and the matrix is ​​the first... Line number Column elements The magnitude reflects the first Each sampling point corresponds to the pixel position in the image space. The signal contribution intensity and its phase reflect the spatial phase propagation delay from the sampling point to the pixel location.

[0069] Step 23: Based on the modulation frequency domain signal and spatial sampling density distribution, perform global signal accumulation calculation along the coordinate direction of the target imaging space to obtain a complex image spatial pixel matrix containing amplitude and phase components, specifically including:

[0070] Based on the generated spatial sampling density distribution, the density weight value corresponding to each pixel position in the target imaging space is obtained, and the density weight matrix is ​​set as follows. for: ;

[0071] in pixel position in image space The normalized sampling density value after bilinear interpolation resampling to the image pixel resolution is used to... The reciprocal of the value is used as a density compensation factor at that pixel location to balance the non-uniform density deviation caused by undersampling.

[0072] Based on the generated modulation frequency domain signal matrix corresponding to each sampling point For each pixel position in the image space , will all The contribution components of the modulated frequency domain signal at each sampling point to that pixel location are accumulated and multiplied by the corresponding density compensation factor to obtain the complex reconstruction result for that pixel location. : ;

[0073] The real part of the above accumulation operation With the imaginary part Write them out separately: ; ;

[0074] The amplitude component of the pixel location reconstruction result and phase components They are respectively: ; ;

[0075] For all pixel positions in the target imaging space ( , Perform the above cumulative calculation and density compensation pixel by pixel, and convert the complex values ​​at each pixel position... Arranged by row and column position The complex matrix is ​​used to obtain the complex image space pixel matrix containing the magnitude and phase components. : ;

[0076] This matrix, in complex form, represents the magnetization vector distribution at each pixel location reconstructed in the image space after density compensation and phase rotation transformation from a non-uniformly sampled k-space frequency domain signal.

[0077] In this embodiment of the invention, by using the spatial frequency correspondence between physical acquisition trajectory parameters and image spatial pixel coordinates to calculate the spatial phase rotation factor of each sampling point, the interpolation error problem introduced by the traditional method of forcibly mapping the sampling points to a uniform grid is overcome. This achieves the technical effect of completing the non-uniform Fourier transformation from the frequency domain to the image domain based on the non-uniform sampling trajectory and generating a complex image spatial pixel matrix containing complete amplitude and phase information.

[0078] In a preferred embodiment of the present invention, step 3 above, which involves calling a preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix to obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background, includes:

[0079] Step 31: Based on the complex image spatial pixel matrix, extract the amplitude components of each pixel to construct a real-valued image matrix. Then, call a preset discrete wavelet transform basis to perform spatial domain low-pass and high-pass filtering operations on the real-valued image matrix to separate the first low-frequency approximate sub-band and the first high-frequency detail sub-band. Specifically, this includes:

[0080] Based on complex image space pixel matrix Each element is , For amplitude components, For the phase component, the amplitude component is extracted pixel by pixel: ;

[0081] in , , will all Arranged by row and column position, constructing... Real-valued image matrix of the same size Each element in the matrix represents the signal amplitude at the corresponding pixel location, reflecting the magnetization vector intensity generated by the tissue at that location after radio frequency excitation.

[0082] The preset discrete wavelet transform basis is selected from the Daubechies-4 wavelet (db4), whose low-pass decomposition filter coefficients The possible values ​​are: ;

[0083] Corresponding high-pass decomposition filter coefficients Determined by orthogonal mirror relationship: , , , For real-valued image matrices Perform filtering and downsampling sequentially in the row and column directions:

[0084] Using low-pass filter coefficients along Each row slides column by column from left to right, the filter window covering the amplitude of 4 consecutive pixels. The weighted sum of the filter coefficients and their corresponding amplitudes is calculated. Every other pixel (with a step size of 2), a row of low-pass filtered results is output, using the high-pass filter coefficients... Slide along each row in the same manner, with a step size of 2, and output the row high-pass filter result. After this step... The matrix is ​​compressed into two intermediate matrices with half the width in the row direction: a row low-pass matrix and a row high-pass matrix. For each column of the row low-pass matrix, the low-pass filter coefficients and high-pass filter coefficients are slid along the column direction with a step size of 2 to obtain the first low-frequency approximate sub-band and the first horizontal high-frequency detail sub-band. For each column of the row high-pass matrix, the low-pass filter coefficients and high-pass filter coefficients are slid along the column direction with a step size of 2 to obtain the first vertical high-frequency detail sub-band and the first diagonal high-frequency detail sub-band.

[0085] The matrix size of the above four sub-bands is all The first low-frequency approximation sub-band retains the main structural information of the original image after low-pass filtering and 1 / 2 downsampling; the first horizontal, vertical and diagonal high-frequency detail sub-bands capture the gray-level change information of the original image in the direction of horizontal edge, vertical edge and diagonal edge respectively. The three are collectively called the first high-frequency detail sub-band. The preset decomposition layer is 2 layers. This step completes the first layer decomposition.

[0086] Step 32: Based on the first low-frequency approximate sub-band, a spatial domain multi-scale filtering operation is performed using a preset discrete wavelet transform basis to separate the second low-frequency approximate sub-band and the second high-frequency detail sub-band, specifically including:

[0087] Based on the generated first low-frequency approximate sub-band, its matrix size is This is the low-frequency profile obtained after the original image has undergone the first layer of low-pass filtering and 1 / 2 downsampling. The spatial resolution has been reduced to one-quarter of the original image. According to the preset 2-layer decomposition scheme, the Daubechies-4 low-pass filter coefficients, which are exactly the same as those in step 31, are used again. and high-pass filter coefficients Repeatedly perform row-direction and column-direction filtering and downsampling with a step size of 2 on the first low-frequency approximate subband:

[0088] First, low-pass and high-pass filtering and downsampling are performed along the row direction, and then low-pass and high-pass filtering and downsampling are performed along the column direction, resulting in four sub-bands. In the decomposition result, the second low-frequency approximate sub-band is obtained by low-pass filtering along both the row and column directions, and the matrix size is further reduced to... This preserves the overall brightness distribution and global structural contours of the original image at a coarser scale. The second horizontal high-frequency detail sub-band, the second vertical high-frequency detail sub-band, and the second diagonal high-frequency detail sub-band, obtained by combining low-pass / high-pass filtering along the row and column directions respectively, are collectively referred to as the second high-frequency detail sub-band, with a matrix size of [missing information]. It captures the edge directionality information of the first low-frequency approximate subband at a coarser scale.

[0089] Step 33: Based on the first high-frequency detail sub-band, the second high-frequency detail sub-band, and the second low-frequency approximation sub-band, align them according to their spatial coordinate positions and arrange them hierarchically to obtain a sparse transform coefficient set containing high-frequency detail coefficients and low-frequency background coefficients, specifically including:

[0090] Based on the generated first high-frequency detail subband, second high-frequency detail subband, and second low-frequency approximate subband, they are arranged according to the hierarchical nested structure of wavelet decomposition:

[0091] The second low-frequency approximation subband is placed in the upper left corner, with the second horizontal high-frequency detail subband and the second diagonal high-frequency detail subband placed adjacent to it on its right. Below it, the second vertical high-frequency detail subband and the three directional subbands of the first high-frequency detail subband are placed adjacent to it. When arranged, the relative spatial position of the coefficient elements in each subband corresponds to their spatial coordinates in the original image. The element in the upper left corner of the matrix corresponds to the upper left region of the image, and the element in the lower right corner corresponds to the lower right region of the image. All the coefficients after this nested arrangement constitute a sparse transform coefficient set. The low-frequency background coefficients contained in the second low-frequency approximation subband correspond to the global brightness distribution and overall anatomical structure contour of the original image. The values ​​are concentrated near zero, and a few large-value coefficients carry the main energy of the image. The high-frequency detail coefficients contained in the six high-frequency detail subbands correspond to the edge and texture information of the original image along the horizontal, vertical and diagonal directions at different spatial scales. Most of the coefficients are close to zero, and only at tissue boundaries and structural abrupt changes do obvious non-zero values ​​appear.

[0092] In a preferred embodiment of the present invention, step 4 above, which involves calculating the grayscale spatial change rate between adjacent pixels point by point based on the complex image spatial pixel matrix to generate an edge response distribution map reflecting the sharpness of tissue edges, includes:

[0093] Step 41: Based on the complex image spatial pixel matrix, extract the grayscale values ​​of the amplitude components corresponding to each pixel to construct a two-dimensional grayscale distribution matrix, specifically including:

[0094] Based on a complex image spatial pixel matrix, each pixel location stores a complex value containing real and imaginary parts. The real part corresponds to the magnetization vector component intensity received by the radio frequency coil in one of the orthogonal directions, and the imaginary part corresponds to the magnetization vector component intensity in the other orthogonal direction. The square root of the sum of the squares of the real and imaginary parts is taken pixel by pixel to obtain the grayscale value of the amplitude component corresponding to that pixel location. This value reflects the total intensity of the magnetization vector generated by the human tissue at that pixel location after radio frequency excitation, and is positively correlated with the proton density and relaxation characteristics of the tissue.

[0095] The grayscale values ​​of the amplitude components at each pixel location are arranged in a one-to-one correspondence with their row and column positions in the pixel matrix of the complex image space, constructing a two-dimensional grayscale distribution matrix of the same size as the reconstructed image. This matrix records the tissue signal intensity at each pixel location in the complex image space in real numerical form.

[0096] Step 42: Based on the two-dimensional grayscale distribution matrix, extract the grayscale values ​​of adjacent pixels along the horizontal and vertical coordinate directions respectively, and calculate the difference to obtain the set of grayscale differences in the horizontal direction and the set of grayscale differences in the vertical direction, specifically including:

[0097] Based on the generated two-dimensional grayscale distribution matrix, where each element is a scalar grayscale value, the grayscale difference between adjacent pixels is extracted along both spatial coordinate directions. The processing along the horizontal coordinate direction is as follows:

[0098] Keeping the row numbers unchanged, starting from the first pixel column of the row, extract the grayscale values ​​of that pixel column and its immediately right-hand neighboring pixel column sequentially. Extract a pixel pair consisting of two adjacent pixels in each row. For each pixel pair, subtract the grayscale value of the left pixel from the grayscale value of the right pixel to obtain the horizontal grayscale difference corresponding to that pixel position. A positive value indicates an increase in grayscale (signal brightens) from left to right horizontally, while a negative value indicates a decrease in grayscale (signal darkens) from left to right. The absolute value reflects the steepness of the grayscale drop or rise at that position horizontally. Arrange the horizontal grayscale differences of all pixel positions according to row and column positions to form a set of horizontal grayscale differences. The processing along the vertical coordinate direction is as follows:

[0099] Keeping the column numbers unchanged, starting from the first pixel row of the column, extract the grayscale values ​​of that pixel row and its immediately adjacent pixel row below it sequentially. Extract two pixels with adjacent row numbers in each column to form a pixel pair. For each pixel pair, subtract the grayscale value of the pixel above from the grayscale value of the pixel below it to obtain the vertical grayscale difference corresponding to that pixel position. A positive value of this difference indicates an increase in grayscale from top to bottom along the vertical direction, while a negative value indicates a decrease in grayscale from top to bottom. The magnitude of the absolute value reflects the steepness of the grayscale change along the vertical direction at that position. Arrange the vertical grayscale differences of all pixel positions according to their row and column positions to form a set of vertical grayscale differences. The aforementioned sets of horizontal and vertical grayscale differences together describe the grayscale spatial distribution of the two-dimensional grayscale distribution matrix in two orthogonal directions within the imaging plane.

[0100] Step 43: Based on the set of gray-level differences in the horizontal direction and the set of gray-level differences in the vertical direction, calculate the energy aggregation value of the gray-level differences in the horizontal and vertical directions for each pixel, and obtain the spatial variation intensity value of each pixel. Specifically, this includes:

[0101] Based on the horizontal and vertical gray-level difference sets, both sets have the same matrix size, and elements at the same row and column positions correspond to the gray-level differences of the same pixel position in the horizontal and vertical directions, respectively. For each pixel position, its horizontal and vertical gray-level differences are extracted. The squares of the horizontal and vertical gray-level differences are then summed and the square root is taken to obtain the spatial change intensity value of that pixel position. The gray-level change components in the two orthogonal directions are combined into a comprehensive change intensity index, regardless of whether the gray-level abrupt change occurs in the horizontal, vertical, or vertical direction. The intensity value can be sensitively reflected in the oblique direction. In areas with uniform tissue, the gray-level difference in both the horizontal and vertical directions is close to zero, and the spatial change intensity value obtained by aggregation is also close to zero. At the tissue boundary (such as the boundary between gray matter and white matter, or the boundary between lesion and normal tissue), the gray-level difference in at least one direction changes dramatically, and the spatial change intensity value obtained by aggregation also increases accordingly. This effectively distinguishes tissue edges from uniform areas. The spatial change intensity value of each pixel position is stored pixel by pixel, and its value range is a non-negative real number. The larger the value, the more drastic the spatial change of tissue gray-level at that pixel position, and the higher the probability of corresponding tissue edge.

[0102] Step 44: Based on the spatial variation intensity values ​​of each pixel, arrange the pixels according to their spatial coordinates in the two-dimensional grayscale distribution matrix to obtain an edge response distribution map reflecting the sharpness of human tissue edges. This specifically includes:

[0103] Based on the spatial variation intensity values ​​of each generated pixel, they are arranged one-to-one according to their original row and column positions in the two-dimensional grayscale distribution matrix to construct an edge response distribution map of the same size as the reconstructed image. The value of each pixel in this distribution map is the spatial variation intensity value at that location. This edge response distribution map quantifies the edge distribution characteristics of human tissue in complex image space on a pixel-by-pixel basis: in uniform tissue regions (such as cerebrospinal fluid regions and large white matter regions), the intensity value approaches zero, indicating a smooth grayscale transition; in tissue boundary regions (such as the boundary between cortex and white matter, the boundary between blood vessel wall and lumen, and lesion contours), the intensity value exhibits local peaks. The height of the peak directly reflects the sharpness of the tissue edge at that location; the sharper the edge and the steeper the grayscale drop, the higher the peak. This edge response distribution map completely records the positional and sharpness information of various tissue boundaries in the reconstructed image.

[0104] In a preferred embodiment of the present invention, step 5 above, which involves calculating the cumulative phase deviation of each pixel based on the complex image spatial pixel matrix, and generating a global phase compensation field reflecting the global error distribution, includes:

[0105] Step 51: Based on the complex image space pixel matrix, extract the phase components corresponding to each pixel to obtain the initial phase distribution matrix, specifically including:

[0106] Based on the generated complex image space pixel matrix, each pixel location stores a complex value consisting of a real part and an imaginary part. The real part corresponds to the magnetization vector component intensity received by the RF coil in one of the orthogonal directions, and the imaginary part corresponds to the magnetization vector component intensity in the other direction. For each pixel location, the arctangent is calculated based on the ratio of its imaginary part to its real part. The quadrant in which the pixel's phase lies is determined based on the signs of the real and imaginary parts. Quadrant correction is applied to the principal value of the arctangent to obtain the phase component of the pixel. The corrected phase value lies between negative π and positive π. A positive value indicates that the precession phase of the magnetization vector at that pixel location leads the reference phase, while a negative value indicates lag.

[0107] The phase components of each pixel are arranged one-to-one according to their original row and column positions in the complex image space pixel matrix, forming an initial phase distribution matrix of the same size as the reconstructed image. This matrix records the phase spatial distribution of the magnetization vector in the complex image space pixel by pixel: under ideal imaging conditions, the phase of adjacent pixels should transition smoothly and continuously; however, under actual undersampling reconstruction conditions, due to the aliasing effect introduced by the lack of frequency domain information and the non-uniformity of the physical field of the RF coil, this phase distribution is superimposed with non-uniform phase deviations that vary with spatial position.

[0108] Step 52: Based on the initial phase distribution matrix, select pixels in the central region of the imaging field of view to calculate the average phase value, thus obtaining the average phase value of the central region. This specifically includes:

[0109] Based on the generated initial phase distribution matrix, in magnetic resonance imaging (MRI), the central region of the imaging field of view has better main magnetic field uniformity due to its proximity to the sensitivity core area of ​​the radio frequency coil and the center point of the gradient field. The influence of radio frequency field inhomogeneity is also relatively small. Therefore, the phase measurement value in the central region can serve as a unified reference benchmark for measuring the phase deviation of the entire image. Taking the geometric center of the image as the center, a circular region with a radius equal to one-eighth of the smaller row and column size of the image is selected as the central region. The phase values ​​of all pixels within this circular region are traversed. First, the mean and standard deviation of these phase values ​​are calculated. Then, values ​​exceeding the mean plus or minus twice the standard deviation are considered outliers and removed. These outliers occur at abrupt changes in magnetic susceptibility at tissue boundaries or at locations of blood flow artifacts; their phase values ​​do not reflect systematic bias, and retaining them would contaminate the accuracy of the reference benchmark. The arithmetic mean of the remaining effective phase values ​​after outlier removal is calculated, and the resulting value is the phase mean of the central region, serving as the global phase reference benchmark.

[0110] Step 53: Based on the initial phase distribution matrix and the average phase value of the central region, calculate the phase difference between the phase component of each pixel and the average phase value of the central region point by point to obtain the local phase offset set, which specifically includes:

[0111] Based on the phase values ​​of each pixel in the generated initial phase distribution matrix and the calculated average phase value of the central region, for each pixel, the original phase difference is obtained by subtracting the average phase value of the central region from the initial phase value. Since both the initial phase value and the average phase value of the central region are located in the range of -π to +π, the original difference value after subtracting the two ranges from -2π to +2π. The portion exceeding the range of -π to +π belongs to the false jumps caused by the phase value circling around the positive and negative π boundaries. For example, a pixel with an actual phase deviation of +0.1π. Because the original phase value is close to +π while the center mean is close to -π, a direct subtraction would yield a spurious difference of approximately +2π, which is completely inconsistent with the true minute deviation. To eliminate this spurious jump, a phase wrap correction is performed on the original phase difference: if the original difference is greater than positive π, it indicates that it has crossed the positive π boundary and wrapped around, so the original difference is subtracted by 2π to pull it back into the true deviation range; if the original difference is less than or equal to negative π, it indicates that it has crossed the negative π boundary and wrapped around, so the original difference is added by 2π to pull it back into the true deviation range; otherwise, the original difference remains unchanged. The phase difference after this correction falls back into the negative π to positive π range, reflecting the true physical phase deviation magnitude.

[0112] The phase difference values ​​of all pixel positions after wrapping correction are arranged in their original row and column positions to form a set of local phase offsets. Positive values ​​in this set indicate that the magnetization vector precession phase of the pixel is ahead of the central region, while negative values ​​indicate phase lag. These offsets show a continuous trend of slowly diverging or converging from the center to the periphery in space, and are superimposed with random high-frequency disturbances caused by local magnetic susceptibility differences and thermal noise.

[0113] Step 54: Based on the local phase offset set, the continuous trend of phase offset variation is deduced pixel by pixel along the spatial coordinate direction, and spatial smoothing filtering is performed for global summarization and calculation to obtain a global phase compensation field reflecting the global error distribution, specifically including:

[0114] Based on the generated set of local phase offsets, which contains two types of components: systematic global phase errors caused by the inhomogeneity of the RF coil physical field and gradient nonlinearity, exhibiting a slow, smooth, and continuous trend in space; and random local disturbances caused by differences in local tissue magnetic susceptibility, blood flow artifacts, or thermal noise, exhibiting rapid, discrete, and high-frequency fluctuations in space. The former is the target that needs to be compensated, while the latter should be filtered out; otherwise, false local corrections will be introduced into the compensation field, resulting in new phase artifacts in the reconstructed image. To separate the systematic global error component, a two-dimensional Gaussian low-pass spatial smoothing filter is applied to the set of local phase offsets. The window size of the Gaussian filter kernel is set to 5 rows × 5 columns of pixels. The filter coefficients at each position within the kernel are determined by sampling the two-dimensional Gaussian function on a discrete grid and then normalizing it. With the kernel center as the origin, the coefficients at each position are proportional to the negative exponent of the quotient obtained by dividing the square of the distance from the center by twice the square of the standard deviation. The standard deviation of the Gaussian filter kernel... Take as the phase standard deviation of the whole image 0.3 times: ;

[0115] in, It is obtained by calculating the standard deviation of the phase values ​​of all pixels in the initial phase distribution matrix. This value is based on the following:

[0116] The thermal noise level in MRI brain images typically corresponds to approximately The phase standard deviation, take It can effectively filter out thermal noise while retaining The global phase gradient variation trend of pixel width under different imaging conditions (such as field strength 3T vs 1.5T, different number of coil channels). The absolute value may vary, but a scaling factor of 0.3 can maintain consistent filtering characteristics. Alternatively, the window size of the Gaussian filter kernel can be fixed at [value missing]. Pixels The set of parameters is set to 1.0 pixels. Simulation results show that this set of parameters can effectively separate systematic phase errors and local random disturbances under typical cranial imaging conditions. The aim is to effectively filter out high-frequency noise while retaining the global phase change trend at low frequencies. The sum of all coefficients is normalized to equal 1 to ensure that the overall signal amplitude remains unchanged before and after filtering.

[0117] During filtering, the 5×5 Gaussian kernel is slid pixel by pixel across the local phase offset set, with the kernel center aligned with the currently processed pixel. The phase offset at each position within the kernel window coverage area is weighted and accumulated with the corresponding kernel coefficient, and the result is used as the smoothed phase offset value for that pixel position. For pixels at the image edge, when the filtering window overflows the boundary, the kernel coefficients on the overflow side are mirrored and folded back into the corresponding position within the window with the image boundary as the axis, so that the smoothed value at the edge naturally decays to near zero. The negative of the smoothed phase offset value for each pixel position after filtering is taken as the global phase compensation value for that pixel position. The global phase compensation values ​​of all pixels are arranged according to row and column positions to form a global phase compensation field. In this compensation field, positive values ​​correspond to regions where the original phase lags, and negative values ​​correspond to regions where the original phase leads.

[0118] In this embodiment of the invention, a complex image is generated through non-uniform Fourier transform, preserving the true frequency distribution of undersampled data and avoiding interpolation errors. High-frequency details and low-frequency background coefficients are obtained through wavelet multi-level decomposition, and an edge response distribution map is generated using the gray-scale spatial change rate to accurately locate tissue boundaries, overcoming the problem of reduced sharpness of true boundaries caused by uniform constraint strength. A global phase compensation field is generated based on the accumulation of global phase deviation, overcoming the problem that local correction cannot evaluate global errors, leading to phase cancellation and banded artifacts. The edge protection region is expanded through morphological dilation, ensuring smooth transition of adjacent edge regions. After point-by-point fusion of edge protection and phase compensation in the spatial domain, weighted coefficients are generated by wavelet sub-band matching. Weighted inverse transformation is performed on high and low frequency coefficients, enabling the two types of constraints to work synergistically in the reconstruction process. This achieves joint optimization of edge sharpness preservation and global phase error compensation, effectively suppressing Gibbs ring artifacts and banded brightness unevenness, and improving the reconstruction fidelity of medical images under undersampled conditions.

[0119] In a preferred embodiment of the present invention, step 6 above, which involves performing geodesic distance-based region expansion processing on the edge response distribution map to generate an expanded edge protection region, includes:

[0120] Step 61: Based on the edge response distribution map, extract the edge response values ​​of each pixel and filter out pixels that exceed a preset response threshold to obtain a set of edge seed points, specifically including:

[0121] Based on the generated edge response distribution map, the value of each pixel position in the distribution map is the calculated spatial change intensity value, which takes the range of non-negative real numbers. The larger the value, the more drastic the spatial change of the gray level of the tissue at that position, and the higher the probability that it belongs to the tissue boundary. The edge response value of the distribution map is extracted pixel by pixel.

[0122] The preset response threshold is taken as the 85th percentile of all pixel response values ​​in the edge response distribution map. That is, after sorting the response values ​​of all pixels in the entire image in ascending order, the value at the 85th percentile is taken as the response threshold. Using the 85th percentile as the preset response threshold can be adaptively determined based on the statistical distribution characteristics of the current image itself. This effectively eliminates weak false edge responses while covering the main tissue boundaries, avoiding misidentification of noise points as edge seeds. The reason for using the percentile method instead of a fixed absolute value is that the overall response level of the edge response distribution map may differ under different scanning parameters and different anatomical structures of the subjects. The percentile threshold can be adjusted accordingly. Based on the statistical distribution characteristics of the current image itself, the edge response distribution map is traversed. All pixels are marked as edge seed points if the edge response value is not less than the preset response threshold. Their pixel coordinates are recorded. All edge seed points constitute an edge seed point set. Each element in the set is a pixel coordinate, which identifies the location of the tissue boundary in the reconstructed image after preliminary determination. These edge seed points are the wavefront propagation starting points for subsequent morphological dilation and geodesic distance calculation. The dilation process starts from these seed points and spreads outward in the direction of gentle gray-scale change, stopping where the gray-scale change is drastic, thus forming a protective transition zone with controllable width on both sides of the true tissue boundary.

[0123] Step 62: Based on the magnitude components of the complex image spatial pixel matrix, calculate the gray-level difference between adjacent pixels point by point, and construct an anisotropic propagation cost map reflecting the cost of traversing the organization boundary. Specifically, this includes:

[0124] Based on the generated complex image spatial pixel matrix, the amplitude components of each pixel position are extracted to construct a two-dimensional gray-level matrix. For this two-dimensional gray-level matrix, the relationship between each pixel and its four neighboring pixels (i.e., the pixels above it) is calculated pixel by pixel. -1, ), below ( +1, ), left side ( , -1), Right side ( , The absolute value of the grayscale difference between (+1) and (+1) is used as the local propagation cost of the pixel along each neighborhood direction.

[0125] Let the pixel position in this two-dimensional grayscale matrix be ( The amplitude component at v is A( If v), then the propagation costs of this pixel along the four directions are as follows: ; ;

[0126] in , , , These correspond to the crossing costs in the upward, downward, leftward, and rightward directions, respectively. In regions where the absolute value of the grayscale difference between adjacent pixels is small, the grayscale transition is smooth, and the crossing costs in all four directions are low, allowing the wavefront to cross smoothly. In regions with tissue boundaries, the absolute value of the grayscale difference in at least one direction is large, increasing the crossing cost in that direction and hindering the propagation of the wavefront in that direction.

[0127] The propagation cost of each pixel in four directions is recorded according to the direction to form an anisotropic propagation cost map of the same size as the reconstructed image. The cost map has a lower overall cost in the uniformly organized area and a higher overall cost at the organization boundary. Its directional difference makes it easier for the wavefront to propagate along the boundary direction than to propagate across the boundary. The gray-level change of adjacent pixels along the boundary direction is small and the cost is low, while the gray-level change of the direction across the boundary is large and the cost is high.

[0128] Step 63: Using each edge seed point in the edge seed point set as the wavefront source point, perform geodesic distance fast-travel calculation based on the anisotropic propagation cost map to obtain the geodesic distance value from each pixel to the nearest edge seed point, specifically including:

[0129] Based on the generated set of edge seed points and the anisotropic propagation cost map, the geodesic distance from each pixel in the entire map to the nearest edge seed point is calculated using the fast travel method, with the edge seed points as the starting point of the wavefront and the local cost in each direction in the propagation cost map as the cumulative cost required for the wavefront to cross the neighboring pixels in that direction.

[0130] Anisotropic rapid travel method is used to solve geodesic distances. The discretized form of its discrete function equation is as follows: ;

[0131] in and The anisotropic propagation costs are respectively for the horizontal and vertical directions, and are defined as a function of the absolute value of the gray-level difference between adjacent pixels: ; ;

[0132] in, The magnitude component grayscale value at the pixel location. For extremely small positive numbers (such as To prevent division by zero errors and ensure numerical stability, the specific solution steps of the anisotropic FMM algorithm are as follows:

[0133] Initialize the geodesic distance of all pixels in the edge seed point set to... The first pixel is marked as a known state. All four neighboring pixels of a known pixel are included in the narrowband state, and their tentative geodesic distances are calculated from the distance values ​​of the known neighbors according to the discrete process function equation. The remaining pixels are marked as far-off states, with their geodesic distances temporarily set to infinity. The pixel with the smallest tentative geodesic distance is selected from the narrowband set, and its state is changed to known; this tentative value is the final geodesic distance. The pixels in the far-off states of the pixel's four neighbors are included in the narrowband, and their tentative distances are calculated. For neighboring pixels already in the narrowband state, if the tentative distance calculated via the path of the newly known pixel is less than its current value, it is updated to the smaller value. The main loop is repeated until all pixels are in the known state or the narrowband set is empty. This equation uses a two-way difference scheme, consistent with the windward discrete scheme of the standard FMM, and can be directly solved by the standard anisotropic FMM algorithm without iteration. This equation weights the propagation costs in the four directions, prioritizing propagation in low-cost directions and hindering propagation in high-cost directions, calculating the anisotropic cost. Figure 1 The geodesic distance, wavefront arrival pixel The propagation time is determined by the wavefront propagating from the pixel with the smallest distance in its neighborhood along the direction of lowest cost. The propagation rate is constrained by the local cost at that location. The specific calculation process is as follows:

[0134] Initialize the geodesic distance of all seed points in the edge seed point set to zero and mark them as known states. Include non-known pixels that are adjacent to the four neighbors of any known pixel in the narrowband state. For each narrowband pixel, use the above discrete process function equation to solve for its tentative geodesic distance from the distance value of the known neighbors. Mark the remaining pixels as far-instant states, and temporarily set the geodesic distance to infinity. Enter the main loop:

[0135] The pixel with the smallest tentative geodesic distance is selected from the narrowband set, and its state is changed to known. This tentative value is the final geodesic distance. Pixels in the four neighboring regions of this pixel that are in an unknown state are included in the narrowband, and their tentative distances are calculated. For neighboring pixels already in the narrowband state, if the tentative distance calculated through the path of the newly known pixel is smaller, it is updated to the smaller value. The main loop is repeated until all pixels are changed to known or the narrowband set is empty. After the above calculation, each pixel in the entire image is assigned a geodesic distance value, representing the cumulative cost of grayscale changes required to reach the nearest edge seed point from that pixel along the path with the lowest cumulative cost. Within a uniform organization, the grayscale difference between adjacent pixels is small, the propagation cost is low, and the geodesic distance increases slowly with spatial location, allowing subsequent expansion of the protected area to extend in this direction. At the organization boundary, the grayscale difference between adjacent pixels is large, the propagation cost is high, the geodesic distance increases sharply, and the expansion of the protected area naturally stops shrinking at this point.

[0136] Step 64: Based on the geodesic distance values ​​of each pixel, select pixels with geodesic distance values ​​lower than a preset distance threshold to include in the extended edge protection region. The boundary of the extended edge protection region extends along the direction of gradual grayscale change and automatically contracts at grayscale abrupt changes, specifically including:

[0137] Based on the calculated geodesic distance values ​​of all pixels in the image, the smaller the geodesic distance value, the closer the pixel is to the tissue edge seed point in the grayscale space, meaning the grayscale transition path between the pixel and the edge seed point is smoother and the cost is lower; the larger the geodesic distance value, the more obvious the tissue boundary is crossed between them, and the higher the propagation cost. The preset distance threshold is taken as the median of the non-zero geodesic distance values ​​of the entire image. The physical basis for taking the median is:

[0138] In MRI images of the brain, the geodesic distance from pixels within tissue regions to the nearest edge accounts for approximately half of the total number of pixels. Setting the threshold at the median allows these homogeneous tissue regions to be included in the extended protection range, while excluding deep tissue core areas and external regions that cross obvious boundaries. By traversing all pixels in the image, pixels with geodesic distance values ​​below the preset distance threshold are assigned to the extended edge protection region and marked with an edge protection label; those with values ​​at or above the threshold are assigned to the unprotected region, thus obtaining the extended edge protection region. Each pixel in this region stores its geodesic distance value. The smaller the geodesic distance value, the closer the pixel is to the tissue edge seed point in grayscale space; the larger the value, the farther the distance. Within homogeneous tissue regions, the cost is low and the geodesic distance increases slowly, allowing for a large-scale extension along this direction after filtering with the median threshold. At tissue boundaries with abrupt grayscale changes, the cost increases sharply, the geodesic distance rapidly climbs above the threshold, and the protection region automatically shrinks and stops at the boundary.

[0139] In a preferred embodiment of the present invention, step 7 above involves jointly solving the extended edge protection region and the global phase compensation field in the spatial domain to generate a multi-constraint reconstruction parameter set, including:

[0140] Step 71: Based on the extended edge protection region, extract the dilation response value of each pixel and perform proportional scaling based on the maximum dilation response value in the set of dilation response values ​​to obtain the edge protection weight factor set, which specifically includes:

[0141] Based on the geodesic distance values ​​stored in each pixel within the extended edge protection region, the reciprocal of the geodesic distance value for each pixel is taken and normalized to obtain the dilation response value, thus constructing a continuous value map of the same size as the reconstructed image. For pixel locations... The expansion response value is calculated as follows: ;

[0142] in For all pixels in the entire image The maximum value for geodetic distance For seed pixel points with a value of 0, set them to... Set value to full image To avoid division by zero errors, the normalized expansion response value ranges from the maximum value of the value. The dilation response value reaches its maximum at the seed point at the tissue edge and decays with increasing geodesic distance. This dilation response value continuously quantifies the gray-level affinity between each pixel and the nearest tissue edge. The larger the value, the more likely the pixel is at the tissue edge or in the transition zone near the edge, and the more necessary it is to preserve its structural information in the high-frequency detail coefficient weighted reconstruction. The smaller the value, the more likely the pixel is far from the edge, and stronger constraints can be used to suppress aliasing artifacts.

[0143] The dilation response value of each pixel is extracted from the continuous value map. All pixels are traversed to find the maximum value in the dilation response value set. This maximum value is usually located at the original edge seed point. The dilation response value is close to 1. Based on this maximum dilation response value, the dilation response value of each pixel is scaled proportionally. That is, the dilation response value of each pixel is divided by the maximum value. After scaling, the dilation response values ​​of all pixels fall in the range of [0,1], with a maximum value of 1. The rest are scaled down proportionally. This scaling process eliminates the difference in the absolute amplitude of the dilation response value under different imaging conditions, so that the edge protection weight factor maintains a consistent value range and constraint strength in the application across cases and scanning parameters. The scaled dilation response values ​​of each pixel are arranged according to their original spatial coordinate positions to form an edge protection weight factor set. The pixel positions in this set with a weight factor close to 1 correspond to the core area of ​​the tissue edge.

[0144] Step 72: Based on the global phase compensation field, extract the global phase compensation value of each pixel and perform phase direction inversion processing to obtain a set of phase compensation factors, specifically including:

[0145] Based on a global phase compensation field, the global phase compensation value stored at each pixel location reflects the direction and amount of phase correction to be applied at the corresponding pixel location, with its positive and negative signs and magnitudes indicating that the original phase lags behind the central reference and requires positive compensation to advance the phase; negative values ​​indicate that the original phase leads the central reference and requires negative compensation to pull the phase back. The larger the absolute value of the compensation value, the more severe the systematic phase deviation at that location. The global phase compensation value is extracted pixel by pixel from the global phase compensation field. The phase direction reversal processing refers to using the global phase compensation value of each pixel as the phase angle parameter in the complex phase factor, constructing a complex factor with a modulus of 1 according to Euler's formula. The phase angle of this factor is the negative of the global phase compensation value. In the subsequent weighted inverse transformation step, this complex factor is multiplied by the wavelet transform coefficient of the corresponding pixel location, which is equivalent to applying a phase rotation to the transform coefficient that is equal in magnitude and opposite in direction to the original phase error, thus canceling the systematic phase deviation accumulated at that pixel location during the undersampling reconstruction process. Let the pixel... The global phase compensation value at that location is Then the phase compensation factor corresponding to that pixel for: ;

[0146] The modulus of this factor is 1, and the phase angle is... In the subsequent weighted inverse transform, this factor is applied multiplicatively to the corresponding position of the wavelet low-frequency approximation coefficients, causing the phase of the coefficient at that position to rotate. radians, position of phase lag ( Its coefficient phase is rotated positively to catch up with the central reference, and the phase leads the position ( Its coefficient phase is rotated negatively to return to the center reference, thus affecting all pixels. Arranged according to spatial coordinates, they form a set of phase compensation factors of the same size as the global phase compensation field.

[0147] Step 73: Based on the edge protection weight factor set and the phase compensation factor set, perform point-by-point multiplication operations for the same spatial coordinate positions to obtain the local constraint weight values ​​corresponding to each spatial coordinate position, specifically including:

[0148] Based on the generated set of edge protection weight factors and the generated set of phase compensation factors, the two sets have the same matrix size, and elements at the same row and column positions correspond to the same pixel position. For each pixel position... Extract the edge protection weight factor corresponding to this location. (Real number, ranging from [0,1]) and phase compensation factor (Complex number, modulus 1) Perform a complex multiplication of the edge protection weight factor and the phase compensation factor at this pixel location. Since the modulus of the phase compensation factor is 1, this multiplication is equivalent to using the edge protection weight factor as the amplitude and the phase angle of the phase compensation factor as the phase, thus forming a complete complex constraint weight: ;

[0149] The complex number constraint weight It simultaneously carries two types of constraint information, namely the modulus. Edge protection constraints are encoded; a larger modulus indicates that the tissue edge features at that pixel location need to be preserved in the reconstruction; a smaller modulus indicates that the transform coefficients at that pixel location can be more strongly suppressed, and the phase angle... The phase compensation constraint is encoded. The phase angle indicates the direction and magnitude of the phase rotation to be applied at that pixel location. The process iterates through all pixel locations in the image, performing the product operation point by point, and then converts each pixel's... Arranged according to spatial coordinates, they form a set of local constraint weight values, which unifies the two types of constraints, edge protection and phase compensation, in a compact complex number form.

[0150] Step 74: Based on the local constraint weight values ​​corresponding to each spatial coordinate position, channel splicing is performed according to the spatial coordinates to obtain a multi-constraint reconstruction parameter set including edge protection weight factors and phase compensation factors, specifically including:

[0151] Based on the generated set of local constraint weight values Each of them For a complex number, The modulus and phase angle of each complex number are split into two independent channels, and the splitting is performed pixel by pixel:

[0152] Extract complex constraint weights for pixel position (u,v) model Store in the corresponding position of the edge protection constraint channel; extract phase angle Store the corresponding position in the phase compensation constraint channel, and the split edge protection constraint channel. and phase compensation constraint channel They are respectively: ; ;

[0153] in It is a real number matrix, where each element represents the spatial distribution of the edge protection strength at the corresponding pixel location; The first matrix is ​​a real-valued matrix, where each element represents the spatial distribution of the phase compensation angle to be applied at the corresponding pixel location. These two matrices are concatenated along the channel dimension, logically maintaining independence in two-dimensional matrix form while physically serving as parallel subsets of the same parameter set, thus forming a multi-constraint reconstruction parameter set. This parameter set fully carries the edge protection weights and phase compensation factors for each pixel location in the image space.

[0154] like Figure 3 As shown, in a preferred embodiment of the present invention, step 8 above, based on the multi-constraint reconstruction parameter set and the sub-band spatial partitioning method of the wavelet transform basis, generates a weighted coefficient matrix matching the size of the high-frequency detail coefficients and the low-frequency background coefficients, and performs a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets based on the weighted coefficient matrix to obtain the target medical image, including:

[0155] Step 81: Extract the edge protection weight factor set and phase compensation factor set from the multi-constraint reconstruction parameter set, specifically including:

[0156] Based on the generated multi-constraint reconstruction parameter set, which consists of two parallel subsets—edge protection constraint channels and phase compensation constraint channels—both have the same size as the reconstructed image, the following parameters are extracted channel by channel: All values ​​of the edge protection constraint channels are extracted from this parameter set to obtain the edge protection weight factor set, which is a real number matrix, where each element... The strength of tissue edge protection at the corresponding pixel location is represented by a value. A larger value indicates that the high-frequency detail coefficients at that location should be preserved with higher weight in subsequent weighted reconstruction, while a smaller value indicates that aliasing artifacts should be suppressed more strongly. All values ​​of the phase compensation constraint channel are used to obtain the phase compensation factor set, which is a real number matrix, where each element... This represents the phase compensation angle that needs to be applied to the corresponding pixel position.

[0157] Step 82: Based on the subband space partitioning method of the preset discrete wavelet transform basis, downsample the edge protection weight factor set to obtain a high-frequency weighting matrix matching the size of each high-frequency detail subband; downsample the phase compensation factor set to obtain a low-frequency weighting matrix matching the size of each low-frequency background subband, specifically including:

[0158] Based on the extracted set of edge protection weight factors and the set of phase compensation factors The wavelet decomposition uses the Daubechies-4 wavelet basis for two-level decomposition, generating three first-level high-frequency detail subbands, three second-level high-frequency detail subbands, and one second low-frequency approximate subband. In order to ensure that the constraint weights are precisely aligned with the coefficients of each subband in spatial position, the elements at a given row and column position in the weight matrix correspond to the original image spatial regions covered by the coefficients at the same position in the subband. Therefore, the full-resolution weight matrix needs to be downsampled using the same scale and filtering method as the wavelet decomposition.

[0159] Perform the same Daubechies-4 low-pass and high-pass filtering and step-size 2 downsampling on the edge protection weight factor set as described above for wavelet decomposition, to obtain a weighting matrix that matches the size of the first high-frequency detail subband. Repeat the above operation on the result obtained after the first downsampling to obtain a weighted matrix that matches the size of the second high-frequency detail subband. Each element in this weighting matrix is ​​associated with one of the three high-frequency detail subbands of the first (or second) level. , , The coefficients at the same row and column positions in the set of phase compensation factors correspond one-to-one in space. :

[0160] Similarly, using the Daubechies-4 low-pass decomposition filter, for By performing two consecutive low-pass filters and downsampling steps of 2, a size of [size missing] is obtained. Low-frequency phase compensation matrix The spatial resolution of this matrix is ​​similar to that of the second low-frequency approximate subband. Completely identical, each element in the matrix records the corresponding The average phase compensation angle to be applied within the original pixel block, after the above processing... Used for edge-protection weighting of the first-level high-frequency detail coefficients. Used for edge-protection weighting of the second-level high-frequency detail coefficients. Used for phase correction of the second low-frequency approximate subband coefficients, with the three levels aligned and sizes matched.

[0161] Step 83: Perform element-wise multiplication of the high-frequency weighting matrix with the corresponding high-frequency detail coefficients to obtain a weighted high-frequency detail coefficient set; perform element-wise multiplication of the low-frequency weighting matrix with the corresponding low-frequency background coefficients to obtain a weighted low-frequency background coefficient set, specifically including:

[0162] Based on the generated two-level high-frequency weighting matrix , and low-frequency phase compensation matrix And the six high-frequency detail subbands and the second low-frequency approximate subband in the sparse transform coefficient set, with weighted processing for the high-frequency detail coefficients:

[0163] For the three high-frequency detail subbands of the first level ( , , ), and the coefficient elements in each sub-band are compared with The weights at the same row and column positions are multiplied element-wise to obtain the weighted first-level high-frequency detail subband, which is expressed as follows: ;

[0164] in, In the first-order wavelet decomposition, the position is... The edge protection weight coefficient at the location, with a value of The value is calculated from the extended edge protection region and controls the degree of high-frequency detail preservation at this location. In the first-order wavelet decomposition, the direction is... (LH horizontal, HL vertical, HH diagonal) in position The original high-frequency detail subband coefficients at the location reflect the image edge intensity. For the three high-frequency detail subbands of the second level, each subband coefficient is compared with... Performing element-wise multiplication yields the weighted second-level high-frequency detail subband: ;

[0165] in, In the second-level wavelet decomposition, the position... The edge protection weight coefficient at the location, with a value of Similar to the first level but for a coarser scale. In the second-order wavelet decomposition, the direction is... In position The original high-frequency detail subband coefficients are almost completely preserved when the weight value is close to 1, and the sharpness of the tissue edges is maintained after reconstruction; when the weight value is close to 0, the original coefficients are significantly suppressed, and aliasing artifacts and noise are effectively removed. The six weighted high-frequency detail subbands constitute the weighted high-frequency detail coefficient set. For phase correction of low-frequency background coefficients:

[0166] The coefficients in the second low-frequency approximation subband LL2 remain unchanged as the original real coefficients, serving as a set of weighted low-frequency background coefficients. The low-frequency background coefficients belong to the real domain amplitude information and do not carry phase components. The systematic phase error introduced by undersampling has been eliminated in the spatial domain by applying phase correction to the spatial pixel matrix of the complex image. There is no need to perform phase compensation on the low-frequency coefficients in the wavelet domain. The weighted second low-frequency approximation subband and the six weighted high-frequency detail subbands together constitute the set of weighted coefficients to be inversely transformed.

[0167] Step 84: Based on the weighted high-frequency detail coefficient set and the weighted low-frequency background coefficient set, the inverse wavelet transform kernel matching the preset discrete wavelet transform basis is invoked to perform layer-by-layer inverse reconstruction synthesis to obtain the target medical image, specifically including:

[0168] Based on the generated weighted high-frequency detail coefficient set and weighted low-frequency background coefficient The inverse wavelet transform reconstruction uses reconstruction filter coefficients orthogonally matched with the Daubechies-4 wavelet basis used in the forward decomposition of step 31. The low-pass reconstruction filter coefficients are... The coefficients of the high-pass reconstruction filter are The two sets of coefficients are obtained by reversing the order of the decomposition filter coefficients according to an orthogonal mirror relationship and alternating their signs. The reconstruction process is performed layer by layer in a coarse-to-fine hierarchical order, first performing the inverse transform from level 2 to level 1, and then reconstructing the weighted second low-frequency approximate subband. The reverse operation is performed along the column direction. For each column, upsampling with a stride of 2 is performed first, then convolving the coefficients of the low-pass reconstruction filter along the column direction to obtain the low-frequency component of the column reconstruction; convolving the coefficients of the high-pass reconstruction filter along the column direction to obtain the high-frequency component of the column reconstruction. The row and column reconstruction results are then reversed along the row direction. For each row, upsampling with a stride of 2 is performed, and convolving with both the low-pass and high-pass reconstruction filters respectively. The final synthesized size is... Level 1 low-frequency approximate sub-band The weighted second-level three high-frequency detail sub-bands , , After sampling and filtering in the same column and row directions, and with Add the corresponding sub-band positions to complete the full reconstruction from level 2 to level 1.

[0169] The reconstructed With the weighted first-level three high-frequency detail subbands , , Perform upsampling and reconstruction filtering again with a step size of 2 in both the column and row directions, resulting in a final synthesized size of... Full-resolution complex image matrix ,right Pixel-by-pixel modulo operation involves taking the square root of the sum of the squares of the real and imaginary parts of the complex value at each pixel location to obtain the reconstructed grayscale value for that pixel location. The reconstructed grayscale value is then arranged by row and column position. The real grayscale image is the target medical image. While suppressing undersampling aliasing artifacts, the sharpness of tissue edges is maintained and the global phase deviation is effectively compensated. The overall image uniformity is improved compared with the conventional compressed sensing results without multi-constraint weighted reconstruction.

[0170] In this embodiment of the invention, by fusing edge protection weights and phase compensation factors point by point into a complex set of multi-constraint reconstruction parameters, the technical problems of isolated processing of edge preservation and phase compensation constraints in the spatial domain, the mismatch between full-resolution weights and wavelet subband coefficient resolution, and the need to convert real coefficients into complex form for low-frequency phase correction are overcome. This achieves the synergistic effect of edge sharpness preservation and global phase error compensation in the weighted inverse transformation, effectively suppressing Gibbs ring artifacts and banded brightness unevenness, and improving the geometric fidelity and grayscale uniformity of medical image reconstruction under undersampling conditions.

[0171] like Figure 2 As shown, embodiments of the present invention also provide a compressed sensing reconstruction system for medical images, comprising:

[0172] The acquisition module is used to acquire the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal.

[0173] The conversion module is used to perform non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image space pixel matrix;

[0174] The decomposition module is used to call a preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix, and obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background.

[0175] The calculation module is used to calculate the gray-scale spatial change rate between adjacent pixels point by point based on the complex image spatial pixel matrix, and generate an edge response distribution map that reflects the sharpness of tissue edges; based on the complex image spatial pixel matrix, it calculates the cumulative phase deviation of each pixel point across the entire domain, and generates a global phase compensation field that reflects the global error distribution.

[0176] The reconstruction module is used to perform geodesic distance-based regional expansion processing on the edge response distribution map to generate an expanded edge protection area; and to perform spatial domain point-by-point weighted fusion operation on the expanded edge protection area and the global phase compensation field to generate a multi-constraint reconstruction parameter set.

[0177] The matching module is used to generate a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis. Then, it performs a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets based on the weighted coefficient matrix to obtain the target medical image.

[0178] 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.

[0179] Embodiments of the present invention also provide a computing device, including: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0180] 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.

[0181] 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 compressed sensing reconstruction method for medical images, characterized in that, The method includes: Step 1: Acquire the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal; Step 2: Perform non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image spatial pixel matrix; Step 3: Call the preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix to obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background. Step 4: Based on the complex image spatial pixel matrix, calculate the grayscale spatial change rate between adjacent pixels point by point to generate an edge response distribution map that reflects the sharpness of tissue edges. Step 5: Based on the complex image spatial pixel matrix, calculate the cumulative phase deviation of each pixel across the entire domain to generate a global phase compensation field that reflects the global error distribution. Step 6: Perform geodesic distance-based regional expansion processing on the edge response distribution map to generate an expanded edge protection area; Step 7: Perform spatial domain point-by-point weighted fusion calculation between the extended edge protection region and the global phase compensation field to generate a multi-constraint reconstruction parameter set; Step 8: Based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis, generate a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients. Then, perform a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets according to the weighted coefficient matrix to obtain the target medical image.

2. The compressed sensing reconstruction method for medical images according to claim 1, characterized in that, Acquire undersampled frequency domain signals from the radio frequency coil of medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signals, including: The raw undersampled frequency domain signal received by the radio frequency coil of the medical imaging equipment is acquired, and the radio frequency pulse timing control information carried by the raw undersampled frequency domain signal is analyzed to obtain the acquisition sequence containing spatial coordinates and timestamps. Based on the spatial coordinates and timestamps contained in the acquisition sequence, the spatial coordinates of each sampling point are arranged in ascending order of timestamps. The spatial displacement vector between sampling points corresponding to adjacent timestamps is calculated to obtain the physical acquisition trajectory parameters. Based on the spatial displacement vector contained in the physical acquisition trajectory parameters, the target imaging space is divided into multiple uniform grids, and the number of sampling points falling in each uniform grid is calculated to obtain the initial sampling density distribution. Based on the number of sampling points included in the initial sampling density distribution, the proportion of the number of sampling points in each uniform grid to the total number of global sampling points is calculated to obtain the spatial sampling density distribution.

3. The compressed sensing reconstruction method for medical images according to claim 2, characterized in that, Based on the physical acquisition trajectory parameters and sampling density distribution, a non-uniform Fourier transform is performed on the undersampled frequency domain signal to generate a complex image spatial pixel matrix, including: Based on the spatial frequency correspondence between the spatial coordinates contained in the physical acquisition trajectory parameters and the pixel coordinates of the target imaging space, the spatial phase rotation factor corresponding to each sampling point is calculated to obtain the set of phase rotation factors. Based on the set of phase rotation factors and the undersampled frequency domain signal, perform point-by-point complex multiplication to obtain the modulation frequency domain signal corresponding to each sampling point; Based on the modulation frequency domain signal and spatial sampling density distribution, a global signal accumulation calculation is performed along the coordinate direction of the target imaging space to obtain a complex image spatial pixel matrix containing amplitude and phase components.

4. The compressed sensing reconstruction method for medical images according to claim 3, characterized in that, The complex image spatial pixel matrix is ​​decomposed into multiple coefficients by invoking a preset wavelet transform basis, obtaining a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background, including: Based on the complex image spatial pixel matrix, the amplitude components of each pixel are extracted to construct a real-valued image matrix. The preset discrete wavelet transform basis is called to perform spatial domain low-pass filtering and high-pass filtering operations on the real-valued image matrix to separate the first low-frequency approximate sub-band and the first high-frequency detail sub-band. Based on the first low-frequency approximate sub-band, a spatial domain multi-scale filtering operation is performed by calling a preset discrete wavelet transform basis to separate the second low-frequency approximate sub-band and the second high-frequency detail sub-band. Based on the first high-frequency detail subband, the second high-frequency detail subband, and the second low-frequency approximation subband, they are aligned according to their spatial coordinate positions and arranged in a hierarchical nested manner to obtain a sparse transform coefficient set containing high-frequency detail coefficients and low-frequency background coefficients.

5. The compressed sensing reconstruction method for medical images according to claim 4, characterized in that, Based on the complex image spatial pixel matrix, the gray-level spatial change rate between adjacent pixels is calculated point by point to generate an edge response distribution map reflecting the sharpness of tissue edges, including: Based on the complex image space pixel matrix, extract the gray value of the amplitude component corresponding to each pixel and construct a two-dimensional gray distribution matrix. Based on the two-dimensional gray-level distribution matrix, the gray-level values ​​of adjacent pixels are extracted along the horizontal and vertical coordinate directions respectively, and the difference is calculated to obtain the set of gray-level differences in the horizontal direction and the set of gray-level differences in the vertical direction. Based on the set of grayscale differences in the horizontal direction and the set of grayscale differences in the vertical direction, the energy aggregation value of the grayscale difference in the horizontal and vertical directions of the same pixel is calculated point by point to obtain the spatial change intensity value of each pixel. Based on the spatial variation intensity value of each pixel, the pixels are arranged according to their spatial coordinates in the two-dimensional grayscale distribution matrix to obtain an edge response distribution map that reflects the edge sharpness of human tissue.

6. The compressed sensing reconstruction method for medical images according to claim 5, characterized in that, Based on the complex image spatial pixel matrix, the cumulative phase deviation of each pixel is calculated globally, generating a global phase compensation field that reflects the global error distribution, including: Based on the complex image space pixel matrix, the phase components corresponding to each pixel are extracted to obtain the initial phase distribution matrix; Based on the initial phase distribution matrix, the average phase value of the pixels in the central region of the imaging field of view is calculated to obtain the average phase value of the central region. Based on the initial phase distribution matrix and the average phase value of the central region, the phase difference between the phase component of each pixel and the average phase value of the central region is calculated point by point to obtain the set of local phase offsets. Based on the set of local phase offsets, the continuous trend of phase offset change is deduced pixel by pixel along the spatial coordinate direction, and spatial smoothing filtering is performed to perform global summary calculation to obtain a global phase compensation field that reflects the global error distribution.

7. The compressed sensing reconstruction method for medical images according to claim 6, characterized in that, Perform geodesic distance-based region extension processing on the edge response distribution map to generate an extended edge protection area, including: Based on the edge response distribution map, the edge response values ​​of each pixel are extracted and pixels exceeding the preset response threshold are filtered to obtain the edge seed point set; Based on the magnitude components of the complex image spatial pixel matrix, the gray-level difference between adjacent pixels is calculated point by point to construct an anisotropic propagation cost map that reflects the cost of crossing the organizational boundary. Using each edge seed point in the edge seed point set as the wavefront source point, geodesic distance fast travel calculation is performed based on the anisotropic propagation cost map to obtain the geodesic distance value from each pixel to the nearest edge seed point. Based on the geodesic distance value of each pixel, pixels with geodesic distance values ​​lower than a preset distance threshold are selected and included in the extended edge protection area. The boundary of the extended edge protection area extends along the direction of gradual grayscale change and automatically shrinks when encountering abrupt grayscale changes.

8. The compressed sensing reconstruction method for medical images according to claim 7, characterized in that, The extended edge protection region and the global phase compensation field are subjected to spatial domain point-by-point weighted fusion calculation to generate a multi-constraint reconstruction parameter set, including: Based on the extended edge protection region, the dilation response value of each pixel is extracted and the maximum dilation response value in the dilation response value set is used as the benchmark to perform proportional scaling to obtain the edge protection weight factor set. Based on the global phase compensation field, the global phase compensation value of each pixel is extracted and phase direction inversion processing is performed to obtain a set of phase compensation factors. Based on the set of edge protection weight factors and the set of phase compensation factors, point-by-point multiplication is performed for the same spatial coordinate positions to obtain the local constraint weight values ​​corresponding to each spatial coordinate position. Based on the local constraint weight values ​​corresponding to each spatial coordinate position, channels are spliced ​​according to spatial coordinates to obtain a multi-constraint reconstruction parameter set that includes edge protection weight factors and phase compensation factors.

9. The compressed sensing reconstruction method for medical images according to claim 8, characterized in that, Based on the multi-constraint reconstruction parameter set and the sub-band spatial partitioning method of the wavelet transform basis, a weighted coefficient matrix matching the size of the high-frequency detail coefficients and the low-frequency background coefficients is generated. Then, a weighted inverse transform is performed on the high-frequency and low-frequency sparse transform coefficient sets based on the weighted coefficient matrix to obtain the target medical image, including: Extract the edge protection weight factor set and the phase compensation factor set from the multi-constraint reconstruction parameter set respectively; Based on the subband space partitioning method of the preset discrete wavelet transform basis, the edge protection weight factor set is downsampled to obtain a high-frequency weighting matrix that matches the size of each high-frequency detail subband; the phase compensation factor set is downsampled to obtain a low-frequency weighting matrix that matches the size of each low-frequency background subband. Perform element-wise multiplication of the high-frequency weighting matrix with the corresponding high-frequency detail coefficients to obtain the weighted high-frequency detail coefficient set; perform element-wise multiplication of the low-frequency weighting matrix with the corresponding low-frequency background coefficients to obtain the weighted low-frequency background coefficient set. Based on the weighted high-frequency detail coefficient set and the weighted low-frequency background coefficient set, the inverse wavelet transform kernel that matches the preset discrete wavelet transform basis is invoked to perform layer-by-layer inverse reconstruction synthesis to obtain the target medical image.

10. A compressed sensing reconstruction system for medical images, wherein the system implements the method as described in any one of claims 1 to 9, characterized in that, include: The acquisition module is used to acquire the undersampled frequency domain signal collected by the radio frequency coil of the medical imaging equipment, and extract the physical acquisition trajectory parameters and sampling density distribution corresponding to the undersampled frequency domain signal. The conversion module is used to perform non-uniform Fourier transform on the undersampled frequency domain signal based on the physical acquisition trajectory parameters and sampling density distribution to generate a complex image space pixel matrix; The decomposition module is used to call a preset wavelet transform basis to perform multi-level coefficient decomposition on the complex image spatial pixel matrix, and obtain a set of high-frequency and low-frequency sparse transform coefficients containing high-frequency details and low-frequency background. The calculation module is used to calculate the gray-scale spatial change rate between adjacent pixels point by point based on the complex image spatial pixel matrix, and generate an edge response distribution map that reflects the sharpness of tissue edges; based on the complex image spatial pixel matrix, it calculates the cumulative phase deviation of each pixel point across the entire domain, and generates a global phase compensation field that reflects the global error distribution. The reconstruction module is used to perform geodesic distance-based regional expansion processing on the edge response distribution map to generate an expanded edge protection area; The extended edge protection region and the global phase compensation field are subjected to spatial domain point-by-point weighted fusion calculation to generate a multi-constraint reconstruction parameter set; The matching module is used to generate a weighted coefficient matrix that matches the size of the high-frequency detail coefficients and the low-frequency background coefficients based on the multi-constraint reconstruction parameter set and the sub-band spatial division method of the wavelet transform basis. Then, it performs a weighted inverse transform on the high-frequency and low-frequency sparse transform coefficient sets based on the weighted coefficient matrix to obtain the target medical image.

Images