A digital x-ray image grid artifact suppression method, system and device

Abstract

The application discloses a digital X-ray image grid artifact suppression method, system and device, relates to the technical field of image data processing, and comprises the following steps: performing Laplacian pyramid decomposition on a digital X-ray image to be processed to obtain Laplacian pyramid layer images of different scales; for each Laplacian pyramid layer image, performing artifact suppression processing, which comprises the following steps: determining target frequency positions and conjugate frequency positions of grid artifacts according to a two-dimensional spectrum image; constructing a local analysis area and performing spectrum projection to obtain a one-dimensional energy distribution curve; filtering the two-dimensional spectrum image, performing inverse Fourier transform, and obtaining a grid artifact removal result of a current scale layer; and performing Laplacian pyramid reconstruction on all grid artifact removal results to obtain a digital X-ray image in which grid artifacts are suppressed. The application improves the image quality of the digital X-ray image and the accuracy of medical diagnosis.

Description

Technical Field

[0001] This application relates to the field of image data processing technology, and in particular to a method, system and device for suppressing artifacts in digital X-ray image filter grids. Background Technology

[0002] In digital X-ray imaging, a grid is typically placed in front of the detector to filter out Compton scatter rays and improve image contrast. As a periodic spatial modulation structure, the grid introduces stable high-frequency components into the image. When these high-frequency components are discretely sampled by the detector with finite spatial resolution, and their frequency approaches or exceeds the system's Nyquist frequency, spatial aliasing occurs. This causes the originally invisible grid lines to fold into visible low-frequency stripe artifacts, reducing image quality and the accuracy of medical diagnosis.

[0003] Currently, most methods for suppressing grid artifacts in digital X-ray images are based on frequency domain analysis. This involves performing a two-dimensional Fourier transform on the entire digital X-ray image, detecting the target frequency position corresponding to the grid artifact in the spectrum, suppressing the target frequency component based on a notch filter, and finally restoring it to the spatial domain through an inverse Fourier transform.

[0004] However, the above method directly performs spectral analysis and filtering on the entire original digital X-ray image. The frequency components such as anatomical structures, high-frequency textures, and random noise in the image are superimposed on the frequency of the filter grid artifacts in the same spectrum. This can easily cause the target frequency peak and the image detail frequency to overlap, reducing the stability and accuracy of target frequency localization. Furthermore, the overall filtering process loses image edge and texture details, affecting image quality and the accuracy of medical diagnosis. Therefore, there is room for improvement. Summary of the Invention

[0005] In order to overcome the shortcomings of the existing technology and improve the image quality and medical diagnostic accuracy of digital X-ray images, this application provides a method, system and device for suppressing digital X-ray image filter grid artifacts.

[0006] Firstly, the objective of this invention is achieved through the following technical solution: A method for suppressing filter grid artifacts in digital X-ray images includes: Acquire the digital X-ray image to be processed; The digitized X-ray image was decomposed into a Laplacian pyramid to obtain multiple Laplacian pyramid layer images at different scales; For each of the Laplacian pyramid layer images, artifact suppression processing is performed, which includes: A two-dimensional Fourier transform is performed on the Laplace pyramid layer image to obtain a two-dimensional spectral image; Based on the two-dimensional spectrum image, the target frequency position and conjugate frequency position of the filter grid artifact are adaptively determined. A local analysis region is constructed with the target frequency position as the center, and spectral projection is performed along the frequency peak width direction and the frequency peak length direction to obtain a one-dimensional energy distribution curve. Based on the one-dimensional energy distribution curve, the spectral width parameter and spectral length parameter are calculated. A two-dimensional Gaussian notch filter is adaptively generated based on the target frequency position, the conjugate frequency position, the spectral width parameter, and the spectral length parameter. The two-dimensional spectral image is filtered using the two-dimensional Gaussian notch filter, and an inverse Fourier transform is performed on the filtered spectrum to obtain the artifact removal result of the current scale layer. Laplacian pyramid reconstruction was performed on the artifact removal results of all scale layers to obtain a digital X-ray image with suppressed grid artifacts.

[0007] By employing the above technical solution, multiple Laplacian pyramid layer images of different scales are used to characterize different spatial frequency components in the image, with a two-dimensional Fourier transform converting the image from the spatial domain to the frequency domain. After artifact suppression processing is performed on the image, Laplacian pyramid reconstruction is then carried out. The two-dimensional Gaussian notch filter forms a frequency suppression region at the target frequency position and the conjugate frequency position, adapting to the actual frequency distribution range of the digital X-ray image to be processed in the frequency domain. This allows the two-dimensional Gaussian notch filter to adaptively adjust its filtering parameters according to the digital X-ray image to be processed, thereby optimizing the image data processing effect and improving the parameter determination accuracy of the two-dimensional Gaussian notch filter.

[0008] Specifically, firstly, by introducing Laplacian pyramid decomposition, the digitized X-ray image to be processed is decomposed into multiple layers of different scales according to spatial frequency. This separates the grid artifacts, which were originally overlapping and spanned a wide frequency band in the original image, into layers of different scales corresponding to their stripe thicknesses. Since each scale layer only carries frequency components with a limited bandwidth, the grid artifacts appear as more concentrated and easily identifiable frequency peaks in the spectrum of the corresponding layer, thus overcoming the deficiency of single-scale overall frequency domain processing in the previous technique, which could not take into account stripes of different thicknesses. On this basis, artifact suppression is performed on each layer separately, and then the images are fused by pyramid reconstruction. This allows artifact components of different scales to be removed layer by layer in a targeted manner, while the normal image structure of each layer is complementaryly restored during reconstruction, significantly reducing the damage of filtering to normal anatomical structures and edge details.

[0009] Secondly, by adaptively determining the target frequency position and its conjugate frequency position for each layer of the spectrum, the filtering position no longer relies on preset values ​​or manual experience, but is dynamically determined according to the image content, the grid specifications, and the actual spectral distribution of the scale layer. Consistency verification is then performed using the conjugate symmetry characteristics of the Fourier spectrum, ensuring the accuracy and robustness of frequency positioning. This overcomes the shortcomings of previous technologies where fixed notch filters had fixed filtering positions and could not adapt to different acquisition conditions.

[0010] Furthermore, by constructing a local analysis region centered on the target frequency location, and projecting the spectrum along both the width and length directions of the frequency peak, a one-dimensional energy distribution curve is obtained. Based on this, the spectral width and length parameters are calculated, allowing for a quantitative characterization of the actual bandwidth and extension range of the artifact frequency peak. Then, driven by these parameters and the target and conjugate frequency locations, a two-dimensional Gaussian notch filter is adaptively generated, ensuring that the center, width, and length of the notch filter precisely match the true shape of the artifact frequency peak. Compared to a fixed-size notch template, this adaptive notch filter achieves sufficient suppression at the artifact frequency while limiting the filtering range to the area actually occupied by the artifact energy. This avoids the problems of accidentally damaging adjacent normal frequency components due to an excessively large filtering window, or incomplete suppression due to an excessively small window, thus achieving a better balance between "artifact removal" and "detail preservation." In addition, the two-dimensional Gaussian notch filter employs a smooth attenuation transition.

