DBT reconstruction method, medium and electronic device based on three-dimensional coordinate displacement

Abstract

The present invention discloses a DBT reconstruction method, medium and electronic device based on three-dimensional coordinate displacement, the method comprising: firstly obtaining multi-angle X-ray projection data and preprocessing the data; then correcting the preprocessed projection data in a three-dimensional coordinate system; then interpolating the projection data and processing it using an autoregressive source moving average separation method to obtain a tomographic image of any depth; finally enhancing the tomographic image. The present invention utilizes multi-angle X-ray projection to perform DBT reconstruction through three-dimensional displacement, and effectively focuses the structure within the reconstruction layer through multi-angle projection displacement correction and magnification standardization, making the structure within the reconstruction layer clearer, reducing artifacts in the reconstruction layer, and improving the contrast of the image of the reconstruction layer, which is more beneficial for observing breast lesions.

Description

Technical Field

[0001] The present invention belongs to the technical field of medical image processing, and in particular relates to a method for DBT reconstruction by using multi-angle X-ray projections through three-dimensional coordinate displacement. Background Art

[0002] Breast cancer is the second leading cause of cancer death among women worldwide, posing a great threat to the lives and health of all women. However, breast cancer screening has been proven to reduce breast cancer mortality, so in order to reduce breast cancer mortality, many medical companies around the world have developed many X-ray-based technologies for screening and lesion detection. Through early screening of breast cancer and timely diagnosis and treatment, 90% of patients can prolong their lives.

[0003] Digital Breast Tomosynthesis (DBT) is one of these new imaging methods, which is characterized by its ability to overcome the limitations of traditional mammography systems due to tissue overlap. In DBT, a limited number of low-dose projection images can be obtained, so the total dose to the patient is still comparable to traditional two-view mammography. It is a three-dimensional mammography technology. Instead of obtaining only two views of images like mammography, it obtains scans at limited angles to provide better depth resolution. Taking the Hologic system as an example, the X-ray tube moves in an isocentric arc of 15 degrees for a total of 15 views. As with traditional mammography, the breast is slightly compressed and the detector remains stationary. Objects at different heights in the breast are projected onto the detector. The entire data acquisition scanning process is short, and the limited-angle projection data acquired by the scan contains all the anatomical structure information of the breast tissue at different angles. Compared with the full-angle data acquisition method of the currently widely used computed tomography (CT), DBT has a smaller acquisition angle range and collects less data, so the total radiation damage is far less than conventional CT.

[0004] However, most of the current DBT imaging methods use the parallel operation mode of the X-ray tube relative to the stationary detection surface, that is, the two-dimensional characteristics, to derive the displacement formula, ignoring the three-dimensional characteristics of DBT, resulting in low contrast, poor focusing effect, blurred boundaries, and obvious metal artifacts. Summary of the invention

[0005] The present invention aims to solve one of the technical problems existing in the related art at least to a certain extent.

[0006] One object of the present invention is to provide a DBT reconstruction method based on three-dimensional coordinate displacement, which uses multi-angle X-ray projection to perform DBT reconstruction through three-dimensional displacement, so that the reconstructed focus layer structure is clearer, metal artifacts are reduced, contrast is improved, and the observation of breast lesions is convenient.

[0007] In order to achieve the above-mentioned object, the present invention provides a DBT reconstruction method based on three-dimensional coordinate displacement, comprising:

[0008] Acquire multi-angle X-ray projection data and pre-process them;

[0009] The pre-processed projection data is subjected to coordinate transformation in three-dimensional space, and then image interpolation processing is performed to obtain corrected projection data;

[0010] The corrected projection data are processed by an autoregressive source moving average separation method to obtain a tomographic image of any depth layer;

[0011] The tomographic images are quantitatively evaluated and adaptively histogram equalized to output enhanced images.

[0012] A further preferred technical solution of the present invention is to obtain multi-angle X-ray projection data and pre-process it; the specific steps include:

[0013] Performing negative log-gamma mapping on the acquired projection data;

[0014] The projection data is filtered in the Fourier domain by adjusting the window through ramp filtering.

[0015] Preferably, coordinate transformation is performed on the preprocessed projection data in three-dimensional space; the steps include:

[0016] Taking the central projection as the reference, calculate the displacement of the projection data along the x-axis and y-axis at different projection angles;

[0017] Convert the displacement of the projection data along the x-axis and the y-axis into pixel displacement in pixel units;

[0018] Calculate the change in projection size according to the displacement of the projection data along the x-axis and the y-axis respectively;

[0019] The preprocessed projection data is subjected to coordinate and size transformation and image gradient interpolation processing to obtain the displaced projection data.

[0020] Preferably, coordinate transformation is performed on the preprocessed projection data in three-dimensional space, and the specific method is as follows:

[0021] The projection angle The projection data of is expressed as follows after translation:

[0022]

[0023] in, Represents the projection data after translation; Represents the x-axis coordinate of the translated projection data, , Represents the x-axis coordinate of the original projection; Represents the y-axis coordinate of the translated projection data, , Represents the y-axis coordinate of the original projection; represents the reconstruction layer depth;