[0011] In summary, this invention, through the synergy of "Laplacian pyramid multi-scale layer-by-layer processing and reconstruction" and "adaptive estimation of filter parameters based on frequency peak energy distribution," effectively suppresses filter grid artifacts while maximizing the preservation of effective image information. It combines high suppression accuracy, low detail loss, and strong adaptability, thereby improving the image quality of digital X-ray images and the accuracy of medical diagnosis.

[0012] In a preferred embodiment of this application, the adaptive determination of the target frequency position and conjugate frequency position of the filter grid artifact based on the two-dimensional spectral image includes: The spectrum center region is adaptively determined based on the size information of the two-dimensional spectrum image, and the DC component region within the spectrum center region is suppressed to obtain a spectrum image after suppressing the DC component. Based on the central region of the spectrum, frequency search band regions with adaptive widths are constructed along the horizontal axis and the vertical axis, respectively. The frequency search band regions are used to cover the frequency range in which filter grid artifacts appear. Based on the conjugate symmetry characteristics of the two-dimensional Fourier spectrum, the conjugate consistency of the candidate frequency regions within the frequency search band is verified to obtain the target frequency position and conjugate frequency position of the filter grid artifact.

[0013] By adopting the above technical solution, the adaptive determination method of the target frequency position and conjugate frequency position is further defined. First, the central region of the spectrum is adaptively determined based on the spectrum size information, and the DC component within it is suppressed. This avoids interference from strong DC energy on subsequent frequency peak detection, making the artifact frequency peaks more prominent in the spectrum. Second, an adaptive-width frequency search band is constructed along the horizontal and vertical directions based on the central region of the spectrum. This allows the search range to be dynamically adjusted according to the image size and the frequency distribution of the filter grid. This covers the frequency range where artifacts actually appear while reducing the search overhead in irrelevant areas, improving localization efficiency and robustness. Finally, utilizing the inherent conjugate symmetry of the two-dimensional Fourier spectrum, conjugate consistency verification is performed on the candidate frequencies within the search band. This eliminates misjudgments of random noise peaks and normal structure peaks, retaining only the true artifact frequencies that appear in pairs.

[0014] In a preferred embodiment of this application: the step of constructing a local analysis region centered on the target frequency position, and performing spectral projection along both the frequency peak width direction and the frequency peak length direction to obtain a one-dimensional energy distribution curve includes: Centered on the target frequency position, a preset width d is extended on both sides of the target frequency position to obtain local analysis regions along the horizontal and vertical directions; By projecting the spectral energy onto the local analysis region, we obtain the horizontal and vertical energy distribution curves that constitute the one-dimensional energy distribution curve. The calculation formula is as follows:

[0015] In the formula, S(i,j) represents the spectrum image after suppressing the DC component. This represents the energy distribution curve in the horizontal direction. The curve represents the energy distribution in the vertical direction, (x,y) represents the target frequency position, and d represents the preset width.

[0016] By employing the above technical solution, a local analysis region is constructed by extending a preset width 'd' to both sides of the target frequency location. This focuses the analysis on the artifact frequency peak and its neighborhood, avoiding irrelevant energy interference introduced by full-spectrum statistics. Within this region, the spectral energy is cumulatively projected along both the horizontal and vertical directions, compressing the energy distribution of the two-dimensional spectral peak into two one-dimensional curves. This allows for a clear and quantitative expression of the energy attenuation patterns of the artifact frequency peak in both the width and length directions. Furthermore, the cumulative window width 'd' can be preset as needed, balancing the sufficiency and locality of energy statistics.

[0017] In a preferred embodiment of this application, the step of calculating the spectral width parameter and spectral length parameter based on the one-dimensional energy distribution curve includes: The one-dimensional energy distribution curve is smoothed using a Gaussian smoothing algorithm to obtain a smoothed one-dimensional energy distribution curve. Determine the peak energy of the target frequency peak The median of the smoothed one-dimensional energy distribution curve is taken as the background energy B. The boundary detection threshold is calculated based on the difference between the peak energy and the background energy, using the following formula:

[0018] In the formula, T represents the boundary detection threshold. This represents a preset scaling factor, and satisfies 0 < <1; Starting from the position of the target frequency peak, search point by point along the smoothed one-dimensional energy distribution curve to both sides of the target frequency peak. The position where the energy first drops below the boundary detection threshold is determined as the frequency peak boundary position, thus obtaining the left boundary L and right boundary R in the horizontal direction, and the upper boundary U and lower boundary D in the vertical direction.

[0019] By employing the above technical solution, Gaussian smoothing is used to suppress random fluctuations and spikes on the one-dimensional energy distribution curve. The median of the smoothed curve is used as the background energy B, making the background estimation more representative and less affected by a few high-energy points. The calculation and adjustment of the boundary detection threshold allows it to adaptively change with the actual difference between the peak and the background, rather than using a fixed threshold, thus adapting to energy distribution differences at different scales and with different artifact intensities. Finally, starting from the peak, a point-by-point search is performed to both sides, and the point where the energy first falls below the threshold is determined as the boundary, yielding the left and right horizontal boundaries and the top and bottom vertical boundaries.

[0020] In a preferred embodiment of this application, the step of calculating the spectral width parameter and spectral length parameter based on the one-dimensional energy distribution curve further includes: The spectral width parameter is calculated based on the left boundary L and the right boundary R, ​​and the spectral length parameter is calculated based on the upper boundary U and the lower boundary D. The calculation formula is as follows: In the formula, Indicates the spectrum width parameter. The spectrum length parameter is used to characterize the actual influence range of the filter grid artifact in the frequency domain, along with the spectrum width parameter.

[0021] By employing the above technical solution, the actual extension range of the artifact frequency peak in both the horizontal and vertical directions is transformed into clear quantitative indicators. Since the spectral width and spectral length parameters are directly derived from the adaptive boundary detection results, they can accurately reflect the influence range occupied by the artifact in the frequency domain, without relying on a preset fixed size. Using these parameters to characterize the artifact frequency peak shape ensures that the subsequently constructed two-dimensional Gaussian notch filter precisely matches the actual range of the artifact in both width and length: achieving sufficient suppression at the points of artifact energy concentration and maintaining a small filtering range near normal frequencies to avoid detail loss caused by over-filtering.

[0022] In a preferred embodiment of this application: the Laplacian pyramid decomposition of the digitized X-ray image yields multiple Laplacian pyramid layer images at different scales, including: The Gaussian pyramid image of layer i is obtained by convolution and downsampling based on a Gaussian filter kernel, using the following formula:

[0023] In the formula, represents the Gaussian pyramid image of the i-th layer, represents the Gaussian pyramid image of the (i+1)-th layer, H represents the Gaussian filter kernel, * represents convolution calculation, and Down() represents downsampling operation; Based on the repeated execution of the convolution calculation and the downsampling operation, a Gaussian pyramid structure is obtained. , where n represents the number of pyramid layers; Upsampling is performed on the (i+1)th layer of the Gaussian pyramid image to obtain an approximate image restored to the size of the i-th layer. The calculation formula is as follows:

[0024] In the formula, Expand() represents the upsampling operation. This represents the approximate image; The image of the i-th layer of the Laplacian pyramid is calculated based on the difference between the image of the i-th layer of the Gaussian pyramid and the approximate image. The calculation formula is as follows:

[0025] In the formula, This represents the image of the i-th Laplacian pyramid layer.