[0024] It means that when the reconstruction layer depth is , the projection angle is The displacement of the projection data along the x-axis at time ; It means that when the reconstruction layer depth is , the projection angle is The displacement of the projection data along the y-axis at time ; and By angle , reconstruction layer depth , the distance from the X-ray tube to the center of rotation , and the vertical distance from the rotation center to the detector DDR Sure;

[0025] The expression of the displacement of the projection data along the x-axis is:

[0026]

[0027] The expression of the displacement of the projection data along the y-axis is:

[0028]

[0029] The real values ​​of the displacements along the x-axis and y-axis calculated above are converted into pixel units. The conversion formula is:

[0030]

[0031]

[0032] in, Express Round down, that is, the closest and not greater than integer; and Represents the pixel displacement on the x-axis and y-axis respectively;

[0033] Calculate the change in projection size as:

[0034]

[0035]

[0036]

[0037]

[0038] The projection data varies in image size as follows:

[0039]

[0040]

[0041] After performing coordinate and size transformation on the collected multi-angle projections, the displaced projection data is obtained.

[0042] Preferably, image interpolation processing is performed on the displacement projection to obtain corrected projection data; the specific steps are:

[0043] The shifted projection data is converted into a matrix Indicates that for the matrix Each line of , each column Sum, we get and :

[0044]

[0045]

[0046] in, is a matrix In the Row, No. Elements of the column; the resulting row and column vectors , , is the number of matrix rows, is the number of matrix columns;

[0047] if , then The row has missing values; if , then Column has missing values; Find The row and column indices that are zero, determine the rows and columns corresponding to the missing values. and ,in is a set containing all null value row and column indices;

[0048] The gradient is calculated by the difference between the previous and next data, and then the missing data is filled by gradient interpolation based on the calculated gradient value. The expression of interpolation filling is:

[0049]

[0050]

[0051]

[0052]

[0053] in, , Respectively represent the pixel gradient interpolation of the adjacent upper and lower rows or adjacent left and right columns of the pixel that needs to be filled in the row and column; The index in the input projection image array is Elements of Represents the first Row, No. Elements of a column; is a loop variable used to iterate over the all-zero row index.

[0054] For missing rows or columns at the boundary, fill them with the directly adjacent non-missing data:

[0055] The top missing row ;

[0056] The missing line at the bottom ;

[0057] Leftmost missing column ;

[0058] The rightmost missing column ;

[0059] in Indicates the line; Represents a column; Represents the pixel value of a row or column; , Represents the maximum value of the row and column respectively.

[0060] Preferably, in the image interpolation processing of the displacement projection data, a limiting condition is added:

[0061]

[0062]

[0063] in, , are the rows and columns traversed in the projected image, is the number of columns, is the number of rows.

[0064] Preferably, the corrected projection data is processed by an autoregressive source moving average separation method to obtain a tomographic image of any depth layer, including the steps of:

[0065] Obtain corrected projection data and initialize the autoregressive model;

[0066] Set the number of iterations of the autoregressive model, the condition number that controls the weight matrix, and the parameters used to constrain the roots of the autoregressive model;

[0067] Based on the application of blind source separation method to separate and reconstruct arbitrary depth layers, the mean of the column vector is subtracted from each element in the column vector of the projection image matrix to make the data concentrated near zero, and the sliding average model is added to estimate the dynamic characteristics of the signal. The formula is as follows:

[0068]

[0069] in, The size is The mean vector of is a column vector of all 1s, is a row vector of all 1s;

[0070] Then, given the autoregressive order, the time period length , calculate the autocorrelation matrix :

[0071]

[0072] in, Representation Matrix From column 1 to column The column vector of the concatenation of columns; Representation Matrix No. Column to The column vector of the concatenation of columns; represents the vector outer product; Indicates optional matrix;

[0073] Then calculate the average value of each specified submatrix region and its transpose, and assign this average value back to the atomic matrix; according to the input matrix and the current estimate, using an approximate joint diagonalization algorithm with uniform weights to unmix the matrix , the diagonalized Matrix; Diagonalize using the approximate joint diagonalization algorithm with non-uniform weights The matrix is ​​optimized within the number of iterations to obtain the separated signal matrix , the expression is:

[0074]

[0075] Among them, the matrix is the coefficient matrix for processing the diagonal elements of the matrix; For the matrix The transposed matrix of

[0076] Finally, the tomographic image is reconstructed by separating the signals and processing the data with autoregressive moving average.

[0077] Preferably, the tomographic image is quantitatively evaluated and adaptively histogram equalized to enhance the image, comprising the steps of:

[0078] Perform autoregressive model analysis on the separated signal source and calculate the cumulative value of each coefficient in the correlation matrix; and calculate ISR matrix;

[0079] The calculated ISR Perform heatmap visualization on the matrix and ISR visualization evaluation on the tomographic images;

[0080] For images with poor evaluation, a mask matrix is ​​defined, and for each pixel in the image, the weighted sum of the pixels in its surrounding area and the mask matrix is ​​calculated;

[0081] Calculate the local contrast mean of the image and compare it with the threshold. If the local contrast mean is higher than the threshold, use a smaller area for enhancement, otherwise use a larger area for enhancement.