[0026] By employing the above technical solution, a Gaussian filter kernel is convolved and downsampled on the Gaussian pyramid image to generate a Gaussian pyramid structure with decreasing resolution level by level. The Gaussian filter effectively suppresses high frequencies before downsampling to prevent frequency aliasing and ensure the image quality of each layer. Then, the upper layer image is upsampled to obtain an approximate image restored to the current layer size. The difference between the current layer image and the approximate image is used to construct a Laplacian layer image, thereby decomposing the image into a series of bandpass components according to frequency bands. Simultaneously, the Laplacian pyramid exhibits complete reversibility; the artifact removal results of each layer can be losslessly fused after reconstruction, and the normal structure is complementaryly restored during reconstruction.

[0027] In a preferred embodiment, after acquiring the digital X-ray image to be processed, this application further includes: The digital X-ray image to be processed is preprocessed based on edge cropping and mirror filling.

[0028] By employing the above technical solutions, edge cropping can remove invalid regions, black borders, or abnormal pixels that may exist at the edges of digital X-ray images, preventing these regions from introducing false high-frequency components during frequency domain transformation and interfering with the detection of artifact frequency peaks. Mirror padding extends the image boundaries, ensuring the image size meets the requirements of pyramid decomposition and Fourier transform. Simultaneously, mirroring maintains the continuity at the boundaries, suppressing boundary discontinuities, spectral leakage, and ringing effects caused by direct zero-padding or truncation during Fourier transform.

[0029] Secondly, the objective of this invention is achieved through the following technical solution: A digital X-ray image filter grid artifact suppression system includes: The image acquisition module is used to acquire the digital X-ray image to be processed; The pyramid decomposition module is used to perform Laplacian pyramid decomposition on the digitized X-ray image to obtain multiple Laplacian pyramid layer images at different scales. The artifact suppression processing module is used to perform artifact suppression processing for each of the Laplacian pyramid layer images. The artifact suppression processing includes: A two-dimensional Fourier transform is performed on the Laplace pyramid layer image to obtain a two-dimensional spectral image; Based on the two-dimensional spectrum image, the target frequency position and conjugate frequency position of the filter grid artifact are adaptively determined. A local analysis region is constructed with the target frequency position as the center, and spectral projection is performed along the frequency peak width direction and the frequency peak length direction to obtain a one-dimensional energy distribution curve. Based on the one-dimensional energy distribution curve, the spectral width parameter and spectral length parameter are calculated. A two-dimensional Gaussian notch filter is adaptively generated based on the target frequency position, the conjugate frequency position, the spectral width parameter, and the spectral length parameter. The two-dimensional spectral image is filtered using the two-dimensional Gaussian notch filter, and an inverse Fourier transform is performed on the filtered spectrum to obtain the artifact removal result of the current scale layer. The image reconstruction module is used to perform Laplacian pyramid reconstruction on the artifact removal results of all scale layers to obtain a digital X-ray image with suppressed grid artifacts.

[0030] Thirdly, the objective of this invention is achieved through the following technical solution: A computer device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the above-described method for suppressing artifacts in a digital X-ray image filter grid.

[0031] Fourthly, the objective of this invention is achieved through the following technical solution: A computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the above-described method for suppressing artifacts in a digital X-ray image filter grid.

[0032] In summary, this application includes at least one of the following beneficial technical effects: 1. By employing Laplacian pyramid decomposition, filter grate artifacts spanning a wide frequency band are separated into multiple scale layers according to spatial frequency. This allows the originally overlapping coarse and fine stripes to appear as concentrated and easily identifiable frequency peaks in their respective layers, overcoming the limitation of single-scale overall frequency domain processing in handling different coarse and fine stripes. Within each layer, a two-dimensional Fourier transform is first used to obtain the spectrum. Then, the target frequency position and conjugate frequency position of the artifact are adaptively determined, allowing the filtering position to be dynamically determined according to the image content and scale rather than being manually preset. Subsequently, a one-dimensional energy distribution curve is obtained through local spectral projection, based on which the spectral width and length parameters are calculated. These parameters, along with the frequency position, are used to adaptively generate a two-dimensional Gaussian notch filter, ensuring that the notch filter shape accurately matches the true range of the artifact frequency peaks. This effectively suppresses artifacts while avoiding damage to adjacent normal frequency components. Finally, pyramid reconstruction is used to fuse the artifact removal results from each layer, allowing for complementary restoration of the normal structure. Thus, this method effectively suppresses filter grate artifacts while maximizing the preservation of image details, combining high precision, low damage, and strong adaptability. 2. It achieves accurate frequency positioning that does not require manual preset and adapts to changes in acquisition conditions and filter grid specifications, providing a reliable basis for subsequent notch filter construction and overcoming the defect that fixed notch filter positions cannot adapt to different specifications. Attached Figure Description

[0033] Figure 1 This is a flowchart of a method for suppressing grating artifacts in digital X-ray images according to an embodiment of this application; Figure 2 This is a flowchart of step S3 in a method for suppressing filter grid artifacts in a digital X-ray image according to an embodiment of this application; Figure 3 This is a schematic diagram of a device according to one embodiment of this application. Detailed Implementation

[0034] The present application will be further described in detail below with reference to the accompanying drawings.

[0035] In one embodiment, such as Figure 1 As shown, this application discloses a method for suppressing filter grid artifacts in digital X-ray images, which specifically includes the following steps: S1: Acquire the digital X-ray image to be processed.

[0036] In this embodiment, the digital X-ray image is the original image acquired by the detector of the digital radiography (DR) system and output after conventional readout and correction.

[0037] Specifically, when a digital X-ray image is stored in DICOM format, the pixel data array and row and column size information can be read based on the DICOM parsing library. The pixel values ​​are then converted into a floating-point grayscale matrix I(m,n) for subsequent processing based on parameters such as window width, window level, and photoelectric conversion relationship recorded in the DICOM tag. Here, m and n represent the row and column indices of the pixel, respectively. To avoid numerical overflow during the subsequent two-dimensional Fourier transform, this embodiment normalizes the grayscale matrix to the range of [0,1] or [0,255] after reading. The normalization method is extreme value normalization.

[0038] Furthermore, after acquiring the digital X-ray image to be processed, the process also includes: S11: Preprocess the digital X-ray image to be processed based on edge cropping and mirror filling.

[0039] In this embodiment, after obtaining the digital X-ray image to be processed in step S1 and before performing Laplacian pyramid decomposition in step S2, the image is first preprocessed with edge cropping and mirror filling: edge cropping is used to remove collimator occlusion black edges, invalid filling areas, or abnormal pixels at the detector edges that may exist at the edges of the digital X-ray image. In this embodiment, the effective imaging area is determined according to the grayscale distribution around the image, and the edge rows and columns with grayscale close to zero or obviously abnormal are cropped, or the edge area is cropped according to a preset cropping width, for example, several to dozens of pixels are cropped from each side of the image.

[0040] Specifically, mirror filling is used to extend the boundaries of the cropped image. In this embodiment, the pixels inside the boundary are mirrored outwards using the image boundary as the axis of symmetry. The row and column sizes of the image are filled to meet the requirements of the Laplacian pyramid decomposition, so that the row and column sizes can be divided by 2^n (n is the number of pyramid decomposition levels). This ensures that the image size strictly corresponds during the downsampling and upsampling processes and can be accurately restored during reconstruction.

[0041] The reason for using mirror filling instead of direct zero padding or direct truncation is that Fourier transform implicitly assumes the periodic extension of the image. If there is a sudden change in gray level at the left and right or upper and lower boundaries of the image, high-energy spectral leakage components and ringing effects distributed along the coordinate axes will be introduced into the spectrum, interfering with the detection of filter grid artifact frequency peaks and producing ringing stripes at the image edges after filtering.

[0042] S2: Perform Laplacian pyramid decomposition on the digital X-ray image to obtain multiple Laplacian pyramid layer images at different scales.

[0043] In this embodiment, Laplacian pyramid decomposition is used to decompose the input image into several bandpass components according to spatial frequency, so that filter gratings of different thicknesses are separated into scale layers corresponding to their spatial frequencies. The number of pyramid decomposition layers n is determined according to the image size and the filter grating line density. In this embodiment, n is taken as 3 to 5 layers, for example, n=4, so that the resolution of each layer image is 1, 1 / 2, 1 / 4, and 1 / 8 of the original image, respectively.

[0044] Specifically, step S2 includes: S21: Based on the Gaussian filter kernel, perform convolution calculation and downsampling operation on the i-th layer Gaussian pyramid image to obtain the (i+1)-th layer Gaussian pyramid image. The calculation formula is as follows:

[0045] In the formula, This represents the image of the i-th layer of the Gaussian pyramid. Let H represent the (i+1)th layer of the Gaussian pyramid image, H represent the Gaussian filter kernel, * represent convolution calculation, and Down() represent downsampling operation. In this embodiment, the Gaussian filter kernel H is a 5×5 two-dimensional Gaussian kernel with a standard deviation of 0.8 to 1.5, preferably 1.0. It is normalized so that the sum of the coefficients in the kernel is 1 to keep the overall grayscale of the image unchanged. The downsampling operation uses a downsampling factor of 2, that is, after filtering, one pixel is extracted every other row and column, so that the image size is reduced to 1 / 2 of the original size.

[0046] S22: Based on repeated execution of convolution calculation and downsampling operation, a Gaussian pyramid structure is obtained. , where n represents the number of pyramid layers.

[0047] In this embodiment, the convolution calculation and downsampling operation in step S21 are repeated from the bottom layer to the top layer. Each time it is executed, a Gaussian pyramid image with half the resolution is obtained until the preset number of pyramid decomposition layers is reached, thereby obtaining a Gaussian pyramid structure composed of multiple layers of images at different scales.

[0048] Specifically, after repeated execution, a Gaussian pyramid structure is obtained. In the formula, n represents the number of pyramid layers. For the highest resolution underlying image, The top layer image has the lowest resolution; the number of decomposition layers n is determined based on the image size and the grid line density. In this embodiment, n is taken as 3 to 5, and preferably n=4.

[0049] S23: Upsample the (i+1)th layer Gaussian pyramid image to obtain an approximate image restored to the size of the i-th layer. The calculation formula is as follows:

[0050] In the formula, Expand() represents the upsampling operation. This represents the approximate image recovered from the (i+1)th layer. In this embodiment, the upsampling operation uses an upsampling factor of 2, that is, firstly, a row and a column of zero values ​​are inserted at every other position in the row and column directions to double the image size, and then the same Gaussian filter kernel as in step S21 (multiplied by a factor of 4 to compensate for the energy loss caused by zero insertion) is used to convolve the magnified image to fill the inserted gaps, thereby obtaining a smooth and continuous approximate image.

[0051] S24: The i-th layer Laplacian pyramid image is calculated based on the difference between the i-th layer Gaussian pyramid image and the approximate image. In this embodiment, the i-th layer Gaussian pyramid image is subtracted pixel by pixel from the approximate image recovered from the upper layer. The difference obtained is the bandpass detail component lost by the upsampled approximate image in this layer, which is used as the i-th layer Laplacian pyramid image, thereby separating the components of the corresponding frequency band in the image.

[0052] Specifically, the calculation formula is as follows:

[0053] In the formula, This represents the image of the i-th Laplacian pyramid layer.

[0054] After decomposition, the structure of the Laplace pyramid is obtained. and the Gaussian residual image at the top layer. ,in It carries the finest high-frequency components of the image. Carrying relatively coarse mid-frequency components, top layer This carries the low-frequency background components of the image. For example, a 5×5 Gaussian kernel with a standard deviation of approximately 1.0 is used.

[0055] S3: For each Laplacian pyramid layer image, perform artifact suppression processing, which includes: In this embodiment, each bandpass layer of the Pyramid of Laplace is described. to The artifact suppression processing in steps S31 to S36 is performed sequentially, while the top Gaussian residual image... It only carries low-frequency background components and does not contain periodic grid artifacts, so it does not participate in artifact suppression processing and is used directly in subsequent reconstruction; the processing of each layer is independent of each other.

[0056] Specifically, the current image to be processed is denoted as... (k=0, 1, ..., n-1), such as Figure 2 As shown, the following steps S31 to S36 are executed sequentially to obtain the artifact removal result for this layer. .

[0057] S31: Perform a two-dimensional Fourier transform on the Laplace pyramid layer image to obtain a two-dimensional spectral image.

[0058] In this embodiment, for the current layer image A two-dimensional discrete Fourier transform (2D-DFT) is performed and implemented using the fast Fourier transform (FFT) algorithm to obtain the complex spectrum F(u,v). To facilitate subsequent frequency peak detection near the center of the spectrum, quadrant shifting (fftshift) is performed on F(u,v) to move the zero-frequency component to the center of the spectrum.

[0059] Specifically, the amplitude S(u,v)=|F(u,v)| of the spectrum is used as a two-dimensional spectrum image for subsequent DC suppression, frequency search and energy projection, while the complex spectrum F(u,v) itself is retained and used for the inverse Fourier transform after multiplication with the notch filter in step S36. In this embodiment, in order to take into account both the display and detection of weak artifact frequency peaks, the logarithmic transform log(1+S(u,v)) of the amplitude spectrum can be taken for the visualization and threshold analysis of frequency peaks, but the original complex spectrum is still used for filtering.

[0060] S32: Based on the two-dimensional spectrum image, adaptively determine the target frequency position and conjugate frequency position of the filter grid artifact.

[0061] In this embodiment, after processing in step S32, the target frequency position (x, y) of the filter grid artifact and its conjugate frequency position (x', y') symmetrical about the center of the spectrum are obtained. The two satisfy a central symmetry relationship about the center of the spectrum. If no frequency peak satisfying the condition is detected in the current layer spectrum, it is determined that there is no significant filter grid artifact in the current layer, and the subsequent filtering steps are skipped. This is directly used as the artifact removal result for this layer.

[0062] Specifically, step S32 includes: S321: Adaptively determine the central region of the spectrum based on the size information of the two-dimensional spectrum image, and suppress the DC component region within the central region of the spectrum to obtain the spectrum image after suppressing the DC component.

[0063] In this embodiment, after quadrant switching in step S31, the zero-frequency component (DC component) of the spectrum is located at the geometric center of the spectrum image. Let the size of the spectrum image be M×N, then the location of the spectrum center is... Since the DC component and its nearby low-frequency background energy are much higher than the energy of the filter grid artifact frequency peak, the artifact frequency peak will be overwhelmed if not suppressed. Therefore, it is necessary to suppress the DC component region in the center of the spectrum first.