[0082] The histogram is calculated by counting the number of pixels with a specific grayscale value in the window, and the histogram is accumulated and summed to obtain a cumulative histogram representing the cumulative distribution of grayscale levels. This is used to map each pixel value in the window and output an enhanced image.

[0083] Another aspect of the present invention provides a non-transitory computer-readable storage medium having computer instructions stored thereon, the computer instructions enabling a computer to execute the above-mentioned DBT reconstruction method based on three-dimensional coordinate displacement.

[0084] On the other hand, the present invention provides an electronic device, comprising: a processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory communicate with each other through the communication bus, and the processor calls the logic instructions in the memory to execute the above-mentioned DBT reconstruction method based on three-dimensional coordinate displacement.

[0085] Beneficial effects: The DBT reconstruction method based on three-dimensional coordinate displacement provided by the present invention utilizes multi-angle X-ray projection to perform DBT reconstruction through three-dimensional displacement. Through different projection positions based on the isocenter model, all structures on the reconstruction plane are focused, and multi-angle projection displacement correction and magnification standardization are realized, so that the reconstructed focus layer structure is clearer, metal artifacts are reduced, contrast is improved, and observation of breast lesions is convenient. BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Figure 1 A flowchart of a DBT reconstruction method based on three-dimensional coordinate displacement provided in Example 1 of the present invention;

[0087] Figure 2 A simplified three-dimensional schematic diagram of the DBT data collection scanning process;

[0088] Figure 3 for Figure 2 Front and side views of the 3D scanning process;

[0089] in, Figure 3 (a) is a front view, i.e., a schematic diagram of the XZ plane; Figure 3 (b) is a side view, i.e., a schematic diagram of the YZ plane;

[0090] Figure 4 Three representative projection angles in the collected multi-angle X-ray projection data are the 1st, 8th and 15th projection angles;

[0091] Figure 5 A comparison diagram of breast projections processed by the method of the present invention and the original method;

[0092] in, Figure 5 (a) is the focus layer output of the original method for processing breast projection. Figure 5 (b) is the focus layer output of the breast projection processed by the method of the present invention, Figure 5 (c) is Figure 5 The spectrum output diagram of (a) is Figure 5 (d) is Figure 5 (b) Spectrum output diagram;

[0093] Figure 6 This is a locally enlarged view of the subtraction between the output of the focusing layer containing only metal balls processed by the original method and the output of the focusing layer containing only metal balls processed by the method of the present invention. DETAILED DESCRIPTION

[0094] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the drawings in the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments, and they should not be understood as limitations on the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. In the description of the present invention, it should be understood that the terms used are only for descriptive purposes and cannot be understood as indicating or implying relative importance.

[0095] Combine the following Figure 1-Figure 6 The present invention describes a DBT reconstruction method, medium and electronic device based on three-dimensional coordinate displacement.

[0096] Embodiment 1: This embodiment provides a DBT reconstruction method based on three-dimensional coordinate displacement. In view of the problems of low contrast, poor focusing effect, blurred boundaries, and obvious metal artifacts caused by ignoring the three-dimensional characteristics of DBT in traditional DBT imaging methods. This embodiment uses multi-angle X-ray projection to perform DBT reconstruction through three-dimensional displacement. By transforming the projection data at different angles with the center projection position as the basis, all structures on the reconstruction plane are focused, and multi-angle projection displacement correction and magnification standardization are realized, so that the reconstructed focus layer structure is clearer, metal artifacts are significantly reduced, and contrast is significantly improved, which is more beneficial for the observation of breast lesions.

[0097] like Figure 1 As shown, the specific steps include:

[0098] S1. Acquire multi-angle X-ray projection data (X-ray projection data, such as Figure 4 As shown), negative log gamma mapping is performed on the acquired projection data; the processing expression is:

[0099]

[0100] in, These are the projection images before and after preprocessing respectively; is the gamma index, which takes a value greater than 1. Take 1.4, then the image effect is ideal.

[0101] It should be noted that this processing step takes into account that in the process of generating DBT projection data, when X-rays pass through the breast, their intensity decays exponentially. Therefore, the intensity values ​​in the input raw projection data are usually not linearly distributed, but exponentially related. And due to the particularity of the gland, it is necessary to adjust the grayscale of the projected image in a nonlinear way, and apply different degrees of adjustment in areas with different grayscale levels to meet different visual needs. Therefore, combined with the negative logarithmic transformation and gamma correction commonly used in medical image processing, this nonlinear relationship is converted into a linear or more linear relationship, and the contrast is adjusted, the grayscale is mapped, and the visual effect is enhanced.

[0102] S2, filtering the projection data in the Fourier domain by adjusting the window through a ramp filter;

[0103] In DBT imaging, the central projection is used as the reference and an odd index sequence is generated according to the central position. , where the odd index generator is:

[0104]

[0105] This odd index centers the filter at the center of the index, keeping the filter response in the frequency domain as real. is the window filter length, which is determined by a minimum power of 2 that is not less than the projection data length. The calculation expression is as follows:

[0106]

[0107] in, is the input projection image data length.

[0108] Calculate the filter coefficient corresponding to the index at the position, the calculation formula is as follows:

[0109]

[0110] in, is the filter coefficient, which is responsible for the attenuation and enhancement characteristics of the filter in the frequency domain, thereby achieving the function of generating discrete time ramp filtering;