[0064] Specifically, the size of the spectral center region is adaptively determined based on the spectral image size; in this embodiment, it is taken as... A rectangular region centered at M and with side length min(M, N) of a certain proportion (e.g., 1 / 8 to 1 / 4 of this proportion) is designated as the spectrum center region; within this region, ... With center and radius as The circular region is taken as the DC component region in this embodiment. For a range of 3 to 10 pixels, the spectral amplitude within the circular region is set to zero or multiplied by an attenuation coefficient close to 0, thereby obtaining the spectral image S(i,j) after suppressing the DC component.

[0065] S322: Based on the center region of the spectrum, an adaptive width frequency search band region is constructed along the horizontal axis and the vertical axis respectively. The frequency search band region is used to cover the frequency range in which filter grid artifacts appear.

[0066] Specifically, the frequency search band along the horizontal axis is a band with an adaptive width along the vertical direction and the horizontal axis as its central axis. The rectangular strip has a frequency search band along the vertical axis with an adaptive width along the horizontal direction, centered on the vertical axis. Rectangular strip; width Based on the adaptive determination of the spectrum size, this embodiment takes... The search band is a certain percentage (e.g., 3% to 8%) of the dimension corresponding to the spectrum and rounded up to an odd number of pixels. The search band extends along the axial direction from the boundary of the DC component region to the boundary of the central region of the spectrum, thereby covering the frequency range in which filter grid artifacts may occur.

[0067] S323: Based on the conjugate symmetry characteristics of the two-dimensional Fourier spectrum, the conjugate consistency of the candidate frequency regions within the frequency search band is verified to obtain the target frequency position and conjugate frequency position of the filter grid artifact.

[0068] In this embodiment, for a real-valued input image, its two-dimensional Fourier spectrum has conjugate symmetry characteristics, that is, the amplitude spectrum is centrally symmetric about the center of the spectrum. Therefore, the true filter grid artifact frequency peaks must appear in pairs at two positions symmetric about the center of the spectrum. Accordingly, by first detecting candidate frequency peaks in the frequency search band region, and then using this symmetry to verify the candidate peaks, isolated peaks caused by random noise and frequency components of asymmetrical anatomical structures can be effectively eliminated.

[0069] Specifically, within the frequency search band of the spectrum image S(i,j) after suppressing the DC component, several candidate frequency peaks are detected in descending order of amplitude to obtain the position and amplitude of each candidate peak; for each candidate peak position (i,j), its symmetrical position about the center of the spectrum is calculated. The system then detects whether there are candidate peaks with similar amplitudes within the neighborhood of the symmetrical position (e.g., within ±2 pixels). If such peaks exist and the ratio of their amplitudes is within a preset range, the system is deemed to have passed the conjugate consistency verification. The preset range is, for example, 0.8 to 1.25. The position of the candidate peak that has passed the verification and has the largest amplitude is determined as the target frequency position (x, y), and its symmetrical position is determined as the conjugate frequency position (x', y').

[0070] S33: Construct a local analysis region centered on the target frequency position, and perform spectral projection along the frequency peak width direction and the frequency peak length direction respectively to obtain a one-dimensional energy distribution curve.

[0071] In this embodiment, the frequency peak width direction corresponds to the horizontal axis direction of the spectrum, and the frequency peak length direction corresponds to the vertical axis direction of the spectrum; the two are orthogonal to each other.

[0072] Specifically, step S33 includes: S331: Taking the target frequency position as the center, expand the preset width d on both sides of the target frequency position to obtain the local analysis area along the horizontal and vertical directions.

[0073] Specifically, the local analysis region includes a first analysis region extending horizontally and a second analysis region extending vertically: the first analysis region is used to statistically analyze the energy distribution in the horizontal direction, which is limited to the range [y−d, y+d] in the vertical direction and covers the search range in the target frequency peak width direction in the horizontal direction; the second analysis region is used to statistically analyze the energy distribution in the vertical direction, which is limited to the range [x−d, x+d] in the horizontal direction and covers the search range in the target frequency peak length direction in the vertical direction; the preset width d is used to determine the window width for energy accumulation. For example, d is taken as 3 to 10 pixels.

[0074] S332: Project the spectral energy onto the local analysis region to obtain the horizontal and vertical energy distribution curves that constitute the one-dimensional energy distribution curve.

[0075] In this embodiment, the spectral amplitudes are accumulated column by column along the vertical direction in the first analysis region to obtain a one-dimensional curve reflecting the energy distribution in the horizontal direction. By summing the spectral amplitudes row by row along the horizontal direction in the second analysis region, a one-dimensional curve reflecting the energy distribution in the vertical direction is obtained. By compressing the energy of the two-dimensional spectral peak into two one-dimensional curves along orthogonal directions, the energy attenuation law of the frequency peak in both directions is clearly expressed, and the computational complexity of subsequent boundary detection is significantly reduced while suppressing the influence of single-point noise.

[0076] Specifically, the formula for calculating the spectral energy projection is as follows: In the formula, S(i,j) represents the spectrum image after suppressing the DC component. This represents the energy distribution curve in the horizontal direction. The curve represents the energy distribution in the vertical direction, (x,y) represents the target frequency position, and d represents the preset width.

[0077] S34: Based on the one-dimensional energy distribution curve, the spectral width parameter and spectral length parameter are calculated.

[0078] In this embodiment, the spectral width parameter and the spectral length parameter are used to characterize the actual influence range of the filter grid artifact in the frequency domain.

[0079] Specifically, step S34 includes: S341: The one-dimensional energy distribution curve is smoothed using the Gaussian smoothing algorithm to obtain the smoothed one-dimensional energy distribution curve.

[0080] In this embodiment, since the one-dimensional energy distribution curve obtained in step S332 inevitably has local fluctuations and spikes caused by random noise, if the boundary search is performed directly on the original curve, it is easy to misjudge the boundary due to single-point noise. Therefore, the Gaussian smoothing algorithm is first used to smooth the one-dimensional energy distribution curve, so as to suppress random fluctuations while preserving the overall shape of the frequency peak and the energy decay trend.

[0081] Specifically, Gaussian smoothing is achieved by convolving a one-dimensional energy distribution curve with a one-dimensional Gaussian kernel, and its calculation formula is as follows: In the formula, P(i) represents the original projection curve. Let g(k) represent the smoothed curve, where g(k) is a one-dimensional Gaussian kernel and satisfies The size and standard deviation of the Gaussian kernel are determined based on the noise level of the curve. In this embodiment, the Gaussian kernel length is set to 5 to 11 sampling points, and the standard deviation is set to... The value is between 1.0 and 2.0, for example, taking a core length of 7. =1.5.

[0082] S342: Determine the peak energy of the target frequency peak. The median of the smoothed one-dimensional energy distribution curve is taken as the background energy B.

[0083] In this embodiment, on the smoothed one-dimensional energy distribution curve, taking the curve sampling point corresponding to the target frequency position obtained in step S32 as the center, the maximum value of the curve in its neighborhood is taken as the peak energy of the target frequency peak. In this embodiment, the median of the smoothed one-dimensional energy distribution curve is used as the background energy B. Compared with using the mean, using the median can effectively reduce the impact of the target frequency peak and other local abnormal high values ​​on the background estimation results.