[0111] The pixels with a pixel value of 0 in the projection data are assigned a value of 1, and then the SL filter is generated by convolution with the Singer function in the frequency domain, and then the Hanning window is applied to generate a filter with a window shape;

[0112] Use this filter As a coefficient vector:

[0113]

[0114] The expressions of the sinc function and window filtering are:

[0115]

[0116]

[0117] The bandwidth is obtained by rounding the product of the window filter length and the cutoff frequency;

[0118] Calculate the coefficient for each position in the Hanning window;

[0119] The filtered projection image is inversely Fourier transformed to transform the projection data back to the spatial domain to complete the filtering.

[0120] In this step, odd indexes are easier to process and optimize than even indexes, reducing possible boundary effects in the algorithm. The window filter length is calculated by the input breast projection image data length to ensure more efficient operation during frequency domain processing. It is significantly faster than the directly given window filter length and is adaptive to other types of DBT projection images.

[0121] S3, coordinate transformation of the preprocessed projection in three-dimensional space, aligning all projection data to the same spatial coordinate system; this step aligns the reconstruction plane, shifts the image in three dimensions, and performs fixed transformation on projection images with different magnifications to standardize them to the same magnification. The specific method is:

[0122] Pre-construct a simple DBT 3D scanning diagram, such as Figure 2 As shown, the center projection position projection is the projection of the X-rays of the scanner on the z-axis onto the XY plane through the reconstruction plane. Point, projection angle without 3D correction Projection data The angle with the z-axis The scanner uses X-rays to project onto the XY plane through the reconstruction plane. During the 3D correction process, the difference from the original correction is that the original correction is the projection angle Projection data That is, the displacement difference of the projection point relative to the origin, and the displacement difference of the three-dimensional correction is the projection angle Projection data The displacement difference from the projection of the center projection position should be distinguished. Based on the circular motion model of DBT for three-dimensional breast scanning, at different scanning angles, the center projection position projection is used as the transformation basis to calculate the displacement change of the projection data after logarithmic transformation, and obtain the displacement of the projection data along the x-axis and y-axis determined by the angle and reconstruction layer depth, that is, the distance from the center projection position projection on the x-axis and y-axis respectively.

[0123] First of all Figure 3 (a) is described by the projection on the XZ plane, and the points are For any point on the reconstruction layer, when the X-ray tube is at the center projection angle (i.e. 0°), the position of the X-ray tube is recorded as point , then the ray AP intersects the x-axis at point , the projection of point P on the x-axis is denoted as point , the projection of point P on the z-axis is denoted as point ; When the X-ray tube projection angle is When the X-ray tube is located, the position is recorded as point , then the ray BP intersects the x-axis at point , the projection of point B on the x-axis is denoted as point , the projection of point B on the z-axis is denoted as point ;

[0124] Another Figure 3 (b) is described from the YZ plane perspective, that is For any point on the reconstruction layer, when the X-ray tube is at the center projection angle (i.e. 0°), the position of the X-ray tube is recorded as point , then ray AP intersects the y-axis at point , the projection of point P on the y-axis is denoted as point , the projection of point P on the z-axis is denoted as point ; When the X-ray tube projection angle is When the X-ray tube is located, the position is recorded as point , then the ray BP intersects the y-axis at point , the projection of point B on the y-axis and the origin coincide.

[0125] in, for The distance from the tube to the center of rotation, is the depth of the reconstructed focus layer, that is, the depth of the reconstructed focus layer and The distance between the radiation detector planes, for The angle of the ray tube at the current projection relative to the central projection. The angle of the central projection is .

[0126] exist Figure 3 In the triangle constructed in (a), , according to the similarity relationship:

[0127]

[0128] Substituting the corresponding coordinates into the above two formulas, we get:

[0129]

[0130] The calculation results are:

[0131]

[0132] Figure 3 In (a), NM is when the projection angle is The displacement along the x-axis on the detection surface relative to the central projection (projection angle is 0°) is , then the expression of the displacement used for registration is as follows:

[0133]

[0134] in, Indicates the projection angle; The projection angle is Time to be translated Axis coordinates; represents the reconstruction layer depth; express The distance from the tube to the center of rotation; Indicates the vertical distance from the rotation center to the detector: It means that when the reconstruction layer depth is , the projection angle is When the projection data is The displacement of the axis;

[0135] It should be noted here that the two-dimensional correction in the prior art is the above The relative displacement on the axis, but because DBT scanning is a three-dimensional surround type, the two-dimensional correction and initial correction are not rigorous enough in displacement correction, and the influence of metal artifacts is not considered in detail. Therefore, this paper proposes a three-dimensional correction displacement formula on this basis, that is, add the following part:

[0136] exist Figure 3 In the triangle constructed in (b), , according to the similarity relationship, we can get:

[0137]

[0138] Substituting the corresponding coordinates into the above two formulas, we get:

[0139]

[0140] The calculation results are:

[0141]

[0142] Figure 3 In (b), NM is when the projection angle is The displacement along the y-axis on the detection surface relative to the central projection (projection angle is 0°) is , then the expression of the displacement used for registration is as follows:

[0143]

[0144] in, Indicates the projection angle; The projection angle is Time to be translated Axis coordinates; Represents the reconstruction layer depth: It means that when the reconstruction layer depth is , the projection angle is When the projection data is The displacement of the axis; express The distance from the X-ray tube to the center of rotation, Indicates the vertical distance from the rotation center to the detector.

[0145] The real values ​​of the displacements on the x-axis and y-axis calculated on two different planes are converted into pixel units. The conversion expression is:

[0146]

[0147]

[0148] in, Express Round down, that is, the closest and not greater than integer; and Represents the pixel displacement on the x-axis and y-axis respectively;

[0149] Calculate the change in projection size as:

[0150]

[0151]

[0152]

[0153]

[0154] The overall range of variation required in the projected image data size is:

[0155]

[0156]

[0157] According to the displacement of the projection data along the x-axis and the y-axis, a fixed transformation is performed on the projection images with different magnifications so as to standardize them to the same magnification and obtain the corrected projection data.

[0158] Therefore, the projection angle The projection data of is expressed as follows after translation:

[0159]

[0160] in, Represents the projection data after translation; Represents the x-axis coordinate of the translated projection data, , Represents the x-axis coordinate of the original projection; Represents the y-axis coordinate of the translated projection data, , Represents the y-axis coordinate of the original projection; Represents the reconstruction layer depth.

[0161] S4, performing gradient-based interpolation correction on rows and columns respectively to obtain corrected projection data;

[0162] The null values ​​of the corrected projection image data are interpolated based on gradients in rows and columns. In order to avoid transitional access and the situation where there is no gradient difference in boundary rows and columns, the following restrictions are added:

[0163]

[0164]

[0165] in, , are the rows and columns traversed in the projected image, is the number of columns, is the number of rows.

[0166] The shifted projection data is converted into a matrix Indicates that for the matrix Each line of , each column Sum, we get and :

[0167]

[0168]

[0169] in, is a matrix In the Row, No. Elements of the column; the resulting row and column vectors , , is the number of matrix rows, is the number of matrix columns;

[0170] if , then The row has missing values; if , then Column has missing values; Find The row and column indices that are zero, determine the rows and columns corresponding to the missing values. and ,in is a set containing all null value row and column indices;

[0171] The gradient is calculated by the difference between the previous and next data, and then the missing data is filled by gradient interpolation based on the calculated gradient value. The expression of interpolation filling is:

[0172]

[0173]

[0174]

[0175]

[0176] in, , Respectively represent the pixel gradient interpolation of the adjacent upper and lower rows or adjacent left and right columns of the pixel that needs to be filled in the row and column; The index in the input projection image array is Elements of Represents the first Row, No. Elements of a column; is a loop variable used to iterate over the all-zero row index;

[0177] For missing rows or columns at the boundary, fill them with the directly adjacent non-missing data:

[0178] The top missing row ;

[0179] The missing line at the bottom ;

[0180] Leftmost missing column ;

[0181] The rightmost missing column ;

[0182] in Indicates the line; Represents a column; Represents the pixel value of a row or column; , Represents the maximum value of the row and column respectively.

[0183] S5. The corrected projection data is processed by the autoregressive model (AR) source moving average separation method to obtain a tomographic image of any depth layer; that is, the dynamic characteristics of the signal are estimated by the AR source moving average separation method, the separation matrix is ​​optimized for multiple iterations to process the data, and the reconstructed layer is output. The specific method is:

[0184] Obtain corrected projection data and initialize the autoregressive model;

[0185] First, the number of iterations, the condition number of the control weight matrix, and the parameter used to limit the root of the autoregressive model are given. In this embodiment, the number of iterations is set to [3,5]; the condition number of the control weight matrix is ​​set to 0.000001; and the parameter used to limit the root of the autoregressive model is set to 0.99. Based on the application of Blind Source Separation (BSS) to separate and reconstruct any depth layer, each element in the column vector in the projection image matrix is ​​first subtracted from the mean of the column vector to make the data concentrated near zero, and a sliding average model is added to estimate the dynamic characteristics of the signal. The formula is as follows:

[0186]

[0187] in, The size is The mean vector of is a column vector of all 1s, is a row vector of all 1s;

[0188] Then, given the autoregressive order, the time period length , calculate the autocorrelation matrix :

[0189]

[0190] in, Representation Matrix From column 1 to column The column vector of the concatenation of columns; Representation Matrix No. Column to The column vector of the concatenation of columns; represents the vector outer product; Indicates optional matrix;

[0191] Then calculate the average value of each specified submatrix region and its transpose, and assign this average value back to the atomic matrix; according to the input matrix and the current estimate, using an approximate joint diagonalization algorithm with uniform weights to unmix the matrix , the diagonalized Matrix; Diagonalize using the approximate joint diagonalization algorithm with non-uniform weights The matrix is ​​optimized within the number of iterations to obtain the separated signal matrix , the expression is:

[0192]

[0193] Among them, the matrix is the coefficient matrix for processing the diagonal elements of the matrix; For the matrix The transposed matrix of .

[0194] Finally, the signal is reconstructed by separating the signal and processing the data with autoregressive moving average.