[0084] S343: The boundary detection threshold is calculated based on the difference between the peak energy and the background energy. The calculation formula is as follows:

[0085] In the formula, T represents the boundary detection threshold. This represents a preset scaling factor, and satisfies 0 < <1.

[0086] In this embodiment, in order to enable the boundary detection threshold to adapt to the energy distribution differences under different scale layers and different artifact intensities, rather than using a fixed threshold, the boundary detection threshold is set between the peak energy Pmax and the background energy B, so that it is both higher than the background baseline to exclude background fluctuations and lower than the peak value to accurately define the extension range of the frequency peak.

[0087] S344: Starting from the position of the target frequency peak, search point by point along the smoothed one-dimensional energy distribution curve to both sides of the target frequency peak. The position where the energy first drops below the boundary detection threshold is determined as the frequency peak boundary position, thus obtaining the left boundary L and right boundary R in the horizontal direction, and the upper boundary U and lower boundary D in the vertical direction.

[0088] Specifically, regarding the horizontal energy distribution curve The above bidirectional search yields the left boundary L and right boundary R in the horizontal direction; the energy distribution curve in the vertical direction is then obtained. Perform the same bidirectional search to obtain the upper boundary U and lower boundary D in the vertical direction. To prevent misjudgment caused by single-point drops in the curve due to noise, a continuity constraint can be added, that is, the energy of two or more consecutive sampling points from the boundary position must be lower than the threshold T before it is determined to be a valid boundary. If the search to a certain side reaches the boundary of the local analysis area and still does not fall below the threshold, then the boundary of the analysis area is taken as the boundary position of that side.

[0089] Furthermore, based on the one-dimensional energy distribution curve, the spectral width parameter and spectral length parameter are calculated, including: S3400: The spectral width parameter is calculated based on the left boundary L and the right boundary R, ​​and the spectral length parameter is calculated based on the upper boundary U and the lower boundary D. The calculation formula is as follows: In the formula, Indicates the spectrum width parameter. The spectral length parameter, along with the spectral width parameter, is used to characterize the actual influence range of the filter grid artifact in the frequency domain; where L, R, U, and D are all represented by spectral pixel coordinates. Take the total number of sampling points spanning between the two corresponding boundaries and ensure that it is a positive value.

[0090] S35: Adaptively generate a two-dimensional Gaussian notch filter based on the target frequency position, conjugate frequency position, spectral width parameter, and spectral length parameter.

[0091] In this embodiment, the target frequency position (x, y) and the conjugate frequency position (x', y') are used as the center, and the spectral width parameter is used as the parameter. With spectral length parameter The broadening of the Gaussian notch filter in two directions is determined, generating a two-dimensional Gaussian notch filter that forms frequency suppression regions at the target frequency and conjugate frequency positions. The notch filter is set to a value close to 0 at both centers to fully suppress artifact frequency components, and a value close to 1 at points far from the center to preserve normal frequency components. It also has a Gaussian smooth transition to avoid the ringing effect introduced by the ideal notch filter at the boundary.

[0092] For example, one implementable expression for a two-dimensional Gaussian notch filter is:

[0093] in Let be the standard deviations of the notch filter in the horizontal and vertical directions, respectively. In this embodiment, let be... With spectral width parameter Proportional With spectral length parameter Proportional, for example, take .

[0094] S36: Use a two-dimensional Gaussian notch filter to filter the two-dimensional spectral image, and perform an inverse Fourier transform on the filtered spectrum to obtain the artifact removal result of the current scale layer.

[0095] In this embodiment, the complex spectrum F(u,v) retained in step S31 is compared with the two-dimensional Gaussian notch filter generated in step S35. Multiplying point by point, we obtain the filtered spectrum. This suppresses the artifact frequency components at the target frequency position and the conjugate frequency position.

[0096] Specifically, the filtered spectrum F'(u,v) is first subjected to inverse quadrant shift (ifftshift) to restore the original arrangement of the spectrum, and then a two-dimensional inverse fast Fourier transform (IFFT) is performed. The real part of the result is taken as the artifact removal result of the current layer. .

[0097] S4: Perform Laplacian pyramid reconstruction on the artifact removal results of all scale layers to obtain a digital X-ray image with suppressed grid artifacts.

[0098] In this embodiment, the artifact removal results of each layer obtained after step S3 are... Compared with the unprocessed top-level Gaussian residual image By fusing layers from top to bottom according to the Laplace pyramid reconstruction rule, a digital X-ray image with suppressed grid artifacts is obtained.

[0099] Specifically, the reconstruction process is an iterative process of layer-by-layer upsampling and addition, that is, letting Calculate sequentially from top to bottom. (k=n-1, n-2, ..., 0), where Expand(·) represents the same upsampling operation as during decomposition, and the final result is... This is the digital X-ray image with the same size as the original image and the filter grid artifacts suppressed; since edge cropping and mirror filling preprocessing were performed in step S1, the filling area is cropped here by the reverse operation of the preprocessing and restored to the original image size to obtain the final output image.

[0100] In one embodiment, before filtering the two-dimensional spectral image using a two-dimensional Gaussian notch filter and performing an inverse Fourier transform on the filtered spectrum to obtain the artifact removal result for the current scale layer, the method further includes: S3601: Using the target frequency position and the conjugate frequency position as the center, extract the residual energy region in the filtered spectrum, and calculate the residual energy value of the residual energy region and the residual difference value between the adjacent background energy.

[0101] Specifically, the residual energy region is defined by its size and spectral width, centered at the target frequency and its conjugate frequency, respectively. and spectrum length parameter The corresponding rectangular region; the filtered spectral amplitude within the residual energy region is statistically analyzed, and the average or maximum amplitude within the region is taken as the residual energy value. Simultaneously, the statistical spectral amplitude is measured in the annular or strip neighborhood adjacent to the residual energy region, and the average value is taken as the neighboring background energy. Then, the residual difference between the two is calculated. In this embodiment, the residual difference is taken as the ratio of the two. .

[0102] S3602: Compare the residual difference value with the preset residual threshold. If the residual difference value is outside the residual threshold range, it is determined that there is artifact residue, and the notch filter secondary correction command is triggered.

[0103] For example, the residual threshold is set to 1.1 to 1.5, preferably 1.2. If ΔE is within the residual threshold range, it indicates that the residual energy is close to the background level, and it is determined that there is no artifact residue. The secondary correction process is skipped, and the current filtered spectrum is directly used as the spectrum to be inversely transformed into the inverse Fourier transform in step S36.

[0104] S3603: According to the notch filter secondary correction instruction, the notch depth and notch bandwidth of the two-dimensional Gaussian notch filter are adaptively increased according to the residual difference value to generate a secondary Gaussian notch filter. The secondary Gaussian notch filter is then used to re-filter the two-dimensional spectrum image, and the re-filtered spectrum is used as the spectrum to be inversely transformed.