[0195] S6. Perform quantitative evaluation and adaptive histogram equalization on the tomographic image, and output an enhanced image. In this embodiment, the reconstructed layer is firstly subjected to ISR visual evaluation, and then adaptive dynamic local contrast enhancement is performed after local weighted averaging through mask operation to obtain the final output; the specific method is:

[0196] ISR Visualization Evaluation:

[0197] Acquire a reconstructed breast image;

[0198] The separated signals are derived from the correlation matrix The coefficients of the autoregressive (AR) signal source are calculated and accumulated:

[0199]

[0200] And calculate matrix:

[0201]

[0202] Finally, calculate matrix:

[0203]

[0204] Representation Matrix No. Line The elements of the column, yes The first Line The elements of the column, is a matrix The Line The elements of the column, is a A column vector of all 1s, Representation Matrix Line 1 of The elements of the column, yes The transpose of is a The row vector of .

[0205] Will The matrix is ​​visualized using a heatmap. If all data on the heatmap are minimum values, the reconstruction effect is better.

[0206] Adaptive dynamic local contrast enhancement is performed by local weighted averaging through mask operation:

[0207] Define a mask matrix , perform weighted averaging on the image: for each pixel in the image . The weighted sum of the pixels in the area and the mask matrix is ​​calculated as:

[0208]

[0209] in, and Respectively represent the current pixel The surrounding pixel offset, and Represent the rows and columns of pixels respectively. Stored at the corresponding position in the new array.

[0210] Then, the reconstruction layer is divided into multiple For a small local block, calculate its local contrast mean and compare it with the threshold. The threshold is selected according to the histogram. The calculation formula of the contrast mean is as follows:

[0211]

[0212]

[0213] in, is the local contrast value at each pixel position in the image, is the total number of pixels in the local image, that is ; is the image at position The gray value at ; It is a field The mean of the grayscale values ​​of the inner pixels; Representation field The number of pixels inside.

[0214] If the local contrast mean is higher than the threshold, a smaller area (4x4 block) is used during enhancement to enhance the details of the image; otherwise, a larger block area (8x8 block) is used to avoid over-enhancement. Finally, the grayscale value in the statistical window is equal to The number of pixels to calculate the histogram :

[0215]

[0216] in, is the Kreonecker function, which is 1 when the parameters are equal and 0 otherwise. Calculate the cumulative histogram within the window to represent the gray level The cumulative distribution of . Then the histogram Perform cumulative summation to obtain a cumulative histogram :

[0217]

[0218] Finally, the cumulative distribution function is used To map each pixel value within the window and output the enhanced image.

[0219] Embodiment 2: This embodiment provides a non-transitory computer-readable storage medium having computer instructions stored thereon, the computer instructions causing a computer to execute a DBT reconstruction method based on three-dimensional coordinate displacement, the method comprising the following steps:

[0220] S1, acquiring X-ray projection data at multiple angles, and performing negative log-gamma mapping on the acquired projection data;

[0221] S2, filtering the projection data in the Fourier domain by adjusting the window through a ramp filter;

[0222] S3, performing three-dimensional correction on the projection data, aligning all the projection data to the same spatial coordinate system;

[0223] S4, performing image interpolation processing to obtain corrected projection data;

[0224] S5, processing the corrected projection data by using an autoregressive source moving average separation method to obtain a tomographic image of any depth layer;

[0225] S6. Perform quantitative evaluation and adaptive histogram equalization on the tomographic image and output an enhanced image.

[0226] Embodiment 3: This embodiment provides an electronic device, which may include: a processor, a communications interface, a memory, and a communication bus, wherein the processor, the communications interface, and the memory communicate with each other via the communication bus. The processor may call the logic instructions in the memory to execute a DBT reconstruction method based on three-dimensional coordinate displacement, the method comprising the following steps:

[0227] S1, acquiring X-ray projection data at multiple angles, and performing negative log-gamma mapping on the acquired projection data;

[0228] S2, filtering the projection data in the Fourier domain by adjusting the window through a ramp filter;

[0229] S3, performing three-dimensional correction on the projection data, aligning all the projection data to the same spatial coordinate system;

[0230] S4, performing image interpolation processing to obtain corrected projection data;

[0231] S5, processing the corrected projection data by using an autoregressive source moving average separation method to obtain a tomographic image of any depth layer;

[0232] S6. Perform quantitative evaluation and adaptive histogram equalization on the tomographic image and output an enhanced image.

[0233] In addition, the logic instructions in the above-mentioned memory can be implemented in the form of software functional units and can be stored in a computer-readable storage medium when sold or used as an independent product. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, etc., which can store program code.

[0234] Embodiment 4: This embodiment provides a computer program product, which includes a computer program. The computer program can be stored on a non-transitory computer-readable storage medium. When the computer program is executed by a processor, the computer can perform a DBT reconstruction method based on three-dimensional coordinate displacement, which includes the following steps:

[0235] S1, acquiring X-ray projection data at multiple angles, and performing negative log-gamma mapping on the acquired projection data;

[0236] S2, filtering the projection data in the Fourier domain by adjusting the window through a ramp filter;

[0237] S3, performing three-dimensional correction on the projection data, aligning all the projection data to the same spatial coordinate system;

[0238] S4, performing image interpolation processing to obtain corrected projection data;

[0239] S5, processing the corrected projection data by using an autoregressive source moving average separation method to obtain a tomographic image of any depth layer;

[0240] S6. Perform quantitative evaluation and adaptive histogram equalization on the tomographic image and output an enhanced image.

[0241] The DBT reconstruction of the present invention was compared with the original method, and the validation test results showed that:

[0242] like Figure 5 As shown, we select the output images of layers with better focusing effects for comparison. Figure 5 (a) is the focus layer output of the original method for processing breast projection. Figure 5 (b) is the focus layer output of the breast projection processed by the method of the present invention. It can be seen that the contrast of the reconstructed points obtained using the original displacement formula is too low, and its boundary points become blurred and have large shadows; while the image corrected by the present invention not only has better focus, but also has significantly increased contrast and clearer boundaries. It can be verified that the displacement formula of the present invention is significantly higher than the original displacement formula and the two-dimensional correction formula in terms of registration accuracy, reducing the blur and artifacts of the focus point. This makes the contrast higher, the texture of the focus reconstruction layer clearer, the boundaries clearer, and less prone to artifacts.

[0243] In addition, comparing the spectrum output diagrams of the two, Figure 5 (c) and Figure 5 (d), Figure 5 (c) looks messy and has a lot of high-frequency noise on its vertical axis. In contrast, Figure 5The spectral output diagram corresponding to (d) is not only clear and neat as a whole but also has less noise on the horizontal and vertical axes, indicating that after the three-dimensional displacement correction is performed by the method of the present invention, all structures on the reconstructed plane are focused, multi-angle projection displacement correction and magnification standardization are achieved, making the reconstructed focus layer structure clearer.

[0244] In order to better observe the correctness of the displacement correction of the present invention, 15 projection images containing only metal balls are merged twice and re-output as 15 angles containing only metal balls and merging two layers of projection images. It should be noted that removing breast tissue can achieve better observation results when determining the displacement correction. The output images corrected by the original method and the output images corrected by the three-dimensional displacement correction of the present invention are subtracted as follows Figure 6 As shown, the local metal ball near the central axis is enlarged, and it can be clearly observed that the farther away from the central axis in the horizontal direction, the greater the difference, that is, the greater the difference is to the right. This shows that after the displacement correction is performed by the method of the present invention, the projection registration is more accurate.

[0245] The device embodiments described above are merely illustrative, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. Ordinary technicians in this field can understand and implement it without paying creative labor.