[0105] Specifically, the notch depth and bandwidth of the two-dimensional Gaussian notch filter are adaptively increased according to the residual difference value to generate a secondary Gaussian notch filter: the increase in notch bandwidth is achieved by increasing the standard deviation of the notch filter, making the secondary standard deviation... , , where β is a preset bandwidth adjustment coefficient, for example, 0.2 to 0.5; The residual threshold is used. The notch depth is increased by increasing the suppression intensity at the center of the notch filter, even if the value of the second-order Gaussian notch filter at the center approaches 0. For example, the two Gaussian terms in the notch filter expression in step S35 are multiplied by a depth enhancement coefficient greater than 1, and the result is truncated to 0. For instance, the depth enhancement coefficient is between 1.1 and 1.5. Subsequently, the two-dimensional spectral image is re-filtered using the second-order Gaussian notch filter, and the re-filtered spectrum is output as the spectrum to be inversely transformed in step S36.

[0106] S3604: If the residual difference value is within the residual threshold range, it is determined that there is no artifact residue, and the secondary correction process is skipped.

[0107] In this embodiment, the secondary correction process is executed at most twice. When the upper limit of the number of executions is reached, the correction is stopped regardless of whether the residual difference value is still outside the residual threshold range, and the current spectrum is used as the spectrum to be inversely transformed, thereby achieving a balance between suppressing residual artifacts and protecting normal frequency components.

[0108] In one embodiment, prior to step S3, i.e., before performing artifact suppression processing for each Laplacian pyramid layer image, the method further includes: S301: Perform frequency domain energy statistics on the Laplacian pyramid layer image and detect candidate peak energy states within the frequency search band region. The candidate peak energy states include the first energy state and the second energy state.

[0109] Specifically, after performing a two-dimensional Fourier transform and DC suppression on the current layer image, the candidate frequency peak with the largest amplitude is detected within the constructed frequency search band region, and its amplitude is taken as the candidate peak energy. The median of the curves within the search band or the median of the spectral amplitude within the search band is taken as the background energy. ; Calculate the ratio of candidate peak energy to background energy. and compare it with the preset hierarchical activation threshold. A comparison is made to determine whether the current layer is in the first or second energy state. (Hierarchical activation threshold) In this embodiment, the value is taken as 1.5 to 3.0, preferably 2.0.

[0110] S302: Based on the candidate peak energy state, adaptively determine the processing mode of the current scale layer. The processing mode includes the artifact processing mode corresponding to the first energy state and triggering artifact suppression processing, and the preservation mode corresponding to the second energy state and skipping artifact suppression processing. The first energy state refers to the ratio of candidate peak energy to background energy being greater than a preset level activation threshold, and the second energy state refers to the ratio of candidate peak energy to background energy not being greater than a preset level activation threshold.

[0111] Specifically, the first energy state indicates that there are candidate frequency peaks significantly higher than the background within the frequency search band region of the current layer, meaning that the layer contains obvious grid artifacts. At this time, the current layer is determined to enter the artifact processing mode, triggering and executing the artifact suppression processing in steps S31 to S36. The second energy state indicates that the current layer does not contain significant grid artifact frequency peaks. At this time, the current layer is determined to enter the preservation mode, skipping the artifact suppression processing, and directly using the current layer image as the artifact removal result for that layer in the Laplacian pyramid reconstruction in step S4.

[0112] In this embodiment, the above-mentioned hierarchical adaptive activation determination is executed independently for each bandpass layer of the Laplace pyramid, and each layer independently determines the processing mode according to its own energy state.

[0113] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0114] In one embodiment, a digital X-ray image grid artifact suppression system is provided, which corresponds to a digital X-ray image grid artifact suppression method described in the above embodiment.

[0115] A digital X-ray image filter grid artifact suppression system includes an image acquisition module, a pyramid decomposition module, an artifact suppression processing module, and an image reconstruction module. Detailed descriptions of each functional module are as follows: The image acquisition module is used to acquire the digital X-ray image to be processed; The pyramid decomposition module is used to decompose digital X-ray images into Laplacian pyramids to obtain Laplacian pyramid layer images at multiple different scales. The artifact suppression module performs artifact suppression processing on each Laplacian pyramid layer image. The artifact suppression processing includes: A two-dimensional Fourier transform is performed on the Laplace pyramid layer image to obtain a two-dimensional spectral image; Based on the two-dimensional spectrum image, the target frequency position and conjugate frequency position of the filter grid artifact are adaptively determined. A local analysis region is constructed with the target frequency position as the center, and the spectrum is projected along the frequency peak width direction and the frequency peak length direction to obtain a one-dimensional energy distribution curve. Based on the one-dimensional energy distribution curve, the spectral width parameter and spectral length parameter are calculated. A two-dimensional Gaussian notch filter is adaptively generated based on the target frequency position, conjugate frequency position, spectral width parameter, and spectral length parameter. A two-dimensional Gaussian notch filter is used to filter the two-dimensional spectral image, and an inverse Fourier transform is performed on the filtered spectrum to obtain the artifact removal result of the current scale layer. The image reconstruction module is used to perform Laplacian pyramid reconstruction on the artifact removal results of all scale layers to obtain a digital X-ray image with suppressed grid artifacts.

[0116] For specific limitations regarding a digital X-ray image grid artifact suppression system, please refer to the limitations of a digital X-ray image grid artifact suppression method mentioned above, which will not be repeated here. Each module in the aforementioned digital X-ray image grid artifact suppression system can be implemented entirely or partially through software, hardware, or a combination thereof. Each module can be embedded in or independent of the processor in a computer device in hardware form, or it can be stored in the memory of a computer device in software form, so that the processor can call and execute the corresponding operations of each module.

[0117] In one embodiment, a computer device is provided, which may be a server, and its internal structure diagram may be as follows: Figure 3 As shown, the computer device includes a processor, memory, network interface, and database connected via a system bus. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and database. The internal memory provides the environment for the operation of the operating system and computer programs in the non-volatile storage media. The database stores Laplacian pyramid layer images at different scales, two-dimensional spectral images, etc. The network interface communicates with external terminals via a network connection. When executed by the processor, the computer program implements a method for suppressing artifacts in digital X-ray image grids.

[0118] In one embodiment, a computer device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the digital X-ray image grid artifact suppression method described above.

[0119] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the steps of the digital X-ray image filter grid artifact suppression method described above.

[0120] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in a variety of forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0121] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is used as an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above.

[0122] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A method for suppressing grid artifacts in digital X-ray images, characterized in that, include: Acquire the digital X-ray image to be processed; The digitized X-ray image was decomposed into a Laplacian pyramid to obtain multiple Laplacian pyramid layer images at different scales; For each of the Laplacian pyramid layer images, artifact suppression processing is performed, which includes: A two-dimensional Fourier transform is performed on the Laplace pyramid layer image to obtain a two-dimensional spectrum image, and the DC component region of the two-dimensional spectrum image is suppressed to obtain a spectrum image after suppressing the DC component. Based on the conjugate symmetry characteristics of the two-dimensional Fourier spectrum, the conjugate consistency of the spectrum image after suppressing the DC component is verified to obtain the target frequency position and conjugate frequency position of the paired filter grid artifacts. A local analysis region is constructed by extending a preset width on both sides of the target frequency position as the center. The spectral energy is projected onto the local analysis region along the frequency peak width direction and the frequency peak length direction, respectively, to obtain the horizontal energy distribution curve and the vertical energy distribution curve that constitute a one-dimensional energy distribution curve. For the horizontal and vertical energy distribution curves, a boundary detection threshold determined by the peak energy of the target frequency peak and the background energy is used to search point by point on both sides of the frequency peak. The position where the energy first drops below the boundary detection threshold is determined as the frequency peak boundary. Based on the frequency peak boundary, the spectral width parameter and spectral length parameter used to characterize the actual influence range of the filter grid artifact in the frequency domain are calculated. Based on the target frequency position, the conjugate frequency position, the spectral width parameter, and the spectral length parameter, a two-dimensional Gaussian notch filter is adaptively generated. The two-dimensional Gaussian notch filter forms a frequency suppression region at the target frequency position and the conjugate frequency position, and its broadening in the horizontal and vertical directions is positively correlated with the spectral width parameter and the spectral length parameter, respectively. The two-dimensional spectral image is filtered using the two-dimensional Gaussian notch filter, and an inverse Fourier transform is performed on the filtered spectrum to obtain the artifact removal result of the current scale layer. Laplacian pyramid reconstruction was performed on the artifact removal results of all scale layers to obtain a digital X-ray image with suppressed grid artifacts.