[0246] Through the description of the above implementation methods, those skilled in the art can clearly understand that each implementation method can be implemented by means of software plus a necessary general hardware platform, and of course, can also be implemented by hardware. Based on this understanding, the above technical solution is essentially or the part that contributes to the prior art can be embodied in the form of a software product, and the computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a disk, an optical disk, etc., including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0247] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A DBT reconstruction method based on three-dimensional coordinate displacement, characterized in that: include: Acquire multi-angle X-ray projection data and pre-process it; the specific steps include: Performing negative log-gamma mapping on the acquired projection data; In the Fourier domain, the projection data is filtered by adjusting the window through a ramp filter; The preprocessed projection data is transformed into coordinates in three-dimensional space. The specific method is as follows: The projection data of the projection angle θ is expressed as follows after translation: I' θ =(x′ θ ,y′ θ ,z0) Among them, I′ θ represents the projection data after translation; x′ θ Represents the x-axis coordinate of the projected data after translation, x′ θ =x θ -shift(x θ ,z0),x θ Represents the x-axis coordinate of the original projection; y′ θ Represents the y-axis coordinate of the projected data after translation, y′ θ =y θ -shift(y θ ,z0),y θ Represents the y-axis coordinate of the original projection; z0 represents the depth of the reconstruction layer; shift(x θ ,z0) represents the displacement of the projection data along the x-axis when the reconstruction layer depth is z0 and the projection angle is θ; shift(y θ ,z0) represents the displacement of the projection data along the y-axis when the reconstruction layer depth is z0 and the projection angle is θ; shift(x θ ,z0) and shift(y θ , z0) is determined by the angle θ, the reconstruction layer depth z0, the distance L from the X-ray tube to the rotation center, and the vertical distance DDR from the rotation center to the detector; The expression of the displacement of the projection data along the x-axis is: The expression of the displacement of the projection data along the y-axis is: The calculated real value of the displacement along the x-axis and y-axis is converted into pixel units. The conversion formula is: in, Indicates rounding down; shift(x)_new and shift(y)_new represent the pixel displacement on the x-axis and y-axis respectively; Calculate the change in projection size as: spreadupper(x)=|max(max(shift(x θ ,z0)))| spreadlower(x)=|min(min(shift(x θ ,z0)))| spreadupper(y)=|max(max(shift(y θ ,z0)))| spreadlower(y)=|min(min(shift(y θ ,z0)))| The projection data varies in image size as follows: spreadsum_x=spreadupper(x)+spreadlower(x) spreadsum_y=spreadupper(y)+spreadlower(y) After performing coordinate and size transformation on the collected multi-angle projections, the displaced projection data are obtained; Then, image interpolation processing is performed to obtain corrected projection data; the specific steps are: The projection data after displacement is represented by matrix T, and each row i and column j of matrix T are summed to obtain S i and S j : Among them, T ij is the element in the i-th row and j-th column of matrix T; the resulting row and column vectors S row =[S1, S2, ..., S m ],S col =[S1, S2, ..., S n ], m is the number of matrix rows, n is the number of matrix columns; If S i =0, then the i-th row has missing values; if S j = 0, then the jth column has missing values; find S i , S j The row and column indices that are zero, determine the rows and columns corresponding to the missing values ​​as C row =[i|S i =0} and C col =[j|S j =0}, where C row , C col is a set containing all null value row and column indices; The gradient is calculated by the difference between the previous and next data, and then the missing data is filled by gradient interpolation based on the calculated gradient value. The expression of interpolation filling is: temp(j,C col(k+1) )=temp(j,C col(k+1) -1)+gradient col Among them, gradient row 、gradient col Respectively represent the pixel gradient interpolation of the adjacent upper and lower rows or adjacent left and right columns of the pixel to be filled in the row or column; temp(C ron(k+1) , j) indicates that the index in the input projection image array is C ron(k+1) , element of j, temp(j,C col(k+1) ) represents the jth row and the Cth row in the input projection image array. col(k+1) The elements of the column; k is a loop variable used to iterate over the all-zero row index; For missing rows or columns at the boundary, fill them with the directly adjacent non-missing data: The top missing row (r=1): T(1,j)=T(2,j); The bottom missing row (r=M): T(M,j)=T(M-1,j); Leftmost missing column (c=1): T(i,1)=T(i,2); The missing rightmost column (c=N): T(i,N)=T(i,N-1); Where r represents the row; c represents the column; T(·) represents the pixel value of a row and a column; M and N represent the maximum values ​​of the row and column respectively; The corrected projection data is processed by using an autoregressive source moving average separation method to obtain a tomographic image of any depth; the steps include: Obtain the corrected projection data and initialize the autoregressive model; Set the number of iterations of the autoregressive model, the condition number that controls the weight matrix, and the parameters used to constrain the roots of the autoregressive model; Based on the application of blind source separation method to separate and reconstruct arbitrary depth layers, the mean of the column vector is subtracted from each element in the column vector of the projection image matrix to make the data concentrated near zero, and the sliding average model is added to estimate the dynamic characteristics of the signal. The formula is as follows: in, is the mean vector of size d×1, 1 c is a column vector of all 1s, is a row vector of all 1s; Then, according to the given autoregressive order and time period length T, the autocorrelation matrix R is calculated hij : Among them, x G The column vector representing the concatenation of the 1st column to the Tth column of the matrix x; x H The column vector representing the concatenation of the index-th column to the index+T-1-th column of the matrix x; W represents the vector outer product; pre represents an optional K×K matrix; Then calculate the average value of each specified submatrix region and its transpose, and assign this average value back to the atomic matrix; according to the input matrix M and the current estimate, use the approximate joint diagonalization algorithm with uniform weights to unmix the matrix W hij , the diagonalized M is calculated by adjusting the weights through the correlation between data components i matrix; using the non-uniform weighted approximate joint diagonalization algorithm to diagonalize M i The matrix is ​​optimized within the number of iterations to obtain the separated signal matrix The expression is: Among them, matrix A is the coefficient matrix used to process the diagonal elements of the matrix; W′ hij is the matrix W hij The transposed matrix of Finally, the tomographic image is reconstructed by separating the signals and the data processed by autoregressive moving average; The tomographic image is quantitatively evaluated and adaptively histogram equalized to output an enhanced image; the steps include: Perform autoregressive model analysis on the separated signal source, calculate the cumulative value of each coefficient in the correlation matrix, and calculate the ISR matrix based on this; The calculated ISR matrix is ​​visualized by heatmap, and the ISR visualization evaluation is performed on the tomographic image; For images with poor evaluation, a mask matrix is ​​defined, and for each pixel in the image, the weighted sum of the pixels in its surrounding area and the mask matrix is ​​calculated; Calculate the local contrast mean of the image and compare it with the threshold. If the local contrast mean is higher than the threshold, use a smaller area for enhancement, otherwise use a larger area for enhancement. The histogram is calculated by counting the number of pixels with a specific grayscale value in the window, and the histogram is accumulated and summed to obtain a cumulative histogram representing the cumulative distribution of grayscale levels. This is used to map each pixel value in the window and output an enhanced image.

2. The DBT reconstruction method based on three-dimensional coordinate displacement according to claim 1, characterized in that: In the image interpolation processing of the displacement projection data, the following limiting conditions are added: C LDM(YZG) ∈(1,size(temp,2)-1) c CDE(YZG) ∈(1,size(temp,1)-1) Among them, C LDM(YZG) , C CDE(YZG) are the rows and columns traversed in the projected image, size(temp,2) is the number of columns, and size(temp,1) is the number of rows.

3. A non-transitory computer-readable storage medium having computer instructions stored thereon, the computer instructions enabling a computer to execute the DBT reconstruction method based on three-dimensional coordinate displacement as claimed in one of claims 1 or 2.

4. An electronic device comprising: A processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory communicate with each other via the communication bus, and the processor calls the logic instructions in the memory to execute the DBT reconstruction method based on three-dimensional coordinate displacement as described in one of claims 1 or 2.

Images