2. The method for suppressing grating artifacts in digital X-ray images according to claim 1, characterized in that, The method of verifying the conjugate consistency of the spectrum image after suppressing the DC component based on the conjugate symmetry characteristics of the two-dimensional Fourier spectrum, and obtaining the target frequency position and conjugate frequency position of the paired filter grid artifacts, includes: The spectrum center region is adaptively determined based on the size information of the two-dimensional spectrum image, and the DC component region within the spectrum center region is suppressed to obtain a spectrum image after suppressing the DC component. Based on the central region of the spectrum, frequency search band regions with adaptive widths are constructed along the horizontal axis and the vertical axis, respectively. The frequency search band regions are used to cover the frequency range in which filter grid artifacts appear. Based on the conjugate symmetry characteristics of the two-dimensional Fourier spectrum, the conjugate consistency of the candidate frequency regions within the frequency search band is verified to obtain the target frequency position and conjugate frequency position of the filter grid artifact.

3. The method for suppressing grating artifacts in digital X-ray images according to claim 2, characterized in that, The process involves constructing a local analysis region centered on the target frequency location, and projecting the spectrum along both the frequency peak width and frequency peak length directions to obtain a one-dimensional energy distribution curve, including: Centered on the target frequency position, a preset width d is extended on both sides of the target frequency position to obtain local analysis regions along the horizontal and vertical directions; By projecting the spectral energy onto the local analysis region, we obtain the horizontal and vertical energy distribution curves that constitute the one-dimensional energy distribution curve. The calculation formula is as follows: In the formula, S(i,j) represents the spectrum image after suppressing the DC component. This represents the energy distribution curve in the horizontal direction. The curve represents the energy distribution in the vertical direction, (x,y) represents the target frequency position, and d represents the preset width.

4. The method for suppressing grating artifacts in digital X-ray images according to claim 3, characterized in that, The process of searching point-by-point towards both sides of the frequency peak using the boundary detection threshold determined by the peak energy of the target frequency peak and the background energy for the horizontal and vertical energy distribution curves, and defining the position where the energy first drops below the boundary detection threshold as the frequency peak boundary, includes: The one-dimensional energy distribution curve is smoothed using a Gaussian smoothing algorithm to obtain a smoothed one-dimensional energy distribution curve. Determine the peak energy of the target frequency peak The median of the smoothed one-dimensional energy distribution curve is taken as the background energy B. The boundary detection threshold is calculated based on the difference between the peak energy and the background energy, using the following formula: In the formula, T represents the boundary detection threshold. This represents a preset scaling factor, and satisfies 0 < <1; Starting from the position of the target frequency peak, search point by point along the smoothed one-dimensional energy distribution curve to both sides of the target frequency peak. The position where the energy first drops below the boundary detection threshold is determined as the frequency peak boundary position, thus obtaining the left boundary L and right boundary R in the horizontal direction, and the upper boundary U and lower boundary D in the vertical direction.

5. The method for suppressing grating artifacts in digital X-ray images according to claim 4, characterized in that, The spectral width and spectral length parameters, calculated based on the frequency peak boundary to characterize the actual influence range of the filter grid artifact in the frequency domain, include: The spectral width parameter is calculated based on the left boundary L and the right boundary R, ​​and the spectral length parameter is calculated based on the upper boundary U and the lower boundary D. The calculation formula is as follows: In the formula, Indicates the spectrum width parameter. The spectrum length parameter is used to characterize the actual influence range of the filter grid artifact in the frequency domain, along with the spectrum width parameter.

6. The method for suppressing grating artifacts in digital X-ray images according to claim 1, characterized in that, The process of performing Laplacian pyramid decomposition on the digitized X-ray image yields multiple Laplacian pyramid layer images at different scales, including: The Gaussian pyramid image of layer i is obtained by convolution and downsampling based on a Gaussian filter kernel, using the following formula: In the formula, This represents the image of the i-th layer of the Gaussian pyramid. This represents the (i+1)th layer of the Gaussian pyramid image, H represents the Gaussian filter kernel, * represents convolution calculation, and Down() represents downsampling operation; Based on the repeated execution of the convolution calculation and the downsampling operation, a Gaussian pyramid structure is obtained. , where n represents the number of pyramid layers; Upsampling is performed on the (i+1)th layer of the Gaussian pyramid image to obtain an approximate image restored to the size of the i-th layer. The calculation formula is as follows: In the formula, Expand() represents the upsampling operation. This represents the approximate image; The image of the i-th layer of the Laplacian pyramid is calculated based on the difference between the image of the i-th layer of the Gaussian pyramid and the approximate image. The calculation formula is as follows: In the formula, This represents the image of the i-th Laplacian pyramid layer.

7. The method for suppressing grating artifacts in digital X-ray images according to claim 1, characterized in that, After acquiring the digital X-ray image to be processed, the process further includes: The digital X-ray image to be processed is preprocessed based on edge cropping and mirror filling.

8. A digital X-ray image filter grid artifact suppression system, characterized in that, include: The image acquisition module is used to acquire the digital X-ray image to be processed; The pyramid decomposition module is used to perform Laplacian pyramid decomposition on the digitized X-ray image to obtain multiple Laplacian pyramid layer images at different scales. The artifact suppression processing module is used to perform artifact suppression processing for each of the Laplacian pyramid layer images. The artifact suppression processing includes: A two-dimensional Fourier transform is performed on the Laplace pyramid layer image to obtain a two-dimensional spectral image; Based on the two-dimensional spectrum image, the target frequency position and conjugate frequency position of the filter grid artifact are adaptively determined. A local analysis region is constructed with the target frequency position as the center, and spectral projection is performed along the frequency peak width direction and the frequency peak length direction to obtain a one-dimensional energy distribution curve. Based on the one-dimensional energy distribution curve, the spectral width parameter and spectral length parameter are calculated. A two-dimensional Gaussian notch filter is adaptively generated based on the target frequency position, the conjugate frequency position, the spectral width parameter, and the spectral length parameter. The two-dimensional spectral image is filtered using the two-dimensional Gaussian notch filter, and an inverse Fourier transform is performed on the filtered spectrum to obtain the artifact removal result of the current scale layer. The image reconstruction module is used to perform Laplacian pyramid reconstruction on the artifact removal results of all scale layers to obtain a digital X-ray image with suppressed grid artifacts.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the digital X-ray image filter grid artifact suppression method as described in any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the digital X-ray image filter grid artifact suppression method as described in any one of claims 1 to 7.

Images